PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx

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