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