Adding the possibility to switch off the centrality usage (pp case) + addition of...
authorpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 22 Oct 2011 22:33:55 +0000 (22:33 +0000)
committerpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 22 Oct 2011 22:33:55 +0000 (22:33 +0000)
PWG2/EBYE/AliAnalysisTaskBF.cxx
PWG2/EBYE/AliAnalysisTaskBF.h
PWG2/EBYE/macros/AddTaskBalanceFunctionInpp.C [new file with mode: 0644]
PWG2/EBYE/macros/runBalanceFunctionInpp.C [new file with mode: 0755]

index aa339c9..9bce442 100755 (executable)
@@ -55,6 +55,7 @@ AliAnalysisTaskBF::AliAnalysisTaskBF(const char *name)
   fHistRefTracks(0),\r
   fESDtrackCuts(0),\r
   fCentralityEstimator("V0M"),\r
+  fUseCentrality(kFALSE),\r
   fCentralityPercentileMin(0.), \r
   fCentralityPercentileMax(5.),\r
   fUseOfflineTrigger(kFALSE),\r
@@ -270,97 +271,98 @@ void AliAnalysisTaskBF::UserExec(Option_t *) {
     if(isSelected) {\r
       fHistEventStats->Fill(2); //triggered events\r
 \r
-      //Centrality stuff\r
-      AliCentrality *centrality = gESD->GetCentrality();\r
-      //Int_t nCentrality = 0;\r
-      //nCentrality = centrality->GetCentralityClass5(fCentralityEstimator.Data());\r
-      //cout<<nCentrality<<" "<<centrality->IsEventInCentralityClass(fCentralityPercentileMin,fCentralityPercentileMax,fCentralityEstimator.Data())<<endl;\r
-\r
-      // take only events inside centrality class\r
-      if(centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
-                                             fCentralityPercentileMax,\r
-                                             fCentralityEstimator.Data())){\r
-\r
+      if(fUseCentrality) {\r
+       //Centrality stuff\r
+       AliCentrality *centrality = gESD->GetCentrality();\r
+       //Int_t nCentrality = 0;\r
+       //nCentrality = centrality->GetCentralityClass5(fCentralityEstimator.Data());\r
+       //cout<<nCentrality<<" "<<centrality->IsEventInCentralityClass(fCentralityPercentileMin,fCentralityPercentileMax,fCentralityEstimator.Data())<<endl;\r
+       \r
+       // take only events inside centrality class\r
+       if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
+                                                fCentralityPercentileMax,\r
+                                                fCentralityEstimator.Data()))\r
+         return;\r
+       \r
        // centrality QA (V0M)\r
        fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
+      }\r
        \r
