adding analysis option MCAODrec = AODs from MC with using MCtruth information for...
authormiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 7 Aug 2013 12:06:01 +0000 (12:06 +0000)
committermiweber <miweber@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 7 Aug 2013 12:06:01 +0000 (12:06 +0000)
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx
PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.h
PWGCF/EBYE/macros/AddTaskBalancePsiCentralityTrain.C
PWGCF/EBYE/macros/configBalanceFunctionPsiAnalysis.C
PWGCF/EBYE/macros/runBalanceFunctionPsi.C

index a372843..7254286 100755 (executable)
@@ -54,7 +54,8 @@ ClassImp(AliAnalysisTaskBFPsi)
 \r
 //________________________________________________________________________\r
 AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) \r
-: AliAnalysisTaskSE(name), \r
+: AliAnalysisTaskSE(name),\r
+  fArrayMC(0), //+++++++++++++\r
   fBalance(0),\r
   fRunShuffling(kFALSE),\r
   fShuffledBalance(0),\r
@@ -70,7 +71,7 @@ AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name)
   fHistListPIDQA(0),\r
   fHistEventStats(0),\r
   fHistCentStats(0),\r
-  fHistCentStatsUsed(0),    //++++++++++++++++++++++++++\r
+  fHistCentStatsUsed(0),\r
   fHistTriggerStats(0),\r
   fHistTrackStats(0),\r
   fHistVx(0),\r
@@ -757,7 +758,7 @@ Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){
 \r
     //Centrality stuff \r
     if(fUseCentrality) {\r
-      if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD") { //centrality in AOD header\r
+      if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec") { //centrality in AOD header   //+++++++++++\r
        AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
        if(header){\r
          gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
@@ -961,7 +962,7 @@ Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){
   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 \r
   if(fEventClass == "Centrality"){\r
-    if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD") { //centrality in AOD header\r
+    if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header  //++++++++++++++\r
       AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
       if(header){\r
         gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
@@ -1332,6 +1333,125 @@ TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gC
   }//MCAOD\r
   //==============================================================================================================\r
 \r
+  //==============================================================================================================\r
+  else if(gAnalysisLevel == "MCAODrec") {\r
+    \r
+    /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());\r
+    if (!fAOD) {\r
+      printf("ERROR: fAOD not available\n");\r
+      return;\r
+      }*/\r
+    \r
+    fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName())); \r
+    if (!fArrayMC) { \r
+       AliError("No array of MC particles found !!!");\r
+    }\r
+\r
+    AliMCEvent* mcEvent = MCEvent();\r
+    if (!mcEvent) {\r
+       AliError("ERROR: Could not retrieve MC event");\r
+    }\r
+     \r
+    for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
+      AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
+      if (!aodTrack) {\r
+       AliError(Form("Could not receive track %d", iTracks));\r
+       continue;\r
+      }\r
+      \r
+      for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
+       fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
+      }\r
+      if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
+      \r
+      vCharge = aodTrack->Charge();\r
+      vEta    = aodTrack->Eta();\r
+      vY      = aodTrack->Y();\r
+      vPhi    = aodTrack->Phi();// * TMath::RadToDeg();\r
+      vPt     = aodTrack->Pt();\r
+      \r
+      Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
+      Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
+            \r
+      // Kinematics cuts from ESD track cuts\r
+      if( vPt < fPtMin || vPt > fPtMax)      continue;\r
+      if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
+      \r
+      // Extra DCA cuts (for systematic studies [!= -1])\r
+      if( fDCAxyCut != -1 && fDCAzCut != -1){\r
+       if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
+         continue;  // 2D cut\r
+       }\r
+      }\r
+      \r
+      // Extra TPC cuts (for systematic studies [!= -1])\r
+      if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
+       continue;\r
+      }\r
+      if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
+       continue;\r
+      }\r
+\r
+      //Exclude resonances\r
+      if(fExcludeResonancesInMC) {\r
+       \r
+       Bool_t kExcludeParticle = kFALSE;\r
+\r
+       Int_t label = TMath::Abs(aodTrack->GetLabel());\r
+       AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);\r
+      \r
+        if (AODmcTrack){ \r
+         //if (AODmcTrack->IsPhysicalPrimary()){\r
+         \r
+         Int_t gMotherIndex = AODmcTrack->GetMother();\r
+         if(gMotherIndex != -1) {\r
+           AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
+           if(motherTrack) {\r
+             Int_t pdgCodeOfMother = motherTrack->GetPdgCode();\r
+             if(pdgCodeOfMother == 113  // rho0\r
+                || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
+                // || pdgCodeOfMother == 221  // eta\r
+                // || pdgCodeOfMother == 331  // eta'\r
+                // || pdgCodeOfMother == 223  // omega\r
+                // || pdgCodeOfMother == 333  // phi\r
+                || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0\r
+                // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*\r
+                // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*\r
+                || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r
+                || pdgCodeOfMother == 111  // pi0 Dalitz\r
+                ) {\r
+               kExcludeParticle = kTRUE;\r
+             }\r
+           }\r
+         }\r
+       }       \r
+       //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
+       if(kExcludeParticle) continue;\r
+      }\r
+      \r
+      // fill QA histograms\r
+      fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
+      fHistDCA->Fill(dcaZ,dcaXY);\r
+      fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
+      fHistPt->Fill(vPt,gCentrality);\r
+      fHistEta->Fill(vEta,gCentrality);\r
+      fHistRapidity->Fill(vY,gCentrality);\r
+      if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
+      else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
+      fHistPhi->Fill(vPhi,gCentrality);\r
+      if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);                 \r
+      else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
+      \r
+      //=======================================correction\r
+      Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
+      //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
+      \r
+      // add the track to the TObjArray\r
+      tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
+    }//track loop\r
+  }//MCAODrec\r
+  //==============================================================================================================\r
+\r
   else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD\r
 \r
     AliESDtrack *trackTPC   = NULL;\r
