PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBF.cxx

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