PWG2/EBYE/AliAnalysisTaskBF.cxx

   1 #include "TChain.h"\r
   2 #include "TList.h"\r
   3 #include "TCanvas.h"\r
   4 #include "TLorentzVector.h"\r
   5 #include "TGraphErrors.h"\r
   6 #include "TH1F.h"\r
   7 #include "TH2F.h"\r
   8 #include "TArrayF.h"\r
   9 \r
  10 #include "AliAnalysisTaskSE.h"\r
  11 #include "AliAnalysisManager.h"\r
  12 \r
  13 #include "AliESDVertex.h"\r
  14 #include "AliESDEvent.h"\r
  15 #include "AliESDInputHandler.h"\r
  16 #include "AliAODEvent.h"\r
  17 #include "AliAODTrack.h"\r
  18 #include "AliAODInputHandler.h"\r
  19 #include "AliGenEventHeader.h"\r
  20 #include "AliGenHijingEventHeader.h"\r
  21 #include "AliMCEventHandler.h"\r
  22 #include "AliMCEvent.h"\r
  23 #include "AliStack.h"\r
  24 #include "AliESDtrackCuts.h"\r
  25 \r
  26 #include "AliAnalysisTaskBF.h"\r
  27 #include "AliBalance.h"\r
  28 \r
  29 \r
  30 // Analysis task for the BF code\r
  31 // Authors: Panos.Christakoglou@nikhef.nl\r
  32 \r
  33 ClassImp(AliAnalysisTaskBF)\r
  34 \r
  35 //________________________________________________________________________\r
  36 AliAnalysisTaskBF::AliAnalysisTaskBF(const char *name) \r
  37 : AliAnalysisTaskSE(name), \r
  38   fBalance(0),\r
  39   fRunShuffling(kFALSE),\r
  40   fShuffledBalance(0),\r
  41   fList(0),\r
  42   fListBF(0),\r
  43   fListBFS(0),\r
  44   fHistEventStats(0),\r
  45   fHistCentStats(0),\r
  46   fHistTriggerStats(0),\r
  47   fHistTrackStats(0),\r
  48   fHistVx(0),\r
  49   fHistVy(0),\r
  50   fHistVz(0),\r
  51   fHistClus(0),\r
  52   fHistDCA(0),\r
  53   fHistChi2(0),\r
  54   fHistPt(0),\r
  55   fHistEta(0),\r
  56   fHistPhi(0),\r
  57   fHistV0M(0),\r
  58   fHistRefTracks(0),\r
  59   fESDtrackCuts(0),\r
  60   fCentralityEstimator("V0M"),\r
  61   fUseCentrality(kFALSE),\r
  62   fCentralityPercentileMin(0.), \r
  63   fCentralityPercentileMax(5.),\r
  64   fImpactParameterMin(0.),\r
  65   fImpactParameterMax(20.),\r
  66   fUseMultiplicity(kFALSE),\r
  67   fNumberOfAcceptedTracksMin(0),\r
  68   fNumberOfAcceptedTracksMax(10000),\r
  69   fUseOfflineTrigger(kFALSE),\r
  70   fVxMax(0.3),\r
  71   fVyMax(0.3),\r
  72   fVzMax(10.),\r
  73   nAODtrackCutBit(128),\r
  74   fPtMin(0.3),\r
  75   fPtMax(1.5),\r
  76   fEtaMin(-0.8),\r
  77   fEtaMax(-0.8),\r
  78   fDCAxyCut(-1),\r
  79   fDCAzCut(-1),\r
  80   fTPCchi2Cut(-1),\r
  81   fNClustersTPCCut(-1){\r
  82   // Constructor\r
  83 \r
  84   // Define input and output slots here\r
  85   // Input slot #0 works with a TChain\r
  86   DefineInput(0, TChain::Class());\r
  87   // Output slot #0 writes into a TH1 container\r
  88   DefineOutput(1, TList::Class());\r
  89   DefineOutput(2, TList::Class());\r
  90   DefineOutput(3, TList::Class());\r
  91 }\r
  92 \r
  93 //________________________________________________________________________\r
  94 AliAnalysisTaskBF::~AliAnalysisTaskBF() {\r
  95 \r
  96   // delete fBalance; \r
  97   // delete fShuffledBalance; \r
  98   // delete fList;\r
  99   // delete fListBF; \r
 100   // delete fListBFS;\r
 101 \r
 102   // delete fHistEventStats; \r
 103   // delete fHistTrackStats; \r
 104   // delete fHistVx; \r
 105   // delete fHistVy; \r
 106   // delete fHistVz; \r
 107 \r
 108   // delete fHistClus;\r
 109   // delete fHistDCA;\r
 110   // delete fHistChi2;\r
 111   // delete fHistPt;\r
 112   // delete fHistEta;\r
 113   // delete fHistPhi;\r
 114   // delete fHistV0M;\r
 115 }\r
 116 \r
 117 //________________________________________________________________________\r
 118 void AliAnalysisTaskBF::UserCreateOutputObjects() {\r
 119   // Create histograms\r
 120   // Called once\r
 121   if(!fBalance) {\r
 122     fBalance = new AliBalance();\r
 123     fBalance->SetAnalysisLevel("ESD");\r
 124     //fBalance->SetNumberOfBins(-1,16);\r
 125     fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
 126   }\r
 127   if(fRunShuffling) {\r
 128     if(!fShuffledBalance) {\r
 129       fShuffledBalance = new AliBalance();\r
 130       fShuffledBalance->SetAnalysisLevel("ESD");\r
 131       //fShuffledBalance->SetNumberOfBins(-1,16);\r
 132       fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6);\r
 133     }\r
 134   }\r
 135 \r
 136   //QA list\r
 137   fList = new TList();\r
 138   fList->SetName("listQA");\r
 139   fList->SetOwner();\r
 140 \r
 141   //Balance Function list\r
 142   fListBF = new TList();\r
 143   fListBF->SetName("listBF");\r
 144   fListBF->SetOwner();\r
 145 \r
 146   if(fRunShuffling) {\r
 147     fListBFS = new TList();\r
 148     fListBFS->SetName("listBFShuffled");\r
 149     fListBFS->SetOwner();\r
 150   }\r
 151 \r
 152   //Event stats.\r
 153   TString gCutName[4] = {"Total","Offline trigger",\r
 154                          "Vertex","Analyzed"};\r
 155   fHistEventStats = new TH1F("fHistEventStats",\r
 156                              "Event statistics;;N_{events}",\r
 157                              4,0.5,4.5);\r
 158   for(Int_t i = 1; i <= 4; i++)\r
 159     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
 160   fList->Add(fHistEventStats);\r
 161 \r
 162   TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
 163   fHistCentStats = new TH2F("fHistCentStats",\r
 164                              "Centrality statistics;;Cent percentile",\r
 165                             9,-0.5,8.5,220,-5,105);\r
 166   for(Int_t i = 1; i <= 9; i++)\r
 167     fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
 168   fList->Add(fHistCentStats);\r
 169 \r
 170   fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);\r
 171   fList->Add(fHistTriggerStats);\r
 172 \r
 173   fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",130,0,130);\r
 174   fList->Add(fHistTrackStats);\r
 175 \r
 176   // Vertex distributions\r
 177   fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r
 178   fList->Add(fHistVx);\r
 179   fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r
 180   fList->Add(fHistVy);\r
 181   fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);\r
 182   fList->Add(fHistVz);\r
 183 \r
 184   // QA histograms\r
 185   fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r
 186   fList->Add(fHistClus);\r
 187   fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);\r
 188   fList->Add(fHistChi2);\r
 189   fHistDCA  = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r
 190   fList->Add(fHistDCA);\r
 191   fHistPt   = new TH1F("fHistPt","p_{T} distribution",200,0,10);\r
 192   fList->Add(fHistPt);\r
 193   fHistEta  = new TH1F("fHistEta","#eta distribution",200,-2,2);\r
 194   fList->Add(fHistEta);\r
 195   fHistPhi  = new TH1F("fHistPhi","#phi distribution",200,-20,380);\r
 196   fList->Add(fHistPhi);\r
 197   fHistV0M  = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r
 198   fList->Add(fHistV0M);\r
 199   TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r
 200   fHistRefTracks  = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r
 201   for(Int_t i = 1; i <= 6; i++)\r
 202     fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
 203   fList->Add(fHistRefTracks);\r
 204 \r
 205 \r
 206   // Balance function histograms\r
 207 \r
 208   // Initialize histograms if not done yet\r
 209   if(!fBalance->GetHistNp(0)){\r
 210     AliWarning("Histograms not yet initialized! --> Will be done now");\r
 211     AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
 212     fBalance->InitHistograms();\r
 213   }\r
 214 \r
 215   if(fRunShuffling) {\r
 216     if(!fShuffledBalance->GetHistNp(0)) {\r
 217       AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
 218       AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
 219       fShuffledBalance->InitHistograms();\r
 220     }\r
 221   }\r
 222 \r
 223   for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r
 224     fListBF->Add(fBalance->GetHistNp(a));\r
 225     fListBF->Add(fBalance->GetHistNn(a));\r
 226     fListBF->Add(fBalance->GetHistNpn(a));\r
 227     fListBF->Add(fBalance->GetHistNnn(a));\r
 228     fListBF->Add(fBalance->GetHistNpp(a));\r
 229     fListBF->Add(fBalance->GetHistNnp(a));\r
 230 \r
 231     if(fRunShuffling) {\r
 232       fListBFS->Add(fShuffledBalance->GetHistNp(a));\r
 233       fListBFS->Add(fShuffledBalance->GetHistNn(a));\r
 234       fListBFS->Add(fShuffledBalance->GetHistNpn(a));\r
 235       fListBFS->Add(fShuffledBalance->GetHistNnn(a));\r
 236       fListBFS->Add(fShuffledBalance->GetHistNpp(a));\r
 237       fListBFS->Add(fShuffledBalance->GetHistNnp(a));\r
 238     }  \r
 239   }\r
 240 \r
 241   if(fESDtrackCuts) fList->Add(fESDtrackCuts);\r
 242 \r
 243   // Post output data.\r
 244   PostData(1, fList);\r
 245   PostData(2, fListBF);\r
 246   if(fRunShuffling) PostData(3, fListBFS);\r
 247 }\r
 248 \r
 249 //________________________________________________________________________\r
 250 void AliAnalysisTaskBF::UserExec(Option_t *) {\r
 251   // Main loop\r
 252   // Called for each event\r
 253   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 254 \r
 255   TObjArray *array         = new TObjArray();\r
 256   AliESDtrack *track_TPC   = NULL;\r
 257 \r
 258   Int_t gNumberOfAcceptedTracks = 0;\r
 259 \r
 260   // vector holding the charges of all tracks\r
 261   vector<Int_t> chargeVectorShuffle;   // this will be shuffled\r
 262   vector<Int_t> chargeVector;          // original charge \r
 263   \r
 264   //ESD analysis\r
 265   if(gAnalysisLevel == "ESD") {\r
 266     AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
 267     if (!gESD) {\r
 268       Printf("ERROR: gESD not available");\r
 269       return;\r
 270     }\r
 271 \r
 272     // store offline trigger bits\r
 273     fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
 274 \r
 275     // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
 276     fHistEventStats->Fill(1); //all events\r
 277     Bool_t isSelected = kTRUE;\r
 278     if(fUseOfflineTrigger)\r
 279       isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
 280     if(isSelected) {\r
 281       fHistEventStats->Fill(2); //triggered events\r
 282 \r
 283       if(fUseCentrality) {\r
 284         //Centrality stuff\r
 285         AliCentrality *centrality = gESD->GetCentrality();\r
 286         //Int_t nCentrality = 0;\r
 287         //nCentrality = centrality->GetCentralityClass5(fCentralityEstimator.Data());\r
 288         //cout<<nCentrality<<" "<<centrality->IsEventInCentralityClass(fCentralityPercentileMin,fCentralityPercentileMax,fCentralityEstimator.Data())<<endl;\r
 289         \r
 290         // take only events inside centrality class\r
 291         if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
 292                                                  fCentralityPercentileMax,\r
 293                                                  fCentralityEstimator.Data()))\r
 294           return;\r
 295         \r
 296         // centrality QA (V0M)\r
 297         fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
 298       }\r
 299         \r
 300       const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
 301       if(vertex) {\r
 302         if(vertex->GetNContributors() > 0) {\r
 303           if(vertex->GetZRes() != 0) {\r
 304             fHistEventStats->Fill(3); //events with a proper vertex\r
 305             if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
 306               if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
 307                 if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
 308                   fHistEventStats->Fill(4); //analayzed events\r
 309                   fHistVx->Fill(vertex->GetXv());\r
 310                   fHistVy->Fill(vertex->GetYv());\r
 311                   fHistVz->Fill(vertex->GetZv());\r
 312                   \r
 313                   //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
 314                   for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
 315                     AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
 316                     if (!track) {\r
 317                       Printf("ERROR: Could not receive track %d", iTracks);\r
 318                       continue;\r
 319                     }   \r
 320                     \r
 321                     // take only TPC only tracks\r
 322                     track_TPC   = new AliESDtrack();\r
 323                     if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
 324                     \r
 325                     //ESD track cuts\r
 326                     if(fESDtrackCuts) \r
 327                       if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
 328                     \r
 329                     // fill QA histograms\r
 330                     Float_t b[2];\r
 331                     Float_t bCov[3];\r
 332                     track_TPC->GetImpactParameters(b,bCov);\r
 333                     if (bCov[0]<=0 || bCov[2]<=0) {\r
 334                       AliDebug(1, "Estimated b resolution lower or equal zero!");\r
 335                       bCov[0]=0; bCov[2]=0;\r
 336                     }\r
 337                     \r
 338                     Int_t nClustersTPC = -1;\r
 339                     nClustersTPC = track_TPC->GetTPCNclsIter1();   // TPC standalone\r
 340                     //nClustersTPC = track->GetTPCclusters(0);   // global track\r
 341                     Float_t chi2PerClusterTPC = -1;\r
 342                     if (nClustersTPC!=0) {\r
 343                       chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
 344                       //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
 345                     }\r
 346                     \r
 347                     fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
 348                     fHistDCA->Fill(b[1],b[0]);\r
 349                     fHistChi2->Fill(chi2PerClusterTPC);\r
 350                     fHistPt->Fill(track_TPC->Pt());\r
 351                     fHistEta->Fill(track_TPC->Eta());\r
 352                     fHistPhi->Fill(track_TPC->Phi()*TMath::RadToDeg());\r
 353                     \r
 354                     // fill BF array\r
 355                     array->Add(track_TPC);\r
 356                     \r
 357                     // fill charge vector\r
 358                     chargeVector.push_back(track_TPC->Charge());\r
 359                     if(fRunShuffling){\r
 360                       chargeVectorShuffle.push_back(track_TPC->Charge());\r
 361                     }\r
 362                     \r
 363                     delete track_TPC;\r
 364                     \r
 365                   } //track loop\r
 366                 }//Vz cut\r
 367               }//Vy cut\r
 368             }//Vx cut\r
 369           }//proper vertex resolution\r
 370         }//proper number of contributors\r
 371       }//vertex object valid\r
 372     }//triggered event \r
 373   }//ESD analysis\r
 374   \r
 375   //AOD analysis (vertex and track cuts also here!!!!)\r
 376   else if(gAnalysisLevel == "AOD") {\r
 377     AliAODEvent* gAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE\r
 378     if(!gAOD) {\r
 379       Printf("ERROR: gAOD not available");\r
 380       return;\r
 381     }\r
 382 \r
 383     AliAODHeader *aodHeader = gAOD->GetHeader();\r
 384 \r
 385     // store offline trigger bits\r
 386     fHistTriggerStats->Fill(aodHeader->GetOfflineTrigger());\r
 387     \r
 388     // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
 389     fHistEventStats->Fill(1); //all events\r
 390     Bool_t isSelected = kTRUE;\r
 391     if(fUseOfflineTrigger)\r
 392       isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
 393     if(isSelected) {\r
 394       fHistEventStats->Fill(2); //triggered events\r
 395                   \r
 396       //Centrality stuff (centrality in AOD header)\r
 397       if(fUseCentrality) {\r
 398         Float_t fCentrality     = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
 399         // cout<<fCentralityEstimator.Data()<<" = "<<fCentrality<<" ,  others are V0M =  "\r
 400         //        << aodHeader->GetCentralityP()->GetCentralityPercentile("V0M")\r
 401         //        <<"  FMD = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("FMD")\r
 402         //        <<"  TRK = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("TRK")\r
 403         //        <<"  TKL = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("TKL")\r
 404         //        <<"  CL0 ="<<aodHeader->GetCentralityP()->GetCentralityPercentile("CL0")\r
 405         //        <<"  CL1 ="<<aodHeader->GetCentralityP()->GetCentralityPercentile("CL1")\r
 406         //        <<"  V0MvsFMD = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD")\r
 407         //        <<"  TKLvsV0M = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M")\r
 408         //        <<"  ZEMvsZDC = "<<aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC")\r
 409         //        <<endl;\r
 410         \r
 411         // QA for centrality estimators\r
 412         fHistCentStats->Fill(0.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0M"));\r
 413         fHistCentStats->Fill(1.,aodHeader->GetCentralityP()->GetCentralityPercentile("FMD"));\r
 414         fHistCentStats->Fill(2.,aodHeader->GetCentralityP()->GetCentralityPercentile("TRK"));\r
 415         fHistCentStats->Fill(3.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKL"));\r
 416         fHistCentStats->Fill(4.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL0"));\r
 417         fHistCentStats->Fill(5.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL1"));\r
 418         fHistCentStats->Fill(6.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
 419         fHistCentStats->Fill(7.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
 420         fHistCentStats->Fill(8.,aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
 421         \r
 422         // take only events inside centrality class\r
 423         if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax)) \r
 424           return;\r
 425         \r
 426         // centrality QA (V0M)\r
 427         fHistV0M->Fill(gAOD->GetVZEROData()->GetMTotV0A(), gAOD->GetVZEROData()->GetMTotV0C());\r
 428         \r
 429         // centrality QA (reference tracks)\r
 430         fHistRefTracks->Fill(0.,aodHeader->GetRefMultiplicity());\r
 431         fHistRefTracks->Fill(1.,aodHeader->GetRefMultiplicityPos());\r
 432         fHistRefTracks->Fill(2.,aodHeader->GetRefMultiplicityNeg());\r
 433         fHistRefTracks->Fill(3.,aodHeader->GetTPConlyRefMultiplicity());\r
 434         fHistRefTracks->Fill(4.,aodHeader->GetNumberOfITSClusters(0));\r
 435         fHistRefTracks->Fill(5.,aodHeader->GetNumberOfITSClusters(1));\r
 436         fHistRefTracks->Fill(6.,aodHeader->GetNumberOfITSClusters(2));\r
 437         fHistRefTracks->Fill(7.,aodHeader->GetNumberOfITSClusters(3));\r
 438         fHistRefTracks->Fill(8.,aodHeader->GetNumberOfITSClusters(4));\r
 439       }\r
 440 \r
 441       const AliAODVertex *vertex = gAOD->GetPrimaryVertex();\r
 442       \r
 443       if(vertex) {\r
 444         Double32_t fCov[6];\r
 445         vertex->GetCovarianceMatrix(fCov);\r
 446         \r
 447         if(vertex->GetNContributors() > 0) {\r
 448           if(fCov[5] != 0) {\r
 449             fHistEventStats->Fill(3); //events with a proper vertex\r
 450             if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
 451               if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
 452                 if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
 453                   fHistEventStats->Fill(4); //analyzed events\r
 454                   fHistVx->Fill(vertex->GetX());\r
 455                   fHistVy->Fill(vertex->GetY());\r
 456                   fHistVz->Fill(vertex->GetZ());\r
 457                   \r
 458                   //Printf("There are %d tracks in this event", gAOD->GetNumberOfTracks());\r
 459                   for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {\r
 460                     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));\r
 461                     if (!aodTrack) {\r
 462                       Printf("ERROR: Could not receive track %d", iTracks);\r
 463                       continue;\r
 464                     }\r
 465                     \r
 466                     // AOD track cuts\r
 467                     \r
 468                     // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
 469                     // take only TPC only tracks \r
 470                     fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
 471                     if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
 472                     \r
 473                     Float_t pt  = aodTrack->Pt();\r
 474                     Float_t eta = aodTrack->Eta();\r
 475                     \r
 476                     Float_t DCAxy = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
 477                     Float_t DCAz  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
 478                     \r
 479                     \r
 480                     // Kinematics cuts from ESD track cuts\r
 481                     if( pt < fPtMin || pt > fPtMax)      continue;\r
 482                     if( eta < fEtaMin || eta > fEtaMax)  continue;\r
 483                     \r
 484                     // Extra DCA cuts (for systematic studies [!= -1])\r
 485                     if( fDCAxyCut != -1 && fDCAxyCut != -1){\r
 486                       if(TMath::Sqrt((DCAxy*DCAxy)/(fDCAxyCut*fDCAxyCut)+(DCAz*DCAz)/(fDCAzCut*fDCAzCut)) > 1 ){\r
 487                         continue;  // 2D cut\r
 488                       }\r
 489                     }\r
 490                     \r
 491                     // Extra TPC cuts (for systematic studies [!= -1])\r
 492                     if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
 493                       continue;\r
 494                     }\r
 495                     if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
 496                       continue;\r
 497                     }\r
 498                     \r
 499                     \r
 500                     // fill QA histograms\r
 501                     fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
 502                     fHistDCA->Fill(DCAz,DCAxy);\r
 503                     fHistChi2->Fill(aodTrack->Chi2perNDF());\r
 504                     fHistPt->Fill(pt);\r
 505                     fHistEta->Fill(eta);\r
 506                     fHistPhi->Fill(aodTrack->Phi()*TMath::RadToDeg());\r
 507                     \r
 508                     gNumberOfAcceptedTracks += 1;\r
 509                     // fill BF array\r
 510                     array->Add(aodTrack);\r
 511                     \r
 512                     // fill charge vector\r
 513                     chargeVector.push_back(aodTrack->Charge());\r
 514                     if(fRunShuffling) {\r
 515                       chargeVectorShuffle.push_back(aodTrack->Charge());\r
 516                     }\r
 517                   } //track loop\r
 518                 }//Vz cut\r
 519               }//Vy cut\r
 520             }//Vx cut\r
 521           }//proper vertex resolution\r
 522         }//proper number of contributors\r
 523       }//vertex object valid\r
 524     }//triggered event \r
 525   }//AOD analysis\r
 526 \r
 527   //MC-ESD analysis\r
 528   if(gAnalysisLevel == "MCESD") {\r
 529     AliMCEvent*  mcEvent = MCEvent(); \r
 530     if (!mcEvent) {\r
 531       Printf("ERROR: mcEvent not available");\r
 532       return;\r
 533     }\r
 534 \r
 535     AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
 536     if (!gESD) {\r
 537       Printf("ERROR: gESD not available");\r
 538       return;\r
 539     }\r
 540 \r
 541     // store offline trigger bits\r
 542     fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
 543 \r
 544     // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
 545     fHistEventStats->Fill(1); //all events\r
 546     Bool_t isSelected = kTRUE;\r
 547     if(fUseOfflineTrigger)\r
 548       isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
 549     if(isSelected) {\r
 550       fHistEventStats->Fill(2); //triggered events\r
 551 \r
 552       if(fUseCentrality) {\r
 553         //Centrality stuff\r
 554         AliCentrality *centrality = gESD->GetCentrality();\r
 555         //Int_t nCentrality = 0;\r
 556         //nCentrality = centrality->GetCentralityClass5(fCentralityEstimator.Data());\r
 557         //cout<<nCentrality<<" "<<centrality->IsEventInCentralityClass(fCentralityPercentileMin,fCentralityPercentileMax,fCentralityEstimator.Data())<<endl;\r
 558         \r
 559         // take only events inside centrality class\r
 560         if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
 561                                                  fCentralityPercentileMax,\r
 562                                                  fCentralityEstimator.Data()))\r
 563           return;\r
 564         \r
 565         // centrality QA (V0M)\r
 566         fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
 567       }\r
 568         \r
 569       const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
 570       if(vertex) {\r
 571         if(vertex->GetNContributors() > 0) {\r
 572           if(vertex->GetZRes() != 0) {\r
 573             fHistEventStats->Fill(3); //events with a proper vertex\r
 574             if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
 575               if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
 576                 if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
 577                   fHistEventStats->Fill(4); //analayzed events\r
 578                   fHistVx->Fill(vertex->GetXv());\r
 579                   fHistVy->Fill(vertex->GetYv());\r
 580                   fHistVz->Fill(vertex->GetZv());\r
 581                   \r
 582                   //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
 583                   for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
 584                     AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
 585                     if (!track) {\r
 586                       Printf("ERROR: Could not receive track %d", iTracks);\r
 587                       continue;\r
 588                     }   \r
 589                     \r
 590                     Int_t label = TMath::Abs(track->GetLabel());\r
 591                     if(label > mcEvent->GetNumberOfTracks()) continue;\r
 592                     if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
 593                     \r
 594                     AliMCParticle* mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
 595                     if(!mcTrack) continue;\r
 596 \r
 597                     // take only TPC only tracks\r
 598                     track_TPC   = new AliESDtrack();\r
 599                     if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
 600                     \r
 601                     //ESD track cuts\r
 602                     if(fESDtrackCuts) \r
 603                       if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
 604                     \r
 605                     // fill QA histograms\r
 606                     Float_t b[2];\r
 607                     Float_t bCov[3];\r
 608                     track_TPC->GetImpactParameters(b,bCov);\r
 609                     if (bCov[0]<=0 || bCov[2]<=0) {\r
 610                       AliDebug(1, "Estimated b resolution lower or equal zero!");\r
 611                       bCov[0]=0; bCov[2]=0;\r
 612                     }\r
 613                     \r
 614                     Int_t nClustersTPC = -1;\r
 615                     nClustersTPC = track_TPC->GetTPCNclsIter1();   // TPC standalone\r
 616                     //nClustersTPC = track->GetTPCclusters(0);   // global track\r
 617                     Float_t chi2PerClusterTPC = -1;\r
 618                     if (nClustersTPC!=0) {\r
 619                       chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
 620                       //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
 621                     }\r
 622                     \r
 623                     fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
 624                     fHistDCA->Fill(b[1],b[0]);\r
 625                     fHistChi2->Fill(chi2PerClusterTPC);\r
 626                     fHistPt->Fill(mcTrack->Pt());\r
 627                     fHistEta->Fill(mcTrack->Eta());\r
 628                     fHistPhi->Fill(mcTrack->Phi()*TMath::RadToDeg());\r
 629                     \r
 630                     // fill BF array\r
 631                     array->Add(track_TPC);\r
 632                     \r
 633                     // fill charge vector\r
 634                     chargeVector.push_back(track_TPC->Charge());\r
 635                     if(fRunShuffling){\r
 636                       chargeVectorShuffle.push_back(track_TPC->Charge());\r
 637                     }\r
 638                     \r
 639                     delete track_TPC;\r
 640                     \r
 641                   } //track loop\r
 642                 }//Vz cut\r
 643               }//Vy cut\r
 644             }//Vx cut\r
 645           }//proper vertex resolution\r
 646         }//proper number of contributors\r
 647       }//vertex object valid\r
 648     }//triggered event \r
 649   }//MC-ESD analysis\r
 650 \r
 651   //MC analysis\r
 652   else if(gAnalysisLevel == "MC") {\r
 653     AliMCEvent*  mcEvent = MCEvent(); \r
 654     if (!mcEvent) {\r
 655       Printf("ERROR: mcEvent not available");\r
 656       return;\r
 657     }\r
 658 \r
 659     Double_t gReactionPlane = 0., gImpactParameter = 0.;\r
 660     if(fUseCentrality) {\r
 661       //Get the MC header\r
 662       AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());\r
 663       if (headerH) {\r
 664         //Printf("=====================================================");\r
 665         //Printf("Reaction plane angle: %lf",headerH->ReactionPlaneAngle());\r
 666         //Printf("Impact parameter: %lf",headerH->ImpactParameter());\r
 667         //Printf("=====================================================");\r
 668         gReactionPlane = headerH->ReactionPlaneAngle();\r
 669         gImpactParameter = headerH->ImpactParameter();\r
 670       }\r
 671       // take only events inside centrality class\r
 672       if((fImpactParameterMin > gImpactParameter) || (fImpactParameterMax < gImpactParameter))\r
 673         return;\r
 674     }\r
 675     \r
 676     AliGenEventHeader *header = mcEvent->GenEventHeader();\r
 677     if(!header) return;\r
 678     \r
 679     TArrayF gVertexArray;\r
 680     header->PrimaryVertex(gVertexArray);\r
 681     //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
 682     //gVertexArray.At(0),\r
 683     //gVertexArray.At(1),\r
 684     //gVertexArray.At(2));\r
 685     fHistEventStats->Fill(3); //events with a proper vertex\r
 686     if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
 687       if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
 688         if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
 689           fHistEventStats->Fill(4); //analayzed events\r
 690           fHistVx->Fill(gVertexArray.At(0));\r
 691           fHistVy->Fill(gVertexArray.At(1));\r
 692           fHistVz->Fill(gVertexArray.At(2));\r
 693           \r
 694           Printf("There are %d tracks in this event", mcEvent->GetNumberOfPrimaries());\r
 695           for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfPrimaries(); iTracks++) {\r
 696             AliMCParticle* track = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iTracks));\r
 697             if (!track) {\r
 698               Printf("ERROR: Could not receive particle %d", iTracks);\r
 699               continue;\r
 700             }\r
 701             \r
 702             if( track->Pt() < fPtMin || track->Pt() > fPtMax)      \r
 703               continue;\r
 704             if( track->Eta() < fEtaMin || track->Eta() > fEtaMax)  \r
 705               continue;\r
 706             \r
 707             array->Add(track);\r
 708             \r
 709             // fill charge vector\r
 710             chargeVector.push_back(track->Charge());\r
 711             if(fRunShuffling) {\r
 712               chargeVectorShuffle.push_back(track->Charge());\r
 713             }\r
 714           } //track loop\r
 715         }//Vz cut\r
 716       }//Vy cut\r
 717     }//Vx cut\r
 718   }//MC analysis\r
 719   \r
 720   if(fUseMultiplicity) \r
 721     if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))\r
 722       return;\r
 723   \r
 724   // calculate balance function\r
 725   fBalance->CalculateBalance(array,chargeVector);\r
 726 \r
 727   if(fRunShuffling) {\r
 728     // shuffle charges\r
 729     random_shuffle( chargeVectorShuffle.begin(), chargeVectorShuffle.end() );\r
 730     fShuffledBalance->CalculateBalance(array,chargeVectorShuffle);\r
 731   }\r
 732   \r
 733   delete array;\r
 734 }      \r
 735 \r
 736 //________________________________________________________________________\r
 737 void  AliAnalysisTaskBF::FinishTaskOutput(){\r
 738   //Printf("END BF");\r
 739 \r
 740   if (!fBalance) {\r
 741     Printf("ERROR: fBalance not available");\r
 742     return;\r
 743   }  \r
 744   if(fRunShuffling) {\r
 745     if (!fShuffledBalance) {\r
 746       Printf("ERROR: fShuffledBalance not available");\r
 747       return;\r
 748     }\r
 749   }\r
 750 \r
 751 }\r
 752 \r
 753 //________________________________________________________________________\r
 754 void AliAnalysisTaskBF::Terminate(Option_t *) {\r
 755   // Draw result to the screen\r
 756   // Called once at the end of the query\r
 757 \r
 758   // not implemented ...\r
 759 \r
 760 }\r