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