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 = dynamic_cast<AliAODHeader*>(gAOD->GetHeader());
 674     if(!aodHeader) AliFatal("Not a standard AOD");
 675
 676     // store offline trigger bits
 677     fHistTriggerStats->Fill(aodHeader->GetOfflineTrigger());
 678
 679     if(fUseCentrality) {
 680       fCentrality = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 681     }
 682
 683     //event selection done in AliAnalysisTaskSE::Exec() --> this is not used
 684     fHistEventStats->Fill(1,fCentrality); //all events
 685     Bool_t isSelected = kTRUE;
 686     if(fUseOfflineTrigger)
 687       isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 688     if(isSelected) {
 689       fHistEventStats->Fill(2,fCentrality); //triggered events
 690
 691       //Centrality stuff (centrality in AOD header)
 692       if(fUseCentrality) {
 693         //fCentrality = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
 694         // in OLD AODs (i.e. AOD049) fCentrality can be == 0
 695         if(fCentrality == 0)
 696           return;
 697
 698         // QA for centrality estimators
 699         fHistCentStats->Fill(0.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0M"));
 700         fHistCentStats->Fill(1.,aodHeader->GetCentralityP()->GetCentralityPercentile("FMD"));
 701         fHistCentStats->Fill(2.,aodHeader->GetCentralityP()->GetCentralityPercentile("TRK"));
 702         fHistCentStats->Fill(3.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKL"));
 703         fHistCentStats->Fill(4.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL0"));
 704         fHistCentStats->Fill(5.,aodHeader->GetCentralityP()->GetCentralityPercentile("CL1"));
 705         fHistCentStats->Fill(6.,aodHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
 706         fHistCentStats->Fill(7.,aodHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
 707         fHistCentStats->Fill(8.,aodHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
 708
 709         // take only events inside centrality class
 710         if((fCentrality < fCentralityPercentileMin) || (fCentrality > fCentralityPercentileMax))
 711           return;
 712
 713         // centrality QA (V0M)
 714         fHistV0M->Fill(gAOD->GetVZEROData()->GetMTotV0A(), gAOD->GetVZEROData()->GetMTotV0C());
 715
 716         // centrality QA (reference tracks)
 717         fHistRefTracks->Fill(0.,aodHeader->GetRefMultiplicity());
 718         fHistRefTracks->Fill(1.,aodHeader->GetRefMultiplicityPos());
 719         fHistRefTracks->Fill(2.,aodHeader->GetRefMultiplicityNeg());
 720         fHistRefTracks->Fill(3.,aodHeader->GetTPConlyRefMultiplicity());
 721         fHistRefTracks->Fill(4.,aodHeader->GetNumberOfITSClusters(0));
 722         fHistRefTracks->Fill(5.,aodHeader->GetNumberOfITSClusters(1));
 723         fHistRefTracks->Fill(6.,aodHeader->GetNumberOfITSClusters(2));
 724         fHistRefTracks->Fill(7.,aodHeader->GetNumberOfITSClusters(3));
 725         fHistRefTracks->Fill(8.,aodHeader->GetNumberOfITSClusters(4));
 726       }
 727
 728       const AliAODVertex *vertex = gAOD->GetPrimaryVertex();
 729
 730       if(vertex) {
 731         Double32_t fCov[6];
 732         vertex->GetCovarianceMatrix(fCov);
 733
 734         if(vertex->GetNContributors() > 0) {
 735           if(fCov[5] != 0) {
 736             fHistEventStats->Fill(3,fCentrality); //events with a proper vertex
 737             if(TMath::Abs(vertex->GetX()) < fVxMax) {
 738               if(TMath::Abs(vertex->GetY()) < fVyMax) {
 739                 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
 740                   fHistEventStats->Fill(4,fCentrality); //analyzed events
 741                   fHistVx->Fill(vertex->GetX());
 742                   fHistVy->Fill(vertex->GetY());
 743                   fHistVz->Fill(vertex->GetZ());
 744
 745                   //===========================================//
 746                   TExMap *trackMap = new TExMap();
 747                   for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {
 748                     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));
 749                     if (!aodTrack) {
 750                       AliError(Form("ERROR: Could not receive track %d", iTracks));
 751                       continue;
 752                     }
 753                     Int_t gID = aodTrack->GetID();
 754                     //if (!aodTrack->TestFilterBit(fAODtrackCutBit)) trackMap->Add(gID, iTracks);
 755                     if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, iTracks);
 756                   }
 757                   AliAODTrack* newAodTrack;
 758                   //===========================================//
 759
 760                   //Printf("There are %d tracks in this event", gAOD->GetNumberOfTracks());
 761                   for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {
 762                     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));
 763                     if (!aodTrack) {
 764                       AliError(Form("ERROR: Could not receive track %d", iTracks));
 765                       continue;
 766                     }
 767
 768                     // AOD track cuts
 769                     // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
 770                     //===========================================//
 771                     // take only TPC only tracks
 772                     fHistTrackStats->Fill(aodTrack->GetFilterMap());
 773                     if(!aodTrack->TestFilterBit(fAODtrackCutBit)) continue;
 774
 775                     Int_t gID = aodTrack->GetID();
 776                     newAodTrack = gID >= 0 ? aodTrack : dynamic_cast<AliAODTrack*>(gAOD->GetTrack(trackMap->GetValue(-1-gID)));
 777                     if(!newAodTrack) AliFatal("Not a standard AOD");
 778                     //Printf("Label: %d - Pt: %lf (old) - %d - Pt: %lf(new)",gID,aodTrack->Pt(), newAodTrack->GetID(), newAodTrack->Pt());
 779                     //===========================================//
 780
 781                     //fHistTrackStats->Fill(aodTrack->GetFilterMap());
 782                     //if(!aodTrack->TestFilterBit(fAODtrackCutBit)) continue;
 783
 784                     vCharge = aodTrack->Charge();
 785                     vY      = aodTrack->Y();
 786                     vEta    = aodTrack->Eta();
 787                     vPhi    = aodTrack->Phi() * TMath::RadToDeg();
 788                     vE      = aodTrack->E();
 789                     vPt     = aodTrack->Pt();
 790                     aodTrack->PxPyPz(vP);
 791
 792                     Float_t dcaXY = aodTrack->DCA();      // this is the DCA from global track (not exactly what is cut on)
 793                     Float_t dcaZ  = aodTrack->ZAtDCA();   // this is the DCA from global track (not exactly what is cut on)
 794
 795
 796                     // Kinematics cuts from ESD track cuts
 797                     if( vPt < fPtMin || vPt > fPtMax)      continue;
 798
 799                     if (!fUsePID) {
 800                       if( vEta < fEtaMin || vEta > fEtaMax)  continue;
 801                     }
 802
 803                     else if (fUsePID){
 804                       if( vY < fEtaMin || vY > fEtaMax)  continue;
 805                     }
 806
 807                     // Extra DCA cuts (for systematic studies [!= -1])
 808                     if( fDCAxyCut != -1 && fDCAzCut != -1){
 809                       if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){
 810                         continue;  // 2D cut
 811                       }
 812                     }
 813
 814                     // Extra TPC cuts (for systematic studies [!= -1])
 815                     if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){
 816                       continue;
 817                     }
 818                     if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){
 819                       continue;
 820                     }
 821
 822                    //===============================================PID==================================//
 823                     if(fUsePID) {
 824                       Double_t prob[AliPID::kSPECIES]={0.};
 825                       Double_t probTPC[AliPID::kSPECIES]={0.};
 826                       Double_t probTOF[AliPID::kSPECIES]={0.};
 827                       Double_t probTPCTOF[AliPID::kSPECIES]={0.};
 828
 829                       Double_t nSigma = 0.;
 830                       UInt_t detUsedTPC = 0;
 831                       UInt_t detUsedTOF = 0;
 832                       UInt_t detUsedTPCTOF = 0;
 833
 834                       //Decide what detector configuration we want to use
 835                       switch(fPidDetectorConfig) {
 836                       case kTPCpid:
 837                         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
 838                         nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)fParticleOfInterest));
 839                         //detUsedTPC = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTPC);
 840                         detUsedTPC = (AliPIDResponse::kDetTPC);
 841                         for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
 842                           prob[iSpecies] = probTPC[iSpecies];
 843                         break;
 844                       case kTOFpid:
 845                         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
 846                         nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,(AliPID::EParticleType)fParticleOfInterest));
 847                         //detUsedTOF = fPIDCombined->ComputeProbabilities(aodTrack, fPIDResponse, probTOF);
 848                         detUsedTPC = (AliPIDResponse::kDetTPC);
 849                         for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
 850                           prob[iSpecies] = probTOF[iSpecies];
 851                         break;
 852                       case kTPCTOF:
 853                         fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
 854                         //detUsedTPCTOF = fPIDCombined->ComputeProbabilities(newAodTrack, fPIDResponse, probTPCTOF);
 855                         detUsedTPC = (AliPIDResponse::kDetTPC);
 856                         for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
 857                           prob[iSpecies] = probTPCTOF[iSpecies];
 858                         break;
 859                       default:
 860                         break;
 861                       }//end switch: define detector mask
 862
 863                       //Filling the PID QA
 864                       Double_t tofTime = -999., tof = -999.; //length = 999., tof = -999.;
 865                       //Double_t c = TMath::C()*1.E-9;// m/ns
 866                       //Double_t beta = -999.;
 867                       Double_t  nSigmaTOFForParticleOfInterest = -999.;
 868                       if ( (newAodTrack->IsOn(AliAODTrack::kTOFin)) &&
 869                            (newAodTrack->IsOn(AliAODTrack::kTIME))  ) {
 870                         tofTime = newAodTrack->GetTOFsignal();//in ps
 871                         //length = newAodTrack->GetIntegratedLength();
 872                         tof = tofTime*1E-3; // ns
 873
 874                         if (tof <= 0) {
 875                           //Printf("WARNING: track with negative TOF time found! Skipping this track for PID checks\n");
 876                           continue;
 877                         }
 878                         //if (length <= 0){
 879                           //printf("WARNING: track with negative length found!Skipping this track for PID checks\n");
 880                           //continue;
 881                         //}
 882
 883                         //length = length*0.01; // in meters
 884                         //tof = tof*c;
 885                         //beta = length/tof;
 886
 887                         nSigmaTOFForParticleOfInterest = fPIDResponse->NumberOfSigmasTOF(newAodTrack,(AliPID::EParticleType)fParticleOfInterest);
 888                         //fHistBetavsPTOFbeforePID ->Fill(aodTrack->P()*aodTrack->Charge(),beta);
 889                         fHistProbTOFvsPtbeforePID ->Fill(newAodTrack->Pt(),probTOF[fParticleOfInterest]);
 890                         fHistNSigmaTOFvsPtbeforePID ->Fill(newAodTrack->Pt(),nSigmaTOFForParticleOfInterest);
 891                       }//TOF signal
 892
 893
 894                       Double_t  nSigmaTPCForParticleOfInterest = fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)fParticleOfInterest);
 895                       fHistdEdxVsPTPCbeforePID -> Fill(newAodTrack->P()*newAodTrack->Charge(),newAodTrack->GetTPCsignal());
 896                       fHistProbTPCvsPtbeforePID -> Fill(newAodTrack->Pt(),probTPC[fParticleOfInterest]);
 897                       fHistNSigmaTPCvsPtbeforePID -> Fill(newAodTrack->Pt(),nSigmaTPCForParticleOfInterest);
 898                       fHistProbTPCTOFvsPtbeforePID -> Fill(newAodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
 899                       //end of QA-before pid
 900
 901                       if ((detUsedTPC != 0)||(detUsedTOF != 0)||(detUsedTPCTOF != 0)) {
 902                         //Make the decision based on the n-sigma
 903                         if(fUsePIDnSigma) {
 904                           if(nSigma > fPIDNSigma) continue;}
 905
 906                         //Make the decision based on the bayesian
 907                         else if(fUsePIDPropabilities) {
 908                           if(fParticleOfInterest != TMath::LocMax(AliPID::kSPECIES,prob)) continue;
 909                           if (prob[fParticleOfInterest] < fMinAcceptedPIDProbability) continue;
 910                           }
 911
 912                         //Fill QA after the PID
 913                         //fHistBetavsPTOFafterPID ->Fill(newAodTrack->P()*newAodTrack->Charge(),beta);
 914                         fHistProbTOFvsPtafterPID ->Fill(newAodTrack->Pt(),probTOF[fParticleOfInterest]);
 915                         fHistNSigmaTOFvsPtafterPID ->Fill(newAodTrack->Pt(),nSigmaTOFForParticleOfInterest);
 916
 917                         fHistdEdxVsPTPCafterPID -> Fill(newAodTrack->P()*newAodTrack->Charge(),newAodTrack->GetTPCsignal());
 918                         fHistProbTPCvsPtafterPID -> Fill(newAodTrack->Pt(),probTPC[fParticleOfInterest]);
 919                         fHistProbTPCTOFvsPtafterPID -> Fill(newAodTrack->Pt(),probTPCTOF[fParticleOfInterest]);
 920                         fHistNSigmaTPCvsPtafterPID -> Fill(newAodTrack->Pt(),nSigmaTPCForParticleOfInterest);
 921                       }
 922
 923                       PostData(4, fHistListPIDQA);
 924                     }
 925
 926                     //=========================================================PID=================================================================//
 927
 928                     // fill QA histograms
 929                     fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
 930                     fHistDCA->Fill(dcaZ,dcaXY);
 931                     fHistChi2->Fill(aodTrack->Chi2perNDF());
 932                     fHistPt->Fill(vPt);
 933                     fHistEta->Fill(vEta);
 934                     fHistPhi->Fill(vPhi);
 935                     fHistRapidity->Fill(vY);
 936                     if(vCharge > 0) fHistPhiPos->Fill(vPhi);
 937                     else if(vCharge < 0) fHistPhiNeg->Fill(vPhi);
 938
 939                     // fill charge vector
 940                     chargeVector[0]->push_back(vCharge);
 941                     chargeVector[1]->push_back(vY);
 942                     chargeVector[2]->push_back(vEta);
 943                     chargeVector[3]->push_back(vPhi);
 944                     chargeVector[4]->push_back(vP[0]);
 945                     chargeVector[5]->push_back(vP[1]);
 946                     chargeVector[6]->push_back(vP[2]);
 947                     chargeVector[7]->push_back(vPt);
 948                     chargeVector[8]->push_back(vE);
 949
 950                     if(fRunShuffling) {
 951                       chargeVectorShuffle[0]->push_back(vCharge);
 952                       chargeVectorShuffle[1]->push_back(vY);
 953                       chargeVectorShuffle[2]->push_back(vEta);
 954                       chargeVectorShuffle[3]->push_back(vPhi);
 955                       chargeVectorShuffle[4]->push_back(vP[0]);
 956                       chargeVectorShuffle[5]->push_back(vP[1]);
 957                       chargeVectorShuffle[6]->push_back(vP[2]);
 958                       chargeVectorShuffle[7]->push_back(vPt);
 959                       chargeVectorShuffle[8]->push_back(vE);
 960                     }
 961
 962                     gNumberOfAcceptedTracks += 1;
 963
 964                   } //track loop
 965                 }//Vz cut
 966               }//Vy cut
 967             }//Vx cut
 968           }//proper vertex resolution
 969         }//proper number of contributors
 970       }//vertex object valid
 971     }//triggered event
 972   }//AOD analysis
 973
 974   //MC-ESD analysis
 975   if(gAnalysisLevel == "MCESD") {
 976     AliMCEvent*  mcEvent = MCEvent();
 977     if (!mcEvent) {
 978       AliError("ERROR: mcEvent not available");
 979       return;
 980     }
 981
 982     AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
 983     if (!gESD) {
 984       AliError("ERROR: gESD not available");
 985       return;
 986     }
 987
 988     // store offline trigger bits
 989     fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
 990
 991     AliCentrality *centrality = 0x0;
 992     if(fUseCentrality) {
 993         centrality = gESD->GetCentrality();
 994         fCentrality = centrality->GetCentralityPercentile(fCentralityEstimator.Data());
 995     }
 996
 997     // event selection done in AliAnalysisTaskSE::Exec() --> this is not used
 998     fHistEventStats->Fill(1,fCentrality); //all events
 999     Bool_t isSelected = kTRUE;
1000     if(fUseOfflineTrigger)
1001       isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1002     if(isSelected) {
1003       fHistEventStats->Fill(2,fCentrality); //triggered events
1004
1005       if(fUseCentrality) {
1006         //Centrality stuff
1007         // take only events inside centrality class
1008         if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,
1009                                                  fCentralityPercentileMax,
1010                                                  fCentralityEstimator.Data()))
1011           return;
1012         // centrality QA (V0M)
1013         fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());
1014       }
1015
1016       const AliESDVertex *vertex = gESD->GetPrimaryVertex();
1017       if(vertex) {
1018         if(vertex->GetNContributors() > 0) {
1019           if(vertex->GetZRes() != 0) {
1020             fHistEventStats->Fill(3,fCentrality); //events with a proper vertex
1021             if(TMath::Abs(vertex->GetX()) < fVxMax) {
1022               if(TMath::Abs(vertex->GetY()) < fVyMax) {
1023                 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
1024                   fHistEventStats->Fill(4,fCentrality); //analayzed events
1025                   fHistVx->Fill(vertex->GetX());
1026                   fHistVy->Fill(vertex->GetY());
1027                   fHistVz->Fill(vertex->GetZ());
1028
1029                   //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());
1030                   for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {
1031                     AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));
1032                     if (!track) {
1033                       AliError(Form("ERROR: Could not receive track %d", iTracks));
1034                       continue;
1035                     }
1036
1037                     Int_t label = TMath::Abs(track->GetLabel());
1038                     if(label > mcEvent->GetNumberOfTracks()) continue;
1039                     if(label > mcEvent->GetNumberOfPrimaries()) continue;
1040
1041                     AliMCParticle* mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));
1042                     if(!mcTrack) continue;
1043
1044                     // take only TPC only tracks
1045                     trackTPC   = new AliESDtrack();
1046                     if(!track->FillTPCOnlyTrack(*trackTPC)) continue;
1047
1048                     //ESD track cuts
1049                     if(fESDtrackCuts)
1050                       if(!fESDtrackCuts->AcceptTrack(trackTPC)) continue;
1051
1052                     // fill QA histograms
1053                     Float_t b[2];
1054                     Float_t bCov[3];
1055                     trackTPC->GetImpactParameters(b,bCov);
1056                     if (bCov[0]<=0 || bCov[2]<=0) {
1057                       AliDebug(1, "Estimated b resolution lower or equal zero!");
1058                       bCov[0]=0; bCov[2]=0;
1059                     }
1060
1061                     Int_t nClustersTPC = -1;
1062                     nClustersTPC = trackTPC->GetTPCNclsIter1();   // TPC standalone
1063                     //nClustersTPC = track->GetTPCclusters(0);   // global track
1064                     Float_t chi2PerClusterTPC = -1;
1065                     if (nClustersTPC!=0) {
1066                       chi2PerClusterTPC = trackTPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone
1067                       //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track
1068                     }
1069
1070                     vCharge = trackTPC->Charge();
1071                     vY      = trackTPC->Y();
1072                     vEta    = trackTPC->Eta();
1073                     vPhi    = trackTPC->Phi() * TMath::RadToDeg();
1074                     vE      = trackTPC->E();
1075                     vPt     = trackTPC->Pt();
1076                     trackTPC->PxPyPz(vP);
1077
1078                     fHistClus->Fill(trackTPC->GetITSclusters(0),nClustersTPC);
1079                     fHistDCA->Fill(b[1],b[0]);
1080                     fHistChi2->Fill(chi2PerClusterTPC);
1081                     fHistPt->Fill(vPt);
1082                     fHistEta->Fill(vEta);
1083                     fHistPhi->Fill(vPhi);
1084                     fHistRapidity->Fill(vY);
1085                     if(vCharge > 0) fHistPhiPos->Fill(vPhi);
1086                     else if(vCharge < 0) fHistPhiNeg->Fill(vPhi);
1087
1088                     // fill charge vector
1089                     chargeVector[0]->push_back(vCharge);
1090                     chargeVector[1]->push_back(vY);
1091                     chargeVector[2]->push_back(vEta);
1092                     chargeVector[3]->push_back(vPhi);
1093                     chargeVector[4]->push_back(vP[0]);
1094                     chargeVector[5]->push_back(vP[1]);
1095                     chargeVector[6]->push_back(vP[2]);
1096                     chargeVector[7]->push_back(vPt);
1097                     chargeVector[8]->push_back(vE);
1098
1099                     if(fRunShuffling) {
1100                       chargeVectorShuffle[0]->push_back(vCharge);
1101                       chargeVectorShuffle[1]->push_back(vY);
1102                       chargeVectorShuffle[2]->push_back(vEta);
1103                       chargeVectorShuffle[3]->push_back(vPhi);
1104                       chargeVectorShuffle[4]->push_back(vP[0]);
1105                       chargeVectorShuffle[5]->push_back(vP[1]);
1106                       chargeVectorShuffle[6]->push_back(vP[2]);
1107                       chargeVectorShuffle[7]->push_back(vPt);
1108                       chargeVectorShuffle[8]->push_back(vE);
1109                     }
1110
1111                     delete trackTPC;
1112                     gNumberOfAcceptedTracks += 1;
1113
1114                   } //track loop
1115                   //cout<<"Centrality: "<<fCentrality<<" - Accepted tracks: "<<gNumberOfAcceptedTracks<<endl;
1116                 }//Vz cut
1117               }//Vy cut
1118             }//Vx cut
1119           }//proper vertex resolution
1120         }//proper number of contributors
1121       }//vertex object valid
1122     }//triggered event
1123   }//MC-ESD analysis
1124
1125   //MC analysis
1126   else if(gAnalysisLevel == "MC") {
1127     AliMCEvent*  mcEvent = MCEvent();
1128     if (!mcEvent) {
1129       AliError("ERROR: mcEvent not available");
1130       return;
1131     }
1132
1133     //fHistEventStats->Fill(1,fCentrality); //total events
1134     //fHistEventStats->Fill(2,fCentrality); //offline trigger
1135
1136     Double_t gReactionPlane = 0., gImpactParameter = 0.;
1137     if(fUseCentrality) {
1138       //Get the MC header
1139       AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
1140       if (headerH) {
1141         //Printf("=====================================================");
1142         //Printf("Reaction plane angle: %lf",headerH->ReactionPlaneAngle());
1143         //Printf("Impact parameter: %lf",headerH->ImpactParameter());
1144         //Printf("=====================================================");
1145         gReactionPlane = headerH->ReactionPlaneAngle();
1146         gImpactParameter = headerH->ImpactParameter();
1147         fCentrality = gImpactParameter;
1148       }
1149       fCentrality = gImpactParameter;
1150
1151       // take only events inside centrality class (DIDN'T CHANGE THIS UP TO NOW)
1152       if((fImpactParameterMin > gImpactParameter) || (fImpactParameterMax < gImpactParameter))
1153         return;
1154     }
1155
1156     fHistEventStats->Fill(1,fCentrality); //total events
1157     fHistEventStats->Fill(2,fCentrality); //offline trigger
1158
1159     AliGenEventHeader *header = mcEvent->GenEventHeader();
1160     if(!header) return;
1161
1162     TArrayF gVertexArray;
1163     header->PrimaryVertex(gVertexArray);
1164     //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",
1165     //gVertexArray.At(0),
1166     //gVertexArray.At(1),
1167     //gVertexArray.At(2));
1168     fHistEventStats->Fill(3,fCentrality); //events with a proper vertex
1169     if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {
1170       if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {
1171         if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {
1172           fHistEventStats->Fill(4,fCentrality); //analayzed events
1173           fHistVx->Fill(gVertexArray.At(0));
1174           fHistVy->Fill(gVertexArray.At(1));
1175           fHistVz->Fill(gVertexArray.At(2));
1176
1177           AliInfo(Form("There are %d tracks in this event", mcEvent->GetNumberOfPrimaries()));
1178           for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfPrimaries(); iTracks++) {
1179             AliMCParticle* track = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iTracks));
1180             if (!track) {
1181               AliError(Form("ERROR: Could not receive particle %d", iTracks));
1182               continue;
1183             }
1184
1185             //exclude non stable particles
1186             if(!mcEvent->IsPhysicalPrimary(iTracks)) continue;
1187
1188             vEta    = track->Eta();
1189             vPt     = track->Pt();
1190             vY      = track->Y();
1191
1192             if( vPt < fPtMin || vPt > fPtMax)
1193               continue;
1194             if (!fUsePID) {
1195               if( vEta < fEtaMin || vEta > fEtaMax)  continue;
1196             }
1197             else if (fUsePID){
1198               if( vY < fEtaMin || vY > fEtaMax)  continue;
1199             }
1200
1201             //analyze one set of particles
1202             if(fUseMCPdgCode) {
1203               TParticle *particle = track->Particle();
1204               if(!particle) continue;
1205
1206               Int_t gPdgCode = particle->GetPdgCode();
1207               if(TMath::Abs(fPDGCodeToBeAnalyzed) != TMath::Abs(gPdgCode))
1208                 continue;
1209             }
1210
1211             //Use the acceptance parameterization
1212             if(fAcceptanceParameterization) {
1213               Double_t gRandomNumber = gRandom->Rndm();
1214               if(gRandomNumber > fAcceptanceParameterization->Eval(track->Pt()))
1215                 continue;
1216             }
1217
1218             //Exclude resonances
1219             if(fExcludeResonancesInMC) {
1220               TParticle *particle = track->Particle();
1221               if(!particle) continue;
1222
1223               Bool_t kExcludeParticle = kFALSE;
1224               Int_t gMotherIndex = particle->GetFirstMother();
1225               if(gMotherIndex != -1) {
1226                 AliMCParticle* motherTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(gMotherIndex));
1227                 if(motherTrack) {
1228                   TParticle *motherParticle = motherTrack->Particle();
1229                   if(motherParticle) {
1230                     Int_t pdgCodeOfMother = motherParticle->GetPdgCode();
1231                     //if((pdgCodeOfMother == 113)||(pdgCodeOfMother == 213)||(pdgCodeOfMother == 221)||(pdgCodeOfMother == 223)||(pdgCodeOfMother == 331)||(pdgCodeOfMother == 333)) {
1232                     if(pdgCodeOfMother == 113) {
1233                       kExcludeParticle = kTRUE;
1234                     }
1235                   }
1236                 }
1237               }
1238
1239               //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi
1240               if(kExcludeParticle) continue;
1241             }
1242
1243             vCharge = track->Charge();
1244             vPhi    = track->Phi();
1245             vE      = track->E();
1246             track->PxPyPz(vP);
1247             //Printf("phi (before): %lf",vPhi);
1248
1249             fHistPt->Fill(vPt);
1250             fHistEta->Fill(vEta);
1251             fHistPhi->Fill(vPhi);
1252             fHistRapidity->Fill(vY);
1253             if(vCharge > 0) fHistPhiPos->Fill(vPhi);
1254             else if(vCharge < 0) fHistPhiNeg->Fill(vPhi);
1255
1256             //Flow after burner
1257             if(fUseFlowAfterBurner) {
1258               Double_t precisionPhi = 0.001;
1259               Int_t maxNumberOfIterations = 100;
1260
1261               Double_t phi0 = vPhi;
1262               Double_t gV2 = fDifferentialV2->Eval(vPt);
1263
1264               for (Int_t j = 0; j < maxNumberOfIterations; j++) {
1265                 Double_t phiprev = vPhi;
1266                 Double_t fl = vPhi - phi0 + gV2*TMath::Sin(2.*(vPhi - gReactionPlane));
1267                 Double_t fp = 1.0 + 2.0*gV2*TMath::Cos(2.*(vPhi - gReactionPlane));
1268                 vPhi -= fl/fp;
1269                 if (TMath::AreEqualAbs(phiprev,vPhi,precisionPhi)) break;
1270               }
1271               //Printf("phi (after): %lf\n",vPhi);
1272               Double_t vDeltaphiBefore = phi0 - gReactionPlane;
1273               if(vDeltaphiBefore < 0) vDeltaphiBefore += 2*TMath::Pi();
1274               fHistPhiBefore->Fill(vDeltaphiBefore);
1275
1276               Double_t vDeltaphiAfter = vPhi - gReactionPlane;
1277               if(vDeltaphiAfter < 0) vDeltaphiAfter += 2*TMath::Pi();
1278               fHistPhiAfter->Fill(vDeltaphiAfter);
1279             }
1280
1281             vPhi *= TMath::RadToDeg();
1282
1283             // fill charge vector
1284             chargeVector[0]->push_back(vCharge);
1285             chargeVector[1]->push_back(vY);
1286             chargeVector[2]->push_back(vEta);
1287             chargeVector[3]->push_back(vPhi);
1288             chargeVector[4]->push_back(vP[0]);
1289             chargeVector[5]->push_back(vP[1]);
1290             chargeVector[6]->push_back(vP[2]);
1291             chargeVector[7]->push_back(vPt);
1292             chargeVector[8]->push_back(vE);
1293
1294             if(fRunShuffling) {
1295               chargeVectorShuffle[0]->push_back(vCharge);
1296               chargeVectorShuffle[1]->push_back(vY);
1297               chargeVectorShuffle[2]->push_back(vEta);
1298               chargeVectorShuffle[3]->push_back(vPhi);
1299               chargeVectorShuffle[4]->push_back(vP[0]);
1300               chargeVectorShuffle[5]->push_back(vP[1]);
1301               chargeVectorShuffle[6]->push_back(vP[2]);
1302               chargeVectorShuffle[7]->push_back(vPt);
1303               chargeVectorShuffle[8]->push_back(vE);
1304             }
1305             gNumberOfAcceptedTracks += 1;
1306
1307           } //track loop
1308         }//Vz cut
1309       }//Vy cut
1310     }//Vx cut
1311   }//MC analysis
1312
1313   //multiplicity cut (used in pp)
1314   if(fUseMultiplicity) {
1315     if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))
1316       return;
1317   }
1318   fHistNumberOfAcceptedTracks->Fill(gNumberOfAcceptedTracks, fCentrality);
1319
1320   // calculate balance function
1321   if(fUseMultiplicity)
1322     fBalance->CalculateBalance(gNumberOfAcceptedTracks,chargeVector,bSign);
1323   else
1324     fBalance->CalculateBalance(fCentrality,chargeVector,bSign);
1325
1326   if(fRunShuffling) {
1327     // shuffle charges
1328     random_shuffle( chargeVectorShuffle[0]->begin(), chargeVectorShuffle[0]->end() );
1329     if(fUseMultiplicity)
1330       fShuffledBalance->CalculateBalance(gNumberOfAcceptedTracks,chargeVectorShuffle,bSign);
1331     else
1332       fShuffledBalance->CalculateBalance(fCentrality,chargeVectorShuffle,bSign);
1333   }
1334 }
1335
1336 //________________________________________________________________________
1337 void  AliAnalysisTaskBF::FinishTaskOutput(){
1338   //Printf("END BF");
1339
1340   if (!fBalance) {
1341     AliError("ERROR: fBalance not available");
1342     return;
1343   }
1344   if(fRunShuffling) {
1345     if (!fShuffledBalance) {
1346       AliError("ERROR: fShuffledBalance not available");
1347       return;
1348     }
1349   }
1350
1351 }
1352
1353 //________________________________________________________________________
1354 void AliAnalysisTaskBF::Terminate(Option_t *) {
1355   // Draw result to the screen
1356   // Called once at the end of the query
1357
1358   // not implemented ...
1359
1360 }