PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx

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