PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBF.cxx

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