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