Browse Source

Added basic Event/Track cuts

PeterParfenov 8 months ago
parent
commit
3324cf9f01
2 changed files with 169 additions and 11 deletions
  1. 161 6
      macro/FemtoDstAnalyzer.C
  2. 8 5
      scripts/GenerateLists.sh

+ 161 - 6
macro/FemtoDstAnalyzer.C

@@ -40,6 +40,30 @@ R__LOAD_LIBRARY(/mnt/pool/rhic/1/nigmatkulov/soft/StFemtoEvent/libStFemtoDst)
 //          of a name.lis(t) files, that contains a list of
 //          name1.FemtoDst.root, name2.FemtoDst.root, ... files
 
+// Used functions (see them below)
+Bool_t isGoodEvent(StFemtoEvent *const &event);
+Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
+Bool_t isGoodTrackFlow(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
+Float_t BBC_GetPhi(const Int_t eastWest, const Int_t tileId);
+
+// Used constants (global)
+const Int_t energy = 200.;
+
+// Used constants for the event cut
+const std::map<Float_t, Double_t> cutVtxZEnergy = {{7.7, 70.}, {11.5, 50.}, {19.6, 70}, {27., 70.}, {39., 40.}, {62., 40.}, {200., 30.}};
+const Double_t cutVtxR = 2.;
+const Double_t cutVpdVz = 3.;
+
+// Used constants for the track cut
+const std::map<Float_t, Double_t> cutDCA = {{7.7, 1.}, {11.5, 1.}, {19.6, 1.}, {27., 1.}, {39., 1.}, {62., 1.}, {200., 3.}};
+const Double_t cutEta = 1.;
+const Int_t    cutNhits = 15;
+const Double_t cutNhitsPoss = 0.;
+const Double_t cutNhitsRatio = 0.51;
+const std::map<Float_t, Double_t> cutPtMin = {{7.7, 0.2}, {11.5, 0.2}, {19.6, 0.2}, {27., 0.2}, {39., 0.2}, {62., 0.2}, {200., 0.15}};
+const Double_t cutPtMax = 2.;
+const Double_t cutPMax = 10.;
+
 //_________________
 void FemtoDstAnalyzer(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",const Char_t *oFileName = "oTest.root") {
 
@@ -92,11 +116,29 @@ void FemtoDstAnalyzer(const Char_t *inFile = "st_physics_12150008_raw_4030001.fe
       break;
     }
 
+    // Event selection
+    if ( !isGoodEvent(event) ) continue;
+
     // Return primary vertex position
     TVector3 pVtx = event->primaryVertex();
 
-    // Reject vertices that are far from the central membrane along the beam
-    if( TMath::Abs( pVtx.Z() ) > 200. ) continue;
+    // Going through the inner BBC rings
+    Double_t phiBBC_East, phiBBC_West;
+    Double_t adcBBC_East, adcBBC_West;
+    for (int iTile=0; iTile<16; iTile++)
+    {
+      // Get BBC East phi angle
+      phiBBC_East = BBC_GetPhi(0,iTile);
+
+      // Get BBC West phi angle
+      phiBBC_West = BBC_GetPhi(1,iTile);
+
+      // Get energy deposition from BBC East iTile module
+      adcBBC_East = event->bbcAdcEast(iTile);
+
+      // Get energy deposition from BBC West iTile module
+      adcBBC_West = event->bbcAdcWest(iTile);
+    }
 
     // Track analysis
     Int_t nTracks = dst->numberOfTracks();
@@ -112,10 +154,13 @@ void FemtoDstAnalyzer(const Char_t *inFile = "st_physics_12150008_raw_4030001.fe
       // Must be a primary track
       if ( !femtoTrack->isPrimary() ) continue;
 
-      // Simple single-track cut
-      if( femtoTrack->gMom().Mag() < 0.1 || femtoTrack->gDCA(pVtx).Mag() > 3. ) {
-	       continue;
-      }
+
+      //Track selection for EP determination in TPC
+      if (!isGoodTrack(femtoTrack, energy, pVtx)) continue;
+
+      //Track selection for flow calculation in TPC
+      if (!isGoodTrackFlow(femtoTrack, energy, pVtx)) continue;
+
 
       // Check if track has TOF signal
       if ( femtoTrack->isTofTrack() ) {
@@ -131,3 +176,113 @@ void FemtoDstAnalyzer(const Char_t *inFile = "st_physics_12150008_raw_4030001.fe
   std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
 	          << std::endl;
 }
+
+Bool_t isGoodEvent(StFemtoEvent *const &event)
+{
+  if (!event) return false;
+  if (event == nullptr) return false;
+
+  if (event->primaryVertex().Perp() > cutVtxR) return false;
+  if (TMath::Abs(event->primaryVertex().Z()) > cutVtxZEnergy.at(energy)) return false;
+
+  if ((energy == 200.) && TMath::Abs(event->primaryVertex().Z() - event->vpdVz()) > cutVpdVz) return false;
+
+  return true;
+}
+
+Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
+{
+  if (!track) return false;
+  // if (!track->isPrimary()) return false;
+  if (track->nHitsFit() < cutNhits) return false;
+  if (track->nHitsPoss() <= cutNhitsPoss) return false;
+  if ((Double_t) track->nHitsFit()/track->nHitsPoss() < cutNhitsRatio) return false;
+  if (TMath::Abs(track->eta()) >= cutEta) return false;
+  
+  if (track->pt() <= cutPtMin.at(_energy)) return false;
+  if (track->pt() > cutPtMax) return false;
+  if (track->ptot() > cutPMax) return false;
+  if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy)) return false;
+  return true;
+}
+
+Bool_t isGoodTrackFlow(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
+{
+  if (!track) return false;
+  // if (!track->isPrimary()) return false;
+  if (track->nHitsFit() < cutNhits) return false;
+  if (track->nHitsPoss() <= cutNhitsPoss) return false;
+  if ((Double_t) track->nHitsFit()/track->nHitsPoss() < cutNhitsRatio) return false;
+  if (TMath::Abs(track->eta()) >= cutEta) return false;
+  
+  if (track->pt() <= cutPtMin.at(_energy)) return false;
+  //if (track->pt() > cutPtMax) return false;
+  if (track->ptot() > cutPMax) return false;
+  if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy)) return false;
+  return true;
+}
+
+//--------------------------------------------------------------------------------------------------------------------------//
+//this function simply connects the gain values read in to the BBC azimuthal distribution
+//since tiles 7 and 9 (+ 13 and 15) share a gain value it is ambiguous how to assign the geometry here
+//I prefer assigning the angle between the tiles thus "greying out" the adcs. 
+//Others have assigned all of the adc to one (exclusive) or the the other. 
+Float_t BBC_GetPhi(const Int_t eastWest, const Int_t tileId)
+{
+  //float GetPhiInBBC(int eastWest, int bbcN) { //tileId=0 to 23
+  const float Pi = TMath::Pi();
+  const float Phi_div = Pi / 6;
+  float bbc_phi = Phi_div;
+  switch(tileId) {
+    case 0: bbc_phi = 3.*Phi_div;
+  break;
+    case 1: bbc_phi = Phi_div;
+  break;
+    case 2: bbc_phi = -1.*Phi_div;
+  break;
+    case 3: bbc_phi = -3.*Phi_div;
+  break;
+    case 4: bbc_phi = -5.*Phi_div;
+  break;
+    case 5: bbc_phi = 5.*Phi_div;
+  break;
+    //case 6: bbc_phi= (mRndm.Rndm() > 0.5) ? 2.*Phi_div:4.*Phi_div;	//tiles 7 and 9 are gained together we randomly assign the gain to one XOR the other
+    case 6: bbc_phi = 3.*Phi_div;
+  break;
+    case 7: bbc_phi = 3.*Phi_div;
+  break;
+    case 8: bbc_phi = Phi_div;
+  break;
+    case 9: bbc_phi = 0.;
+  break;
+    case 10: bbc_phi = -1.*Phi_div;
+  break;
+    //case 11: bbc_phi = (mRndm.Rndm() > 0.5) ? -2.*Phi_div:-4.*Phi_div;	//tiles 13 and 15 are gained together
+    case 11: bbc_phi = -3.*Phi_div;
+  break;
+    case 12: bbc_phi = -3.*Phi_div;
+  break;
+    case 13: bbc_phi = -5.*Phi_div;
+  break;
+    case 14: bbc_phi = Pi;
+  break;
+    case 15: bbc_phi = 5.*Phi_div;
+  break;
+  }
+
+  //if we're looking at the east BBC we need to flip around x in the STAR coordinates, 
+  //a line parallel to the beam would go through tile 1 on the W BBC and tile 3 on the 
+  if(0 == eastWest){
+    if (bbc_phi > -0.001){ //this is not a >= since we are talking about finite adcs -- not to important
+      bbc_phi = Pi - bbc_phi;
+    }
+    else {
+      bbc_phi= -Pi - bbc_phi;
+    }
+  }
+
+  if(bbc_phi < 0.0) bbc_phi += 2.*Pi;
+  if(bbc_phi > 2.*Pi) bbc_phi -= 2.*Pi;
+
+  return bbc_phi;
+}

+ 8 - 5
scripts/GenerateLists.sh

@@ -4,6 +4,8 @@
 INPUT_DIR=$1
 #-Set the maximum number of files in 1 list
 NUM_DIVIDES=$2
+#-Name pattern for output lists
+NAME_PATTERN=$3
 
 #-Example usage
 # . GenerateLists.sh /mnt/pool/rhic/2/nigmatkulov/femtoDst/auau/200gev/12135/ 100
@@ -11,7 +13,8 @@ NUM_DIVIDES=$2
 CURRENT_DIR=${PWD}
 OUTPUT_DIR="../lists"
 
-TOTAL_NUM_FILES=`ls $INPUT_DIR/*.femtoDst.root | wc -l`
+#TOTAL_NUM_FILES=`ls $INPUT_DIR/*.femtoDst.root | wc -l`
+TOTAL_NUM_FILES=`find $INPUT_DIR -type f -name '*.femtoDst.root' | wc -l`
 echo "Total number of DST files: ${TOTAL_NUM_FILES}"
 
 MULT=$((TOTAL_NUM_FILES / NUM_DIVIDES))
@@ -21,12 +24,12 @@ echo "Mult: $MULT, Residue: $RESIDUE"
 
 for i in `seq 1 $MULT`
 do
-  OUTPUT_FILE=$OUTPUT_DIR/StRuns${i}.list
+  OUTPUT_FILE=$OUTPUT_DIR/${NAME_PATTERN}${i}.list
   echo $OUTPUT_FILE
-  ls $INPUT_DIR/*.femtoDst.root | head -n$((i*NUM_DIVIDES)) | tail -n$NUM_DIVIDES &> $OUTPUT_FILE
+  find $INPUT_DIR -type f -name '*.femtoDst.root' | head -n$((i*NUM_DIVIDES)) | tail -n$NUM_DIVIDES &> $OUTPUT_FILE
 done
 
 #-Generate last filelist for residual files
-OUTPUT_FILE=$OUTPUT_DIR/StRuns$((i+1)).list
+OUTPUT_FILE=$OUTPUT_DIR/${NAME_PATTERN}$((i+1)).list
 echo $OUTPUT_FILE
-ls $INPUT_DIR/*.femtoDst.root |  tail -n$RESIDUE &> $OUTPUT_FILE
+find $INPUT_DIR -type f -name '*.femtoDst.root' |  tail -n$RESIDUE &> $OUTPUT_FILE