-       const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
-       if(vertex) {\r
-         if(vertex->GetNContributors() > 0) {\r
-           if(vertex->GetZRes() != 0) {\r
-             fHistEventStats->Fill(3); //events with a proper vertex\r
-             if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
-               if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
-                 if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
-                   fHistEventStats->Fill(4); //analayzed events\r
-                   fHistVx->Fill(vertex->GetXv());\r
-                   fHistVy->Fill(vertex->GetYv());\r
-                   fHistVz->Fill(vertex->GetZv());\r
+      const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
+      if(vertex) {\r
+       if(vertex->GetNContributors() > 0) {\r
+         if(vertex->GetZRes() != 0) {\r
+           fHistEventStats->Fill(3); //events with a proper vertex\r
+           if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
+             if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
+               if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
+                 fHistEventStats->Fill(4); //analayzed events\r
+                 fHistVx->Fill(vertex->GetXv());\r
+                 fHistVy->Fill(vertex->GetYv());\r
+                 fHistVz->Fill(vertex->GetZv());\r
+                 \r
+                 //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
+                 for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
+                   AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
+                   if (!track) {\r
+                     Printf("ERROR: Could not receive track %d", iTracks);\r
+                     continue;\r
+                   }   \r
                    \r
-                   //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
-                   for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
-                     AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
-                     if (!track) {\r
-                       Printf("ERROR: Could not receive track %d", iTracks);\r
-                       continue;\r
-                     } \r
-\r
-                     // take only TPC only tracks\r
-                     track_TPC   = new AliESDtrack();\r
-                     if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
-                     \r
-                     //ESD track cuts\r
-                     if(fESDtrackCuts) \r
-                       if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
-                     \r
-                     // fill QA histograms\r
-                     Float_t b[2];\r
-                     Float_t bCov[3];\r
-                     track_TPC->GetImpactParameters(b,bCov);\r
-                     if (bCov[0]<=0 || bCov[2]<=0) {\r
-                       AliDebug(1, "Estimated b resolution lower or equal zero!");\r
-                       bCov[0]=0; bCov[2]=0;\r
-                     }\r
-                     \r
-                     Int_t nClustersTPC = -1;\r
-                     nClustersTPC = track_TPC->GetTPCNclsIter1();   // TPC standalone\r
-                     //nClustersTPC = track->GetTPCclusters(0);   // global track\r
-                     Float_t chi2PerClusterTPC = -1;\r
-                     if (nClustersTPC!=0) {\r
-                       chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
-                       //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
-                     }\r
-                     \r
-                     fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
-                     fHistDCA->Fill(b[1],b[0]);\r
-                     fHistChi2->Fill(chi2PerClusterTPC);\r
-                     fHistPt->Fill(track_TPC->Pt());\r
-                     fHistEta->Fill(track_TPC->Eta());\r
-                     fHistPhi->Fill(track_TPC->Phi()*TMath::RadToDeg());\r
-                     \r
-                     // fill BF array\r
-                     array->Add(track_TPC);\r
-\r
-                     // fill charge vector\r
-                     chargeVector.push_back(track_TPC->Charge());\r
-                     if(fRunShuffling){\r
-                       chargeVectorShuffle.push_back(track_TPC->Charge());\r
-                     }\r
-      \r
-                     delete track_TPC;\r
-\r
-                   } //track loop\r
-                 }//Vz cut\r
-               }//Vy cut\r
-             }//Vx cut\r
-           }//proper vertex resolution\r
-         }//proper number of contributors\r
-       }//vertex object valid\r
-      }//centrality\r
+                   // take only TPC only tracks\r
+                   track_TPC   = new AliESDtrack();\r
+                   if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
+                   \r
+                   //ESD track cuts\r
+                   if(fESDtrackCuts) \r
+                     if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
+                   \r
+                   // fill QA histograms\r
+                   Float_t b[2];\r
+                   Float_t bCov[3];\r
+                   track_TPC->GetImpactParameters(b,bCov);\r
+                   if (bCov[0]<=0 || bCov[2]<=0) {\r
+                     AliDebug(1, "Estimated b resolution lower or equal zero!");\r
+                     bCov[0]=0; bCov[2]=0;\r
+                   }\r
+                   \r
+                   Int_t nClustersTPC = -1;\r
+                   nClustersTPC = track_TPC->GetTPCNclsIter1();   // TPC standalone\r
+                   //nClustersTPC = track->GetTPCclusters(0);   // global track\r
+                   Float_t chi2PerClusterTPC = -1;\r
+                   if (nClustersTPC!=0) {\r
+                     chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
+                     //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
+                   }\r
+                   \r
+                   fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
+                   fHistDCA->Fill(b[1],b[0]);\r
+                   fHistChi2->Fill(chi2PerClusterTPC);\r
+                   fHistPt->Fill(track_TPC->Pt());\r
+                   fHistEta->Fill(track_TPC->Eta());\r
+                   fHistPhi->Fill(track_TPC->Phi()*TMath::RadToDeg());\r
+                   \r
+                   // fill BF array\r
+                   array->Add(track_TPC);\r
+                   \r
+                   // fill charge vector\r
+                   chargeVector.push_back(track_TPC->Charge());\r
+                   if(fRunShuffling){\r
+                     chargeVectorShuffle.push_back(track_TPC->Charge());\r
+                   }\r
+                   \r
+                   delete track_TPC;\r
+                   \r
+                 } //track loop\r
+               }//Vz cut\r
+             }//Vy cut\r
+           }//Vx cut\r
+         }//proper vertex resolution\r
+       }//proper number of contributors\r
+      }//vertex object valid\r
     }//triggered event \r
   }//ESD analysis\r
   \r
