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