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 "AliAODMCParticle.h" \r
  25 #include "AliCollisionGeometry.h"\r
  26 #include "AliGenEventHeader.h"\r
  27 #include "AliMCEventHandler.h"\r
  28 #include "AliMCEvent.h"\r
  29 #include "AliStack.h"\r
  30 #include "AliESDtrackCuts.h"\r
  31 #include "AliEventplane.h"\r
  32 #include "AliTHn.h"    \r
  33 #include "AliLog.h"\r
  34 #include "AliAnalysisUtils.h"\r
  35 \r
  36 #include "AliEventPoolManager.h"           \r
  37 \r
  38 #include "AliPID.h"                \r
  39 #include "AliPIDResponse.h"        \r
  40 #include "AliPIDCombined.h"        \r
  41 \r
  42 #include "AliAnalysisTaskBFPsi.h"\r
  43 #include "AliBalancePsi.h"\r
  44 #include "AliAnalysisTaskTriggeredBF.h"\r
  45 \r
  46 \r
  47 // Analysis task for the BF vs Psi code\r
  48 // Authors: Panos.Christakoglou@nikhef.nl\r
  49 \r
  50 using std::cout;\r
  51 using std::endl;\r
  52 \r
  53 ClassImp(AliAnalysisTaskBFPsi)\r
  54 \r
  55 //________________________________________________________________________\r
  56 AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) \r
  57 : AliAnalysisTaskSE(name),\r
  58   fArrayMC(0), //+++++++++++++\r
  59   fBalance(0),\r
  60   fRunShuffling(kFALSE),\r
  61   fShuffledBalance(0),\r
  62   fRunMixing(kFALSE),\r
  63   fRunMixingEventPlane(kFALSE),\r
  64   fMixingTracks(50000),\r
  65   fMixedBalance(0),\r
  66   fPoolMgr(0),\r
  67   fList(0),\r
  68   fListBF(0),\r
  69   fListBFS(0),\r
  70   fListBFM(0),\r
  71   fHistListPIDQA(0),\r
  72   fHistEventStats(0),\r
  73   fHistCentStats(0),\r
  74   fHistCentStatsUsed(0),\r
  75   fHistTriggerStats(0),\r
  76   fHistTrackStats(0),\r
  77   fHistVx(0),\r
  78   fHistVy(0),\r
  79   fHistVz(0),\r
  80   fHistEventPlane(0),\r
  81   fHistClus(0),\r
  82   fHistDCA(0),\r
  83   fHistChi2(0),\r
  84   fHistPt(0),\r
  85   fHistEta(0),\r
  86   fHistRapidity(0),\r
  87   fHistPhi(0),\r
  88   fHistEtaPhiPos(0),             \r
  89   fHistEtaPhiNeg(0), \r
  90   fHistPhiBefore(0),\r
  91   fHistPhiAfter(0),\r
  92   fHistPhiPos(0),\r
  93   fHistPhiNeg(0),\r
  94   fHistV0M(0),\r
  95   fHistRefTracks(0),\r
  96   fHistdEdxVsPTPCbeforePID(NULL),\r
  97   fHistBetavsPTOFbeforePID(NULL), \r
  98   fHistProbTPCvsPtbeforePID(NULL), \r
  99   fHistProbTOFvsPtbeforePID(NULL), \r
 100   fHistProbTPCTOFvsPtbeforePID(NULL),\r
 101   fHistNSigmaTPCvsPtbeforePID(NULL), \r
 102   fHistNSigmaTOFvsPtbeforePID(NULL), \r
 103   fHistdEdxVsPTPCafterPID(NULL),\r
 104   fHistBetavsPTOFafterPID(NULL), \r
 105   fHistProbTPCvsPtafterPID(NULL), \r
 106   fHistProbTOFvsPtafterPID(NULL), \r
 107   fHistProbTPCTOFvsPtafterPID(NULL),\r
 108   fHistNSigmaTPCvsPtafterPID(NULL), \r
 109   fHistNSigmaTOFvsPtafterPID(NULL),  \r
 110   fCentralityArrayBinsForCorrections(kCENTRALITY),\r
 111   fPIDResponse(0x0),\r
 112   fPIDCombined(0x0),\r
 113   fParticleOfInterest(kPion),\r
 114   fPidDetectorConfig(kTPCTOF),\r
 115   fUsePID(kFALSE),\r
 116   fUsePIDnSigma(kTRUE),\r
 117   fUsePIDPropabilities(kFALSE), \r
 118   fPIDNSigma(3.),\r
 119   fMinAcceptedPIDProbability(0.8),\r
 120   fElectronRejection(kFALSE),\r
 121   fElectronRejectionNSigma(-1.),\r
 122   fESDtrackCuts(0),\r
 123   fCentralityEstimator("V0M"),\r
 124   fUseCentrality(kFALSE),\r
 125   fCentralityPercentileMin(0.), \r
 126   fCentralityPercentileMax(5.),\r
 127   fImpactParameterMin(0.),\r
 128   fImpactParameterMax(20.),\r
 129   fUseMultiplicity(kFALSE),\r
 130   fNumberOfAcceptedTracksMin(0),\r
 131   fNumberOfAcceptedTracksMax(10000),\r
 132   fHistNumberOfAcceptedTracks(0),\r
 133   fUseOfflineTrigger(kFALSE),\r
 134   fCheckFirstEventInChunk(kFALSE),\r
 135   fCheckPileUp(kFALSE),\r
 136   fVxMax(0.3),\r
 137   fVyMax(0.3),\r
 138   fVzMax(10.),\r
 139   nAODtrackCutBit(128),\r
 140   fPtMin(0.3),\r
 141   fPtMax(1.5),\r
 142   fPtMinForCorrections(0.3),//=================================correction\r
 143   fPtMaxForCorrections(1.5),//=================================correction\r
 144   fPtBinForCorrections(36), //=================================correction\r
 145   fEtaMin(-0.8),\r
 146   fEtaMax(-0.8),\r
 147   fEtaMinForCorrections(-0.8),//=================================correction\r
 148   fEtaMaxForCorrections(0.8),//=================================correction\r
 149   fEtaBinForCorrections(16), //=================================correction\r
 150   fPhiMin(0.),\r
 151   fPhiMax(360.),\r
 152   fPhiMinForCorrections(0.),//=================================correction\r
 153   fPhiMaxForCorrections(360.),//=================================correction\r
 154   fPhiBinForCorrections(100), //=================================correction\r
 155   fDCAxyCut(-1),\r
 156   fDCAzCut(-1),\r
 157   fTPCchi2Cut(-1),\r
 158   fNClustersTPCCut(-1),\r
 159   fAcceptanceParameterization(0),\r
 160   fDifferentialV2(0),\r
 161   fUseFlowAfterBurner(kFALSE),\r
 162   fExcludeResonancesInMC(kFALSE),\r
 163   fUseMCPdgCode(kFALSE),\r
 164   fPDGCodeToBeAnalyzed(-1),\r
 165   fEventClass("EventPlane") \r
 166 {\r
 167   // Constructor\r
 168   // Define input and output slots here\r
 169   // Input slot #0 works with a TChain\r
 170 \r
 171   //======================================================correction\r
 172   for (Int_t i=0; i<kCENTRALITY; i++){\r
 173     fHistCorrectionPlus[i] = NULL; \r
 174     fHistCorrectionMinus[i] = NULL; \r
 175     fCentralityArrayForCorrections[i] = -1.;\r
 176   }\r
 177   //=====================================================correction\r
 178 \r
 179   DefineInput(0, TChain::Class());\r
 180   // Output slot #0 writes into a TH1 container\r
 181   DefineOutput(1, TList::Class());\r
 182   DefineOutput(2, TList::Class());\r
 183   DefineOutput(3, TList::Class());\r
 184   DefineOutput(4, TList::Class());\r
 185   DefineOutput(5, TList::Class());\r
 186 }\r
 187 \r
 188 //________________________________________________________________________\r
 189 AliAnalysisTaskBFPsi::~AliAnalysisTaskBFPsi() {\r
 190 \r
 191   // delete fBalance; \r
 192   // delete fShuffledBalance; \r
 193   // delete fList;\r
 194   // delete fListBF; \r
 195   // delete fListBFS;\r
 196 \r
 197   // delete fHistEventStats; \r
 198   // delete fHistTrackStats; \r
 199   // delete fHistVx; \r
 200   // delete fHistVy; \r
 201   // delete fHistVz; \r
 202 \r
 203   // delete fHistClus;\r
 204   // delete fHistDCA;\r
 205   // delete fHistChi2;\r
 206   // delete fHistPt;\r
 207   // delete fHistEta;\r
 208   // delete fHistPhi;\r
 209   // delete fHistEtaPhiPos;                      \r
 210   // delete fHistEtaPhiNeg;\r
 211   // delete fHistV0M;\r
 212 }\r
 213 \r
 214 //________________________________________________________________________\r
 215 void AliAnalysisTaskBFPsi::UserCreateOutputObjects() {\r
 216   // Create histograms\r
 217   // Called once\r
 218 \r
 219   // global switch disabling the reference \r
 220   // (to avoid "Replacing existing TH1" if several wagons are created in train)\r
 221   Bool_t oldStatus = TH1::AddDirectoryStatus();\r
 222   TH1::AddDirectory(kFALSE);\r
 223 \r
 224   if(!fBalance) {\r
 225     fBalance = new AliBalancePsi();\r
 226     fBalance->SetAnalysisLevel("ESD");\r
 227     fBalance->SetEventClass(fEventClass);\r
 228     //fBalance->SetNumberOfBins(-1,16);\r
 229     //fBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
 230   }\r
 231   if(fRunShuffling) {\r
 232     if(!fShuffledBalance) {\r
 233       fShuffledBalance = new AliBalancePsi();\r
 234       fShuffledBalance->SetAnalysisLevel("ESD");\r
 235       fShuffledBalance->SetEventClass(fEventClass);\r
 236       //fShuffledBalance->SetNumberOfBins(-1,16);\r
 237       //fShuffledBalance->SetInterval(-1,-0.8,0.8,16,0.,1.6,15.);\r
 238     }\r
 239   }\r
 240   if(fRunMixing) {\r
 241     if(!fMixedBalance) {\r
 242       fMixedBalance = new AliBalancePsi();\r
 243       fMixedBalance->SetAnalysisLevel("ESD");\r
 244       fMixedBalance->SetEventClass(fEventClass);\r
 245     }\r
 246   }\r
 247 \r
 248   //QA list\r
 249   fList = new TList();\r
 250   fList->SetName("listQA");\r
 251   fList->SetOwner();\r
 252 \r
 253   //Balance Function list\r
 254   fListBF = new TList();\r
 255   fListBF->SetName("listBF");\r
 256   fListBF->SetOwner();\r
 257 \r
 258   if(fRunShuffling) {\r
 259     fListBFS = new TList();\r
 260     fListBFS->SetName("listBFShuffled");\r
 261     fListBFS->SetOwner();\r
 262   }\r
 263 \r
 264   if(fRunMixing) {\r
 265     fListBFM = new TList();\r
 266     fListBFM->SetName("listTriggeredBFMixed");\r
 267     fListBFM->SetOwner();\r
 268   }\r
 269 \r
 270   //PID QA list\r
 271   if(fUsePID || fElectronRejection) {\r
 272     fHistListPIDQA = new TList();\r
 273     fHistListPIDQA->SetName("listQAPID");\r
 274     fHistListPIDQA->SetOwner();\r
 275   }\r
 276 \r
 277   //Event stats.\r
 278   TString gCutName[7] = {"Total","Offline trigger",\r
 279                          "Vertex","Analyzed","sel. Centrality","Not1stEvInChunk","No Pile-Up"};\r
 280   fHistEventStats = new TH2F("fHistEventStats",\r
 281                              "Event statistics;;Centrality percentile;N_{events}",\r
 282                              7,0.5,7.5,220,-5,105);\r
 283   for(Int_t i = 1; i <= 7; i++)\r
 284     fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());\r
 285   fList->Add(fHistEventStats);\r
 286 \r
 287   TString gCentName[13] = {"V0M","V0A","V0C","FMD","TRK","TKL","CL0","CL1","ZNA","ZPA","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};\r
 288   fHistCentStats = new TH2F("fHistCentStats",\r
 289                              "Centrality statistics;;Cent percentile",\r
 290                             13,-0.5,12.5,220,-5,105);\r
 291   for(Int_t i = 1; i <= 13; i++){\r
 292     fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());\r
 293     //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());  //++++++++++++++++++++++\r
 294   }\r
 295   fList->Add(fHistCentStats);\r
 296 \r
 297   fHistCentStatsUsed = new TH2F("fHistCentStatsUsed","Centrality statistics;;Cent percentile", 1,-0.5,0.5,220,-5,105); //++++++++++++++++++++++\r
 298   fHistCentStatsUsed->GetXaxis()->SetBinLabel(1,fCentralityEstimator.Data());  //++++++++++++++++++++++\r
 299   fList->Add(fHistCentStatsUsed); //++++++++++++++++++++++\r
 300 \r
 301   fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",1025,0,1025);\r
 302   fList->Add(fHistTriggerStats);\r
 303 \r
 304   fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TrackFilterBit;N_{events}",16,0,16);\r
 305   fList->Add(fHistTrackStats);\r
 306 \r
 307   fHistNumberOfAcceptedTracks = new TH2F("fHistNumberOfAcceptedTracks",";N_{acc.};Centrality percentile;Entries",4001,-0.5,4000.5,220,-5,105);\r
 308   fList->Add(fHistNumberOfAcceptedTracks);\r
 309 \r
 310   // Vertex distributions\r
 311   fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);\r
 312   fList->Add(fHistVx);\r
 313   fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);\r
 314   fList->Add(fHistVy);\r
 315   fHistVz = new TH2F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Centrality percentile;Entries",100,-20.,20.,220,-5,105);\r
 316   fList->Add(fHistVz);\r
 317 \r
 318   //Event plane\r
 319   fHistEventPlane = new TH2F("fHistEventPlane",";#Psi_{2} [deg.];Centrality percentile;Counts",100,0,360.,220,-5,105);\r
 320   fList->Add(fHistEventPlane);\r
 321 \r
 322   // QA histograms\r
 323   fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);\r
 324   fList->Add(fHistClus);\r
 325   fHistChi2 = new TH2F("fHistChi2","Chi2/NDF distribution;#chi^{2}/ndf;Centrality percentile",200,0,10,220,-5,105);\r
 326   fList->Add(fHistChi2);\r
 327   fHistDCA  = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5); \r
 328   fList->Add(fHistDCA);\r
 329   fHistPt   = new TH2F("fHistPt","p_{T} distribution;p_{T} (GeV/c);Centrality percentile",200,0,10,220,-5,105);\r
 330   fList->Add(fHistPt);\r
 331   fHistEta  = new TH2F("fHistEta","#eta distribution;#eta;Centrality percentile",200,-2,2,220,-5,105);\r
 332   fList->Add(fHistEta);\r
 333   fHistRapidity  = new TH2F("fHistRapidity","y distribution;y;Centrality percentile",200,-2,2,220,-5,105);\r
 334   fList->Add(fHistRapidity);\r
 335   fHistPhi  = new TH2F("fHistPhi","#phi distribution;#phi (rad);Centrality percentile",200,0.0,2.*TMath::Pi(),220,-5,105);\r
 336   fList->Add(fHistPhi);\r
 337   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
 338   fList->Add(fHistEtaPhiPos);                    \r
 339   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
 340   fList->Add(fHistEtaPhiNeg);\r
 341   fHistPhiBefore  = new TH2F("fHistPhiBefore","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
 342   fList->Add(fHistPhiBefore);\r
 343   fHistPhiAfter  = new TH2F("fHistPhiAfter","#phi distribution;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
 344   fList->Add(fHistPhiAfter);\r
 345   fHistPhiPos  = new TH2F("fHistPhiPos","#phi distribution for positive particles;#phi;Centrality percentile",200,0.,2*TMath::Pi(),220,-5,105);\r
 346   fList->Add(fHistPhiPos);\r
 347   fHistPhiNeg  = new TH2F("fHistPhiNeg","#phi distribution for negative particles;#phi;Centrality percentile",200,0.,2.*TMath::Pi(),220,-5,105);\r
 348   fList->Add(fHistPhiNeg);\r
 349   fHistV0M  = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);\r
 350   fList->Add(fHistV0M);\r
 351   TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};\r
 352   fHistRefTracks  = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 400, 0, 20000);\r
 353   for(Int_t i = 1; i <= 6; i++)\r
 354     fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());\r
 355   fList->Add(fHistRefTracks);\r
 356 \r
 357   // QA histograms for different cuts\r
 358   fList->Add(fBalance->GetQAHistHBTbefore());\r
 359   fList->Add(fBalance->GetQAHistHBTafter());\r
 360   fList->Add(fBalance->GetQAHistConversionbefore());\r
 361   fList->Add(fBalance->GetQAHistConversionafter());\r
 362   fList->Add(fBalance->GetQAHistPsiMinusPhi());\r
 363   fList->Add(fBalance->GetQAHistResonancesBefore());\r
 364   fList->Add(fBalance->GetQAHistResonancesRho());\r
 365   fList->Add(fBalance->GetQAHistResonancesK0());\r
 366   fList->Add(fBalance->GetQAHistResonancesLambda());\r
 367   fList->Add(fBalance->GetQAHistQbefore());\r
 368   fList->Add(fBalance->GetQAHistQafter());\r
 369 \r
 370   // Balance function histograms\r
 371   // Initialize histograms if not done yet\r
 372   if(!fBalance->GetHistNp()){\r
 373     AliWarning("Histograms not yet initialized! --> Will be done now");\r
 374     AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
 375     fBalance->InitHistograms();\r
 376   }\r
 377 \r
 378   if(fRunShuffling) {\r
 379     if(!fShuffledBalance->GetHistNp()) {\r
 380       AliWarning("Histograms (shuffling) not yet initialized! --> Will be done now");\r
 381       AliWarning("--> Add 'gBalance->InitHistograms()' in your configBalanceFunction");\r
 382       fShuffledBalance->InitHistograms();\r
 383     }\r
 384   }\r
 385 \r
 386   //for(Int_t a = 0; a < ANALYSIS_TYPES; a++){\r
 387   fListBF->Add(fBalance->GetHistNp());\r
 388   fListBF->Add(fBalance->GetHistNn());\r
 389   fListBF->Add(fBalance->GetHistNpn());\r
 390   fListBF->Add(fBalance->GetHistNnn());\r
 391   fListBF->Add(fBalance->GetHistNpp());\r
 392   fListBF->Add(fBalance->GetHistNnp());\r
 393 \r
 394   if(fRunShuffling) {\r
 395     fListBFS->Add(fShuffledBalance->GetHistNp());\r
 396     fListBFS->Add(fShuffledBalance->GetHistNn());\r
 397     fListBFS->Add(fShuffledBalance->GetHistNpn());\r
 398     fListBFS->Add(fShuffledBalance->GetHistNnn());\r
 399     fListBFS->Add(fShuffledBalance->GetHistNpp());\r
 400     fListBFS->Add(fShuffledBalance->GetHistNnp());\r
 401   }  \r
 402 \r
 403   if(fRunMixing) {\r
 404     fListBFM->Add(fMixedBalance->GetHistNp());\r
 405     fListBFM->Add(fMixedBalance->GetHistNn());\r
 406     fListBFM->Add(fMixedBalance->GetHistNpn());\r
 407     fListBFM->Add(fMixedBalance->GetHistNnn());\r
 408     fListBFM->Add(fMixedBalance->GetHistNpp());\r
 409     fListBFM->Add(fMixedBalance->GetHistNnp());\r
 410   }\r
 411   //}\r
 412 \r
 413 \r
 414   // Event Mixing\r
 415   if(fRunMixing){\r
 416     Int_t trackDepth = fMixingTracks; \r
 417     Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager\r
 418     \r
 419     // centrality bins\r
 420     Double_t centralityBins[] = {0.,1.,2.,3.,4.,5.,7.,10.,20.,30.,40.,50.,60.,70.,80.,100.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 421     Double_t* centbins        = centralityBins;\r
 422     Int_t nCentralityBins     = sizeof(centralityBins) / sizeof(Double_t) - 1;\r
 423 \r
 424     // multiplicity bins\r
 425     Double_t multiplicityBins[] = {0,10,20,30,40,50,60,70,80,100,100000}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 426     Double_t* multbins        = multiplicityBins;\r
 427     Int_t nMultiplicityBins     = sizeof(multiplicityBins) / sizeof(Double_t) - 1;\r
 428     \r
 429     // Zvtx bins\r
 430     Double_t vertexBins[] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 431     Double_t* vtxbins     = vertexBins;\r
 432     Int_t nVertexBins     = sizeof(vertexBins) / sizeof(Double_t) - 1;\r
 433     \r
 434     // Event plane angle (Psi) bins\r
 435     Double_t psiBins[] = {0.,45.,135.,215.,305.,360.}; // SHOULD BE DEDUCED FROM CREATED ALITHN!!!\r
 436     Double_t* psibins     = psiBins;\r
 437     Int_t nPsiBins     = sizeof(psiBins) / sizeof(Double_t) - 1;\r
 438     \r
 439     // run the event mixing also in bins of event plane (statistics!)\r
 440     if(fRunMixingEventPlane){\r
 441       if(fEventClass=="Multiplicity"){\r
 442         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
 443       }\r
 444       else{\r
 445         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins, nPsiBins, psibins);\r
 446       }\r
 447     }\r
 448     else{\r
 449       if(fEventClass=="Multiplicity"){\r
 450         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nMultiplicityBins, multbins, nVertexBins, vtxbins);\r
 451       }\r
 452       else{\r
 453         fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centbins, nVertexBins, vtxbins);\r
 454       }\r
 455     }\r
 456   }\r
 457 \r
 458   if(fESDtrackCuts) fList->Add(fESDtrackCuts);\r
 459 \r
 460   //====================PID========================//\r
 461   if(fUsePID) {\r
 462     fPIDCombined = new AliPIDCombined();\r
 463     fPIDCombined->SetDefaultTPCPriors();\r
 464 \r
 465     fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
 466     fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
 467     \r
 468     fHistBetavsPTOFbeforePID = new TH2D ("BetavsPTOFbefore","BetavsPTOFbefore", 1000, -10.0, 10., 1000, 0, 1.2); \r
 469     fHistListPIDQA->Add(fHistBetavsPTOFbeforePID); //addition\r
 470     \r
 471     fHistProbTPCvsPtbeforePID = new TH2D ("ProbTPCvsPtbefore","ProbTPCvsPtbefore", 1000, -10.0,10.0, 1000, 0, 2.0); \r
 472     fHistListPIDQA->Add(fHistProbTPCvsPtbeforePID); //addition \r
 473     \r
 474     fHistProbTOFvsPtbeforePID = new TH2D ("ProbTOFvsPtbefore","ProbTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
 475     fHistListPIDQA->Add(fHistProbTOFvsPtbeforePID); //addition \r
 476 \r
 477     fHistProbTPCTOFvsPtbeforePID =new TH2D ("ProbTPCTOFvsPtbefore","ProbTPCTOFvsPtbefore", 1000, -50, 50, 1000, 0, 2.0); \r
 478     fHistListPIDQA->Add(fHistProbTPCTOFvsPtbeforePID); //addition \r
 479     \r
 480     fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
 481     fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
 482     \r
 483     fHistNSigmaTOFvsPtbeforePID = new TH2D ("NSigmaTOFvsPtbefore","NSigmaTOFvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
 484     fHistListPIDQA->Add(fHistNSigmaTOFvsPtbeforePID); //addition \r
 485     \r
 486     fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
 487     fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
 488     \r
 489     fHistBetavsPTOFafterPID = new TH2D ("BetavsPTOFafter","BetavsPTOFafter", 1000, -10, 10, 1000, 0, 1.2); \r
 490     fHistListPIDQA->Add(fHistBetavsPTOFafterPID); //addition \r
 491     \r
 492     fHistProbTPCvsPtafterPID = new TH2D ("ProbTPCvsPtafter","ProbTPCvsPtafter", 1000, -10, 10, 1000, 0, 2); \r
 493     fHistListPIDQA->Add(fHistProbTPCvsPtafterPID); //addition \r
 494   \r
 495     fHistProbTOFvsPtafterPID = new TH2D ("ProbTOFvsPtafter","ProbTOFvsPtafter", 1000,  -10, 10, 1000, 0, 2); \r
 496     fHistListPIDQA->Add(fHistProbTOFvsPtafterPID); //addition  \r
 497     \r
 498     fHistProbTPCTOFvsPtafterPID =new TH2D ("ProbTPCTOFvsPtafter","ProbTPCTOFvsPtafter", 1000, -50, 50, 1000, 0, 2.0); \r
 499     fHistListPIDQA->Add(fHistProbTPCTOFvsPtafterPID); //addition \r
 500 \r
 501     fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
 502     fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition  \r
 503     \r
 504     fHistNSigmaTOFvsPtafterPID = new TH2D ("NSigmaTOFvsPtafter","NSigmaTOFvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
 505     fHistListPIDQA->Add(fHistNSigmaTOFvsPtafterPID); //addition \r
 506   }\r
 507 \r
 508   // for electron rejection only TPC nsigma histograms\r
 509   if(!fUsePID && fElectronRejection) {\r
 510  \r
 511     fHistdEdxVsPTPCbeforePID = new TH2D ("dEdxVsPTPCbefore","dEdxVsPTPCbefore", 1000, -10.0, 10.0, 1000, 0, 1000); \r
 512     fHistListPIDQA->Add(fHistdEdxVsPTPCbeforePID); //addition \r
 513     \r
 514     fHistNSigmaTPCvsPtbeforePID = new TH2D ("NSigmaTPCvsPtbefore","NSigmaTPCvsPtbefore", 1000, -10, 10, 1000, 0, 500); \r
 515     fHistListPIDQA->Add(fHistNSigmaTPCvsPtbeforePID); //addition \r
 516     \r
 517     fHistdEdxVsPTPCafterPID = new TH2D ("dEdxVsPTPCafter","dEdxVsPTPCafter", 1000, -10, 10, 1000, 0, 1000); \r
 518     fHistListPIDQA->Add(fHistdEdxVsPTPCafterPID); //addition \r
 519 \r
 520     fHistNSigmaTPCvsPtafterPID = new TH2D ("NSigmaTPCvsPtafter","NSigmaTPCvsPtafter", 1000, -10, 10, 1000, 0, 500); \r
 521     fHistListPIDQA->Add(fHistNSigmaTPCvsPtafterPID); //addition  \r
 522   }\r
 523   //====================PID========================//\r
 524 \r
 525   // Post output data.\r
 526   PostData(1, fList);\r
 527   PostData(2, fListBF);\r
 528   if(fRunShuffling) PostData(3, fListBFS);\r
 529   if(fRunMixing) PostData(4, fListBFM);\r
 530   if(fUsePID || fElectronRejection) PostData(5, fHistListPIDQA);       //PID\r
 531 \r
 532   TH1::AddDirectory(oldStatus);\r
 533 }\r
 534 \r
 535 \r
 536 //________________________________________________________________________\r
 537 void AliAnalysisTaskBFPsi::SetInputCorrection(TString filename, \r
 538                                               Int_t nCentralityBins, \r
 539                                               Double_t *centralityArrayForCorrections) {\r
 540   //Open files that will be used for correction\r
 541   fCentralityArrayBinsForCorrections = nCentralityBins;\r
 542   for (Int_t i=0; i<nCentralityBins; i++)\r
 543     fCentralityArrayForCorrections[i] = centralityArrayForCorrections[i];\r
 544 \r
 545   // No file specified -> run without corrections\r
 546   if(!filename.Contains(".root")) {\r
 547     AliInfo(Form("No correction file specified (= %s) --> run without corrections",filename.Data()));\r
 548     return;\r
 549   }\r
 550 \r
 551   //Open the input file\r
 552   TFile *f = TFile::Open(filename);\r
 553   if(!f->IsOpen()) {\r
 554     AliInfo(Form("File %s not found --> run without corrections",filename.Data()));\r
 555     return;\r
 556   }\r
 557     \r
 558   //TString listEffName = "";\r
 559   for (Int_t iCent = 0; iCent < fCentralityArrayBinsForCorrections-1; iCent++) {    \r
 560     //Printf("iCent %d:",iCent);    \r
 561     TString histoName = "fHistCorrectionPlus";\r
 562     histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));\r
 563     fHistCorrectionPlus[iCent]= dynamic_cast<TH3F *>(f->Get(histoName.Data()));\r
 564     if(!fHistCorrectionPlus[iCent]) {\r
 565       AliError(Form("fHist %s not found!!!",histoName.Data()));\r
 566       return;\r
 567     }\r
 568     \r
 569     histoName = "fHistCorrectionMinus";\r
 570     histoName += Form("%d-%d",(Int_t)(fCentralityArrayForCorrections[iCent]),(Int_t)(fCentralityArrayForCorrections[iCent+1]));\r
 571     fHistCorrectionMinus[iCent] = dynamic_cast<TH3F *>(f->Get(histoName.Data())); \r
 572     if(!fHistCorrectionMinus[iCent]) {\r
 573       AliError(Form("fHist %s not found!!!",histoName.Data()));\r
 574       return; \r
 575     }\r
 576   }//loop over centralities: ONLY the PbPb case is covered\r
 577 \r
 578   if(fHistCorrectionPlus[0]){\r
 579     fEtaMinForCorrections = fHistCorrectionPlus[0]->GetXaxis()->GetXmin();\r
 580     fEtaMaxForCorrections = fHistCorrectionPlus[0]->GetXaxis()->GetXmax();\r
 581     fEtaBinForCorrections = fHistCorrectionPlus[0]->GetNbinsX();\r
 582     \r
 583     fPtMinForCorrections = fHistCorrectionPlus[0]->GetYaxis()->GetXmin();\r
 584     fPtMaxForCorrections = fHistCorrectionPlus[0]->GetYaxis()->GetXmax();\r
 585     fPtBinForCorrections = fHistCorrectionPlus[0]->GetNbinsY();\r
 586     \r
 587     fPhiMinForCorrections = fHistCorrectionPlus[0]->GetZaxis()->GetXmin();\r
 588     fPhiMaxForCorrections = fHistCorrectionPlus[0]->GetZaxis()->GetXmax();\r
 589     fPhiBinForCorrections = fHistCorrectionPlus[0]->GetNbinsZ();\r
 590   }\r
 591 }\r
 592 \r
 593 //________________________________________________________________________\r
 594 void AliAnalysisTaskBFPsi::UserExec(Option_t *) {\r
 595   // Main loop\r
 596   // Called for each event\r
 597 \r
 598   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 599   Int_t gNumberOfAcceptedTracks = 0;\r
 600   Double_t gCentrality          = -1.;\r
 601   Double_t gReactionPlane       = -1.; \r
 602   Float_t bSign = 0.;\r
 603   \r
 604   // get the event (for generator level: MCEvent())\r
 605   AliVEvent* eventMain = NULL;\r
 606   if(gAnalysisLevel == "MC") {\r
 607     eventMain = dynamic_cast<AliVEvent*>(MCEvent()); \r
 608   }\r
 609   else{\r
 610     eventMain = dynamic_cast<AliVEvent*>(InputEvent()); \r
 611     \r
 612     // for HBT like cuts need magnetic field sign\r
 613     bSign = (eventMain->GetMagneticField() > 0) ? 1 : -1;\r
 614   }\r
 615   if(!eventMain) {\r
 616     AliError("eventMain not available");\r
 617     return;\r
 618   }\r
 619   \r
 620   // PID Response task active?\r
 621   if(fUsePID || fElectronRejection) {\r
 622     fPIDResponse = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetPIDResponse();\r
 623     if (!fPIDResponse) AliFatal("This Task needs the PID response attached to the inputHandler");\r
 624   }\r
 625   \r
 626   // check event cuts and fill event histograms\r
 627   if((gCentrality = IsEventAccepted(eventMain)) < 0){\r
 628     return;\r
 629   }\r
 630   \r
 631   //Compute Multiplicity or Centrality variable\r
 632   Double_t lMultiplicityVar = GetRefMultiOrCentrality( eventMain );\r
 633 \r
 634   // get the reaction plane\r
 635   if(fEventClass != "Multiplicity") {\r
 636     gReactionPlane = GetEventPlane(eventMain);\r
 637     fHistEventPlane->Fill(gReactionPlane,lMultiplicityVar);\r
 638     if(gReactionPlane < 0){\r
 639       return;\r
 640     }\r
 641   }\r
 642   \r
 643   // get the accepted tracks in main event\r
 644   TObjArray *tracksMain = GetAcceptedTracks(eventMain,lMultiplicityVar,gReactionPlane);\r
 645   gNumberOfAcceptedTracks = tracksMain->GetEntriesFast();\r
 646 \r
 647   //multiplicity cut (used in pp)\r
 648   fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks,lMultiplicityVar);\r
 649 \r
 650   // store charges of all accepted tracks,shuffle and reassign(two extra loops)\r
 651   TObjArray* tracksShuffled = NULL;\r
 652   if(fRunShuffling){\r
 653     tracksShuffled = GetShuffledTracks(tracksMain,lMultiplicityVar);\r
 654   }\r
 655   \r
 656   // Event mixing \r
 657   if (fRunMixing)\r
 658     {\r
 659       // 1. First get an event pool corresponding in mult (cent) and\r
 660       //    zvertex to the current event. Once initialized, the pool\r
 661       //    should contain nMix (reduced) events. This routine does not\r
 662       //    pre-scan the chain. The first several events of every chain\r
 663       //    will be skipped until the needed pools are filled to the\r
 664       //    specified depth. If the pool categories are not too rare, this\r
 665       //    should not be a problem. If they are rare, you could lose`\r
 666       //    statistics.\r
 667       \r
 668       // 2. Collect the whole pool's content of tracks into one TObjArray\r
 669       //    (bgTracks), which is effectively a single background super-event.\r
 670       \r
 671       // 3. The reduced and bgTracks arrays must both be passed into\r
 672       //    FillCorrelations(). Also nMix should be passed in, so a weight\r
 673       //    of 1./nMix can be applied.\r
 674       \r
 675       AliEventPool* pool = fPoolMgr->GetEventPool(lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane);\r
 676       \r
 677       if (!pool){\r
 678         AliFatal(Form("No pool found for centrality = %f, zVtx = %f, psi = %f", lMultiplicityVar, eventMain->GetPrimaryVertex()->GetZ(),gReactionPlane));\r
 679       }\r
 680       else{\r
 681         \r
 682         //pool->SetDebug(1);\r
 683 \r
 684         if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5){ \r
 685           \r
 686           \r
 687           Int_t nMix = pool->GetCurrentNEvents();\r
 688           //cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;\r
 689           \r
 690           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);\r
 691           //((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());\r
 692           //if (pool->IsReady())\r
 693           //((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);\r
 694           \r
 695           // Fill mixed-event histos here  \r
 696           for (Int_t jMix=0; jMix<nMix; jMix++) \r
 697             {\r
 698               TObjArray* tracksMixed = pool->GetEvent(jMix);\r
 699               fMixedBalance->CalculateBalance(gReactionPlane,tracksMain,tracksMixed,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
 700             }\r
 701         }\r
 702         \r
 703         // Update the Event pool\r
 704         pool->UpdatePool(tracksMain);\r
 705         //pool->PrintInfo();\r
 706         \r
 707       }//pool NULL check  \r
 708     }//run mixing\r
 709   \r
 710   // calculate balance function\r
 711   fBalance->CalculateBalance(gReactionPlane,tracksMain,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
 712   \r
 713   // calculate shuffled balance function\r
 714   if(fRunShuffling && tracksShuffled != NULL) {\r
 715     fShuffledBalance->CalculateBalance(gReactionPlane,tracksShuffled,NULL,bSign,lMultiplicityVar,eventMain->GetPrimaryVertex()->GetZ());\r
 716   }\r
 717 }      \r
 718 \r
 719 //________________________________________________________________________\r
 720 Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){\r
 721   // Checks the Event cuts\r
 722   // Fills Event statistics histograms\r
 723   \r
 724   Bool_t isSelectedMain = kTRUE;\r
 725   Float_t gCentrality = -1.;\r
 726   Float_t gRefMultiplicity = -1.;\r
 727   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 728 \r
 729   AliMCEvent *mcevent = dynamic_cast<AliMCEvent*>(event);\r
 730 \r
 731   fHistEventStats->Fill(1,gCentrality); //all events\r
 732 \r
 733   // check first event in chunk (is not needed for new reconstructions)\r
 734   if(fCheckFirstEventInChunk){\r
 735     AliAnalysisUtils ut;\r
 736     if(ut.IsFirstEventInChunk(event)) \r
 737       return -1.;\r
 738     fHistEventStats->Fill(6,gCentrality); \r
 739   }\r
 740 \r
 741   // check for pile-up event\r
 742   if(fCheckPileUp){\r
 743     AliAnalysisUtils ut;\r
 744     ut.SetUseMVPlpSelection(kTRUE);\r
 745     ut.SetUseOutOfBunchPileUp(kTRUE);\r
 746     if(ut.IsPileUpEvent(event))\r
 747       return -1.;\r
 748     fHistEventStats->Fill(7,gCentrality); \r
 749   }\r
 750 \r
 751   // Event trigger bits\r
 752   fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
 753   if(fUseOfflineTrigger)\r
 754     isSelectedMain = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
 755   \r
 756   if(isSelectedMain) {\r
 757     fHistEventStats->Fill(2,gCentrality); //triggered events\r
 758 \r
 759     //Centrality stuff \r
 760     if(fUseCentrality) {\r
 761       if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec") { //centrality in AOD header   //+++++++++++\r
 762         AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
 763         if(header){\r
 764           gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
 765 \r
 766           // QA for centrality estimators\r
 767           fHistCentStats->Fill(0.,header->GetCentralityP()->GetCentralityPercentile("V0M"));\r
 768           fHistCentStats->Fill(1.,header->GetCentralityP()->GetCentralityPercentile("V0A"));\r
 769           fHistCentStats->Fill(2.,header->GetCentralityP()->GetCentralityPercentile("V0C"));\r
 770           fHistCentStats->Fill(3.,header->GetCentralityP()->GetCentralityPercentile("FMD"));\r
 771           fHistCentStats->Fill(4.,header->GetCentralityP()->GetCentralityPercentile("TRK"));\r
 772           fHistCentStats->Fill(5.,header->GetCentralityP()->GetCentralityPercentile("TKL")); \r
 773           fHistCentStats->Fill(6.,header->GetCentralityP()->GetCentralityPercentile("CL0"));\r
 774           fHistCentStats->Fill(7.,header->GetCentralityP()->GetCentralityPercentile("CL1"));\r
 775           fHistCentStats->Fill(8.,header->GetCentralityP()->GetCentralityPercentile("ZNA"));\r
 776           fHistCentStats->Fill(9.,header->GetCentralityP()->GetCentralityPercentile("ZPA"));\r
 777           fHistCentStats->Fill(10.,header->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));\r
 778           fHistCentStats->Fill(11.,header->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));\r
 779           fHistCentStats->Fill(12.,header->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));\r
 780 \r
 781           // Centrality estimator USED   ++++++++++++++++++++++++++++++\r
 782           fHistCentStatsUsed->Fill(0.,header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()));\r
 783           \r
 784           // centrality QA (V0M)\r
 785           fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
 786           \r
 787           // centrality QA (reference tracks)\r
 788           fHistRefTracks->Fill(0.,header->GetRefMultiplicity());\r
 789           fHistRefTracks->Fill(1.,header->GetRefMultiplicityPos());\r
 790           fHistRefTracks->Fill(2.,header->GetRefMultiplicityNeg());\r
 791           fHistRefTracks->Fill(3.,header->GetTPConlyRefMultiplicity());\r
 792           fHistRefTracks->Fill(4.,header->GetNumberOfITSClusters(0));\r
 793           fHistRefTracks->Fill(5.,header->GetNumberOfITSClusters(1));\r
 794           fHistRefTracks->Fill(6.,header->GetNumberOfITSClusters(2));\r
 795           fHistRefTracks->Fill(7.,header->GetNumberOfITSClusters(3));\r
 796           fHistRefTracks->Fill(8.,header->GetNumberOfITSClusters(4));\r
 797         }//AOD header\r
 798       }//AOD\r
 799 \r
 800       else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
 801         AliCentrality *centrality = event->GetCentrality();\r
 802         gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
 803 \r
 804         // QA for centrality estimators\r
 805         fHistCentStats->Fill(0.,centrality->GetCentralityPercentile("V0M"));\r
 806         fHistCentStats->Fill(1.,centrality->GetCentralityPercentile("V0A"));\r
 807         fHistCentStats->Fill(2.,centrality->GetCentralityPercentile("V0C"));\r
 808         fHistCentStats->Fill(3.,centrality->GetCentralityPercentile("FMD"));\r
 809         fHistCentStats->Fill(4.,centrality->GetCentralityPercentile("TRK"));\r
 810         fHistCentStats->Fill(5.,centrality->GetCentralityPercentile("TKL"));\r
 811         fHistCentStats->Fill(6.,centrality->GetCentralityPercentile("CL0"));\r
 812         fHistCentStats->Fill(7.,centrality->GetCentralityPercentile("CL1"));\r
 813         fHistCentStats->Fill(8.,centrality->GetCentralityPercentile("ZNA"));\r
 814         fHistCentStats->Fill(9.,centrality->GetCentralityPercentile("ZPA"));\r
 815         fHistCentStats->Fill(10.,centrality->GetCentralityPercentile("V0MvsFMD"));\r
 816         fHistCentStats->Fill(11.,centrality->GetCentralityPercentile("TKLvsV0M"));\r
 817         fHistCentStats->Fill(12.,centrality->GetCentralityPercentile("ZEMvsZDC"));\r
 818 \r
 819         // Centrality estimator USED   ++++++++++++++++++++++++++++++\r
 820         fHistCentStatsUsed->Fill(0.,centrality->GetCentralityPercentile(fCentralityEstimator.Data()));\r
 821 \r
 822         // centrality QA (V0M)\r
 823         fHistV0M->Fill(event->GetVZEROData()->GetMTotV0A(), event->GetVZEROData()->GetMTotV0C());\r
 824       }//ESD\r
 825       else if(gAnalysisLevel == "MC"){\r
 826         Double_t gImpactParameter = 0.;\r
 827         if(mcevent) {\r
 828           AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(dynamic_cast<AliMCEvent*>(mcevent)->GenEventHeader());\r
 829           if(headerH){\r
 830             gImpactParameter = headerH->ImpactParameter();\r
 831             gCentrality      = gImpactParameter;\r
 832           }//MC header\r
 833         }//MC event cast\r
 834       }//MC\r
 835       else{\r
 836         gCentrality = -1.;\r
 837       }\r
 838     }//centrality\r
 839 \r
 840     //Multiplicity stuff \r
 841     if(fUseMultiplicity)\r
 842       gRefMultiplicity = GetRefMultiOrCentrality(event);\r
 843 \r
 844     // Event Vertex MC\r
 845     if(gAnalysisLevel == "MC") {\r
 846       if(!event) {\r
 847         AliError("mcEvent not available");\r
 848         return 0x0;\r
 849       }\r
 850       \r
 851       if(mcevent){\r
 852         AliGenEventHeader *header = dynamic_cast<AliGenEventHeader*>(mcevent->GenEventHeader());\r
 853         if(header){  \r
 854           TArrayF gVertexArray;\r
 855           header->PrimaryVertex(gVertexArray);\r
 856           //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
 857           //gVertexArray.At(0),\r
 858           //gVertexArray.At(1),\r
 859           //gVertexArray.At(2));\r
 860           if(fUseMultiplicity) \r
 861             fHistEventStats->Fill(3,gRefMultiplicity); //events with a proper vertex\r
 862           else \r
 863             fHistEventStats->Fill(3,gCentrality); //events with a proper vertex\r
 864           if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
 865             if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
 866               if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
 867                 if(fUseMultiplicity) \r
 868                   fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
 869                 else \r
 870                   fHistEventStats->Fill(4,gCentrality); //analyzed events\r
 871                 fHistVx->Fill(gVertexArray.At(0));\r
 872                 fHistVy->Fill(gVertexArray.At(1));\r
 873                 fHistVz->Fill(gVertexArray.At(2),gCentrality);\r
 874                 \r
 875                 // take only events inside centrality class\r
 876                 if(fUseCentrality) {\r
 877                   if((fImpactParameterMin < gCentrality) && (fImpactParameterMax > gCentrality)){\r
 878                     fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
 879                     return gCentrality;     \r
 880                   }//centrality class\r
 881                 }\r
 882                 // take events only within the same multiplicity class\r
 883                 else if(fUseMultiplicity) {\r
 884                   if((gRefMultiplicity > fNumberOfAcceptedTracksMin)||(gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
 885                     fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r
 886                     return gRefMultiplicity;\r
 887                   }\r
 888                 }//multiplicity range\r
 889               }//Vz cut\r
 890             }//Vy cut\r
 891           }//Vx cut\r
 892         }//header    \r
 893       }//MC event object\r
 894     }//MC\r
 895     \r
 896     // Event Vertex AOD, ESD, ESDMC\r
 897     else{\r
 898       const AliVVertex *vertex = event->GetPrimaryVertex();\r
 899       \r
 900       if(vertex) {\r
 901         Double32_t fCov[6];\r
 902         vertex->GetCovarianceMatrix(fCov);\r
 903         if(vertex->GetNContributors() > 0) {\r
 904           if(fCov[5] != 0) {\r
 905             if(fUseMultiplicity) \r
 906               fHistEventStats->Fill(3,gRefMultiplicity); //proper vertex\r
 907             else if(fUseCentrality)\r
 908               fHistEventStats->Fill(3,gCentrality); //proper vertex\r
 909             if(TMath::Abs(vertex->GetX()) < fVxMax) {\r
 910               if(TMath::Abs(vertex->GetY()) < fVyMax) {\r
 911                 if(TMath::Abs(vertex->GetZ()) < fVzMax) {\r
 912                   if(fUseMultiplicity) {\r
 913                     //cout<<"Filling 4 for multiplicity..."<<endl;\r
 914                     fHistEventStats->Fill(4,gRefMultiplicity);//analyzed events\r
 915                   }\r
 916                   else if(fUseCentrality) {\r
 917                     //cout<<"Filling 4 for centrality..."<<endl;\r
 918                     fHistEventStats->Fill(4,gCentrality); //analyzed events\r
 919                   }\r
 920                   fHistVx->Fill(vertex->GetX());\r
 921                   fHistVy->Fill(vertex->GetY());\r
 922                   fHistVz->Fill(vertex->GetZ(),gCentrality);\r
 923                   \r
 924                   // take only events inside centrality class\r
 925                   if(fUseCentrality) {\r
 926                     //cout<<"Centrality..."<<endl;\r
 927                     if((gCentrality > fCentralityPercentileMin) && (gCentrality < fCentralityPercentileMax)){\r
 928                       fHistEventStats->Fill(5,gCentrality); //events with correct centrality\r
 929                       return gCentrality;               \r
 930                     }//centrality class\r
 931                   }\r
 932                   // take events only within the same multiplicity class\r
 933                   else if(fUseMultiplicity) {\r
 934                     //cout<<"N(min): "<<fNumberOfAcceptedTracksMin<<" - N(max): "<<fNumberOfAcceptedTracksMax<<" - Nref: "<<gRefMultiplicity<<endl;\r
 935 \r
 936                     if((gRefMultiplicity > fNumberOfAcceptedTracksMin)||(gRefMultiplicity < fNumberOfAcceptedTracksMax)) {\r
 937                       fHistEventStats->Fill(5,gRefMultiplicity); //events with correct multiplicity\r
 938                       return gRefMultiplicity;\r
 939                     }\r
 940                   }//multiplicity range\r
 941                 }//Vz cut\r
 942               }//Vy cut\r
 943             }//Vx cut\r
 944           }//proper vertex resolution\r
 945         }//proper number of contributors\r
 946       }//vertex object valid\r
 947     }//triggered event \r
 948   }//AOD,ESD,ESDMC\r
 949   \r
 950   // in all other cases return -1 (event not accepted)\r
 951   return -1;\r
 952 }\r
 953 \r
 954 \r
 955 //________________________________________________________________________\r
 956 Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){\r
 957     // Checks the Event cuts\r
 958     // Fills Event statistics histograms\r
 959   \r
 960   Float_t gCentrality = -1.;\r
 961   Double_t gMultiplicity = 0.;\r
 962   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
 963 \r
 964   if(fEventClass == "Centrality"){\r
 965     if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header  //++++++++++++++\r
 966       AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
 967       if(header){\r
 968         gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());\r
 969       }//AOD header\r
 970     }//AOD\r
 971     \r
 972     else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD"){ // centrality class for ESDs or MC-ESDs\r
 973       AliCentrality *centrality = event->GetCentrality();\r
 974       gCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());\r
 975     }//ESD\r
 976     else if(gAnalysisLevel == "MC"){\r
 977       Double_t gImpactParameter = 0.;\r
 978       AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
 979       if(gMCEvent){\r
 980         AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());      \r
 981         if(headerH){\r
 982           gImpactParameter = headerH->ImpactParameter();\r
 983           gCentrality      = gImpactParameter;\r
 984         }//MC header\r
 985       }//MC event cast\r
 986     }//MC\r
 987     else{\r
 988       gCentrality = -1.;\r
 989     }\r
 990   }//End if "Centrality"\r
 991   if(fEventClass=="Multiplicity"&&gAnalysisLevel == "ESD"){\r
 992     AliESDEvent* gESDEvent = dynamic_cast<AliESDEvent*>(event);\r
 993     if(gESDEvent){\r
 994       gMultiplicity = fESDtrackCuts->GetReferenceMultiplicity(gESDEvent, AliESDtrackCuts::kTrackletsITSTPC,0.5);\r
 995     }//AliESDevent cast\r
 996   }\r
 997   else if(fEventClass=="Multiplicity"&&gAnalysisLevel == "AOD"){\r
 998     AliAODHeader *header = (AliAODHeader*) event->GetHeader();\r
 999     if(header){\r
1000       gMultiplicity = header->GetRefMultiplicity();\r
1001     }//AOD header\r
1002   }\r
1003   else if(fEventClass=="Multiplicity"&&gAnalysisLevel == "MC") {\r
1004     AliMCEvent* gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
1005     //Calculating the multiplicity as the number of charged primaries\r
1006     //within \pm 0.8 in eta and pT > 0.1 GeV/c\r
1007     for(Int_t iParticle = 0; iParticle < gMCEvent->GetNumberOfPrimaries(); iParticle++) {\r
1008       AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iParticle));\r
1009       if (!track) {\r
1010         AliError(Form("Could not receive particle %d", iParticle));\r
1011         continue;\r
1012       }\r
1013       \r
1014       //exclude non stable particles\r
1015       if(!(gMCEvent->IsPhysicalPrimary(iParticle))) continue;\r
1016       if(track->Pt() < 0.1)  continue;\r
1017       if(track->Eta() < fEtaMin || track->Eta() > fEtaMax)  continue;\r
1018       if(track->Charge() == 0) continue;\r
1019 \r
1020       gMultiplicity += 1;\r
1021     }//loop over primaries\r
1022   }//MC mode & multiplicity class\r
1023 \r
1024   Double_t lReturnVal = -100;\r
1025   if(fEventClass=="Multiplicity"){\r
1026     lReturnVal = gMultiplicity;\r
1027   }else if(fEventClass=="Centrality"){\r
1028     lReturnVal = gCentrality;\r
1029   }\r
1030   return lReturnVal;\r
1031 }\r
1032 \r
1033 //________________________________________________________________________\r
1034 Double_t AliAnalysisTaskBFPsi::GetEventPlane(AliVEvent *event){\r
1035   // Get the event plane\r
1036 \r
1037   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
1038 \r
1039   Float_t gVZEROEventPlane    = -10.;\r
1040   Float_t gReactionPlane      = -10.;\r
1041   Double_t qxTot = 0.0, qyTot = 0.0;\r
1042 \r
1043   //MC: from reaction plane\r
1044   if(gAnalysisLevel == "MC"){\r
1045     if(!event) {\r
1046       AliError("mcEvent not available");\r
1047       return 0x0;\r
1048     }\r
1049 \r
1050     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
1051     if(gMCEvent){\r
1052       AliCollisionGeometry* headerH = dynamic_cast<AliCollisionGeometry*>(gMCEvent->GenEventHeader());    \r
1053       if (headerH) {\r
1054         gReactionPlane = headerH->ReactionPlaneAngle();\r
1055         //gReactionPlane *= TMath::RadToDeg();\r
1056       }//MC header\r
1057     }//MC event cast\r
1058   }//MC\r
1059   \r
1060   // AOD,ESD,ESDMC: from VZERO Event Plane\r
1061   else{\r
1062    \r
1063     AliEventplane *ep = event->GetEventplane();\r
1064     if(ep){ \r
1065       gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);\r
1066       if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();\r
1067       //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();\r
1068       gReactionPlane = gVZEROEventPlane;\r
1069     }\r
1070   }//AOD,ESD,ESDMC\r
1071 \r
1072   return gReactionPlane;\r
1073 }\r
1074 \r
1075 //________________________________________________________________________\r
1076 Double_t AliAnalysisTaskBFPsi::GetTrackbyTrackCorrectionMatrix( Double_t vEta, \r
1077                                                                 Double_t vPhi, \r
1078                                                                 Double_t vPt, \r
1079                                                                 Short_t vCharge, \r
1080                                                                 Double_t gCentrality) {\r
1081   // -- Get efficiency correction of particle dependent on (eta, phi, pt, charge, centrality) \r
1082 \r
1083   Double_t correction = 1.;\r
1084   Int_t binEta = 0, binPt = 0, binPhi = 0;\r
1085 \r
1086   //Printf("EtaMAx: %lf - EtaMin: %lf - EtaBin: %lf", fEtaMaxForCorrections,fEtaMinForCorrections,fEtaBinForCorrections);\r
1087   if(fEtaBinForCorrections != 0) {\r
1088     Double_t widthEta = (fEtaMaxForCorrections - fEtaMinForCorrections)/fEtaBinForCorrections;\r
1089     if(fEtaMaxForCorrections != fEtaMinForCorrections) \r
1090       binEta = (Int_t)((vEta-fEtaMinForCorrections)/widthEta)+1;\r
1091   }\r
1092 \r
1093   if(fPtBinForCorrections != 0) {\r
1094     Double_t widthPt = (fPtMaxForCorrections - fPtMinForCorrections)/fPtBinForCorrections;\r
1095     if(fPtMaxForCorrections != fPtMinForCorrections) \r
1096       binPt = (Int_t)((vPt-fPtMinForCorrections)/widthPt) + 1;\r
1097   }\r
1098  \r
1099   if(fPhiBinForCorrections != 0) {\r
1100     Double_t widthPhi = (fPhiMaxForCorrections - fPhiMinForCorrections)/fPhiBinForCorrections;\r
1101     if(fPhiMaxForCorrections != fPhiMinForCorrections) \r
1102       binPhi = (Int_t)((vPhi-fPhiMinForCorrections)/widthPhi)+ 1;\r
1103   }\r
1104 \r
1105   Int_t gCentralityInt = -1;\r
1106   for (Int_t i=0; i<fCentralityArrayBinsForCorrections-1; i++){\r
1107     if((fCentralityArrayForCorrections[i] <= gCentrality)&&(gCentrality <= fCentralityArrayForCorrections[i+1])){\r
1108       gCentralityInt = i;\r
1109       break;\r
1110     }\r
1111   }  \r
1112 \r
1113   // centrality not in array --> no correction\r
1114   if(gCentralityInt < 0){\r
1115     correction = 1.;\r
1116   }\r
1117   else{\r
1118     \r
1119     //Printf("//=============CENTRALITY=============// %d:",gCentralityInt);\r
1120     \r
1121     if(fHistCorrectionPlus[gCentralityInt]){\r
1122       if (vCharge > 0) {\r
1123         correction = fHistCorrectionPlus[gCentralityInt]->GetBinContent(fHistCorrectionPlus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
1124         //Printf("CORRECTIONplus: %.2f | Centrality %d",correction,gCentralityInt);  \r
1125       }\r
1126       if (vCharge < 0) {\r
1127         correction = fHistCorrectionMinus[gCentralityInt]->GetBinContent(fHistCorrectionMinus[gCentralityInt]->GetBin(binEta, binPt, binPhi));\r
1128         //Printf("CORRECTIONminus: %.2f | Centrality %d",correction,gCentralityInt);\r
1129       }\r
1130     }\r
1131     else {\r
1132       correction = 1.;\r
1133     }\r
1134   }//centrality in array\r
1135   \r
1136   if (correction == 0.) { \r
1137     AliError(Form("Should not happen : bin content = 0. >> eta: %.2f | phi : %.2f | pt : %.2f | cent %d",vEta, vPhi, vPt, gCentralityInt)); \r
1138     correction = 1.; \r
1139   } \r
1140   \r
1141   return correction;\r
1142 }\r
1143 \r
1144 //________________________________________________________________________\r
1145 TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gCentrality, Double_t gReactionPlane){\r
1146   // Returns TObjArray with tracks after all track cuts (only for AOD!)\r
1147   // Fills QA histograms\r
1148 \r
1149   TString gAnalysisLevel = fBalance->GetAnalysisLevel();\r
1150 \r
1151   //output TObjArray holding all good tracks\r
1152   TObjArray* tracksAccepted = new TObjArray;\r
1153   tracksAccepted->SetOwner(kTRUE);\r
1154 \r
1155   Short_t vCharge;\r
1156   Double_t vEta;\r
1157   Double_t vY;\r
1158   Double_t vPhi;\r
1159   Double_t vPt;\r
1160 \r
1161 \r
1162   if(gAnalysisLevel == "AOD") { // handling of TPC only tracks different in AOD and ESD\r
1163     // Loop over tracks in event\r
1164     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1165       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
1166       if (!aodTrack) {\r
1167         AliError(Form("Could not receive track %d", iTracks));\r
1168         continue;\r
1169       }\r
1170       \r
1171       // AOD track cuts\r
1172       \r
1173       // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C\r
1174       // take only TPC only tracks \r
1175       //fHistTrackStats->Fill(aodTrack->GetFilterMap());\r
1176       for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
1177         fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
1178       }\r
1179       if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
1180 \r
1181       //===========================PID (so far only for electron rejection)===============================//                \r
1182       if(fElectronRejection) {\r
1183         \r
1184         // get the electron nsigma\r
1185         Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodTrack,(AliPID::EParticleType)AliPID::kElectron));\r
1186         \r
1187         //Fill QA before the PID\r
1188         fHistdEdxVsPTPCbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
1189         fHistNSigmaTPCvsPtbeforePID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
1190         //end of QA-before pid\r
1191         \r
1192         //Make the decision based on the n-sigma\r
1193         if(nSigma < fElectronRejectionNSigma) continue;\r
1194         \r
1195         //Fill QA after the PID\r
1196         fHistdEdxVsPTPCafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),aodTrack->GetTPCsignal());\r
1197         fHistNSigmaTPCvsPtafterPID -> Fill(aodTrack->P()*aodTrack->Charge(),nSigma); \r
1198       }\r
1199       //===========================end of PID (so far only for electron rejection)===============================//\r
1200       \r
1201       vCharge = aodTrack->Charge();\r
1202       vEta    = aodTrack->Eta();\r
1203       vY      = aodTrack->Y();\r
1204       vPhi    = aodTrack->Phi();// * TMath::RadToDeg();\r
1205       vPt     = aodTrack->Pt();\r
1206       \r
1207       Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
1208       Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
1209       \r
1210       \r
1211       // Kinematics cuts from ESD track cuts\r
1212       if( vPt < fPtMin || vPt > fPtMax)      continue;\r
1213       if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1214       \r
1215       // Extra DCA cuts (for systematic studies [!= -1])\r
1216       if( fDCAxyCut != -1 && fDCAzCut != -1){\r
1217         if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
1218           continue;  // 2D cut\r
1219         }\r
1220       }\r
1221       \r
1222       // Extra TPC cuts (for systematic studies [!= -1])\r
1223       if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
1224         continue;\r
1225       }\r
1226       if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
1227         continue;\r
1228       }\r
1229       \r
1230       // fill QA histograms\r
1231       fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
1232       fHistDCA->Fill(dcaZ,dcaXY);\r
1233       fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
1234       fHistPt->Fill(vPt,gCentrality);\r
1235       fHistEta->Fill(vEta,gCentrality);\r
1236       fHistRapidity->Fill(vY,gCentrality);\r
1237       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1238       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,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       \r
1243       //=======================================correction\r
1244       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1245       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
1246       \r
1247       // add the track to the TObjArray\r
1248       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
1249     }//track loop\r
1250   }// AOD analysis\r
1251 \r
1252   //==============================================================================================================\r
1253   else if(gAnalysisLevel == "MCAOD") {\r
1254     \r
1255     AliMCEvent* mcEvent = MCEvent();\r
1256     if (!mcEvent) {\r
1257       AliError("ERROR: Could not retrieve MC event");\r
1258     }\r
1259     else{\r
1260       \r
1261       for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {\r
1262         AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTracks); \r
1263         if (!aodTrack) {\r
1264           AliError(Form("ERROR: Could not receive track %d (mc loop)", iTracks));\r
1265           continue;\r
1266         }\r
1267         \r
1268         if(!aodTrack->IsPhysicalPrimary()) continue;   \r
1269         \r
1270         vCharge = aodTrack->Charge();\r
1271         vEta    = aodTrack->Eta();\r
1272         vY      = aodTrack->Y();\r
1273         vPhi    = aodTrack->Phi();// * TMath::RadToDeg();\r
1274         vPt     = aodTrack->Pt();\r
1275         \r
1276         // Kinematics cuts from ESD track cuts\r
1277         if( vPt < fPtMin || vPt > fPtMax)      continue;\r
1278         if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1279         \r
1280         // Remove neutral tracks\r
1281         if( vCharge == 0 ) continue;\r
1282         \r
1283         //Exclude resonances\r
1284         if(fExcludeResonancesInMC) {\r
1285           \r
1286           Bool_t kExcludeParticle = kFALSE;\r
1287           Int_t gMotherIndex = aodTrack->GetMother();\r
1288           if(gMotherIndex != -1) {\r
1289             AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
1290             if(motherTrack) {\r
1291               Int_t pdgCodeOfMother = motherTrack->GetPdgCode();\r
1292               //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
1293               //if(pdgCodeOfMother == 113) {\r
1294               if(pdgCodeOfMother == 113  // rho0\r
1295                  || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
1296                  // || pdgCodeOfMother == 221  // eta\r
1297                  // || pdgCodeOfMother == 331  // eta'\r
1298                  // || pdgCodeOfMother == 223  // omega\r
1299                  // || pdgCodeOfMother == 333  // phi\r
1300                  || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0\r
1301                  // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*\r
1302                  // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*\r
1303                  || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r
1304                  || pdgCodeOfMother == 111  // pi0 Dalitz\r
1305                  || pdgCodeOfMother == 22  // photon\r
1306                  ) {\r
1307                 kExcludeParticle = kTRUE;\r
1308               }\r
1309             }\r
1310           }\r
1311           \r
1312           //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1313           if(kExcludeParticle) continue;\r
1314         }\r
1315         \r
1316         // fill QA histograms\r
1317         fHistPt->Fill(vPt,gCentrality);\r
1318         fHistEta->Fill(vEta,gCentrality);\r
1319         fHistRapidity->Fill(vY,gCentrality);\r
1320         if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1321         else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1322         fHistPhi->Fill(vPhi,gCentrality);\r
1323         if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);                \r
1324         else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1325         \r
1326         //=======================================correction\r
1327         Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1328         //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);   \r
1329         \r
1330         // add the track to the TObjArray\r
1331         tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
1332       }//aodTracks\r
1333     }//MC event\r
1334   }//MCAOD\r
1335   //==============================================================================================================\r
1336 \r
1337   //==============================================================================================================\r
1338   else if(gAnalysisLevel == "MCAODrec") {\r
1339     \r
1340     /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());\r
1341     if (!fAOD) {\r
1342       printf("ERROR: fAOD not available\n");\r
1343       return;\r
1344       }*/\r
1345     \r
1346     fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName())); \r
1347     if (!fArrayMC) { \r
1348        AliError("No array of MC particles found !!!");\r
1349     }\r
1350 \r
1351     AliMCEvent* mcEvent = MCEvent();\r
1352     if (!mcEvent) {\r
1353        AliError("ERROR: Could not retrieve MC event");\r
1354     }\r
1355      \r
1356     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1357       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
1358       if (!aodTrack) {\r
1359         AliError(Form("Could not receive track %d", iTracks));\r
1360         continue;\r
1361       }\r
1362       \r
1363       for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
1364         fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
1365       }\r
1366       if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
1367       \r
1368       vCharge = aodTrack->Charge();\r
1369       vEta    = aodTrack->Eta();\r
1370       vY      = aodTrack->Y();\r
1371       vPhi    = aodTrack->Phi();// * TMath::RadToDeg();\r
1372       vPt     = aodTrack->Pt();\r
1373       \r
1374       Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
1375       Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
1376             \r
1377       // Kinematics cuts from ESD track cuts\r
1378       if( vPt < fPtMin || vPt > fPtMax)      continue;\r
1379       if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1380       \r
1381       // Extra DCA cuts (for systematic studies [!= -1])\r
1382       if( fDCAxyCut != -1 && fDCAzCut != -1){\r
1383         if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
1384           continue;  // 2D cut\r
1385         }\r
1386       }\r
1387       \r
1388       // Extra TPC cuts (for systematic studies [!= -1])\r
1389       if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
1390         continue;\r
1391       }\r
1392       if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
1393         continue;\r
1394       }\r
1395 \r
1396       //Exclude resonances\r
1397       if(fExcludeResonancesInMC) {\r
1398         \r
1399         Bool_t kExcludeParticle = kFALSE;\r
1400 \r
1401         Int_t label = TMath::Abs(aodTrack->GetLabel());\r
1402         AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);\r
1403       \r
1404         if (AODmcTrack){ \r
1405           //if (AODmcTrack->IsPhysicalPrimary()){\r
1406           \r
1407           Int_t gMotherIndex = AODmcTrack->GetMother();\r
1408           if(gMotherIndex != -1) {\r
1409             AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
1410             if(motherTrack) {\r
1411               Int_t pdgCodeOfMother = motherTrack->GetPdgCode();\r
1412               if(pdgCodeOfMother == 113  // rho0\r
1413                  || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
1414                  // || pdgCodeOfMother == 221  // eta\r
1415                  // || pdgCodeOfMother == 331  // eta'\r
1416                  // || pdgCodeOfMother == 223  // omega\r
1417                  // || pdgCodeOfMother == 333  // phi\r
1418                  || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0\r
1419                  // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*\r
1420                  // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*\r
1421                  || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r
1422                  || pdgCodeOfMother == 111  // pi0 Dalitz\r
1423                  || pdgCodeOfMother == 22   // photon\r
1424                  ) {\r
1425                 kExcludeParticle = kTRUE;\r
1426               }\r
1427             }\r
1428           }\r
1429         }       \r
1430         //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1431         if(kExcludeParticle) continue;\r
1432       }\r
1433       \r
1434       // fill QA histograms\r
1435       fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
1436       fHistDCA->Fill(dcaZ,dcaXY);\r
1437       fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
1438       fHistPt->Fill(vPt,gCentrality);\r
1439       fHistEta->Fill(vEta,gCentrality);\r
1440       fHistRapidity->Fill(vY,gCentrality);\r
1441       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1442       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1443       fHistPhi->Fill(vPhi,gCentrality);\r
1444       if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);                  \r
1445       else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1446       \r
1447       //=======================================correction\r
1448       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1449       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
1450       \r
1451       // add the track to the TObjArray\r
1452       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
1453     }//track loop\r
1454   }//MCAODrec\r
1455   //==============================================================================================================\r
1456 \r
1457   else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD\r
1458 \r
1459     AliESDtrack *trackTPC   = NULL;\r
1460     AliMCParticle *mcTrack   = NULL;\r
1461 \r
1462     AliMCEvent*  mcEvent     = NULL;\r
1463     \r
1464     // for MC ESDs use also MC information\r
1465     if(gAnalysisLevel == "MCESD"){\r
1466       mcEvent = MCEvent(); \r
1467       if (!mcEvent) {\r
1468         AliError("mcEvent not available");\r
1469         return tracksAccepted;\r
1470       }\r
1471     }\r
1472     \r
1473     // Loop over tracks in event\r
1474     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1475       AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));\r
1476       if (!track) {\r
1477         AliError(Form("Could not receive track %d", iTracks));\r
1478         continue;\r
1479       } \r
1480 \r
1481       // for MC ESDs use also MC information --> MC track not used further???\r
1482       if(gAnalysisLevel == "MCESD"){\r
1483         Int_t label = TMath::Abs(track->GetLabel());\r
1484         if(label > mcEvent->GetNumberOfTracks()) continue;\r
1485         if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
1486         \r
1487         mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
1488         if(!mcTrack) continue;\r
1489       }\r
1490 \r
1491       // take only TPC only tracks\r
1492       trackTPC   = new AliESDtrack();\r
1493       if(!track->FillTPCOnlyTrack(*trackTPC)) continue;\r
1494       \r
1495       //ESD track cuts\r
1496       if(fESDtrackCuts) \r
1497         if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;\r
1498       \r
1499       // fill QA histograms\r
1500       Float_t b[2];\r
1501       Float_t bCov[3];\r
1502       trackTPC->GetImpactParameters(b,bCov);\r
1503       if (bCov[0]<=0 || bCov[2]<=0) {\r
1504         AliDebug(1, "Estimated b resolution lower or equal zero!");\r
1505         bCov[0]=0; bCov[2]=0;\r
1506       }\r
1507       \r
1508       Int_t nClustersTPC = -1;\r
1509       nClustersTPC = trackTPC->GetTPCNclsIter1();   // TPC standalone\r
1510       //nClustersTPC = track->GetTPCclusters(0);   // global track\r
1511       Float_t chi2PerClusterTPC = -1;\r
1512       if (nClustersTPC!=0) {\r
1513         chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
1514         //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
1515       }\r
1516       \r
1517       //===========================PID===============================//             \r
1518       if(fUsePID) {\r
1519         Double_t prob[AliPID::kSPECIES]={0.};\r
1520         Double_t probTPC[AliPID::kSPECIES]={0.};\r
1521         Double_t probTOF[AliPID::kSPECIES]={0.};\r
1522         Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
1523         \r
1524         Double_t nSigma = 0.;\r
1525         UInt_t detUsedTPC = 0;\r
1526         UInt_t detUsedTOF = 0;\r
1527         UInt_t detUsedTPCTOF = 0;\r
1528         \r
1529         //Decide what detector configuration we want to use\r
1530         switch(fPidDetectorConfig) {\r
1531         case kTPCpid:\r
1532           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
1533           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
1534           detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
1535           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1536             prob[iSpecies] = probTPC[iSpecies];\r
1537           break;\r
1538         case kTOFpid:\r
1539           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
1540           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
1541           detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
1542           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1543             prob[iSpecies] = probTOF[iSpecies];\r
1544           break;\r
1545         case kTPCTOF:\r
1546           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
1547           detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
1548           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1549             prob[iSpecies] = probTPCTOF[iSpecies];\r
1550           break;\r
1551         default:\r
1552           break;\r
1553         }//end switch: define detector mask\r
1554         \r
1555         //Filling the PID QA\r
1556         Double_t tofTime = -999., length = 999., tof = -999.;\r
1557         Double_t c = TMath::C()*1.E-9;// m/ns\r
1558         Double_t beta = -999.;\r
1559         Double_t  nSigmaTOFForParticleOfInterest = -999.;\r
1560         if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
1561              (track->IsOn(AliESDtrack::kTIME))  ) { \r
1562           tofTime = track->GetTOFsignal();//in ps\r
1563           length = track->GetIntegratedLength();\r
1564           tof = tofTime*1E-3; // ns     \r
1565           \r
1566           if (tof <= 0) {\r
1567             //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
1568             continue;\r
1569           }\r
1570           if (length <= 0){\r
1571             //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
1572             continue;\r
1573           }\r
1574           \r
1575           length = length*0.01; // in meters\r
1576           tof = tof*c;\r
1577           beta = length/tof;\r
1578           \r
1579           nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
1580           fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
1581           fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1582           fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1583         }//TOF signal \r
1584         \r
1585         \r
1586         Double_t  nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
1587         fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1588         fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1589         fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1590         fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1591         //end of QA-before pid\r
1592         \r
1593         if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
1594           //Make the decision based on the n-sigma\r
1595           if(fUsePIDnSigma) {\r
1596             if(nSigma > fPIDNSigma) continue;}\r
1597           \r
1598           //Make the decision based on the bayesian\r
1599           else if(fUsePIDPropabilities) {\r
1600             if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
1601             if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;      \r
1602           }\r
1603           \r
1604           //Fill QA after the PID\r
1605           fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
1606           fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1607           fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1608           \r
1609           fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1610           fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1611           fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1612           fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1613         }\r
1614       }\r
1615       //===========================PID===============================//\r
1616       vCharge = trackTPC->Charge();\r
1617       vY      = trackTPC->Y();\r
1618       vEta    = trackTPC->Eta();\r
1619       vPhi    = trackTPC->Phi();// * TMath::RadToDeg();\r
1620       vPt     = trackTPC->Pt();\r
1621       fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);\r
1622       fHistDCA->Fill(b[1],b[0]);\r
1623       fHistChi2->Fill(chi2PerClusterTPC,gCentrality);\r
1624       fHistPt->Fill(vPt,gCentrality);\r
1625       fHistEta->Fill(vEta,gCentrality);\r
1626       fHistPhi->Fill(vPhi,gCentrality);\r
1627       if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
1628       else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1629       fHistRapidity->Fill(vY,gCentrality);\r
1630       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1631       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1632       \r
1633       //=======================================correction\r
1634       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1635       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
1636 \r
1637       // add the track to the TObjArray\r
1638       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));   \r
1639 \r
1640       delete trackTPC;\r
1641     }//track loop\r
1642   }// ESD analysis\r
1643 \r
1644   else if(gAnalysisLevel == "MC"){\r
1645     if(!event) {\r
1646       AliError("mcEvent not available");\r
1647       return 0x0;\r
1648     }\r
1649 \r
1650     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
1651     if(gMCEvent) {\r
1652       // Loop over tracks in event\r
1653       for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {\r
1654         AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));\r
1655         if (!track) {\r
1656           AliError(Form("Could not receive particle %d", iTracks));\r
1657           continue;\r
1658         }\r
1659         \r
1660         //exclude non stable particles\r
1661         if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;\r
1662         \r
1663         vCharge = track->Charge();\r
1664         vEta    = track->Eta();\r
1665         vPt     = track->Pt();\r
1666         vY      = track->Y();\r
1667         \r
1668         if( vPt < fPtMin || vPt > fPtMax)      \r
1669           continue;\r
1670         if (!fUsePID) {\r
1671           if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1672         }\r
1673         else if (fUsePID){\r
1674           if( vY < fEtaMin || vY > fEtaMax)  continue;\r
1675         }\r
1676 \r
1677         // Remove neutral tracks\r
1678         if( vCharge == 0 ) continue;\r
1679         \r
1680         //analyze one set of particles\r
1681         if(fUseMCPdgCode) {\r
1682           TParticle *particle = track->Particle();\r
1683           if(!particle) continue;\r
1684           \r
1685           Int_t gPdgCode = particle->GetPdgCode();\r
1686           if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
1687             continue;\r
1688         }\r
1689         \r
1690         //Use the acceptance parameterization\r
1691         if(fAcceptanceParameterization) {\r
1692           Double_t gRandomNumber = gRandom->Rndm();\r
1693           if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
1694             continue;\r
1695         }\r
1696         \r
1697         //Exclude resonances\r
1698         if(fExcludeResonancesInMC) {\r
1699           TParticle *particle = track->Particle();\r
1700           if(!particle) continue;\r
1701           \r
1702           Bool_t kExcludeParticle = kFALSE;\r
1703           Int_t gMotherIndex = particle->GetFirstMother();\r
1704           if(gMotherIndex != -1) {\r
1705             AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));\r
1706             if(motherTrack) {\r
1707               TParticle *motherParticle = motherTrack->Particle();\r
1708               if(motherParticle) {\r
1709                 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
1710                 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
1711                 if(pdgCodeOfMother == 113  // rho0\r
1712                    || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
1713                    // || pdgCodeOfMother == 221  // eta\r
1714                    // || pdgCodeOfMother == 331  // eta'\r
1715                    // || pdgCodeOfMother == 223  // omega\r
1716                    // || pdgCodeOfMother == 333  // phi\r
1717                    || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0\r
1718                    // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*\r
1719                    // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*\r
1720                    || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r
1721                    || pdgCodeOfMother == 111  // pi0 Dalitz\r
1722                    ) {\r
1723                   kExcludeParticle = kTRUE;\r
1724                 }\r
1725               }\r
1726             }\r
1727           }\r
1728           \r
1729           //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1730           if(kExcludeParticle) continue;\r
1731         }\r
1732         \r
1733         vPhi    = track->Phi();\r
1734         //Printf("phi (before): %lf",vPhi);\r
1735         \r
1736         fHistPt->Fill(vPt,gCentrality);\r
1737         fHistEta->Fill(vEta,gCentrality);\r
1738         fHistPhi->Fill(vPhi,gCentrality);\r
1739         if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
1740         else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1741         //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1742         fHistRapidity->Fill(vY,gCentrality);\r
1743         //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1744         //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1745         if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1746         else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1747         \r
1748         //Flow after burner\r
1749         if(fUseFlowAfterBurner) {\r
1750           Double_t precisionPhi = 0.001;\r
1751           Int_t maxNumberOfIterations = 100;\r
1752           \r
1753           Double_t phi0 = vPhi;\r
1754           Double_t gV2 = fDifferentialV2->Eval(vPt);\r
1755           \r
1756           for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
1757             Double_t phiprev = vPhi;\r
1758             Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));\r
1759             Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad())); \r
1760             vPhi -= fl/fp;\r
1761             if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;\r
1762           }\r
1763           //Printf("phi (after): %lf\n",vPhi);\r
1764           Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();\r
1765           if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();\r
1766           fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);\r
1767           \r
1768           Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();\r
1769           if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();\r
1770           fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);\r
1771           \r
1772         }\r
1773         \r
1774         //vPhi *= TMath::RadToDeg();\r
1775         \r
1776       //=======================================correction\r
1777       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1778       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
1779 \r
1780         tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
1781       } //track loop\r
1782     }//MC event object\r
1783   }//MC\r
1784   \r
1785   return tracksAccepted;  \r
1786 }\r
1787 \r
1788 //________________________________________________________________________\r
1789 TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){\r
1790   // Clones TObjArray and returns it with tracks after shuffling the charges\r
1791 \r
1792   TObjArray* tracksShuffled = new TObjArray;\r
1793   tracksShuffled->SetOwner(kTRUE);\r
1794 \r
1795   vector<Short_t> *chargeVector = new vector<Short_t>;   //original charge of accepted tracks \r
1796 \r
1797   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)\r
1798   {\r
1799     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
1800     chargeVector->push_back(track->Charge());\r
1801   }  \r
1802  \r
1803   random_shuffle(chargeVector->begin(), chargeVector->end());\r
1804   \r
1805   for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){\r
1806     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
1807     //==============================correction\r
1808     Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);\r
1809     //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
1810     tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));\r
1811   }\r
1812 \r
1813   delete chargeVector;\r
1814    \r
1815   return tracksShuffled;\r
1816 }\r
1817 \r
1818 \r
1819 //________________________________________________________________________\r
1820 void  AliAnalysisTaskBFPsi::FinishTaskOutput(){\r
1821   //Printf("END BF");\r
1822 \r
1823   if (!fBalance) {\r
1824     AliError("fBalance not available");\r
1825     return;\r
1826   }  \r
1827   if(fRunShuffling) {\r
1828     if (!fShuffledBalance) {\r
1829       AliError("fShuffledBalance not available");\r
1830       return;\r
1831     }\r
1832   }\r
1833 \r
1834 }\r
1835 \r
1836 //________________________________________________________________________\r
1837 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {\r
1838   // Draw result to the screen\r
1839   // Called once at the end of the query\r
1840 \r
1841   // not implemented ...\r
1842 \r
1843 }\r
1844 \r