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