index 716d333..6100b3d 100755 (executable)
@@ -172,7 +172,8 @@ class AliAnalysisTaskBFPsi : public AliAnalysisTaskSE {
   //===============================correction\r
   TObjArray* GetAcceptedTracks(AliVEvent* event, Double_t gCentrality, Double_t gReactionPlane);\r
   TObjArray* GetShuffledTracks(TObjArray* tracks, Double_t gCentrality);\r
-  \r
\r
+  TClonesArray* fArrayMC; //! AOD object  //+++++++++++++++++++++\r
   AliBalancePsi *fBalance; //BF object\r
   Bool_t fRunShuffling;//run shuffling or not\r
   AliBalancePsi *fShuffledBalance; //BF object (shuffled)\r
index 44a66b8..8e81385 100644 (file)
@@ -104,6 +104,11 @@ AliAnalysisTaskBFPsi *AddTaskBalancePsiCentralityTrain(Double_t centrMin=0.,
     if(gRunShuffling) bfs = GetBalanceFunctionObject("MCAOD",centralityEstimator,centrMin,centrMax,kTRUE,bResonancesCut,bHBTcut,bConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);\r
     if(gRunMixing)    bfm = GetBalanceFunctionObject("MCAOD",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,bConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);\r
   }\r
+  else if (analysisType=="MCAODrec"){\r
+    bf  = GetBalanceFunctionObject("MCAODrec",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,bConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);\r
+    if(gRunShuffling) bfs = GetBalanceFunctionObject("MCAODrec",centralityEstimator,centrMin,centrMax,kTRUE,bResonancesCut,bHBTcut,bConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);\r
+    if(gRunMixing)    bfm = GetBalanceFunctionObject("MCAODrec",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,bConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning);\r
+  }\r
   else{\r
     ::Error("AddTaskBF", "analysis type NOT known.");\r
     return NULL;\r
@@ -181,6 +186,17 @@ AliAnalysisTaskBFPsi *AddTaskBalancePsiCentralityTrain(Double_t centrMin=0.,
     taskBF->SetAODtrackCutBit(AODfilterBit);\r
     taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax);    \r
   }\r
+  else if(analysisType == "MCAODrec") {     //++++++++++++++++\r
+    // pt and eta cut (pt_min, pt_max, eta_min, eta_max)\r
+    taskBF->SetAODtrackCutBit(AODfilterBit);\r
+    taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax); \r
+\r
+    // set extra DCA cuts (-1 no extra cut)\r
+    taskBF->SetExtraDCACutsAOD(DCAxy,DCAz);\r
+\r
+    // set extra TPC chi2 / nr of clusters cut\r
+    taskBF->SetExtraTPCCutsAOD(maxTPCchi2, minNClustersTPC);\r
+  }//++++++++++++++++\r
 \r
   // offline trigger selection (AliVEvent.h)\r
   // taskBF->UseOfflineTrigger(); // NOT used (selection is done with the AliAnalysisTaskSE::SelectCollisionCandidates()) \r
index ea3343f..7c3d5a8 100644 (file)
@@ -1,5 +1,5 @@
 //__________________________________________________//\r
-AliBalancePsi *GetBalanceFunctionObject(const char* analysisLevel = "AOD", \r
+AliBalancePsi *GetBalanceFunctionObject(const char* analysisLevel = "MCAOD",   //\r
                                        const char* centralityName = 0x0,\r
                                        Double_t centrMin = 0.,\r
                                        Double_t centrMax = 100.,\r
index b07bf64..d219e38 100644 (file)
@@ -133,7 +133,7 @@ void runBalanceFunctionPsi(
     }//local mode\r
 \r
     // input handler (ESD or AOD)\r
-    AliVEventHandler* inputH = NULL;\r
+    Aliveventhandler* inputH = NULL;\r
     if(!bAOD){\r
       inputH = new AliESDInputHandler();\r
     }\r