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