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