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