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