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                  ) {\r
1306                 kExcludeParticle = kTRUE;\r
1307               }\r
1308             }\r
1309           }\r
1310           \r
1311           //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1312           if(kExcludeParticle) continue;\r
1313         }\r
1314         \r
1315         // fill QA histograms\r
1316         fHistPt->Fill(vPt,gCentrality);\r
1317         fHistEta->Fill(vEta,gCentrality);\r
1318         fHistRapidity->Fill(vY,gCentrality);\r
1319         if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1320         else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1321         fHistPhi->Fill(vPhi,gCentrality);\r
1322         if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);                \r
1323         else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1324         \r
1325         //=======================================correction\r
1326         Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1327         //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);   \r
1328         \r
1329         // add the track to the TObjArray\r
1330         tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
1331       }//aodTracks\r
1332     }//MC event\r
1333   }//MCAOD\r
1334   //==============================================================================================================\r
1335 \r
1336   //==============================================================================================================\r
1337   else if(gAnalysisLevel == "MCAODrec") {\r
1338     \r
1339     /* fAOD = dynamic_cast<AliAODEvent*>(InputEvent());\r
1340     if (!fAOD) {\r
1341       printf("ERROR: fAOD not available\n");\r
1342       return;\r
1343       }*/\r
1344     \r
1345     fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName())); \r
1346     if (!fArrayMC) { \r
1347        AliError("No array of MC particles found !!!");\r
1348     }\r
1349 \r
1350     AliMCEvent* mcEvent = MCEvent();\r
1351     if (!mcEvent) {\r
1352        AliError("ERROR: Could not retrieve MC event");\r
1353     }\r
1354      \r
1355     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1356       AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(event->GetTrack(iTracks));\r
1357       if (!aodTrack) {\r
1358         AliError(Form("Could not receive track %d", iTracks));\r
1359         continue;\r
1360       }\r
1361       \r
1362       for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){\r
1363         fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<<iTrackBit));\r
1364       }\r
1365       if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;\r
1366       \r
1367       vCharge = aodTrack->Charge();\r
1368       vEta    = aodTrack->Eta();\r
1369       vY      = aodTrack->Y();\r
1370       vPhi    = aodTrack->Phi();// * TMath::RadToDeg();\r
1371       vPt     = aodTrack->Pt();\r
1372       \r
1373       Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)\r
1374       Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)\r
1375             \r
1376       // Kinematics cuts from ESD track cuts\r
1377       if( vPt < fPtMin || vPt > fPtMax)      continue;\r
1378       if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1379       \r
1380       // Extra DCA cuts (for systematic studies [!= -1])\r
1381       if( fDCAxyCut != -1 && fDCAzCut != -1){\r
1382         if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){\r
1383           continue;  // 2D cut\r
1384         }\r
1385       }\r
1386       \r
1387       // Extra TPC cuts (for systematic studies [!= -1])\r
1388       if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){\r
1389         continue;\r
1390       }\r
1391       if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){\r
1392         continue;\r
1393       }\r
1394 \r
1395       //Exclude resonances\r
1396       if(fExcludeResonancesInMC) {\r
1397         \r
1398         Bool_t kExcludeParticle = kFALSE;\r
1399 \r
1400         Int_t label = TMath::Abs(aodTrack->GetLabel());\r
1401         AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);\r
1402       \r
1403         if (AODmcTrack){ \r
1404           //if (AODmcTrack->IsPhysicalPrimary()){\r
1405           \r
1406           Int_t gMotherIndex = AODmcTrack->GetMother();\r
1407           if(gMotherIndex != -1) {\r
1408             AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));\r
1409             if(motherTrack) {\r
1410               Int_t pdgCodeOfMother = motherTrack->GetPdgCode();\r
1411               if(pdgCodeOfMother == 113  // rho0\r
1412                  || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
1413                  // || pdgCodeOfMother == 221  // eta\r
1414                  // || pdgCodeOfMother == 331  // eta'\r
1415                  // || pdgCodeOfMother == 223  // omega\r
1416                  // || pdgCodeOfMother == 333  // phi\r
1417                  || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0\r
1418                  // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*\r
1419                  // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*\r
1420                  || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r
1421                  || pdgCodeOfMother == 111  // pi0 Dalitz\r
1422                  ) {\r
1423                 kExcludeParticle = kTRUE;\r
1424               }\r
1425             }\r
1426           }\r
1427         }       \r
1428         //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1429         if(kExcludeParticle) continue;\r
1430       }\r
1431       \r
1432       // fill QA histograms\r
1433       fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());\r
1434       fHistDCA->Fill(dcaZ,dcaXY);\r
1435       fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality);\r
1436       fHistPt->Fill(vPt,gCentrality);\r
1437       fHistEta->Fill(vEta,gCentrality);\r
1438       fHistRapidity->Fill(vY,gCentrality);\r
1439       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1440       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1441       fHistPhi->Fill(vPhi,gCentrality);\r
1442       if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);                  \r
1443       else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1444       \r
1445       //=======================================correction\r
1446       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1447       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
1448       \r
1449       // add the track to the TObjArray\r
1450       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));  \r
1451     }//track loop\r
1452   }//MCAODrec\r
1453   //==============================================================================================================\r
1454 \r
1455   else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD\r
1456 \r
1457     AliESDtrack *trackTPC   = NULL;\r
1458     AliMCParticle *mcTrack   = NULL;\r
1459 \r
1460     AliMCEvent*  mcEvent     = NULL;\r
1461     \r
1462     // for MC ESDs use also MC information\r
1463     if(gAnalysisLevel == "MCESD"){\r
1464       mcEvent = MCEvent(); \r
1465       if (!mcEvent) {\r
1466         AliError("mcEvent not available");\r
1467         return tracksAccepted;\r
1468       }\r
1469     }\r
1470     \r
1471     // Loop over tracks in event\r
1472     for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) {\r
1473       AliESDtrack* track = dynamic_cast<AliESDtrack *>(event->GetTrack(iTracks));\r
1474       if (!track) {\r
1475         AliError(Form("Could not receive track %d", iTracks));\r
1476         continue;\r
1477       } \r
1478 \r
1479       // for MC ESDs use also MC information --> MC track not used further???\r
1480       if(gAnalysisLevel == "MCESD"){\r
1481         Int_t label = TMath::Abs(track->GetLabel());\r
1482         if(label > mcEvent->GetNumberOfTracks()) continue;\r
1483         if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
1484         \r
1485         mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
1486         if(!mcTrack) continue;\r
1487       }\r
1488 \r
1489       // take only TPC only tracks\r
1490       trackTPC   = new AliESDtrack();\r
1491       if(!track->FillTPCOnlyTrack(*trackTPC)) continue;\r
1492       \r
1493       //ESD track cuts\r
1494       if(fESDtrackCuts) \r
1495         if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;\r
1496       \r
1497       // fill QA histograms\r
1498       Float_t b[2];\r
1499       Float_t bCov[3];\r
1500       trackTPC->GetImpactParameters(b,bCov);\r
1501       if (bCov[0]<=0 || bCov[2]<=0) {\r
1502         AliDebug(1, "Estimated b resolution lower or equal zero!");\r
1503         bCov[0]=0; bCov[2]=0;\r
1504       }\r
1505       \r
1506       Int_t nClustersTPC = -1;\r
1507       nClustersTPC = trackTPC->GetTPCNclsIter1();   // TPC standalone\r
1508       //nClustersTPC = track->GetTPCclusters(0);   // global track\r
1509       Float_t chi2PerClusterTPC = -1;\r
1510       if (nClustersTPC!=0) {\r
1511         chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
1512         //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
1513       }\r
1514       \r
1515       //===========================PID===============================//             \r
1516       if(fUsePID) {\r
1517         Double_t prob[AliPID::kSPECIES]={0.};\r
1518         Double_t probTPC[AliPID::kSPECIES]={0.};\r
1519         Double_t probTOF[AliPID::kSPECIES]={0.};\r
1520         Double_t probTPCTOF[AliPID::kSPECIES]={0.};\r
1521         \r
1522         Double_t nSigma = 0.;\r
1523         UInt_t detUsedTPC = 0;\r
1524         UInt_t detUsedTOF = 0;\r
1525         UInt_t detUsedTPCTOF = 0;\r
1526         \r
1527         //Decide what detector configuration we want to use\r
1528         switch(fPidDetectorConfig) {\r
1529         case kTPCpid:\r
1530           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);\r
1531           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest));\r
1532           detUsedTPC = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPC);\r
1533           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1534             prob[iSpecies] = probTPC[iSpecies];\r
1535           break;\r
1536         case kTOFpid:\r
1537           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);\r
1538           nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest));\r
1539           detUsedTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTOF);\r
1540           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1541             prob[iSpecies] = probTOF[iSpecies];\r
1542           break;\r
1543         case kTPCTOF:\r
1544           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);\r
1545           detUsedTPCTOF = fPIDCombined->ComputeProbabilities(track, fPIDResponse, probTPCTOF);\r
1546           for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)\r
1547             prob[iSpecies] = probTPCTOF[iSpecies];\r
1548           break;\r
1549         default:\r
1550           break;\r
1551         }//end switch: define detector mask\r
1552         \r
1553         //Filling the PID QA\r
1554         Double_t tofTime = -999., length = 999., tof = -999.;\r
1555         Double_t c = TMath::C()*1.E-9;// m/ns\r
1556         Double_t beta = -999.;\r
1557         Double_t  nSigmaTOFForParticleOfInterest = -999.;\r
1558         if ( (track->IsOn(AliESDtrack::kTOFin)) &&\r
1559              (track->IsOn(AliESDtrack::kTIME))  ) { \r
1560           tofTime = track->GetTOFsignal();//in ps\r
1561           length = track->GetIntegratedLength();\r
1562           tof = tofTime*1E-3; // ns     \r
1563           \r
1564           if (tof <= 0) {\r
1565             //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");\r
1566             continue;\r
1567           }\r
1568           if (length <= 0){\r
1569             //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");\r
1570             continue;\r
1571           }\r
1572           \r
1573           length = length*0.01; // in meters\r
1574           tof = tof*c;\r
1575           beta = length/tof;\r
1576           \r
1577           nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)fParticleOfInterest);\r
1578           fHistBetavsPTOFbeforePID ->Fill(track->P()*track->Charge(),beta);\r
1579           fHistProbTOFvsPtbeforePID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1580           fHistNSigmaTOFvsPtbeforePID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1581         }//TOF signal \r
1582         \r
1583         \r
1584         Double_t  nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)fParticleOfInterest);\r
1585         fHistdEdxVsPTPCbeforePID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1586         fHistProbTPCvsPtbeforePID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1587         fHistNSigmaTPCvsPtbeforePID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1588         fHistProbTPCTOFvsPtbeforePID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1589         //end of QA-before pid\r
1590         \r
1591         if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {\r
1592           //Make the decision based on the n-sigma\r
1593           if(fUsePIDnSigma) {\r
1594             if(nSigma > fPIDNSigma) continue;}\r
1595           \r
1596           //Make the decision based on the bayesian\r
1597           else if(fUsePIDPropabilities) {\r
1598             if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;\r
1599             if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;      \r
1600           }\r
1601           \r
1602           //Fill QA after the PID\r
1603           fHistBetavsPTOFafterPID ->Fill(track->P()*track->Charge(),beta);\r
1604           fHistProbTOFvsPtafterPID ->Fill(track->Pt(),probTOF[fParticleOfInterest]);\r
1605           fHistNSigmaTOFvsPtafterPID ->Fill(track->Pt(),nSigmaTOFForParticleOfInterest);\r
1606           \r
1607           fHistdEdxVsPTPCafterPID -> Fill(track->P()*track->Charge(),track->GetTPCsignal());\r
1608           fHistProbTPCvsPtafterPID -> Fill(track->Pt(),probTPC[fParticleOfInterest]); \r
1609           fHistProbTPCTOFvsPtafterPID -> Fill(track->Pt(),probTPCTOF[fParticleOfInterest]);\r
1610           fHistNSigmaTPCvsPtafterPID -> Fill(track->Pt(),nSigmaTPCForParticleOfInterest); \r
1611         }\r
1612       }\r
1613       //===========================PID===============================//\r
1614       vCharge = trackTPC->Charge();\r
1615       vY      = trackTPC->Y();\r
1616       vEta    = trackTPC->Eta();\r
1617       vPhi    = trackTPC->Phi();// * TMath::RadToDeg();\r
1618       vPt     = trackTPC->Pt();\r
1619       fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);\r
1620       fHistDCA->Fill(b[1],b[0]);\r
1621       fHistChi2->Fill(chi2PerClusterTPC,gCentrality);\r
1622       fHistPt->Fill(vPt,gCentrality);\r
1623       fHistEta->Fill(vEta,gCentrality);\r
1624       fHistPhi->Fill(vPhi,gCentrality);\r
1625       if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
1626       else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1627       fHistRapidity->Fill(vY,gCentrality);\r
1628       if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1629       else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1630       \r
1631       //=======================================correction\r
1632       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1633       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
1634 \r
1635       // add the track to the TObjArray\r
1636       tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction));   \r
1637 \r
1638       delete trackTPC;\r
1639     }//track loop\r
1640   }// ESD analysis\r
1641 \r
1642   else if(gAnalysisLevel == "MC"){\r
1643     if(!event) {\r
1644       AliError("mcEvent not available");\r
1645       return 0x0;\r
1646     }\r
1647 \r
1648     AliMCEvent *gMCEvent = dynamic_cast<AliMCEvent*>(event);\r
1649     if(gMCEvent) {\r
1650       // Loop over tracks in event\r
1651       for (Int_t iTracks = 0; iTracks < gMCEvent->GetNumberOfPrimaries(); iTracks++) {\r
1652         AliMCParticle* track = dynamic_cast<AliMCParticle *>(gMCEvent->GetTrack(iTracks));\r
1653         if (!track) {\r
1654           AliError(Form("Could not receive particle %d", iTracks));\r
1655           continue;\r
1656         }\r
1657         \r
1658         //exclude non stable particles\r
1659         if(!(gMCEvent->IsPhysicalPrimary(iTracks))) continue;\r
1660         \r
1661         vCharge = track->Charge();\r
1662         vEta    = track->Eta();\r
1663         vPt     = track->Pt();\r
1664         vY      = track->Y();\r
1665         \r
1666         if( vPt < fPtMin || vPt > fPtMax)      \r
1667           continue;\r
1668         if (!fUsePID) {\r
1669           if( vEta < fEtaMin || vEta > fEtaMax)  continue;\r
1670         }\r
1671         else if (fUsePID){\r
1672           if( vY < fEtaMin || vY > fEtaMax)  continue;\r
1673         }\r
1674 \r
1675         // Remove neutral tracks\r
1676         if( vCharge == 0 ) continue;\r
1677         \r
1678         //analyze one set of particles\r
1679         if(fUseMCPdgCode) {\r
1680           TParticle *particle = track->Particle();\r
1681           if(!particle) continue;\r
1682           \r
1683           Int_t gPdgCode = particle->GetPdgCode();\r
1684           if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode)) \r
1685             continue;\r
1686         }\r
1687         \r
1688         //Use the acceptance parameterization\r
1689         if(fAcceptanceParameterization) {\r
1690           Double_t gRandomNumber = gRandom->Rndm();\r
1691           if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt())) \r
1692             continue;\r
1693         }\r
1694         \r
1695         //Exclude resonances\r
1696         if(fExcludeResonancesInMC) {\r
1697           TParticle *particle = track->Particle();\r
1698           if(!particle) continue;\r
1699           \r
1700           Bool_t kExcludeParticle = kFALSE;\r
1701           Int_t gMotherIndex = particle->GetFirstMother();\r
1702           if(gMotherIndex != -1) {\r
1703             AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(event->GetTrack(gMotherIndex));\r
1704             if(motherTrack) {\r
1705               TParticle *motherParticle = motherTrack->Particle();\r
1706               if(motherParticle) {\r
1707                 Int_t pdgCodeOfMother = motherParticle->GetPdgCode();\r
1708                 //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {\r
1709                 if(pdgCodeOfMother == 113  // rho0\r
1710                    || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+\r
1711                    // || pdgCodeOfMother == 221  // eta\r
1712                    // || pdgCodeOfMother == 331  // eta'\r
1713                    // || pdgCodeOfMother == 223  // omega\r
1714                    // || pdgCodeOfMother == 333  // phi\r
1715                    || pdgCodeOfMother == 311  || pdgCodeOfMother == -311 // K0\r
1716                    // || pdgCodeOfMother == 313  || pdgCodeOfMother == -313 // K0*\r
1717                    // || pdgCodeOfMother == 323  || pdgCodeOfMother == -323 // K+*\r
1718                    || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda\r
1719                    || pdgCodeOfMother == 111  // pi0 Dalitz\r
1720                    ) {\r
1721                   kExcludeParticle = kTRUE;\r
1722                 }\r
1723               }\r
1724             }\r
1725           }\r
1726           \r
1727           //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi\r
1728           if(kExcludeParticle) continue;\r
1729         }\r
1730         \r
1731         vPhi    = track->Phi();\r
1732         //Printf("phi (before): %lf",vPhi);\r
1733         \r
1734         fHistPt->Fill(vPt,gCentrality);\r
1735         fHistEta->Fill(vEta,gCentrality);\r
1736         fHistPhi->Fill(vPhi,gCentrality);\r
1737         if(vCharge > 0)      fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality);\r
1738         else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality);\r
1739         //fHistPhi->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1740         fHistRapidity->Fill(vY,gCentrality);\r
1741         //if(vCharge > 0) fHistPhiPos->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1742         //else if(vCharge < 0) fHistPhiNeg->Fill(vPhi*TMath::RadToDeg(),gCentrality);\r
1743         if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality);\r
1744         else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality);\r
1745         \r
1746         //Flow after burner\r
1747         if(fUseFlowAfterBurner) {\r
1748           Double_t precisionPhi = 0.001;\r
1749           Int_t maxNumberOfIterations = 100;\r
1750           \r
1751           Double_t phi0 = vPhi;\r
1752           Double_t gV2 = fDifferentialV2->Eval(vPt);\r
1753           \r
1754           for (Int_t j = 0; j < maxNumberOfIterations; j++) {\r
1755             Double_t phiprev = vPhi;\r
1756             Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane*TMath::DegToRad()));\r
1757             Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane*TMath::DegToRad())); \r
1758             vPhi -= fl/fp;\r
1759             if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;\r
1760           }\r
1761           //Printf("phi (after): %lf\n",vPhi);\r
1762           Double_t vDeltaphiBefore = phi0 - gReactionPlane*TMath::DegToRad();\r
1763           if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();\r
1764           fHistPhiBefore->Fill(vDeltaphiBefore,gCentrality);\r
1765           \r
1766           Double_t vDeltaphiAfter = vPhi - gReactionPlane*TMath::DegToRad();\r
1767           if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();\r
1768           fHistPhiAfter->Fill(vDeltaphiAfter,gCentrality);\r
1769           \r
1770         }\r
1771         \r
1772         //vPhi *= TMath::RadToDeg();\r
1773         \r
1774       //=======================================correction\r
1775       Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality);  \r
1776       //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
1777 \r
1778         tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); \r
1779       } //track loop\r
1780     }//MC event object\r
1781   }//MC\r
1782   \r
1783   return tracksAccepted;  \r
1784 }\r
1785 \r
1786 //________________________________________________________________________\r
1787 TObjArray* AliAnalysisTaskBFPsi::GetShuffledTracks(TObjArray *tracks, Double_t gCentrality){\r
1788   // Clones TObjArray and returns it with tracks after shuffling the charges\r
1789 \r
1790   TObjArray* tracksShuffled = new TObjArray;\r
1791   tracksShuffled->SetOwner(kTRUE);\r
1792 \r
1793   vector<Short_t> *chargeVector = new vector<Short_t>;   //original charge of accepted tracks \r
1794 \r
1795   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)\r
1796   {\r
1797     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
1798     chargeVector->push_back(track->Charge());\r
1799   }  \r
1800  \r
1801   random_shuffle(chargeVector->begin(), chargeVector->end());\r
1802   \r
1803   for(Int_t i = 0; i < tracks->GetEntriesFast(); i++){\r
1804     AliVParticle* track = (AliVParticle*) tracks->At(i);\r
1805     //==============================correction\r
1806     Double_t correction = GetTrackbyTrackCorrectionMatrix(track->Eta(), track->Phi(),track->Pt(), chargeVector->at(i), gCentrality);\r
1807     //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality);\r
1808     tracksShuffled->Add(new AliBFBasicParticle(track->Eta(), track->Phi(), track->Pt(),chargeVector->at(i), correction));\r
1809   }\r
1810 \r
1811   delete chargeVector;\r
1812    \r
1813   return tracksShuffled;\r
1814 }\r
1815 \r
1816 \r
1817 //________________________________________________________________________\r
1818 void  AliAnalysisTaskBFPsi::FinishTaskOutput(){\r
1819   //Printf("END BF");\r
1820 \r
1821   if (!fBalance) {\r
1822     AliError("fBalance not available");\r
1823     return;\r
1824   }  \r
1825   if(fRunShuffling) {\r
1826     if (!fShuffledBalance) {\r
1827       AliError("fShuffledBalance not available");\r
1828       return;\r
1829     }\r
1830   }\r
1831 \r
1832 }\r
1833 \r
1834 //________________________________________________________________________\r
1835 void AliAnalysisTaskBFPsi::Terminate(Option_t *) {\r
1836   // Draw result to the screen\r
1837   // Called once at the end of the query\r
1838 \r
1839   // not implemented ...\r
1840 \r
1841 }\r
1842 \r