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