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