-\r
   //AOD analysis (vertex and track cuts also here!!!!)\r
   else if(gAnalysisLevel == "AOD") {\r
     AliAODEvent* gAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE\r
@@ -383,36 +385,38 @@ void AliAnalysisTaskBF::UserExec(Option_t *) {
       fHistEventStats->Fill(2); //triggered events\r
                  \r
       //Centrality stuff (centrality in AOD header)\r
-      Float_t fCentrality     = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
-      // cout<<fCentralityEstimator.Data()<<" = "<<fCentrality<<" ,  others are V0M =  "\r
-      //         << aodHeader->GetCentralityP()->GetCentralityPercentile("V0M")\r
-      //         <<"  FMD = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("FMD")\r
-      //         <<"  TRK = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("TRK")\r
-      //         <<"  TKL = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("TKL")\r
-      //         <<"  CL0 ="<<aodHeader->GetCentralityP()->GetCentralityPercentile("CL0")\r
-      //         <<"  CL1 ="<<aodHeader->GetCentralityP()->GetCentralityPercentile("CL1")\r
-      //         <<"  V0MvsFMD = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD")\r
-      //         <<"  TKLvsV0M = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M")\r
-      //         <<"  ZEMvsZDC = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC")\r
-      //         <<endl;\r
-\r
-      // QA for centrality estimators\r
-      fHistCentStats->Fill(0.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0M"));\r
-      fHistCentStats->Fill(1.,aodHeader->GetCentralityP()->GetCentralityPercentile("FMD"));\r
-      fHistCentStats->Fill(2.,aodHeader->GetCentralityP()->GetCentralityPercentile("TRK"));\r
-      fHistCentStats->Fill(3.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKL"));\r
-      fHistCentStats->Fill(4.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL0"));\r
-      fHistCentStats->Fill(5.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL1"));\r
-      fHistCentStats->Fill(6.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
-      fHistCentStats->Fill(7.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
-      fHistCentStats->Fill(8.,aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
-\r
-      // take only events inside centrality class\r
-      if(fCentrality > fCentralityPercentileMin && fCentrality < fCentralityPercentileMax){\r
-\r
+      if(fUseCentrality) {\r
+       Float_t fCentrality     = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
+       // cout<<fCentralityEstimator.Data()<<" = "<<fCentrality<<" ,  others are V0M =  "\r
+       //        << aodHeader->GetCentralityP()->GetCentralityPercentile("V0M")\r
+       //        <<"  FMD = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("FMD")\r
+       //        <<"  TRK = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("TRK")\r
+       //        <<"  TKL = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("TKL")\r
+       //        <<"  CL0 ="<<aodHeader->GetCentralityP()->GetCentralityPercentile("CL0")\r
+       //        <<"  CL1 ="<<aodHeader->GetCentralityP()->GetCentralityPercentile("CL1")\r
+       //        <<"  V0MvsFMD = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD")\r
+       //        <<"  TKLvsV0M = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M")\r
+       //        <<"  ZEMvsZDC = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC")\r
+       //        <<endl;\r
+       \r
+       // QA for centrality estimators\r
+       fHistCentStats->Fill(0.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0M"));\r
+       fHistCentStats->Fill(1.,aodHeader->GetCentralityP()->GetCentralityPercentile("FMD"));\r
+       fHistCentStats->Fill(2.,aodHeader->GetCentralityP()->GetCentralityPercentile("TRK"));\r
+       fHistCentStats->Fill(3.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKL"));\r
+       fHistCentStats->Fill(4.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL0"));\r
+       fHistCentStats->Fill(5.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL1"));\r
+       fHistCentStats->Fill(6.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
+       fHistCentStats->Fill(7.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
+       fHistCentStats->Fill(8.,aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
+       \r
+       // take only events inside centrality class\r
+       if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) \r
+         return;\r
+       \r
        // centrality QA (V0M)\r
        fHistV0M->Fill(gAOD->GetVZEROData()->GetMTotV0A(), gAOD->GetVZEROData()->GetMTotV0C());\r
-      \r
+       \r
        // centrality QA (reference tracks)\r
        fHistRefTracks->Fill(0.,aodHeader->GetRefMultiplicity());\r
        fHistRefTracks->Fill(1.,aodHeader->GetRefMultiplicityPos());\r
@@ -423,91 +427,90 @@ void AliAnalysisTaskBF::UserExec(Option_t *) {
        fHistRefTracks->Fill(6.,aodHeader->GetNumberOfITSClusters(2));\r
        fHistRefTracks->Fill(7.,aodHeader->GetNumberOfITSClusters(3));\r
        fHistRefTracks->Fill(8.,aodHeader->GetNumberOfITSClusters(4));\r
+      }\r
 \r
-\r
-       const AliAODVertex *vertex = gAOD->GetPrimaryVertex();\r
-\r
-       if(vertex) {\r
-         Double32_t fCov[6];\r
-         vertex->GetCovarianceMatrix(fCov);\r
-                 \r
-         if(vertex->GetNContributors() > 0) {\r
-           if(fCov[5] != 0) {\r
-             fHistEventStats->Fill(3); //events with a proper vertex\r
-             if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
-               if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
-                 if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
-                   fHistEventStats->Fill(4); //analyzed events\r
-                   fHistVx->Fill(vertex->GetX());\r
-                   fHistVy->Fill(vertex->GetY());\r
-                   fHistVz->Fill(vertex->GetZ());\r
-\r
-                   //Printf("There are %d tracks in this event", gAOD->GetNumberOfTracks());\r
-                   for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {\r
-                     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));\r
-                     if (!aodTrack) {\r
-                       Printf("ERROR: Could not receive track %d", iTracks);\r
-                       continue;\r
-                     }\r
-\r
-                     // AOD track cuts\r
-                     \r
-                     // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
-                     // take only TPC only tracks \r
-                     fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
-                     if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
-\r
-                     Float_t pt  = aodTrack->Pt();\r
-                     Float_t eta = aodTrack->Eta();\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
-                     \r
-                     // Kinematics cuts from ESD track cuts\r
-                     if( pt < fPtMin || pt > fPtMax)      continue;\r
-                     if( eta < fEtaMin || eta > fEtaMax)  continue;\r
-\r
-                     // Extra DCA cuts (for systematic studies [!= -1])\r
-                     if( fDCAxyCut != -1 && fDCAxyCut != -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
-\r
-                     // fill QA histograms\r
-                     fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
-                     fHistDCA->Fill(DCAz,DCAxy);\r
-                     fHistChi2->Fill(aodTrack->Chi2perNDF());\r
-                     fHistPt->Fill(pt);\r
-                     fHistEta->Fill(eta);\r
-                     fHistPhi->Fill(aodTrack->Phi()*TMath::RadToDeg());\r
-                     \r
-                     // fill BF array\r
-                     array->Add(aodTrack);\r
-\r
-                     // fill charge vector\r
-                     chargeVector.push_back(aodTrack->Charge());\r
-                     if(fRunShuffling) {\r
-                       chargeVectorShuffle.push_back(aodTrack->Charge());\r
+      const AliAODVertex *vertex = gAOD->GetPrimaryVertex();\r
+      \r
+      if(vertex) {\r
+       Double32_t fCov[6];\r
+       vertex->GetCovarianceMatrix(fCov);\r
+       \r
+       if(vertex->GetNContributors() > 0) {\r
+         if(fCov[5] != 0) {\r
+           fHistEventStats->Fill(3); //events with a proper vertex\r
+           if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
+             if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
+               if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
+                 fHistEventStats->Fill(4); //analyzed events\r
+                 fHistVx->Fill(vertex->GetX());\r
+                 fHistVy->Fill(vertex->GetY());\r
+                 fHistVz->Fill(vertex->GetZ());\r
+                 \r
+                 //Printf("There are %d tracks in this event", gAOD->GetNumberOfTracks());\r
+                 for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {\r
+                   AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));\r
+                   if (!aodTrack) {\r
+                     Printf("ERROR: Could not receive track %d", iTracks);\r
+                     continue;\r
+                   }\r
+                   \r
+                   // AOD track cuts\r
+                   \r
+                   // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
+                   // take only TPC only tracks \r
+                   fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
+                   if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
+                   \r
+                   Float_t pt  = aodTrack->Pt();\r
+                   Float_t eta = aodTrack->Eta();\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
+                   \r
+                   // Kinematics cuts from ESD track cuts\r
+                   if( pt < fPtMin || pt > fPtMax)      continue;\r
+                   if( eta < fEtaMin || eta > fEtaMax)  continue;\r
+                   \r
+                   // Extra DCA cuts (for systematic studies [!= -1])\r
+                   if( fDCAxyCut != -1 && fDCAxyCut != -1){\r
+                     if(TMath::Sqrt((DCAxy*DCAxy)/(fDCAxyCut*fDCAxyCut)+(DCAz*DCAz)/(fDCAzCut*fDCAzCut)) > 1 ){\r
+                       continue;  // 2D cut\r
                      }\r
-                   } //track loop\r
-                 }//Vz cut\r
-               }//Vy cut\r
-             }//Vx cut\r
-           }//proper vertex resolution\r
-         }//proper number of contributors\r
-       }//vertex object valid\r
-      }//centrality\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
+                   \r
+                   // fill QA histograms\r
+                   fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
+                   fHistDCA->Fill(DCAz,DCAxy);\r
+                   fHistChi2->Fill(aodTrack->Chi2perNDF());\r
+                   fHistPt->Fill(pt);\r
+                   fHistEta->Fill(eta);\r
+                   fHistPhi->Fill(aodTrack->Phi()*TMath::RadToDeg());\r
+                   \r
+                   // fill BF array\r
+                   array->Add(aodTrack);\r
+                   \r
+                   // fill charge vector\r
+                   chargeVector.push_back(aodTrack->Charge());\r
+                   if(fRunShuffling) {\r
+                     chargeVectorShuffle.push_back(aodTrack->Charge());\r
+                   }\r
+                 } //track loop\r
+               }//Vz cut\r
+             }//Vy cut\r
+           }//Vx cut\r
+         }//proper vertex resolution\r
+       }//proper number of contributors\r
+      }//vertex object valid\r
     }//triggered event \r
   }//AOD analysis\r
 \r
index 3271e87..df68359 100755 (executable)
@@ -68,6 +68,7 @@ class AliAnalysisTaskBF : public AliAnalysisTaskSE {
   void SetCentralityEstimator(const char* centralityEstimator) {fCentralityEstimator = centralityEstimator;}\r
   const char* GetCentralityEstimator(void)                     {return fCentralityEstimator;}\r
   void SetCentralityPercentileRange(Double_t min, Double_t max) { \r
+    fUseCentrality = kTRUE;\r
     fCentralityPercentileMin=min;\r
     fCentralityPercentileMax=max;\r
   }\r
@@ -102,6 +103,7 @@ class AliAnalysisTaskBF : public AliAnalysisTaskSE {
   AliESDtrackCuts *fESDtrackCuts; //ESD track cuts\r
 \r
   TString fCentralityEstimator;      //"V0M","TRK","TKL","ZDC","FMD"\r
+  Bool_t fUseCentrality;//use the centrality (PbPb) or not (pp)\r
   Double_t fCentralityPercentileMin;\r
   Double_t fCentralityPercentileMax;\r
 \r
diff --git a/PWG2/EBYE/macros/AddTaskBalanceFunctionInpp.C b/PWG2/EBYE/macros/AddTaskBalanceFunctionInpp.C
new file mode 100644 (file)
index 0000000..f14eef9
--- /dev/null
@@ -0,0 +1,85 @@
+//_________________________________________________________//\r
+AliAnalysisTaskBF *AddTaskBalanceFunctionInpp(Double_t vertexZ=10.,\r
+                                             Double_t DCAxy=2.4,\r
+                                             Double_t DCAz=3.2,\r
+                                             Double_t ptMin=0.3,\r
+                                             Double_t ptMax=1.5,\r
+                                             Double_t etaMin=-0.8,\r
+                                             Double_t etaMax=0.8,\r
+                                             TString fileNameBase="AnalysisResults") {\r
+  // Creates a balance function analysis task and adds it to the analysis manager.\r
+  // Get the pointer to the existing analysis manager via the static access method.\r
+  TString outputFileName(fileNameBase);\r
+  outputFileName.Append(".root");\r
+\r
+  //===========================================================================\r
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
+  if (!mgr) {\r
+    ::Error("AddTaskBF", "No analysis manager to connect to.");\r
+    return NULL;\r
+  }\r
+\r
+  // Check the analysis type using the event handlers connected to the analysis manager.\r
+  //===========================================================================\r
+  if (!mgr->GetInputEventHandler()) {\r
+    ::Error("AddTaskBF", "This task requires an input event handler");\r
+    return NULL;\r
+  }\r
+  TString analysisType = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"\r
+  gROOT->LoadMacro("./configBalanceFunctionAnalysis.C");\r
+  AliBalance *bf  = 0;  // Balance Function object\r
+  AliBalance *bfs = 0;  // shuffled Balance function object\r
+\r
+  if (analysisType=="ESD"){\r
+    bf  = GetBalanceFunctionObject("ESD");\r
+    bfs = GetBalanceFunctionObject("ESD",kTRUE);\r
+  }\r
+  else if (analysisType=="AOD"){\r
+    bf  = GetBalanceFunctionObject("AOD");\r
+    bfs = GetBalanceFunctionObject("AOD",kTRUE);\r
+  }\r
+  else{\r
+    bf  = GetBalanceFunctionObject("MC");\r
+    bfs = GetBalanceFunctionObject("MC",kTRUE);\r
+  }\r
+\r
+  // Create the task, add it to manager and configure it.\r
+  //===========================================================================\r
+  AliAnalysisTaskBF *taskBF = new AliAnalysisTaskBF("TaskBF");\r
+  taskBF->SetAnalysisObject(bf);\r
+  taskBF->SetShufflingObject(bfs);\r
+  if(analysisType == "ESD") {\r
+    AliESDtrackCuts *trackCuts = GetTrackCutsObject();\r
+    taskBF->SetAnalysisCutObject(trackCuts);\r
+    \r
+    // offline trigger selection (AliVEvent.h)\r
+    // taskBF->UseOfflineTrigger(); // NOT used (selection is done with the AliAnalysisTaskSE::SelectCollisionCandidates()) \r
+    // with this only selected events are analyzed (first 2 bins in event QA histogram are the same))\r
+    taskBF->SelectCollisionCandidates(AliVEvent::kMB);\r
+  }\r
+  else if(analysisType == "AOD") {\r
+    // pt and eta cut (pt_min, pt_max, eta_min, eta_max)\r
+    taskBF->SetAODtrackCutBit(128);\r
+    taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax);\r
+    taskBF->SetExtraDCACutsAOD(DCAxy,DCAz);\r
+  }\r
+  \r
+  // vertex cut (x,y,z)\r
+  taskBF->SetVertexDiamond(.3,.3,vertexZ);\r
+  \r
+  //bf->PrintAnalysisSettings();\r
+  mgr->AddTask(taskBF);\r
+\r
+  // Create ONLY the output containers for the data produced by the task.\r
+  // Get and connect other common input/output containers via the manager as below\r
+  //======================================================================\r
+  AliAnalysisDataContainer *coutQA = mgr->CreateContainer("listQA", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  AliAnalysisDataContainer *coutBF = mgr->CreateContainer("listBF", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  AliAnalysisDataContainer *coutBFS= mgr->CreateContainer("listBFshuffled", TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  mgr->ConnectInput(taskBF, 0, mgr->GetCommonInputContainer());\r
+  mgr->ConnectOutput(taskBF, 1, coutQA);\r
+  mgr->ConnectOutput(taskBF, 2, coutBF);\r
+  mgr->ConnectOutput(taskBF, 3, coutBFS);\r
+\r
+  return taskBF;\r
+}\r
diff --git a/PWG2/EBYE/macros/runBalanceFunctionInpp.C b/PWG2/EBYE/macros/runBalanceFunctionInpp.C
new file mode 100755 (executable)
index 0000000..343bae6
--- /dev/null
@@ -0,0 +1,223 @@
+enum analysisModes {mLocal,mLocalPAR,mPROOF,mGrid,mGridPAR};\r
+enum analysisTypes {mESD,mAOD,mMC,mMCESD};\r
+\r
+//\r
+class AliAnalysisGrid;\r
+class AliAnalysisTaskBF;\r
+class AliBalance;\r
+\r
+//________________________________________________________________________//\r
+void runBalanceFunctionInpp(Int_t mode = mLocal, \r
+                           Int_t type = mAOD,\r
+                           Bool_t DATA = kFALSE) {\r
+  // Time:\r
+  TStopwatch timer;\r
+  timer.Start();\r
+  \r
+  //Check analysis mode\r
+  if((mode < 0) || (mode > 4)) {\r
+    Printf("Analysis mode not recognized!");\r
+    Printf("You can select out of 0: local, 1: local with par files, 2: proof, 3: grid, 4: grid with par files");\r
+    return;\r
+  }\r
+  \r
+  //Check analysis type\r
+  if((type < 0) || (type > 3)) {\r
+    Printf("Analysis type not recognized!");\r
+    Printf("You can select out of 0: ESD, 1: AOD, 2: MC (stack), 3: MC (from the ESD)");\r
+    return;\r
+  }\r
+  \r
+  // Load needed libraries:\r
+  LoadLibraries(mode);\r
+  \r
+ // Create and configure the AliEn plug-in:\r
+  if(mode == mGrid || mode == mGridPAR) {\r
+    gROOT->LoadMacro("CreateAlienHandler.C");\r
+    AliAnalysisGrid *alienHandler = CreateAlienHandler(runListFileName);  \r
+    if (!alienHandler) return;\r
+  }\r
+  // Chains:   \r
+  if(mode == mLocal || mode == mLocalPAR) {\r
+    TChain* chain = 0x0;\r
+    if((type == mESD)||(type == mMCESD))  \r
+      chain = new TChain("esdTree");\r
+    else if(type == mAOD)\r
+      chain = new TChain("aodTree");\r
+    else if(type == mMC)\r
+      chain = new TChain("TE");\r
+\r
+    TString filename;\r
+    for(Int_t i = 1; i < 20; i++) {\r
+      filename = "/data/alice2/pchrist/pp/LHC10c/0.9TeV/Data/";\r
+      filename += "Set"; filename += i; \r
+      if((type == mESD)||(type == mMCESD))  \r
+       filename += "/AliESDs.root";\r
+      else if(type == mAOD)\r
+       filename += "/AliAOD.root";\r
+      else if(type == mMC)\r
+       filename += "/galice.root";\r
+\r
+      chain->Add(filename.Data());\r
+    }\r
+  }\r
+  //Proof\r
+  if(mode == mPROOF) {\r
+    gROOT->ProcessLine(Form(".include %s/include", gSystem->ExpandPathName("$ALICE_ROOT")));\r
+  }\r
+  \r
+  // analysis manager\r
+  AliAnalysisManager* mgr = new AliAnalysisManager("balanceFunctionManager");\r
+  if(mode == mGrid || mode == mGridPAR)\r
+    mgr->SetGridHandler(alienHandler);\r
+    \r
+  // input handler (ESD or AOD)\r
+  AliVEventHandler* inputH = NULL;\r
+  if((type == mESD)||(type == mMCESD))  \r
+    inputH = new AliESDInputHandler();\r
+  else if(type == mAOD)\r
+    inputH = new AliAODInputHandler();\r
+  mgr->SetInputEventHandler(inputH);\r
+    \r
+  // mc event handler\r
+  if((type == mMC) || (type == mMCESD)) {\r
+    AliMCEventHandler* mchandler = new AliMCEventHandler();\r
+    // Not reading track references\r
+    mchandler->SetReadTR(kFALSE);\r
+    mgr->SetMCtruthEventHandler(mchandler);\r
+  }   \r
+\r
+  if(type != mAOD){\r
+    gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");\r
+    AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection();\r
+  }\r
+\r
+  //Add the BF task (all centralities)\r
+  gROOT->LoadMacro("AddTaskBalanceFunctionInpp.C"); \r
+  AddTaskBalanceFunctionInpp();\r
+\r
+  // enable debug printouts\r
+  mgr->SetDebugLevel(2);\r
+  mgr->SetUseProgressBar(1,100);\r
+  if (!mgr->InitAnalysis()) return;\r
+  mgr->PrintStatus();\r
+  \r
+  // start analysis\r
+  if(mode == mLocal || mode == mLocalPAR) \r
+    mgr->StartAnalysis("local",chain);\r
+  else if(mode == mPROOF) \r
+    mgr->StartAnalysis("proof",dataDir,nRuns,offset);\r
+  else if(mode == mGrid || mode == mGridPAR) \r
+    mgr->StartAnalysis("grid");\r
+\r
+  // Print real and CPU time used for analysis:  \r
+  timer.Stop();\r
+  timer.Print();\r
+}\r
+\r
+//=============================================================//\r
+void LoadLibraries(const analysisModes mode) {  \r
+  //--------------------------------------\r
+  // Load the needed libraries most of them already loaded by aliroot\r
+  //--------------------------------------\r
+  gSystem->Load("libCore.so");        \r
+  gSystem->Load("libGeom.so");\r
+  gSystem->Load("libVMC.so");\r
+  gSystem->Load("libPhysics.so");\r
+  gSystem->Load("libTree.so");\r
+\r
+  //----------------------------------------------------------\r
+  // >>>>>>>>>>> Local mode <<<<<<<<<<<<<< \r
+  //----------------------------------------------------------\r
+  if (mode==mLocal || mode==mGrid || mode == mGridPAR) {\r
+    //--------------------------------------------------------\r
+    // If you want to use already compiled libraries \r
+    // in the aliroot distribution\r
+    //--------------------------------------------------------\r
+    gSystem->Load("libSTEERBase.so");\r
+    gSystem->Load("libESD.so");\r
+    gSystem->Load("libAOD.so");\r
+    gSystem->Load("libANALYSIS.so");\r
+    gSystem->Load("libANALYSISalice.so");\r
+    gSystem->Load("libPWG2ebye.so");\r
+    // Use AliRoot includes to compile our task\r
+    gROOT->ProcessLine(".include $ALICE_ROOT/include");\r
+  }\r
+  \r
+  else if (mode == mLocalPAR) {\r
+    //--------------------------------------------------------\r
+    //If you want to use root and par files from aliroot\r
+    //--------------------------------------------------------  \r
+    SetupPar("STEERBase");\r
+    SetupPar("ESD");\r
+    SetupPar("AOD");\r
+    SetupPar("ANALYSIS");\r
+    SetupPar("ANALYSISalice");\r
+    SetupPar("PWG2ebye");\r
+}\r
+  \r
+  //---------------------------------------------------------\r
+  // <<<<<<<<<< PROOF mode >>>>>>>>>>>>\r
+  //---------------------------------------------------------\r
+  else if (mode==mPROOF) {\r
+    // Connect to proof\r
+    printf("*** Connect to PROOF ***\n");\r
+    gEnv->SetValue("XSec.GSI.DelegProxy","2");\r
+    // Put appropriate username here\r
+    TProof::Open("alice-caf.cern.ch");\r
+    //TProof::Open("skaf.saske.sk");\r
+    //TProof::Open("prf000-iep-grid.saske.sk");\r
+\r
+    gProof->EnablePackage("VO_ALICE@AliRoot::v4-21-12-AN");\r
+  }  \r
+  \r
+} // end of void LoadLibraries(const anaModes mode)\r
+\r
+//======================================================================//\r
+void SetupPar(char* pararchivename) {\r
+  //Load par files, create analysis libraries\r
+  //For testing, if par file already decompressed and modified\r
+  //classes then do not decompress.\r
+  \r
+  TString cdir(Form("%s", gSystem->WorkingDirectory() )) ; \r
+  TString parpar(Form("%s.par", pararchivename)) ; \r
+  if ( gSystem->AccessPathName(parpar.Data()) ) {\r
+    gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;\r
+    TString processline(Form(".! make %s", parpar.Data())) ; \r
+    gROOT->ProcessLine(processline.Data()) ;\r
+    gSystem->ChangeDirectory(cdir) ; \r
+    processline = Form(".! mv /tmp/%s .", parpar.Data()) ;\r
+    gROOT->ProcessLine(processline.Data()) ;\r
+  } \r
+  if ( gSystem->AccessPathName(pararchivename) ) {  \r
+    TString processline = Form(".! tar xvzf %s",parpar.Data()) ;\r
+    gROOT->ProcessLine(processline.Data());\r
+  }\r
+  \r
+  TString ocwd = gSystem->WorkingDirectory();\r
+  gSystem->ChangeDirectory(pararchivename);\r
+  \r
+  // check for BUILD.sh and execute\r
+  if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {\r
+    printf("*******************************\n");\r
+    printf("*** Building PAR archive    ***\n");\r
+    cout<<pararchivename<<endl;\r
+    printf("*******************************\n");\r
+    if (gSystem->Exec("PROOF-INF/BUILD.sh")) {\r
+      Error("runProcess","Cannot Build the PAR Archive! - Abort!");\r
+      return -1;\r
+    }\r
+  }\r
+  // check for SETUP.C and execute\r
+  if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {\r
+    printf("*******************************\n");\r
+    printf("*** Setup PAR archive       ***\n");\r
+    cout<<pararchivename<<endl;\r
+    printf("*******************************\n");\r
+    gROOT->Macro("PROOF-INF/SETUP.C");\r
+  }\r
+  \r
+  gSystem->ChangeDirectory(ocwd.Data());\r
+  printf("Current dir: %s\n", ocwd.Data());\r
+\r
+} // end of void SetupPar(char* pararchivename) \r