PWGCF/Correlations/DPhi/DiHadronPID/AliAODTrackCutsDiHadronPID.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 // -----------------------------------------------------------------------
  17 //  Track Cut class for the DiHadronPID analysis.
  18 // -----------------------------------------------------------------------
  19 //  Author: Misha Veldhoen (misha.veldhoen@cern.ch)
  20
  21 #include "AliAODTrackCutsDiHadronPID.h"
  22
  23 #include <iostream>
  24 using namespace std;
  25
  26 #include "TMath.h"
  27 #include "TH1F.h"
  28 #include "TH2F.h"
  29 #include "TH3F.h"
  30 #include "TNamed.h"
  31 #include "TFormula.h"
  32 #include "TIterator.h"
  33
  34 // AOD includes.
  35 #include "AliAODTrack.h"
  36 #include "AliAODEvent.h"
  37 #include "AliAODVertex.h"
  38 #include "AliAODMCParticle.h"
  39
  40 // PID includes.
  41 #include "AliPID.h"
  42 #include "AliPIDResponse.h"
  43 #include "AliTPCPIDResponse.h"
  44
  45 #include "AliTrackDiHadronPID.h"
  46 #include "AliHistToolsDiHadronPID.h"
  47 #include "AliFunctionsDiHadronPID.h"
  48
  49 ClassImp(AliAODTrackCutsDiHadronPID);
  50
  51 // -----------------------------------------------------------------------
  52 AliAODTrackCutsDiHadronPID::AliAODTrackCutsDiHadronPID():
  53         TNamed(),
  54         fMinPt(0.),
  55         fMaxPt(10.),
  56         fFilterMask(0),
  57         fMaxEta(-999.),
  58         fMaxRapidity(-999.),
  59         fMinimumNumberOfTPCClusters(-999),
  60         fDemandedFlags(0),
  61         fMinSPDHitsForPtDeptDCAcut(0),
  62         fPtDeptDCAxyCutFormula(0x0),
  63         fDCAzCut(999.),
  64         fIsMC(kFALSE),
  65         fLowPtNSigmaTOFOnly(kFALSE),
  66         fUseNSigmaOnPIDAxes(kFALSE),
  67         fTestPt(kFALSE),
  68         fTestFilterMask(kFALSE),
  69         fTestMaxEta(kFALSE),
  70         fTestMaxRapidity(kFALSE),
  71         fTestFlags(kFALSE),
  72         fTestNumberOfTPCClusters(kFALSE),
  73         fTestSPDAny(kFALSE),
  74         fTestTOFmismatch(kFALSE),
  75         fTestPtDeptDCAcut(kFALSE),
  76         fDataTrackQAHistos(0x0),
  77         fHistAcceptedFilterBits(0x0),
  78         fRelevantBitsArray(0x0),
  79         fTOFMatchingStat(0x0),
  80         fPrimRecMCTrackQAHistos(0x0),
  81         fPrimGenMCTrackQAHistos(0x0),
  82         fSecRecMCTrackQAHistos(0x0),
  83         fSecGenMCTrackQAHistos(0x0),
  84         fNPtBins(0),
  85         fNEtaBins(32),
  86         fNPhiBins(32),
  87         fDebug(0)
  88
  89  {
  90
  91         //
  92         // Default constructor
  93         //
  94
  95         cout<<"AliAODTrackCutsDiHadronPID Default Constructor Called."<<endl;
  96         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
  97
  98         for (Int_t iHistoName = 0; iHistoName < 12; iHistoName++) {
  99
 100                 // Initialize the data histograms.
 101                 if (iHistoName < 3) {
 102                         fHistDataPt[iHistoName] = 0x0;
 103                         fHistDataPhiEtaPt[iHistoName] = 0x0;
 104                         fHistDataNTracks[iHistoName] = 0x0;
 105
 106                         // DCA histograms.
 107                         fHistDataDCAxy[iHistoName] = 0x0;
 108                         fHistDataDCAz[iHistoName] = 0x0;
 109
 110                         // Initialize PID histograms.
 111                         for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
 112                                 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
 113                                         fHistPrimRecPID[iHistoName][iSpecies][iPtClass] = 0x0;
 114                                         fHistPrimRecMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
 115                                         fHistDataPID[iHistoName][iSpecies][iPtClass] = 0x0;
 116                                         fHistTOFMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
 117                                         fHistTPCTOFMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
 118                                 }
 119                         }
 120                 }
 121
 122                 fNTracks[iHistoName] = 0;
 123
 124                 // Initialize data DCA for one sigma.
 125                 fHistDataDCAxyOneSigma[iHistoName] = 0x0;
 126
 127                 // Initialize MC histograms.
 128                 fHistPrimGenMCPt[iHistoName] = 0x0;
 129                 fHistPrimRecPtGenPt[iHistoName] = 0x0;
 130                 fHistPrimGenMCPhiEtaPt[iHistoName] = 0x0;
 131                 fHistPrimRecMCPt[iHistoName] = 0x0;
 132                 fHistPrimRecMCPhiEtaPt[iHistoName] = 0x0;
 133                 fHistPrimRecNTracks[iHistoName] = 0x0;
 134
 135                 fHistSecGenMCPt[iHistoName] = 0x0;
 136                 fHistSecGenMCPhiEtaPt[iHistoName] = 0x0;
 137                 fHistSecRecMCPt[iHistoName] = 0x0;
 138                 fHistSecRecMCPhiEtaPt[iHistoName] = 0x0;
 139
 140                 // Initialze MC DCA histograms
 141                 fHistPrimRecMCDCA[iHistoName] = 0x0;
 142                 fHistSecRecMCDCAMat[iHistoName] = 0x0;
 143                 fHistSecRecMCDCAWeak[iHistoName] = 0x0;
 144
 145                 // Initialize other stuff.
 146                 fHistRequested[iHistoName] = kFALSE;
 147                 f3DSpectraEnabeled[iHistoName] = kFALSE;
 148                 fPIDHistosEnabeled[iHistoName] = kFALSE;
 149         }
 150
 151         // TODO: It would be a bit neater to initialize all values to zero
 152         // instead of the default values...
 153         InitializeDefaultHistoNamesAndAxes();
 154
 155 }
 156
 157 // -----------------------------------------------------------------------
 158 AliAODTrackCutsDiHadronPID::AliAODTrackCutsDiHadronPID(const char* name):
 159         TNamed(name,"AOD Track Cuts"),
 160         fMinPt(0.),
 161         fMaxPt(10.),
 162         fFilterMask(0),
 163         fMaxEta(-999.),
 164         fMaxRapidity(-999.),
 165         fMinimumNumberOfTPCClusters(-999),
 166         fDemandedFlags(0),
 167         fMinSPDHitsForPtDeptDCAcut(0),
 168         fPtDeptDCAxyCutFormula(0x0),
 169         fDCAzCut(999.),
 170         fIsMC(kFALSE),
 171         fLowPtNSigmaTOFOnly(kFALSE),
 172         fUseNSigmaOnPIDAxes(kFALSE),
 173         fTestPt(kFALSE),
 174         fTestFilterMask(kFALSE),
 175         fTestMaxEta(kFALSE),
 176         fTestMaxRapidity(kFALSE),
 177         fTestFlags(kFALSE),
 178         fTestNumberOfTPCClusters(kFALSE),
 179         fTestSPDAny(kFALSE),
 180         fTestTOFmismatch(kFALSE),
 181         fTestPtDeptDCAcut(kFALSE),
 182         fDataTrackQAHistos(0x0),
 183         fHistAcceptedFilterBits(0x0),
 184         fRelevantBitsArray(0x0),
 185         fTOFMatchingStat(0x0),
 186         fPrimRecMCTrackQAHistos(0x0),
 187         fPrimGenMCTrackQAHistos(0x0),
 188         fSecRecMCTrackQAHistos(0x0),
 189         fSecGenMCTrackQAHistos(0x0),
 190         fNPtBins(0),
 191         fNEtaBins(32),
 192         fNPhiBins(32),
 193         fDebug(0)
 194
 195 {
 196
 197         //
 198         // Named constructor
 199         //
 200
 201         cout<<"AliAODTrackCutsDiHadronPID Named Constructor Called."<<endl;
 202         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 203
 204         for (Int_t iHistoName = 0; iHistoName < 12; iHistoName++) {
 205
 206                 // Initialize the data histograms.
 207                 if (iHistoName < 3) {
 208                         fHistDataPt[iHistoName] = 0x0;
 209                         fHistDataPhiEtaPt[iHistoName] = 0x0;
 210                         fHistDataNTracks[iHistoName] = 0x0;
 211
 212                         // DCA histograms.
 213                         fHistDataDCAxy[iHistoName] = 0x0;
 214                         fHistDataDCAz[iHistoName] = 0x0;
 215
 216                         // Initialize PID histograms.
 217                         for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
 218                                 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
 219                                         fHistPrimRecPID[iHistoName][iSpecies][iPtClass] = 0x0;
 220                                         fHistPrimRecMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
 221                                         fHistDataPID[iHistoName][iSpecies][iPtClass] = 0x0;
 222                                         fHistTOFMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
 223                                         fHistTPCTOFMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
 224
 225                                 }
 226                         }
 227                 }
 228
 229                 fNTracks[iHistoName] = 0;
 230
 231                 // Initialize data DCA for one sigma.
 232                 fHistDataDCAxyOneSigma[iHistoName] = 0x0;
 233
 234                 // Initialize MC histograms.
 235                 fHistPrimGenMCPt[iHistoName] = 0x0;
 236                 fHistPrimRecPtGenPt[iHistoName] = 0x0;
 237                 fHistPrimGenMCPhiEtaPt[iHistoName] = 0x0;
 238                 fHistPrimRecMCPt[iHistoName] = 0x0;
 239                 fHistPrimRecMCPhiEtaPt[iHistoName] = 0x0;
 240                 fHistPrimRecNTracks[iHistoName] = 0x0;
 241
 242                 fHistSecGenMCPt[iHistoName] = 0x0;
 243                 fHistSecGenMCPhiEtaPt[iHistoName] = 0x0;
 244                 fHistSecRecMCPt[iHistoName] = 0x0;
 245                 fHistSecRecMCPhiEtaPt[iHistoName] = 0x0;
 246
 247                 // Initialze MC DCA histograms
 248                 fHistPrimRecMCDCA[iHistoName] = 0x0;
 249                 fHistSecRecMCDCAMat[iHistoName] = 0x0;
 250                 fHistSecRecMCDCAWeak[iHistoName] = 0x0;
 251
 252                 // Initialize other stuff.
 253                 fHistRequested[iHistoName] = kFALSE;
 254                 f3DSpectraEnabeled[iHistoName] = kFALSE;
 255                 fPIDHistosEnabeled[iHistoName] = kFALSE;
 256         }
 257
 258         InitializeDefaultHistoNamesAndAxes();
 259
 260 }
 261
 262 // -----------------------------------------------------------------------
 263 void AliAODTrackCutsDiHadronPID::InitializeDefaultHistoNamesAndAxes() {
 264
 265         // Initializes the histogram name conventions and the ranges of all the histograms
 266         // by their default values. This method should only be called by the
 267         // (default) constructor.
 268
 269         // TODO: User should be able to change these standard values at initialization.
 270         // TODO: User should be able to retrieve all these variables with appropriate Getters.
 271
 272         cout<<"AliAODTrackCutsDiHadronPID - Initializing Default Histogram Names and axes..."<<endl;
 273         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 274
 275         // Setting the Pt axis for all histograms except the PID and Mismatch histograms -> Now the same for
 276         // PID histograms.
 277         Double_t ptaxis[52] = {0.500,0.525,0.550,0.575,0.600,0.625,0.650,0.675,0.700,
 278                                                         0.75,0.80,0.85,0.90,0.95,1.00,
 279                                                         1.05,1.10,1.15,1.20,1.25,1.30,1.35,1.40,1.45,1.50,1.55,1.60,1.65,1.70,
 280                                                         1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,
 281                                                         3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0};
 282         for (Int_t iPtBins = 0; iPtBins < 52; iPtBins++) {fPtAxis[iPtBins] = ptaxis[iPtBins];}
 283         fNPtBins = 51;
 284 /*
 285         Double_t ptaxis[57] = {0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.00,
 286                                                         1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,
 287                                                         3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,9.5,10.0};
 288         for (Int_t iPtBins = 0; iPtBins < 57; iPtBins++) {fPtAxis[iPtBins] = ptaxis[iPtBins];}
 289         fNPtBins = 56;
 290 */
 291
 292         // Setting the pt range of the five PID histogram pt classes.
 293         Double_t ptboundarypid[6] = {0.5,0.7,1.0,1.7,3.0,5.0};
 294         for (Int_t iPtBoundary = 0; iPtBoundary < 6; iPtBoundary++) {fPtBoundaryPID[iPtBoundary] = ptboundarypid[iPtBoundary];}
 295
 296         // Setting the number of bins in pt for the five PID histogram pt classes.
 297         Int_t nptbinspid[5] = {8,6,14,13,10};
 298         for (Int_t iPtBins = 0; iPtBins < 5; iPtBins++) {fNPtBinsPID[iPtBins] = nptbinspid[iPtBins];}
 299
 300         // Setting the TOF axes for the PID histograms.
 301         //Double_t tofsigmaapprox = 80.;
 302         Double_t toflowerbound[5][3] = {{-2000.,-6000.,-10000.},{-2000.,-4000.,-10000.},{-1000.,-2000.,-5000.},{-1000.,-1000.,-2500.},{-500.,-500.,-1000.}};
 303         Double_t tofupperbound[5][3] = {{10000.,10000.,6000.},{10000.,8000.,6000.},{6000.,6000.,6000.},{6000.,6000.,6000.},{4000.,4000.,6000.}};
 304         Int_t tofbins[5][3] = {{120,160,160},{120,120,160},{140,140,165},{140,140,170},{90,90,140}};
 305         for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
 306                 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
 307                         fTOFLowerBound[iPtClass][iSpecies] = toflowerbound[iPtClass][iSpecies];
 308                         //if (fUseNSigmaOnPIDAxes) {toflowerbound[iPtClass][iSpecies] /= tofsigmaapprox;}
 309                         fTOFUpperBound[iPtClass][iSpecies] = tofupperbound[iPtClass][iSpecies];
 310                         //if (fUseNSigmaOnPIDAxes) {tofupperbound[iPtClass][iSpecies] /= tofsigmaapprox;}
 311                         fTOFbins[iPtClass][iSpecies] = tofbins[iPtClass][iSpecies];
 312                 }
 313         }
 314
 315         // Setting the TPC axes for the PID histograms.
 316         //Double_t tpcsigmaaxpprox = 3.5;
 317         Double_t tpclowerbound[5][3] = {{-20.,-50.,-100.},{-20.,-30.,-80.},{-25.,-25.,-45.},{-25.,-25.,-45.},{-25.,-20.,-20.}};
 318         Double_t tpcupperbound[5][3] = {{60.,30.,20.},{60.,40.,20.},{50.,50.,25.},{45.,45.,25.},{25.,30.,30.}}; // Check highest pT bin boundaries for K,p
 319         Int_t tpcbins[5][3] = {{80,80,120},{80,70,100},{75,75,70},{70,70,70},{50,50,50}};
 320         for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
 321                 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
 322                         fTPCLowerBound[iPtClass][iSpecies] = tpclowerbound[iPtClass][iSpecies];
 323                         //if (fUseNSigmaOnPIDAxes) {tpclowerbound[iPtClass][iSpecies] /= tpcsigmaaxpprox;}
 324                         fTPCUpperBound[iPtClass][iSpecies] = tpcupperbound[iPtClass][iSpecies];
 325                         //if (fUseNSigmaOnPIDAxes) {tpcupperbound[iPtClass][iSpecies] /= tpcsigmaaxpprox;}
 326                         fTPCbins[iPtClass][iSpecies] = tpcbins[iPtClass][iSpecies];
 327                 }
 328         }
 329
 330         // Names for the 12 (species,charge) combinations considered.
 331         const char* histoname[12] = {"AllCharged","Pos","Neg","AllPion","PosPion","NegPion","AllKaon","PosKaon","NegKaon","AllProton","PosProton","NegProton"};
 332         const char* histolatex[12] = {"ch","ch^{+}","ch^{-}","#pi^{+/-}","#pi^{+}","#pi^{-}","K^{+/-}","K^{+}","K^{-}","p/#bar{p}","p","#bar{p}"};
 333         for (Int_t iHistoName = 0; iHistoName < 12; iHistoName++) {
 334                 fHistoName[iHistoName] = histoname[iHistoName];
 335                 fHistoLatex[iHistoName] = histolatex[iHistoName];
 336         }
 337
 338         // Names of the 3 species considered.
 339         const char* particlename[3] = {"Pion","Kaon","Proton"};
 340         for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {fParticleName[iSpecies] = particlename[iSpecies];}
 341
 342         // Names of the 5 pt classes considered for the PID histograms.
 343         const char* ptclassname[5] = {"VeryLowPt","LowPt","MedPt","MedHighPt","HighPt"};
 344         for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {fPtClassName[iPtClass] = ptclassname[iPtClass];}
 345
 346 }
 347
 348 // -----------------------------------------------------------------------
 349 AliAODTrackCutsDiHadronPID::~AliAODTrackCutsDiHadronPID() {
 350
 351         //
 352         // Destructor
 353         //
 354
 355         cout<<"AliAODTrackCutsDiHadronPID Destructor Called."<<endl;
 356         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 357
 358         if (fDataTrackQAHistos) delete fDataTrackQAHistos;
 359         fDataTrackQAHistos = 0x0;
 360         if (fPrimRecMCTrackQAHistos) delete fPrimRecMCTrackQAHistos;
 361         fPrimRecMCTrackQAHistos = 0x0;
 362         if (fPrimGenMCTrackQAHistos) delete fPrimGenMCTrackQAHistos;
 363         fPrimGenMCTrackQAHistos = 0x0;
 364         if (fSecRecMCTrackQAHistos) delete fSecRecMCTrackQAHistos;
 365         fSecRecMCTrackQAHistos = 0x0;
 366         if (fSecGenMCTrackQAHistos) delete fSecGenMCTrackQAHistos;
 367         fSecGenMCTrackQAHistos = 0x0;
 368
 369 }
 370
 371 // -----------------------------------------------------------------------
 372 Long64_t AliAODTrackCutsDiHadronPID::Merge(TCollection* list) {
 373
 374         //
 375         // Merger.
 376         //
 377
 378         cout<<"AliAODTrackCutsDiHadronPID Merger Called."<<endl;
 379         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 380
 381         Bool_t HistosOK = kTRUE;
 382
 383         if (fIsMC) {
 384                 if (!fPrimGenMCTrackQAHistos||!fSecGenMCTrackQAHistos||!fPrimRecMCTrackQAHistos||!fSecRecMCTrackQAHistos) {HistosOK = kFALSE;}
 385         } else {
 386                 if (!fDataTrackQAHistos) {HistosOK = kFALSE;}
 387         }
 388
 389         if (!HistosOK) {
 390                 cout<<"AliAODTrackCutsDiHadronPID::Merge() - Warning, current object's histograms are missing... Generating."<<endl;
 391                 CreateHistos();
 392         }
 393
 394         if (!list) {
 395                 //cout<<"No list available..."<<endl;
 396                 return 0;
 397         }
 398         if (list->IsEmpty()) {
 399                 //cout<<"List is empty..."<<endl;
 400                 return 1;
 401         }
 402
 403         //Int_t NEntries = list->GetEntries();
 404         //cout<<"Supplied list has "<<NEntries<<" entries."<<endl;
 405
 406         TIterator* iter = list->MakeIterator();
 407         TObject* obj;
 408
 409         // List of collections
 410         TList collection_fDataTrackQAHistos;
 411         TList collection_fPrimRecMCTrackQAHistos;
 412         TList collection_fPrimGenMCTrackQAHistos;
 413         TList collection_fSecRecMCTrackQAHistos;
 414         TList collection_fSecGenMCTrackQAHistos;
 415
 416         Int_t count = 0;
 417
 418         while ((obj = iter->Next())) {
 419         AliAODTrackCutsDiHadronPID* entry = dynamic_cast<AliAODTrackCutsDiHadronPID*> (obj);
 420         if (entry == 0) continue;
 421
 422         // Check if the object to be merged really has the same name! (FIXME!)
 423
 424         // Getting the lists from obj.
 425         TList* list_fDataTrackQAHistos = entry->GetListOfDataQAHistos();
 426         TList* list_fPrimRecMCTrackQAHistos = entry->GetListOfPrimRecMCTrackQAHistos();
 427         TList* list_fPrimGenMCTrackQAHistos = entry->GetListOfPrimGenMCTrackQAHistos();
 428         TList* list_fSecRecMCTrackQAHistos = entry->GetListOfSecRecMCTrackQAHistos();
 429         TList* list_fSecGenMCTrackQAHistos = entry->GetListOfSecGenMCTrackQAHistos();
 430
 431         // Adding the retrieved lists to the collection.
 432         if (list_fDataTrackQAHistos) collection_fDataTrackQAHistos.Add(list_fDataTrackQAHistos);
 433         if (list_fPrimRecMCTrackQAHistos) collection_fPrimRecMCTrackQAHistos.Add(list_fPrimRecMCTrackQAHistos);
 434         if (list_fPrimGenMCTrackQAHistos) collection_fPrimGenMCTrackQAHistos.Add(list_fPrimGenMCTrackQAHistos);
 435         if (list_fSecRecMCTrackQAHistos) collection_fSecRecMCTrackQAHistos.Add(list_fSecRecMCTrackQAHistos);
 436         if (list_fSecGenMCTrackQAHistos) collection_fSecGenMCTrackQAHistos.Add(list_fSecGenMCTrackQAHistos);
 437
 438         //cout<<"Entries intermediate list Data: "<<collection_fDataTrackQAHistos.GetEntries()<<endl;
 439         //cout<<"Entries intermediate list RecMC: "<<collection_fPrimRecMCTrackQAHistos.GetEntries()<<endl;
 440
 441         count++;
 442     }
 443
 444     // Merging. Note that we require the original list to exist.
 445     //  * Assume that if the collection happens to be empty, then nothing will happen.
 446     //  |-> This of course leads to some problems, since if the first file does not
 447     //  |   have such a list, then the merged file will have no results...
 448     //  |   IDEA: if fDataTrackQAHistos does not exist, then create a new one. (TO DO)
 449     //  * All other variables are taken from the original object.
 450     if (fDataTrackQAHistos) fDataTrackQAHistos->Merge(&collection_fDataTrackQAHistos);
 451     if (fPrimRecMCTrackQAHistos) fPrimRecMCTrackQAHistos->Merge(&collection_fPrimRecMCTrackQAHistos);
 452     if (fPrimGenMCTrackQAHistos) fPrimGenMCTrackQAHistos->Merge(&collection_fPrimGenMCTrackQAHistos);
 453     if (fSecRecMCTrackQAHistos) fSecRecMCTrackQAHistos->Merge(&collection_fSecRecMCTrackQAHistos);
 454     if (fSecGenMCTrackQAHistos) fSecGenMCTrackQAHistos->Merge(&collection_fSecGenMCTrackQAHistos);
 455
 456     delete iter;
 457
 458         return count+1;
 459
 460 }
 461
 462 // -----------------------------------------------------------------------
 463 void AliAODTrackCutsDiHadronPID::PrintCuts() const { /* NOT IMPLEMENTED */ }
 464
 465 // -----------------------------------------------------------------------
 466 TList* AliAODTrackCutsDiHadronPID::GetListOfDataQAHistos() const {
 467
 468         // Returns the list of data histograms.
 469         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 470
 471         if (fDataTrackQAHistos) {
 472                 return fDataTrackQAHistos;
 473         } else {
 474                 return 0x0;
 475         }
 476
 477 }
 478
 479 // -----------------------------------------------------------------------
 480 TList* AliAODTrackCutsDiHadronPID::GetListOfPrimRecMCTrackQAHistos() const {
 481
 482         // Returns the list of histograms of reconstructed primary tracks.
 483         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 484
 485         if (fPrimRecMCTrackQAHistos) {
 486                 return fPrimRecMCTrackQAHistos;
 487         } else {
 488                 return 0x0;
 489         }
 490
 491 }
 492
 493 // -----------------------------------------------------------------------
 494 TList* AliAODTrackCutsDiHadronPID::GetListOfPrimGenMCTrackQAHistos() const {
 495
 496         // Returns the list of histograms of generator level primary particles.
 497         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 498
 499         if (fPrimGenMCTrackQAHistos) {
 500                 return fPrimGenMCTrackQAHistos;
 501         } else {
 502                 return 0x0;
 503         }
 504
 505 }
 506
 507 // -----------------------------------------------------------------------
 508 TList* AliAODTrackCutsDiHadronPID::GetListOfSecRecMCTrackQAHistos() const  {
 509
 510         // Returns the list of histograms of reconstructed secondary tracks.
 511         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 512
 513         if (fSecRecMCTrackQAHistos) {
 514                 return fSecRecMCTrackQAHistos;
 515         } else {
 516                 return 0x0;
 517         }
 518
 519 }
 520
 521 // -----------------------------------------------------------------------
 522 TList* AliAODTrackCutsDiHadronPID::GetListOfSecGenMCTrackQAHistos() const {
 523
 524         // Returns the list of histograms of generator level secondary particles.
 525         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 526
 527         if (fSecGenMCTrackQAHistos) {
 528                 return fSecGenMCTrackQAHistos;
 529         } else {
 530                 return 0x0;
 531         }
 532
 533 }
 534
 535 // -----------------------------------------------------------------------
 536 TH1F* AliAODTrackCutsDiHadronPID::GetHistDataTOFProjection(Int_t charge, Int_t species, Int_t ptbin) {
 537
 538         // Returns a projection in TOF of the data histogram.
 539         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 540
 541         Int_t ptclass = GetPtClass(ptbin);
 542         if (ptclass == -1) return 0x0;
 543         Int_t bininptclass = GetBinInPtClass(ptbin);
 544
 545         // Retrieve original the 3D histogram (we don't own it).
 546         TH3F* htmp = (TH3F*)GetHistData(Form("fHistDataPID%s%s%s",fHistoName[charge].Data(),fParticleName[species].Data(),fPtClassName[ptclass].Data()));
 547
 548         // Make the projection.
 549         TH1F* htmp_proj = (TH1F*)htmp->ProjectionY(Form("TOFprojection_%i_%i_%i",charge,species,ptbin),bininptclass,bininptclass);
 550
 551         // Some settings of the output histogram.
 552         htmp_proj->SetDirectory(0);
 553         htmp_proj->SetTitle(Form("%5.3f < p_{T} < %5.3f",GetPtMinPID(ptbin),GetPtMaxPID(ptbin)));
 554         htmp_proj->Sumw2();
 555
 556         return htmp_proj;
 557
 558 }
 559
 560 // -----------------------------------------------------------------------
 561 TObjArray* AliAODTrackCutsDiHadronPID::GetDataTOFProjection(Int_t charge, Int_t species) {
 562
 563         // Returns a TObjArray with all TOF histograms (as needed by AliSpectrumDiHadronPID)
 564         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 565
 566         TObjArray* aout = new TObjArray(GetNPtBinsPID());
 567         aout->SetOwner(kTRUE);
 568         for (Int_t iPtBin = 1; iPtBin < (GetNPtBinsPID() + 1); iPtBin++) {
 569                 aout->AddLast((TH1F*)GetHistDataTOFProjection(charge, species, iPtBin));
 570         }
 571
 572         return aout;
 573
 574 }
 575
 576 // -----------------------------------------------------------------------
 577 TH1F* AliAODTrackCutsDiHadronPID::GetHistDataTOFMismatch(Int_t charge, Int_t species, Int_t ptbin) {
 578
 579         // Returns a projection in TOF of the mismatch histogram.
 580         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 581
 582         Int_t ptclass = GetPtClass(ptbin);
 583         if (ptclass == -1) return 0x0;
 584         Int_t bininptclass = GetBinInPtClass(ptbin);
 585
 586         // Retrieve original the 2D histogram (we don't own it).
 587         TH2F* htmp = (TH2F*)GetHistData(Form("fHistTOFMismatch%s%s%s",fHistoName[charge].Data(),fParticleName[species].Data(),fPtClassName[ptclass].Data()));
 588
 589         // Make the projection.
 590         TH1F* htmp_proj = (TH1F*)htmp->ProjectionY(Form("TOFprojectionMismatch_%i_%i_%i",charge,species,ptbin),bininptclass,bininptclass);
 591
 592         // Some settings of the output histogram.
 593         htmp_proj->SetDirectory(0);
 594         htmp_proj->SetTitle(Form("%5.3f < p_{T} < %5.3f",GetPtMinPID(ptbin),GetPtMaxPID(ptbin)));
 595         htmp_proj->Sumw2();
 596
 597         return htmp_proj;
 598
 599 }
 600
 601 // -----------------------------------------------------------------------
 602 TObjArray* AliAODTrackCutsDiHadronPID::GetDataTOFMismatch(Int_t charge, Int_t species) {
 603
 604         // Returns a TObjArray with all TOF mismatch histograms (as needed by AliSpectrumDiHadronPID)
 605         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 606
 607         TObjArray* aout = new TObjArray(GetNPtBinsPID());
 608         aout->SetOwner(kTRUE);
 609         for (Int_t iPtBin = 1; iPtBin < (GetNPtBinsPID() + 1); iPtBin++) {
 610                 aout->AddLast((TH1F*)GetHistDataTOFMismatch(charge, species, iPtBin));
 611         }
 612
 613         return aout;
 614
 615 }
 616
 617 // -----------------------------------------------------------------------
 618 TH2F* AliAODTrackCutsDiHadronPID::GetHistDataTPCTOFProjection(Int_t charge, Int_t species, Int_t ptbin) {
 619
 620         // Returns a projection in TOF of the data histogram.
 621         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 622
 623         Int_t ptclass = GetPtClass(ptbin);
 624         if (ptclass == -1) return 0x0;
 625         Int_t bininptclass = GetBinInPtClass(ptbin);
 626
 627         // Retrieve original the 3D histogram (we don't own it).
 628         TH3F* htmp = (TH3F*)GetHistData(Form("fHistDataPID%s%s%s",fHistoName[charge].Data(),fParticleName[species].Data(),fPtClassName[ptclass].Data()));
 629
 630         // Make the projection.
 631         TAxis* ptaxis = htmp->GetXaxis();
 632         Int_t NbinsPt = ptaxis->GetNbins();
 633         //cout<<"bin in pt class: "<<bininptclass<<endl;
 634         ptaxis->SetRange(bininptclass, bininptclass);
 635         TH2F* htmp_proj = (TH2F*)htmp->Project3D("zy");
 636         htmp_proj->SetName(Form("TPCTOFprojection_%i_%i_%i",charge,species,ptbin));
 637
 638         //cout<<"projection: "<<htmp_proj<<endl;
 639
 640         // Some settings of the output histogram.
 641         htmp_proj->SetDirectory(0);
 642         htmp_proj->SetTitle(Form("%5.3f < p_{T} < %5.3f",GetPtMinPID(ptbin),GetPtMaxPID(ptbin)));
 643         htmp_proj->Sumw2();
 644
 645         // Putting back the range on the p_T axis.
 646         ptaxis->SetRange(1, NbinsPt);
 647
 648         return htmp_proj;
 649
 650 }
 651
 652 // -----------------------------------------------------------------------
 653 TObjArray* AliAODTrackCutsDiHadronPID::GetDataTPCTOFProjection(Int_t charge, Int_t species) {
 654
 655         // TODO: This is basically a copy of GetDataTOFProjection -> Make into one function.
 656
 657         // Returns a TObjArray with all TPC/TOF histograms (as needed by AliSpectrumDiHadronPID)
 658         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 659
 660         TObjArray* aout = new TObjArray(GetNPtBinsPID());
 661         aout->SetOwner(kTRUE);
 662         for (Int_t iPtBin = 1; iPtBin < (GetNPtBinsPID() + 1); iPtBin++) {
 663
 664                 aout->AddLast((TH2F*)GetHistDataTPCTOFProjection(charge, species, iPtBin));
 665                 //cout<<"aout entries: "<<aout->GetEntriesFast()<<endl;
 666         }
 667
 668         return aout;
 669
 670 }
 671
 672 // -----------------------------------------------------------------------
 673 TH2F* AliAODTrackCutsDiHadronPID::GetHistDataTPCTOFMismatch(Int_t charge, Int_t species, Int_t ptbin) {
 674
 675         // Returns a projection in TPC/TOF of the mismatch histogram.
 676         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 677
 678         Int_t ptclass = GetPtClass(ptbin);
 679         if (ptclass == -1) return 0x0;
 680         Int_t bininptclass = GetBinInPtClass(ptbin);
 681
 682         // Retrieve original the 3D histogram (we don't own it).
 683         TH3F* htmp = (TH3F*)GetHistData(Form("fHistTPCTOFMismatch%s%s%s",fHistoName[charge].Data(),fParticleName[species].Data(),fPtClassName[ptclass].Data()));
 684
 685         // Make the projection.
 686         TAxis* ptaxis = htmp->GetXaxis();
 687         Int_t NbinsPt = ptaxis->GetNbins();
 688         ptaxis->SetRange(bininptclass, bininptclass);
 689         TH2F* htmp_proj = (TH2F*)htmp->Project3D("zy");
 690         htmp_proj->SetName(Form("TPCTOFprojectionMismatch_%i_%i_%i",charge,species,ptbin));
 691
 692         // Some settings of the output histogram.
 693         htmp_proj->SetDirectory(0);
 694         htmp_proj->SetTitle(Form("%5.3f < p_{T} < %5.3f",GetPtMinPID(ptbin),GetPtMaxPID(ptbin)));
 695         htmp_proj->Sumw2();
 696
 697         // Putting back the range on the p_T axis.
 698         ptaxis->SetRange(1, NbinsPt);
 699
 700         return htmp_proj;
 701
 702 }
 703
 704 // -----------------------------------------------------------------------
 705 TObjArray* AliAODTrackCutsDiHadronPID::GetDataTPCTOFMismatch(Int_t charge, Int_t species) {
 706
 707         // Returns a TObjArray with all TPC/TOF mismatch histograms (as needed by AliSpectrumDiHadronPID)
 708         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 709
 710         TObjArray* aout = new TObjArray(GetNPtBinsPID());
 711         aout->SetOwner(kTRUE);
 712         for (Int_t iPtBin = 1; iPtBin < (GetNPtBinsPID() + 1); iPtBin++) {
 713                 aout->AddLast((TH2F*)GetHistDataTPCTOFMismatch(charge, species, iPtBin));
 714         }
 715
 716         return aout;
 717
 718 }
 719
 720 // -----------------------------------------------------------------------
 721 Double_t AliAODTrackCutsDiHadronPID::GetPtMin(Int_t bin) const {
 722
 723         // Same as: TAxis::GetBinLowEdge()
 724
 725         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 726
 727         if ((bin < 1) || (bin > fNPtBins + 1)) {
 728                 cout<<"Bin is out of range..."<<endl;
 729                 return -999.;
 730         } else {
 731                 return fPtAxis[bin - 1];
 732         }
 733
 734 }
 735
 736 // -----------------------------------------------------------------------
 737 Double_t AliAODTrackCutsDiHadronPID::GetPtMax(Int_t bin) const {
 738
 739         // Same as: TAxis::GetBinUpEdge()
 740
 741         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 742
 743         if ((bin < 1) || (bin > fNPtBins + 1)) {
 744                 cout<<"Bin is out of range..."<<endl;
 745                 return -999.;
 746         } else {
 747                 return fPtAxis[bin];
 748         }
 749
 750 }
 751
 752 // -----------------------------------------------------------------------
 753 Int_t AliAODTrackCutsDiHadronPID::GetNPtBinsPID(Int_t ptclass) const {
 754
 755         // Returns the number of pt bins that are used in every "pt class",
 756         // where a "pt class" is a range in pT which have the same range in
 757         // TOF/ TPC. If class == -1, then the retrun value is the total number
 758         // of bins in all classes together.
 759
 760         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 761
 762         if (ptclass == -1) {
 763                 Int_t nptbinspid = 0;
 764                 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
 765                         nptbinspid += fNPtBinsPID[iPtClass];
 766                 }
 767                 return nptbinspid;
 768         } else if (ptclass >= 0 && ptclass < 5) {
 769                 return fNPtBinsPID[ptclass];
 770         } else {
 771                 return -999;
 772         }
 773
 774 }
 775
 776 // -----------------------------------------------------------------------
 777 Double_t* AliAODTrackCutsDiHadronPID::GetPtAxisPID() const {
 778
 779         // Returns an array representing the pT axis for the PID
 780         // histograms.
 781
 782         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 783
 784         const Int_t nptbinspid = GetNPtBinsPID();
 785         Double_t* ptaxis = new Double_t[nptbinspid + 1];
 786
 787         for (Int_t iPtBin = 0; iPtBin < nptbinspid; iPtBin++) {
 788                 ptaxis[iPtBin] = GetPtMinPID(iPtBin + 1);
 789         }
 790
 791         ptaxis[nptbinspid] = GetPtMaxPID(nptbinspid);
 792         return ptaxis;
 793
 794 }
 795
 796 // -----------------------------------------------------------------------
 797 Double_t AliAODTrackCutsDiHadronPID::GetPtMinPID(Int_t bin) const {
 798
 799         // Same as: TAxis::GetBinLowEdge()
 800
 801         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 802
 803         Int_t ptclass = GetPtClass(bin);
 804
 805         if (ptclass == -1) {return -999.;}
 806
 807         Int_t bininptclass = GetBinInPtClass(bin);
 808
 809         Double_t minpt = fPtBoundaryPID[ptclass];
 810         Double_t maxpt = fPtBoundaryPID[ptclass+1];
 811         Double_t ptres = (maxpt - minpt)/((Double_t)fNPtBinsPID[ptclass]);
 812
 813         return (minpt + ptres * ((Double_t)(bininptclass - 1)) );
 814
 815 }
 816
 817 // -----------------------------------------------------------------------
 818 Double_t AliAODTrackCutsDiHadronPID::GetPtMaxPID(Int_t bin) const {
 819
 820         // Same as: TAxis::GetBinUpEdge()
 821
 822         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 823
 824         Int_t ptclass = GetPtClass(bin);
 825
 826         if (ptclass == -1) {return -999.;}
 827
 828         Int_t bininptclass = GetBinInPtClass(bin);
 829
 830         Double_t minpt = fPtBoundaryPID[ptclass];
 831         Double_t maxpt = fPtBoundaryPID[ptclass+1];
 832         Double_t ptres = (maxpt - minpt)/((Double_t)fNPtBinsPID[ptclass]);
 833
 834         return (minpt + ptres * ((Double_t)(bininptclass)) );
 835
 836 }
 837
 838 // -----------------------------------------------------------------------
 839 Double_t AliAODTrackCutsDiHadronPID::GetPtClassMin(Int_t ptclass) const {
 840
 841         // Returns the minimum p_T of a p_T class.
 842
 843         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 844
 845         if (ptclass >= 0 && ptclass < 5) {
 846                 return fPtBoundaryPID[ptclass];
 847         } else {
 848                 return -999;
 849         }
 850
 851 }
 852
 853 // -----------------------------------------------------------------------
 854 Double_t AliAODTrackCutsDiHadronPID::GetPtClassMax(Int_t ptclass) const {
 855
 856         // Returns the maximum p_T of a p_T class.
 857
 858         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 859
 860         if (ptclass >= 0 && ptclass < 5) {
 861                 return fPtBoundaryPID[ptclass+1];
 862         } else {
 863                 return -999;
 864         }
 865
 866 }
 867
 868 // -----------------------------------------------------------------------
 869 Bool_t AliAODTrackCutsDiHadronPID::RequestQAHistos(Int_t histoclass, Bool_t Enable3DSpectra, Bool_t EnablePIDHistos) {
 870
 871         // Request certain histograms to be filled.
 872
 873         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 874
 875         if ((histoclass > -1) && (histoclass < 12)) {
 876                 fHistRequested[histoclass] = kTRUE;
 877                 f3DSpectraEnabeled[histoclass] = Enable3DSpectra;
 878                 fPIDHistosEnabeled[histoclass] = EnablePIDHistos;
 879                 //cout<<"histoclass: "<<histoclass<<" requested: "<<fHistRequested[histoclass]<<endl;
 880                 return kTRUE;
 881         } else {
 882                 return kFALSE;
 883         }
 884 }
 885
 886 // -----------------------------------------------------------------------
 887 void AliAODTrackCutsDiHadronPID::SetPtRange(Double_t minpt, Double_t maxpt) {
 888
 889         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 890
 891         fMinPt = minpt;
 892         fMaxPt = maxpt;
 893         fTestPt = kTRUE;
 894 }
 895
 896 // -----------------------------------------------------------------------
 897 void AliAODTrackCutsDiHadronPID::SetFilterMask(UInt_t filtermask) {
 898
 899         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 900
 901         fFilterMask = filtermask;
 902         fTestFilterMask = kTRUE;
 903 }
 904
 905 // -----------------------------------------------------------------------
 906 void AliAODTrackCutsDiHadronPID::SetMaxEta(Double_t maxeta) {
 907
 908         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 909
 910         fMaxEta = maxeta;
 911         fTestMaxEta = kTRUE;
 912 }
 913
 914 // -----------------------------------------------------------------------
 915 void AliAODTrackCutsDiHadronPID::SetMaxRapidity(Double_t maxrapidity) {
 916
 917         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 918
 919         fMaxRapidity = maxrapidity;
 920         fTestMaxRapidity = kTRUE;
 921 }
 922
 923 // -----------------------------------------------------------------------
 924 void AliAODTrackCutsDiHadronPID::SetDemandNoMismatch() {
 925
 926         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 927
 928         fTestTOFmismatch = kTRUE;
 929 }
 930
 931 // -----------------------------------------------------------------------
 932 void AliAODTrackCutsDiHadronPID::SetDemandFlags(ULong_t demandedflags) {
 933
 934         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 935
 936         fDemandedFlags = demandedflags;
 937         fTestFlags = kTRUE;
 938 }
 939
 940 // -----------------------------------------------------------------------
 941 void AliAODTrackCutsDiHadronPID::SetMinimumNumberOfTPCClusters(Int_t minimumnumberoftpcclusters) {
 942
 943         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 944
 945         fMinimumNumberOfTPCClusters = minimumnumberoftpcclusters;
 946         fTestNumberOfTPCClusters = kTRUE;
 947 }
 948
 949 // -----------------------------------------------------------------------
 950 void AliAODTrackCutsDiHadronPID::SetDemandSPDCluster() {
 951
 952         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 953
 954         fTestSPDAny = kTRUE;
 955 }
 956
 957 // -----------------------------------------------------------------------
 958 void AliAODTrackCutsDiHadronPID::SetPtDeptDCACut(TFormula* DCAxyCutFormula, Double_t DCAzCut, UInt_t MinSPDHits) {
 959
 960         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 961
 962         fPtDeptDCAxyCutFormula = DCAxyCutFormula;
 963         fDCAzCut = DCAzCut;
 964         fMinSPDHitsForPtDeptDCAcut = MinSPDHits;
 965         fTestPtDeptDCAcut = kTRUE;
 966 }
 967
 968 // -----------------------------------------------------------------------
 969 void AliAODTrackCutsDiHadronPID::StartNewEvent() {
 970
 971         // FIXME: This method is now only suited for Data (3 histo classes only.)
 972         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 973
 974         // Resetting the counters.
 975         for (Int_t iHistoClass = 0; iHistoClass < 3; iHistoClass++) {
 976                 fNTracks[iHistoClass] = 0;
 977         }
 978 }
 979
 980 // -----------------------------------------------------------------------
 981 void AliAODTrackCutsDiHadronPID::EventIsDone(Bool_t IsMC) {
 982
 983         // FIXME: This method is now only suited for Data (3 histo classes only.)
 984         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 985
 986         // Fill NTracks histos.
 987
 988         for (Int_t iHistoClass = 0; iHistoClass < 3; iHistoClass++) {
 989
 990                 if (fHistRequested[iHistoClass]) {
 991
 992                         if (IsMC) {
 993                                 // THE FOLLOWING SHOULD NEVER HAPPEN.
 994                                 //if (!fHistPrimRecMCPt[iHistoClass]) InitializeRecMCHistos(iHistoClass);
 995                                 fHistPrimRecNTracks[iHistoClass]->Fill(fNTracks[iHistoClass]);
 996                         } else {
 997                                 // THE FOLLOWING SHOULD NEVER HAPPEN.
 998                                 //if (!fHistDataPt[iHistoClass]) InitializeDataHistos(iHistoClass);
 999                                 fHistDataNTracks[iHistoClass]->Fill(fNTracks[iHistoClass]);
1000                         }
1001                 }
1002         }
1003 }
1004
1005 // -----------------------------------------------------------------------
1006 void AliAODTrackCutsDiHadronPID::CreateHistos() {
1007
1008         // Function should be called by the UserCreateOutput() function of the
1009         // analysis task. This function then generates all the histograms that
1010         // were requested locally on the workernode. Even if case that the
1011         // histograms are not filled, it is imperative that they are still
1012         // created, and in the same order, otherwise problems with merging will
1013         // arise.
1014
1015         // TODO: In principle it should never happen that if this method is called,
1016         //       that the lists of QA histograms already exist, so this may become a fatal error
1017         //       instead of a warning.
1018
1019         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1020
1021         if (fIsMC) {
1022
1023                 fTOFMatchingStat = new TH1F("fTOFMatchingStat","fTOFMatchingStat",3,0,3);
1024                 (fTOFMatchingStat->GetXaxis())->SetBinLabel(1,"Match");
1025                 (fTOFMatchingStat->GetXaxis())->SetBinLabel(2,"Mismatch");
1026                 (fTOFMatchingStat->GetXaxis())->SetBinLabel(3,"No TOF hit");
1027
1028                 if (!fPrimGenMCTrackQAHistos) {
1029                         cout<<"AliAODTrackCutsDiHadronPID - Creating Prim. Gen. MC Track QA TList..."<<endl;
1030                         fPrimGenMCTrackQAHistos = new TList();
1031                         fPrimGenMCTrackQAHistos->SetName("PrimGenMCTrackQAHistos");
1032                         fPrimGenMCTrackQAHistos->SetOwner(kTRUE);
1033                 } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Prim. Gen. MC Track QA TList was already created..."<<endl;}
1034
1035                 if (!fSecGenMCTrackQAHistos) {
1036                         cout<<"AliAODTrackCutsDiHadronPID - Creating Sec. Gen. MC Track QA TList..."<<endl;
1037                         fSecGenMCTrackQAHistos = new TList();
1038                         fSecGenMCTrackQAHistos->SetName("SecGenMCTrackQAHistos");
1039                         fSecGenMCTrackQAHistos->SetOwner(kTRUE);
1040                         for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {if (fHistRequested[iHistoClass]) {InitializeGenMCHistos(iHistoClass);}}
1041                 } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Sec. Gen. MC Track QA TList was already created..."<<endl;}
1042
1043                 if (!fPrimRecMCTrackQAHistos) {
1044                         cout<<"AliAODTrackCutsDiHadronPID - Creating Prim. Rec. MC Track QA TList..."<<endl;
1045                         fPrimRecMCTrackQAHistos = new TList();
1046                         fPrimRecMCTrackQAHistos->SetName("PrimRecMCTrackQAHistos");
1047                         fPrimRecMCTrackQAHistos->SetOwner(kTRUE);
1048                 } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Prim. Rec. MC Track QA TList was already created..."<<endl;}
1049
1050                 if (!fSecRecMCTrackQAHistos) {
1051                         cout<<"AliAODTrackCutsDiHadronPID - Creating Sec. Rec. MC Track QA TList..."<<endl;
1052                         fSecRecMCTrackQAHistos = new TList();
1053                         fSecRecMCTrackQAHistos->SetName("SecRecMCTrackQAHistos");
1054                         fSecRecMCTrackQAHistos->SetOwner(kTRUE);
1055                         for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
1056                                 if (fHistRequested[iHistoClass]) {
1057                                         InitializeRecMCHistos(iHistoClass);
1058                                 }
1059                         }
1060                 } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Sec. Rec. MC Track QA TList was already created..."<<endl;}
1061
1062         } else {
1063
1064                 if (!fDataTrackQAHistos) {
1065                         cout<<"AliAODTrackCutsDiHadronPID - Creating Data Track QA TList..."<<endl;
1066                         fDataTrackQAHistos = new TList();
1067                         fDataTrackQAHistos->SetName("DataTrackQAHistos");
1068                         fDataTrackQAHistos->SetOwner(kTRUE);
1069
1070                         // Add general histograms.
1071                         fHistAcceptedFilterBits = InitializeAcceptedFilterBits("fHistAcceptedFilterBits");
1072                         fDataTrackQAHistos->Add(fHistAcceptedFilterBits);
1073
1074                         // Add histograms per class.
1075                         for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
1076                                 fHistDataDCAxyOneSigma[iHistoClass] = InitializeDCASpectrum("fHistDataDCAxyOneSigma",iHistoClass);
1077                                 fDataTrackQAHistos->Add(fHistDataDCAxyOneSigma[iHistoClass]);
1078                                 if (fHistRequested[iHistoClass]) {
1079                                         if (iHistoClass < 3) {InitializeDataHistos(iHistoClass);}
1080                                 }
1081                         }
1082                 } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Data QA TList was already created..."<<endl;}
1083
1084         }
1085 }
1086
1087 // -----------------------------------------------------------------------
1088 Bool_t AliAODTrackCutsDiHadronPID::IsSelectedData(AliTrackDiHadronPID* track, Double_t randomhittime) {
1089
1090         //
1091         // Checks performed on a data track.
1092         //
1093
1094         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1095
1096         if (!track) return kFALSE;
1097
1098         if (!fDataTrackQAHistos) {cout<<"AliAODTrackCutsDiHadronPID - Histograms were not created, you should have called CreateHistos()..."<<endl;}
1099
1100         // Check the track cuts (NB. check function return kTRUE if track is accepted.)
1101         if (!CheckPt(track->Pt())) return kFALSE;
1102         if (!CheckMaxEta(track->Eta())) return kFALSE;
1103         if (!CheckFilterMask(track->GetFilterMap())) return kFALSE;
1104         if (!CheckFlags(track->GetFlags())) return kFALSE;
1105         if (!CheckNclsTPC(track->GetNclsTPC())) return kFALSE;
1106         if (!CheckTOFmismatch(track->IsTOFMismatch())) return kFALSE;
1107
1108         Int_t NSPDhits = 0;
1109         if (track->HasPointOnITSLayer(0)) NSPDhits++;
1110         if (track->HasPointOnITSLayer(1)) NSPDhits++;
1111         if (!CheckPtDeptDCACut(track->GetZAtDCA(),track->GetXYAtDCA(),track->Pt(),NSPDhits)) return kFALSE;
1112         if (fTestSPDAny) {if (NSPDhits < 1) return kFALSE;}
1113
1114         // Fill the filterbit histogram.
1115         for (Int_t iBin = (fRelevantBitsArray->GetSize() - 1); iBin >= 0; --iBin) {
1116                 if ( (((Int_t)track->GetFilterMap())&(1<<fRelevantBitsArray->At(iBin))) == (1<<fRelevantBitsArray->At(iBin))) {
1117                         fHistAcceptedFilterBits->Fill(iBin);
1118                         break;
1119                 }
1120         }
1121
1122         // Track has passed the cuts, fill QA histograms.
1123         for (Int_t iHistoClass = 0; iHistoClass < 3; iHistoClass++) {
1124
1125                 if (fHistRequested[iHistoClass]) {
1126
1127                         // Check the charge (could be neater).
1128                         if ((iHistoClass == 0) && (track->Charge() == 0)) continue;
1129                         if ((iHistoClass == 1) && (track->Charge() <= 0)) continue;
1130                         if ((iHistoClass == 2) && (track->Charge() >= 0)) continue;
1131
1132                         FillDataHistos(iHistoClass, track);
1133
1134                         // Ignore if random hit is < -1.e20.
1135                         if (randomhittime > -1.e20) FillTOFMismatchHistos(iHistoClass,track,randomhittime);
1136
1137                         fNTracks[iHistoClass]++;
1138
1139                 }
1140         }
1141
1142         // Track Has Passed.
1143         return kTRUE;
1144
1145 }
1146
1147 // -----------------------------------------------------------------------
1148 Bool_t AliAODTrackCutsDiHadronPID::IsSelectedGeneratedMC(AliAODMCParticle* particle) {
1149
1150         //
1151         // Checks performed on a generated MC particle.
1152         //
1153
1154         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1155
1156         // PDG codes for particles:
1157         //  pi+ 211, K+ 321, p 2212
1158
1159         if (!particle) return kFALSE;
1160
1161         if (!fPrimGenMCTrackQAHistos||!fSecGenMCTrackQAHistos) {cout<<"AliAODTrackCutsDiHadronPID - Histograms were not created, you should have called CreateHistos()..."<<endl;}
1162
1163         // Check the track cuts.
1164         if (!CheckPt(particle->Pt())) return kFALSE;
1165         if (!CheckMaxEta(particle->Eta())) return kFALSE;
1166         if (!CheckRapidity(particle->Y())) return kFALSE;               // NEW: rapidity cut. (not done on data)
1167
1168         // NOT YET IMPLEMENTED! - Histoclasses for specific particles.
1169         for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
1170
1171                 if (fHistRequested[iHistoClass]) {
1172
1173                         // Check the charge (could be neater).
1174                         if ((iHistoClass == 0) && (particle->Charge() == 0)) continue;
1175                         if ((iHistoClass == 1) && (particle->Charge() <= 0)) continue;
1176                         if ((iHistoClass == 2) && (particle->Charge() >= 0)) continue;
1177                         if ((iHistoClass == 3) && (TMath::Abs(particle->GetPdgCode()) != 211)) continue;
1178                         if ((iHistoClass == 4) && (particle->GetPdgCode()) != 211) continue;
1179                         if ((iHistoClass == 5) && (particle->GetPdgCode()) != -211) continue;
1180                         if ((iHistoClass == 6) && (TMath::Abs(particle->GetPdgCode()) != 321)) continue;
1181                         if ((iHistoClass == 7) && (particle->GetPdgCode()) != 321) continue;
1182                         if ((iHistoClass == 8) && (particle->GetPdgCode()) != -321) continue;
1183                         if ((iHistoClass == 9) && (TMath::Abs(particle->GetPdgCode()) != 2212)) continue;
1184                         if ((iHistoClass == 10) && (particle->GetPdgCode()) != 2212) continue;
1185                         if ((iHistoClass == 11) && (particle->GetPdgCode()) != -2212) continue;
1186
1187                         // Secondary specification not set.
1188                         //cout<<"Charge: "<<particle->Charge()<<" PDG code: "<<particle->GetPdgCode()<<endl;
1189                         //cout<<particle->IsPhysicalPrimary()<<" "<<particle->IsSecondaryFromWeakDecay()<<" "<<particle->IsSecondaryFromMaterial()<<endl;
1190
1191                         // These two functions are not implemented...
1192                         //if (particle->IsSecondaryFromWeakDecay()) cout<<"Secondary From Weak Decay!"<<endl;
1193                         //if (particle->IsSecondaryFromMaterial()) cout<<"Secondary From Material!"<<endl;
1194
1195                         FillGenMCHistos(iHistoClass, particle);
1196
1197                 }
1198         }
1199
1200         // Particle has Passed.
1201         return kTRUE;
1202
1203 }
1204
1205 // -----------------------------------------------------------------------
1206 Bool_t AliAODTrackCutsDiHadronPID::IsSelectedReconstructedMC(AliTrackDiHadronPID* track) {
1207
1208         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1209
1210         if (!track) return kFALSE;
1211
1212         if (!fPrimRecMCTrackQAHistos||!fSecRecMCTrackQAHistos) {cout<<"AliAODTrackCutsDiHadronPID - Histograms were not created, you should have called CreateHistos()..."<<endl;}
1213
1214         // Check the track cuts.
1215         if (!CheckPt(track->MCPt())) return kFALSE;                                     // Kinematic cuts are done on the MC particles
1216         if (!CheckMaxEta(track->MCEta())) return kFALSE;                        // and not on the data.
1217         if (!CheckRapidity(track->MCY())) return kFALSE;                        // NEW: rapidity cut.
1218         if (!CheckFilterMask(track->GetFilterMap())) return kFALSE;
1219         if (!CheckFlags(track->GetFlags())) return kFALSE;
1220         if (!CheckNclsTPC(track->GetNclsTPC())) return kFALSE;
1221
1222         Int_t NSPDhits = 0;
1223         if (track->HasPointOnITSLayer(0)) NSPDhits++;
1224         if (track->HasPointOnITSLayer(1)) NSPDhits++;
1225         if (!CheckPtDeptDCACut(track->GetZAtDCA(),track->GetXYAtDCA(),track->Pt(),NSPDhits)) return kFALSE;
1226         if (fTestSPDAny) {if (NSPDhits < 1) return kFALSE;}
1227
1228         // Track has passed the cuts, fill QA histograms.
1229         for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
1230
1231                 if (fHistRequested[iHistoClass]) {
1232
1233                         // Check the charge (could be neater).
1234                         if ((iHistoClass == 0) && (track->Charge() == 0)) continue;
1235                         if ((iHistoClass == 1) && (track->Charge() <= 0)) continue;
1236                         if ((iHistoClass == 2) && (track->Charge() >= 0)) continue;
1237                         if ((iHistoClass == 3) && (TMath::Abs(track->GetPdgCode()) != 211)) continue;
1238                         if ((iHistoClass == 4) && (track->GetPdgCode()) != 211) continue;
1239                         if ((iHistoClass == 5) && (track->GetPdgCode()) != -211) continue;
1240                         if ((iHistoClass == 6) && (TMath::Abs(track->GetPdgCode()) != 321)) continue;
1241                         if ((iHistoClass == 7) && (track->GetPdgCode()) != 321) continue;
1242                         if ((iHistoClass == 8) && (track->GetPdgCode()) != -321) continue;
1243                         if ((iHistoClass == 9) && (TMath::Abs(track->GetPdgCode()) != 2212)) continue;
1244                         if ((iHistoClass == 10) && (track->GetPdgCode()) != 2212) continue;
1245                         if ((iHistoClass == 11) && (track->GetPdgCode()) != -2212) continue;
1246
1247                         FillRecMCHistos(iHistoClass, track);
1248
1249                         fNTracks[iHistoClass]++;
1250                 }
1251         }
1252
1253         // Track Has Passed.
1254         return kTRUE;
1255
1256 }
1257
1258 // -----------------------------------------------------------------------
1259 Int_t AliAODTrackCutsDiHadronPID::GetPtClass(Int_t ptbin) const {
1260
1261         // Returns a p_T class as a function of bin (PID histos)
1262         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1263
1264         Int_t currentptclass = 0;
1265         Int_t currentptbin = fNPtBinsPID[0];
1266
1267         while (currentptbin < ptbin) {
1268                 currentptclass++;
1269                 if (currentptclass == 5) {break;}
1270                 currentptbin += fNPtBinsPID[currentptclass];
1271         }
1272
1273         if (currentptclass == 5) {
1274                 cout<<"GetPtClass -> ptbin out of range!"<<endl;
1275                 return -1;
1276         }
1277
1278         return currentptclass;
1279 }
1280
1281 // -----------------------------------------------------------------------
1282 Int_t AliAODTrackCutsDiHadronPID::GetBinInPtClass(Int_t ptbin) const {
1283
1284         // Returns the bin withing the p_T class (i.e., 1..Nbins)
1285         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1286
1287         Int_t ptclass = GetPtClass(ptbin);
1288         if (ptclass == -1) {return -1;}
1289
1290         Int_t ptbinout = ptbin;
1291         for (Int_t iPtClass = 0; iPtClass < ptclass; iPtClass++) {ptbinout -= fNPtBinsPID[iPtClass];}
1292
1293         return ptbinout;
1294
1295 }
1296
1297 // -----------------------------------------------------------------------
1298 Bool_t AliAODTrackCutsDiHadronPID::CheckPt(Double_t pt) const {
1299
1300         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1301
1302         if (!fTestPt) return kTRUE;
1303         if ((pt > fMinPt) && (pt < fMaxPt)) return kTRUE;
1304         return kFALSE;
1305         }
1306
1307 // -----------------------------------------------------------------------
1308 Bool_t AliAODTrackCutsDiHadronPID::CheckMaxEta(Double_t eta) const {
1309
1310         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1311
1312         if (!fTestMaxEta) return kTRUE;                         // Accepted if there is no check on this parameter.
1313         if (TMath::Abs(eta) < fMaxEta) return kTRUE;
1314         return kFALSE;
1315         }
1316
1317 // -----------------------------------------------------------------------
1318 Bool_t AliAODTrackCutsDiHadronPID::CheckRapidity(Double_t rap) const {
1319
1320         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1321
1322         if (!fTestMaxRapidity) return kTRUE;
1323         if (TMath::Abs(rap) < fMaxRapidity) return kTRUE;
1324         return kFALSE;
1325         }
1326
1327 // -----------------------------------------------------------------------
1328 Bool_t AliAODTrackCutsDiHadronPID::CheckFilterMask(UInt_t filtermap) const {
1329
1330         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1331
1332         if (!fTestFilterMask) return kTRUE;
1333         if (fFilterMask & filtermap) return kTRUE;
1334         return kFALSE;
1335         }
1336
1337 // -----------------------------------------------------------------------
1338 Bool_t AliAODTrackCutsDiHadronPID::CheckFlags(ULong_t flags) const {
1339
1340         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1341
1342         if (!fTestFlags) return kTRUE;
1343         if ((flags & fDemandedFlags) == fDemandedFlags) return kTRUE;
1344         return kFALSE;
1345         }
1346
1347 // -----------------------------------------------------------------------
1348 Bool_t AliAODTrackCutsDiHadronPID::CheckNclsTPC(Int_t ncls) const {
1349
1350         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1351
1352         if (!fTestNumberOfTPCClusters) return kTRUE;
1353         if (ncls > fMinimumNumberOfTPCClusters) return kTRUE;
1354         return kFALSE;
1355         }
1356
1357 // -----------------------------------------------------------------------
1358 Bool_t AliAODTrackCutsDiHadronPID::CheckTOFmismatch(Bool_t ismismatch) const {
1359
1360         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1361
1362         if (!fTestTOFmismatch) return kTRUE; // if we're not cutting on mismatch, then it's accepted.
1363         if (!ismismatch) return kTRUE;          // so if the track is not a mismatch, then it is accepted.
1364         return kFALSE;                                          // if it is a mismatch, then it's not accepted.
1365         }
1366
1367 // -----------------------------------------------------------------------
1368 Bool_t AliAODTrackCutsDiHadronPID::CheckPtDeptDCACut(Double_t dcaz, Double_t dcaxy, Double_t pt, UInt_t SPDhits) const {
1369
1370         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1371
1372         if (!fTestPtDeptDCAcut) return kTRUE;
1373         if (SPDhits < fMinSPDHitsForPtDeptDCAcut) return kTRUE; // If there are not enough SPD hits to do the cut.
1374         if ((dcaz < fDCAzCut) && (dcaxy < fPtDeptDCAxyCutFormula->Eval(pt))) return kTRUE;
1375         return kFALSE;
1376 }
1377
1378 // -----------------------------------------------------------------------
1379 Bool_t AliAODTrackCutsDiHadronPID::FillDataHistos(Int_t histoclass, AliTrackDiHadronPID* track) {
1380
1381         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1382
1383         // Fill the histograms.
1384         fHistDataPt[histoclass]->Fill(track->Pt());
1385         fHistDataPhiEtaPt[histoclass]->Fill(track->Phi(),track->Eta(),track->Pt());
1386
1387         // Fill DCA histograms.
1388         fHistDataDCAxy[histoclass]->Fill(track->GetXYAtDCA());
1389         fHistDataDCAz[histoclass]->Fill(track->GetZAtDCA());
1390
1391         Int_t checkSum = 0;
1392
1393         /* Philip:
1394                 Introduced a separate selection mechanism here as the TPC is off for low pt for
1395                 in particular protons and kaons. Therefor use only the TOF for identifying under 1.8 GeV.
1396                 Setting this with fLowPtNSigmaTOFOnly flag.
1397         */
1398         // histoclass = 0: all charges; histoclass = 1: positive; histoclass = 2: negative;
1399         if (TMath::Sqrt(track->GetNumberOfSigmasTOF(0) * track->GetNumberOfSigmasTOF(0) +
1400                                         track->GetNumberOfSigmasTPC(0) * track->GetNumberOfSigmasTPC(0)) < 1.) {
1401                 fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA());  // All species.
1402                 fHistDataDCAxyOneSigma[3 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA());  // Pions.
1403         //      checkSum++;     cout<<"Pion found: nSigTOF: "<<track->GetNumberOfSigmasTOF(0)<<"; nSigTPC: "<<track->GetNumberOfSigmasTPC(0)<<endl;
1404         }
1405
1406         // for protons and low pt, only when fLowPtNSigmaTOFOnly set to kTRUE:
1407         if ((TMath::Abs(track->Pt()) < 1.8) && fLowPtNSigmaTOFOnly) {
1408
1409                 if (TMath::Abs(track->GetNumberOfSigmasTOF(1)) < 1.) {
1410                         fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA());  // All species.
1411                         fHistDataDCAxyOneSigma[6 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA());  // Kaons.
1412                 }
1413                 if (TMath::Abs(track->GetNumberOfSigmasTOF(2)) < 1.) {
1414                         fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA());  // All species.
1415                         fHistDataDCAxyOneSigma[9 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA());  // Protons.
1416                 //      checkSum++;     cout<<"Proton found: nSigTOF: "<<track->GetNumberOfSigmasTOF(2)<<"; nSigTPC: "<<track->GetNumberOfSigmasTPC(2)<<endl;
1417                 }
1418         } else {
1419                 if (TMath::Sqrt(track->GetNumberOfSigmasTOF(1) * track->GetNumberOfSigmasTOF(1) +
1420                                                 track->GetNumberOfSigmasTPC(1) * track->GetNumberOfSigmasTPC(1)) < 1.) {
1421                         fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA());  // All species.
1422                         fHistDataDCAxyOneSigma[6 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA());  // Kaons.
1423                 //      checkSum++;     cout<<"Kaon found: nSigTOF: "<<track->GetNumberOfSigmasTOF(1)<<"; nSigTPC: "<<track->GetNumberOfSigmasTPC(1)<<endl;
1424                 }
1425                 if (TMath::Sqrt(track->GetNumberOfSigmasTOF(2) * track->GetNumberOfSigmasTOF(2) +
1426                                                 track->GetNumberOfSigmasTPC(2) * track->GetNumberOfSigmasTPC(2)) < 1.) {
1427                         fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA());  // All species.
1428                         fHistDataDCAxyOneSigma[9 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA());  // Protons.
1429                 //      checkSum++;     cout<<"Proton found: nSigTOF: "<<track->GetNumberOfSigmasTOF(2)<<"; nSigTPC: "<<track->GetNumberOfSigmasTPC(2)<<endl;
1430                 }
1431         }
1432
1433         // check for double identification (or triple..)
1434         if(checkSum > 1) {AliError("More than one particle identified for the same track!"); }
1435
1436         // Fill PID histos.
1437         for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
1438
1439                 // Note that a possible Y cut is only done on the PID cuts!
1440                 if (!CheckRapidity(track->Y(iSpecies))) continue;
1441
1442                 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
1443                         fHistDataPID[histoclass][iSpecies][iPtClass]->Fill(track->Pt(),
1444                                 (fUseNSigmaOnPIDAxes ? track->GetNumberOfSigmasTOF(iSpecies) : track->GetTOFsignalMinusExpected(iSpecies)),
1445                                 (fUseNSigmaOnPIDAxes ? track->GetNumberOfSigmasTPC(iSpecies) : track->GetTPCsignalMinusExpected(iSpecies)));
1446                 }
1447         }
1448
1449         return kTRUE;
1450
1451 }
1452
1453 // -----------------------------------------------------------------------
1454 Bool_t AliAODTrackCutsDiHadronPID::FillTOFMismatchHistos(Int_t histoclass, AliTrackDiHadronPID* track, Double_t randomhittime) {
1455
1456         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1457
1458         // Fill TOF mismatch.
1459         for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
1460
1461                 // Note that a possible Y cut is only done on the PID cuts!
1462                 if (!CheckRapidity(track->Y(iSpecies))) continue;
1463 /*
1464                 if (!track->GetNumberOfSigmasTOF(iSpecies)) {
1465                         cout<<"NSigma TOF: "<<track->GetNumberOfSigmasTOF(iSpecies)<<endl;
1466                         cout<<"TOFsignal exp: "<<track->GetTOFsignalExpected(iSpecies)<<endl;
1467                         cout<<"TOFsigma exp: "<<track->GetTOFsigmaExpected(iSpecies)<<endl;
1468                         cout<<endl;
1469                 }
1470 */
1471
1472                 Double_t TOFmismatchSigma = randomhittime - track->GetTOFsignalExpected(iSpecies);
1473                 if (fUseNSigmaOnPIDAxes) {
1474                         if (track->GetTOFsigmaExpected(iSpecies) < 10e-30) {
1475                                 /*
1476                                 cout << "ERROR: division through (almost) zero on the mismatch signal..."<<endl;
1477                                 cout << "Random time: "<<randomhittime<<endl;
1478                                 cout << "Exp TOF signal track for species "<<iSpecies<<" "<<track->GetTOFsignalExpected(iSpecies)<<endl;
1479                                 cout << "Exp Mismatch signal assuming species "<< iSpecies<<": "<<TOFmismatchSigma<<endl;
1480                                 cout << "Expected TOF sigma: "<<track->GetTOFsigmaExpected(iSpecies)<<endl;
1481                                 ULong_t flags_kTOFout_kTIME = (UInt_t)(AliAODTrack::kTOFout)|(UInt_t)(AliAODTrack::kTIME);
1482                                 cout << "Random track TOF status (kTOFout, kTIME): "<<(((track->GetFlags())&flags_kTOFout_kTIME)==flags_kTOFout_kTIME)<<endl;
1483                                 */
1484                         } else {
1485                                 TOFmismatchSigma /= track->GetTOFsigmaExpected(iSpecies);
1486                         }
1487                 }
1488
1489                 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
1490                         fHistTOFMismatch[histoclass][iSpecies][iPtClass]->Fill(track->Pt(), TOFmismatchSigma);
1491                         fHistTPCTOFMismatch[histoclass][iSpecies][iPtClass]->Fill(track->Pt(), TOFmismatchSigma,
1492                                 (fUseNSigmaOnPIDAxes ? track->GetNumberOfSigmasTPC(iSpecies) : track->GetTPCsignalMinusExpected(iSpecies)));
1493                 }
1494         }
1495
1496         return kTRUE;
1497
1498 }
1499
1500 // -----------------------------------------------------------------------
1501 Bool_t AliAODTrackCutsDiHadronPID::FillGenMCHistos(Int_t histoclass, AliAODMCParticle* particle) {
1502
1503         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1504
1505         //cout << "histoclass: "<<histoclass<<" particle: "<<particle<<endl;
1506
1507         // Fill the histograms.
1508         if (particle->IsPhysicalPrimary()) {
1509                 //cout<<"Trying to fill: "<<fHistPrimGenMCPt[histoclass]->GetName()<<" of type: " <<fHistPrimGenMCPt[histoclass]->ClassName()<<" at: "<<fHistPrimGenMCPt[histoclass]<<endl;
1510                 fHistPrimGenMCPt[histoclass]->Fill(particle->Pt());
1511                 //cout<< "phi: "<<particle->Phi()<< " eta: "<<particle->Eta() << " pt: "<<particle->Pt()<<" ";
1512                 fHistPrimGenMCPhiEtaPt[histoclass]->Fill(particle->Phi(), particle->Eta(), particle->Pt());
1513                 //cout<<"OK!"<<endl;
1514         } else {
1515                 //cout<<"Trying to fill: "<<fHistPrimGenMCPt[histoclass]->GetName()<<" of type: " <<fHistPrimGenMCPt[histoclass]->ClassName()<<" at: "<<fHistPrimGenMCPt[histoclass]<<endl;
1516                 fHistSecGenMCPt[histoclass]->Fill(particle->Pt());
1517                 //cout<< "phi: "<<particle->Phi()<< " eta: "<<particle->Eta() << " pt: "<<particle->Pt()<<" ";
1518                 fHistSecGenMCPhiEtaPt[histoclass]->Fill(particle->Phi(), particle->Eta(), particle->Pt());
1519                 //cout<<"OK!"<<endl;
1520         }
1521
1522         return kTRUE;
1523
1524 }
1525
1526 // -----------------------------------------------------------------------
1527 Bool_t AliAODTrackCutsDiHadronPID::FillRecMCHistos(Int_t histoclass, AliTrackDiHadronPID* track) {
1528
1529         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1530
1531         // Fill the Pt and acceptance histograms.
1532         if (track->IsPhysicalPrimary()) {
1533                 fHistPrimRecMCPt[histoclass]->Fill(track->MCPt());
1534                 fHistPrimRecPtGenPt[histoclass]->Fill(track->MCPt(), track->Pt());
1535                 fHistPrimRecMCPhiEtaPt[histoclass]->Fill(track->MCPhi(), track->MCEta(), track->MCPt());
1536         } else {
1537                 fHistSecRecMCPt[histoclass]->Fill(track->MCPt());
1538                 fHistSecRecMCPhiEtaPt[histoclass]->Fill(track->MCPhi(), track->MCEta(), track->MCPt());
1539         }
1540
1541         // Fill the DCA histograms.
1542         if (track->IsPhysicalPrimary()) {
1543                 fHistPrimRecMCDCA[histoclass]->Fill(track->Pt(), track->GetXYAtDCA());
1544         }
1545         if (track->IsSecondaryFromMaterial()) {
1546                 fHistSecRecMCDCAMat[histoclass]->Fill(track->Pt(), track->GetXYAtDCA());
1547         }
1548         if (track->IsSecondaryFromWeakDecay()) {
1549                 fHistSecRecMCDCAWeak[histoclass]->Fill(track->Pt(), track->GetXYAtDCA());
1550         }
1551
1552         // Fill PID histos.
1553         if (histoclass < 3) {
1554
1555                 // Note that the following histogram is only sensible if the fTOFlabel is set properly.
1556                 // If not, the histogram will be filled with "no match".
1557                 fTOFMatchingStat->Fill(((Double_t)track->GetTOFMatchingStatus())+0.5);
1558
1559                 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
1560
1561                         // Note that a possible Y cut is only done on the PID cuts!
1562                         if (!CheckRapidity(track->MCY())) continue;
1563
1564                         for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
1565                                 //cout << "recpt: " << track->Pt() << " mcpt: "<<track->MCPt() << " DTOF: " <<
1566                                 //(fUseNSigmaOnPIDAxes ? track->GetNumberOfSigmasTOF(iSpecies) : track->GetTOFsignalMinusExpected(iSpecies)) << " DTPC: " <<
1567                                 //(fUseNSigmaOnPIDAxes ? track->GetNumberOfSigmasTPC(iSpecies) : track->GetTPCsignalMinusExpected(iSpecies)) << endl;
1568
1569                                 fHistPrimRecPID[histoclass][iSpecies][iPtClass]->Fill(track->MCPt(),
1570                                         (fUseNSigmaOnPIDAxes ? track->GetNumberOfSigmasTOF(iSpecies) : track->GetTOFsignalMinusExpected(iSpecies)));
1571
1572                                 if (track->IsTOFMismatch()) {
1573                                         fHistPrimRecMismatch[histoclass][iSpecies][iPtClass]->Fill(track->MCPt(),
1574                                                 (fUseNSigmaOnPIDAxes ? track->GetNumberOfSigmasTOF(iSpecies) : track->GetTOFsignalMinusExpected(iSpecies)));
1575                                 }
1576                         }
1577                 }
1578         }
1579
1580         return kTRUE;
1581
1582 }
1583
1584 // -----------------------------------------------------------------------
1585 Bool_t AliAODTrackCutsDiHadronPID::InitializeDataHistos(Int_t histoclass) {
1586
1587         cout<<"AliAODTrackCutsDiHadronPID - Creating Data Histograms of Class "<<histoclass<<endl;
1588         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1589
1590         // Add the Pt spectra.
1591         fHistDataPt[histoclass] = InitializePtSpectrum("fHistDataPt",histoclass);
1592         fDataTrackQAHistos->Add(fHistDataPt[histoclass]);
1593         fHistDataPhiEtaPt[histoclass] = InitializePhiEtaPt("fHistDataPhiEtaPt",histoclass);
1594         fDataTrackQAHistos->Add(fHistDataPhiEtaPt[histoclass]);
1595
1596         // Add the NTrack histograms.
1597         fHistDataNTracks[histoclass] = InitializeNTracksHisto("fHistDataNTracks",histoclass);
1598         fDataTrackQAHistos->Add(fHistDataNTracks[histoclass]);
1599
1600         // Add the DCA histograms.
1601         fHistDataDCAxy[histoclass] = InitializeDCAxyHisto("fHistDCAxy",histoclass);
1602         fDataTrackQAHistos->Add(fHistDataDCAxy[histoclass]);
1603         fHistDataDCAz[histoclass] = InitializeDCAzHisto("fHistDCAz",histoclass);
1604         fDataTrackQAHistos->Add(fHistDataDCAz[histoclass]);
1605
1606         // Add the PID histograms. (FIXME shoudl be able to turn the creation of these histos off.)
1607         for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
1608                 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
1609
1610                         fHistDataPID[histoclass][iSpecies][iPtClass] = InitializePIDHisto("fHistDataPID",histoclass,iSpecies,iPtClass);
1611                         fDataTrackQAHistos->Add(fHistDataPID[histoclass][iSpecies][iPtClass]);
1612
1613                         fHistTOFMismatch[histoclass][iSpecies][iPtClass] = InitializeTOFMismatchHisto("fHistTOFMismatch",histoclass,iSpecies,iPtClass);
1614                         fDataTrackQAHistos->Add(fHistTOFMismatch[histoclass][iSpecies][iPtClass]);
1615
1616                         fHistTPCTOFMismatch[histoclass][iSpecies][iPtClass] = InitializePIDHisto("fHistTPCTOFMismatch",histoclass,iSpecies,iPtClass);
1617                         fHistTPCTOFMismatch[histoclass][iSpecies][iPtClass]->SetTitle(Form("PID %s (Exp: %s)",fHistoName[histoclass].Data(),fParticleName[iSpecies].Data()));
1618                         fDataTrackQAHistos->Add(fHistTPCTOFMismatch[histoclass][iSpecies][iPtClass]);
1619
1620                 }
1621         }
1622
1623         return kTRUE;
1624
1625 }
1626
1627 // -----------------------------------------------------------------------
1628 Bool_t AliAODTrackCutsDiHadronPID::InitializeGenMCHistos(Int_t histoclass) {
1629
1630         cout<<"AliAODTrackCutsDiHadronPID - Creating Generated MC Histograms of Class "<<histoclass<<endl;
1631         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1632
1633         // Primary Particles.
1634         fHistPrimGenMCPt[histoclass] = InitializePtSpectrum("fHistPrimGenMCPt",histoclass);
1635         fPrimGenMCTrackQAHistos->Add(fHistPrimGenMCPt[histoclass]);
1636
1637         fHistPrimGenMCPhiEtaPt[histoclass] = InitializePhiEtaPt("fHistPrimGenMCPhiEtaPt",histoclass);
1638         fPrimGenMCTrackQAHistos->Add(fHistPrimGenMCPhiEtaPt[histoclass]);
1639
1640         // Secondary Particles.
1641         fHistSecGenMCPt[histoclass] = InitializePtSpectrum("fHistSecGenMCPt",histoclass);
1642         fSecGenMCTrackQAHistos->Add(fHistSecGenMCPt[histoclass]);
1643
1644         fHistSecGenMCPhiEtaPt[histoclass] = InitializePhiEtaPt("fHistSecGenMCPhiEtaPt",histoclass);
1645         fSecGenMCTrackQAHistos->Add(fHistSecGenMCPhiEtaPt[histoclass]);
1646
1647         return kTRUE;
1648
1649 }
1650
1651 // -----------------------------------------------------------------------
1652 Bool_t AliAODTrackCutsDiHadronPID::InitializeRecMCHistos(Int_t histoclass) {
1653
1654         cout<<"AliAODTrackCutsDiHadronPID - Creating Reconstructed MC Histograms of Class "<<histoclass<<endl;
1655         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1656
1657         // Primary Particles.
1658         fHistPrimRecMCPt[histoclass] = InitializePtSpectrum("fHistPrimRecMCPt",histoclass);
1659         fPrimRecMCTrackQAHistos->Add(fHistPrimRecMCPt[histoclass]);
1660
1661         fHistPrimRecPtGenPt[histoclass] = InitializeRecPtGenPt("fHistPrimRecPtGenPt",histoclass);
1662         fPrimRecMCTrackQAHistos->Add(fHistPrimRecPtGenPt[histoclass]);
1663
1664         fHistPrimRecMCPhiEtaPt[histoclass] = InitializePhiEtaPt("fHistPrimRecMCPhiEtaPt",histoclass);
1665         fPrimRecMCTrackQAHistos->Add(fHistPrimRecMCPhiEtaPt[histoclass]);
1666
1667         fHistPrimRecNTracks[histoclass] = InitializeNTracksHisto("fHistPrimRecNTracks",histoclass);
1668         fPrimRecMCTrackQAHistos->Add(fHistPrimRecNTracks[histoclass]);
1669
1670         fHistPrimRecMCDCA[histoclass] = InitializeDCASpectrum("fHistPrimRecDCA",histoclass);
1671         fPrimRecMCTrackQAHistos->Add(fHistPrimRecMCDCA[histoclass]);
1672
1673         if (histoclass < 3) {
1674                 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
1675                         for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
1676
1677                                 fHistPrimRecPID[histoclass][iSpecies][iPtClass] = InitializeTOFHisto("fHistPrimRecPID",histoclass,iSpecies,iPtClass);
1678                                 fPrimRecMCTrackQAHistos->Add(fHistPrimRecPID[histoclass][iSpecies][iPtClass]);
1679
1680                                 fHistPrimRecMismatch[histoclass][iSpecies][iPtClass] = InitializeTOFMismatchHisto("fHistPrimRecMismatch",histoclass,iSpecies,iPtClass);
1681                                 fPrimRecMCTrackQAHistos->Add(fHistPrimRecMismatch[histoclass][iSpecies][iPtClass]);
1682
1683
1684                         }
1685                 }
1686         }
1687
1688         // Secondary Particles.
1689         fHistSecRecMCPt[histoclass] = InitializePtSpectrum("fHistSecRecMCPt",histoclass);
1690         fSecRecMCTrackQAHistos->Add(fHistSecRecMCPt[histoclass]);
1691
1692         fHistSecRecMCPhiEtaPt[histoclass] = InitializePhiEtaPt("fHistSecRecMCPhiEtaPt",histoclass);
1693         fSecRecMCTrackQAHistos->Add(fHistSecRecMCPhiEtaPt[histoclass]);
1694
1695         fHistSecRecMCDCAMat[histoclass] = InitializeDCASpectrum("fHistSecRecDCAMat",histoclass);
1696         fSecRecMCTrackQAHistos->Add(fHistSecRecMCDCAMat[histoclass]);
1697
1698         fHistSecRecMCDCAWeak[histoclass] = InitializeDCASpectrum("fHistSecRecDCAWeak",histoclass);
1699         fSecRecMCTrackQAHistos->Add(fHistSecRecMCDCAWeak[histoclass]);
1700
1701         return kTRUE;
1702
1703 }
1704
1705 // -----------------------------------------------------------------------
1706 TH1F* AliAODTrackCutsDiHadronPID::InitializeAcceptedFilterBits(const char* name) {
1707
1708         // This histogram keeps track of the filtermask of all accepted tracks, projected
1709         // onto the requested filtermask. For example, we requested mask 2 or 4, then this
1710         // histogram will have three bins, 2, 4 and 6. Suppose now a track is found which has
1711         // mask 12, then since (12 & 2) = 0, (12 & 4) = 4, (12 & 6) = 4, the track will end up
1712         // in bin 4.
1713
1714         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1715         if (!fTestFilterMask || fFilterMask == 0) {cout << Form("%s -> ERROR: No filtermask requested.",__func__) << endl; return 0x0;}
1716
1717         // Determine the labels of the X axis.
1718         SetXaxisAcceptedFilterBits();
1719
1720         // Create histogram.
1721         TH1F* hout = new TH1F(name,"Filtermask of accepted track;Mask;N",fRelevantBitsArray->GetSize(),-0.5,fRelevantBitsArray->GetSize()-0.5);
1722         hout->SetDirectory(0);
1723
1724         // Set bin labels.
1725         TAxis* axistmp = hout->GetXaxis();
1726         for (Int_t iBin = 1; iBin <= axistmp->GetNbins(); ++iBin) {
1727                 axistmp->SetBinLabel(iBin, Form("%i",fRelevantBitsArray->At(iBin-1)));
1728         }
1729         return hout;
1730
1731 }
1732
1733 // -----------------------------------------------------------------------
1734 void AliAODTrackCutsDiHadronPID::SetXaxisAcceptedFilterBits() {
1735
1736         // Creates the axis for the AcceptedFilterBits histogram.
1737         // See exercise: "FindAllCombinations.C"
1738         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1739
1740         // Step 1: Find the largest bit, in the requested filtermask,
1741         Int_t largestBit = 0;
1742         Int_t baseArraySizeTmp = 0;
1743         Int_t fullArraySizeTmp = 0;
1744
1745         while (((Int_t)fFilterMask) >= (1<<(largestBit)) ) {
1746                 if ((((Int_t)fFilterMask)&(1<<largestBit))==(1<<largestBit)) {
1747                         fullArraySizeTmp += AliFunctionsDiHadronPID::Power(2, baseArraySizeTmp);
1748                         baseArraySizeTmp++;
1749                 }
1750                 largestBit++;
1751         }
1752         largestBit--;
1753
1754         // Step 2: Create and fill base array.
1755         const Int_t baseArraySize = baseArraySizeTmp;
1756         Int_t* baseArray = new Int_t[baseArraySize];
1757         for (Int_t ii = 0; ii < baseArraySize; ++ii) {baseArray[ii] = 0;}
1758
1759         Int_t iBaseArray = 0;
1760         for (Int_t iBit = 0; iBit <= largestBit; ++iBit) {
1761
1762                 if ((((Int_t)fFilterMask)&(1<<iBit))==(1<<iBit)) {
1763                         baseArray[iBaseArray] = (1<<iBit);
1764                         iBaseArray++;
1765                 }
1766
1767         }
1768
1769         // Step 3: Create and fill full array.
1770         const Int_t fullArraySize = fullArraySizeTmp;
1771         Int_t* fullArray = new Int_t[fullArraySize];
1772         fullArray[0] = baseArray[0];
1773         Int_t iFullArray = 1;
1774
1775         for (Int_t ii = 1; ii < baseArraySize; ++ii) {
1776                 Int_t range = (iFullArray + AliFunctionsDiHadronPID::Power(2,ii));
1777                 for (Int_t jj = iFullArray; jj < range; ++jj) {
1778
1779                         fullArray[jj] = baseArray[ii];
1780
1781                         // Add beginning part of the array:
1782                         if (jj!=iFullArray) {
1783                                 fullArray[jj] += fullArray[jj - iFullArray - 1];
1784                         }
1785
1786                 }
1787                 iFullArray += AliFunctionsDiHadronPID::Power(2,ii);
1788
1789         }
1790
1791         // Step 4: Convert to TArrayI object.
1792         fRelevantBitsArray = new TArrayI(fullArraySize, fullArray);
1793
1794         // Delete the temporary arrays.
1795         delete baseArray;
1796         delete fullArray;
1797
1798 }
1799
1800 // -----------------------------------------------------------------------
1801 TH1F* AliAODTrackCutsDiHadronPID::InitializePtSpectrum(const char* name, Int_t histoclass) {
1802
1803         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1804
1805         TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
1806                 Form("p_{T} Spectrum (%s);p_{T} (GeV/c);Count",fHistoName[histoclass].Data()),fNPtBins,fPtAxis);
1807
1808         hout->SetDirectory(0);
1809
1810         return hout;
1811
1812 }
1813
1814 // -----------------------------------------------------------------------
1815 TH2F* AliAODTrackCutsDiHadronPID::InitializeRecPtGenPt(const char* name, Int_t histoclass) {
1816
1817         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1818
1819         TH2F* hout = new TH2F(Form("%s%s",name,fHistoName[histoclass].Data()),
1820                 Form("p_{T} Rec vs Gen (%s);p_{T,gen} (GeV/c); p_{T,rec}",fHistoName[histoclass].Data()),
1821                 fNPtBins,fPtAxis,fNPtBins,fPtAxis);
1822
1823         hout->SetDirectory(0);
1824
1825         return hout;
1826
1827 }
1828
1829 // -----------------------------------------------------------------------
1830 TH3F* AliAODTrackCutsDiHadronPID::InitializePhiEtaPt(const char* name, Int_t histoclass) {
1831
1832         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1833
1834         TH3F* hout = AliHistToolsDiHadronPID::MakeHist3D(Form("%s%s",name,fHistoName[histoclass].Data()),
1835                 Form("Spectrum (%s);#phi;#eta;p_{T} (GeV/c)",fHistoName[histoclass].Data()),
1836                 fNPhiBins,0.,2.*TMath::Pi(),
1837                 fNEtaBins,-fMaxEta,fMaxEta,
1838                 fNPtBins, fPtAxis);
1839
1840         hout->SetDirectory(0);
1841
1842         return hout;
1843
1844 }
1845
1846 // -----------------------------------------------------------------------
1847 TH2F* AliAODTrackCutsDiHadronPID::InitializeDCASpectrum(const char* name, Int_t histoclass) {
1848
1849         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1850
1851         TH2F* hout = new TH2F(Form("%s%s",name,fHistoName[histoclass].Data()),
1852                 Form("DCA_{xy} (%s); p_{T} (GeV/c); DCA_{xy} (cm)",fHistoName[histoclass].Data()),
1853                 fNPtBins,fPtAxis,
1854                 300,-3.,3.);
1855
1856         hout->SetDirectory(0);
1857
1858         return hout;
1859
1860 }
1861
1862 // -----------------------------------------------------------------------
1863 TH1F* AliAODTrackCutsDiHadronPID::InitializeNTracksHisto(const char* name, Int_t histoclass) {
1864
1865         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1866
1867         TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
1868                 Form("Number of Accepted Tracks (%s);N%s;N_{event}",fHistoName[histoclass].Data(),fHistoLatex[histoclass].Data()),
1869                 100,0,4000);
1870
1871         hout->SetDirectory(0);
1872
1873         return hout;
1874
1875 }
1876
1877 // -----------------------------------------------------------------------
1878 TH1F* AliAODTrackCutsDiHadronPID::InitializeDCAxyHisto(const char* name, Int_t histoclass) {
1879
1880         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1881
1882         TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
1883                 Form("DCAxy (%s);DCAxy (cm);Count",fHistoName[histoclass].Data()),300,-15.,15.);
1884
1885         hout->SetDirectory(0);
1886
1887         return hout;
1888
1889 }
1890
1891 // -----------------------------------------------------------------------
1892 TH1F* AliAODTrackCutsDiHadronPID::InitializeDCAzHisto(const char* name, Int_t histoclass) {
1893
1894         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1895
1896         TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
1897         Form("DCAz (%s);DCAz (cm);Count",fHistoName[histoclass].Data()),300,-15.,15.);
1898
1899         hout->SetDirectory(0);
1900
1901         return hout;
1902
1903 }
1904
1905 // -----------------------------------------------------------------------
1906 TH3F* AliAODTrackCutsDiHadronPID::InitializeAcceptanceHisto(const char* /*name*/, Int_t /*histoclass*/) {
1907
1908         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1909         return 0x0;
1910
1911 }
1912
1913 // -----------------------------------------------------------------------
1914 TH3F* AliAODTrackCutsDiHadronPID::InitializePIDHisto(const char* name, Int_t histoclass, Int_t expspecies, Int_t ptclass) {
1915
1916         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1917
1918         TString PIDaxeslabel;
1919         if (fUseNSigmaOnPIDAxes) {PIDaxeslabel = ";n#sigma_{TOF};n#sigma_{TPC}";}
1920         else {PIDaxeslabel = ";#Delta t (ps);dE/dx (a.u.)";}
1921
1922         TH3F* hout = new TH3F(Form("%s%s%s%s",name,fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),fPtClassName[ptclass].Data()),
1923                 Form("PID %s (Exp: %s);p_{T} (GeV/c)%s",fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),PIDaxeslabel.Data()),
1924                 fNPtBinsPID[ptclass],fPtBoundaryPID[ptclass],fPtBoundaryPID[ptclass+1],
1925                 fTOFbins[ptclass][expspecies],
1926                 (fUseNSigmaOnPIDAxes ? fTOFLowerBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTOFStd : fTOFLowerBound[ptclass][expspecies]),
1927                 (fUseNSigmaOnPIDAxes ? fTOFUpperBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTOFStd : fTOFUpperBound[ptclass][expspecies]),
1928                 fTPCbins[ptclass][expspecies],
1929                 (fUseNSigmaOnPIDAxes ? fTPCLowerBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTPCStd : fTPCLowerBound[ptclass][expspecies]),
1930                 (fUseNSigmaOnPIDAxes ? fTPCUpperBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTPCStd : fTPCUpperBound[ptclass][expspecies]));
1931
1932         hout->SetDirectory(0);
1933
1934         return hout;
1935
1936 }
1937
1938 // -----------------------------------------------------------------------
1939 TH2F* AliAODTrackCutsDiHadronPID::InitializeTOFMismatchHisto(const char* name, Int_t histoclass, Int_t expspecies, Int_t ptclass) {
1940
1941         // Is basically the same as InitializeTOFHisto -> CAN BE REMOVED!
1942         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1943
1944         TString PIDaxeslabel;
1945         if (fUseNSigmaOnPIDAxes) {PIDaxeslabel = ";n#sigma_{TOF}";}
1946         else {PIDaxeslabel = ";#Delta t (ps)";}
1947
1948         TH2F* hout = new TH2F(Form("%s%s%s%s",name,fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),fPtClassName[ptclass].Data()),
1949                 Form("TOF Mismatch %s (Exp: %s);p_{T} (GeV/c)%s",fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),PIDaxeslabel.Data()),
1950                 fNPtBinsPID[ptclass],fPtBoundaryPID[ptclass],fPtBoundaryPID[ptclass+1],
1951                 fTOFbins[ptclass][expspecies],
1952                 (fUseNSigmaOnPIDAxes ? fTOFLowerBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTOFStd : fTOFLowerBound[ptclass][expspecies]),
1953                 (fUseNSigmaOnPIDAxes ? fTOFUpperBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTOFStd : fTOFUpperBound[ptclass][expspecies]));
1954
1955         hout->SetDirectory(0);
1956
1957         return hout;
1958
1959 }
1960
1961 // -----------------------------------------------------------------------
1962 TH2F* AliAODTrackCutsDiHadronPID::InitializeTOFHisto(const char* name, Int_t histoclass, Int_t expspecies, Int_t ptclass) {
1963
1964         if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1965
1966         TString PIDaxeslabel;
1967         if (fUseNSigmaOnPIDAxes) {PIDaxeslabel = ";n#sigma_{TOF}";}
1968         else {PIDaxeslabel = ";#Delta t (ps)";}
1969
1970         TH2F* hout = new TH2F(Form("%s%s%s%s",name,fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),fPtClassName[ptclass].Data()),
1971                 Form("TOF %s (Exp: %s);p_{T} (GeV/c)%s",fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),PIDaxeslabel.Data()),
1972                 fNPtBinsPID[ptclass],fPtBoundaryPID[ptclass],fPtBoundaryPID[ptclass+1],
1973                 fTOFbins[ptclass][expspecies],
1974                 (fUseNSigmaOnPIDAxes ? fTOFLowerBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTOFStd : fTOFLowerBound[ptclass][expspecies]),
1975                 (fUseNSigmaOnPIDAxes ? fTOFUpperBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTOFStd : fTOFUpperBound[ptclass][expspecies]));
1976
1977         hout->SetDirectory(0);
1978
1979         return hout;
1980
1981 }