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