PWGLF/SPECTRA/ChargedHadrons/multPbPb/AliAnalysisMultPbTrackHistoManager.cxx

   1 #include "AliAnalysisMultPbTrackHistoManager.h"
   2 #include "AliLog.h"
   3 #include "TH1.h"
   4 #include "TH3D.h"
   5 #include "TH1I.h"
   6 #include "TROOT.h"
   7 #include "TMCProcess.h"
   8 #include "AliMCParticle.h"
   9 #include "TMap.h"
  10 #include <iostream>
  11 #include "TH2D.h"
  12 #include "TDatabasePDG.h"
  13 #include "TParameter.h"
  14
  15 using namespace std;
  16
  17 ClassImp(AliAnalysisMultPbTrackHistoManager)
  18
  19 const char * AliAnalysisMultPbTrackHistoManager::kStatStepNames[]     = { "All Events", "After centrality selection",  "After physics Selection", "With Vertex (quality cuts)", "After ZDC cut", "After vertex range cut (10 cm)" };
  20 const char * AliAnalysisMultPbTrackHistoManager::kHistoPtEtaVzNames[] = { "hGenPtEtaVz", "hRecPtEtaVz", "hRecPtEtaVzPrim",
  21                                                                           "hRecPtEtaVzSecWeak", "hRecPtEtaVzSecMaterial", "hRecPtEtaVzFake"};
  22 const char * AliAnalysisMultPbTrackHistoManager::kHistoDCANames[]     = { "hGenDCA", "hRecDCA", "hRecDCAPrim", "hRecDCASecWeak","hRecDCASecMaterial", "hRecDCAFake"};
  23 const char * AliAnalysisMultPbTrackHistoManager::kHistoPrefix[]     = { "hGen", "hRec", "hRecPrim", "hRecSecWeak","hRecSecMaterial", "hRecFake", "hRecHighestMeanPt", "hRecLowestMeanPt"};
  24 const char * AliAnalysisMultPbTrackHistoManager::kSpeciesName[]     = { "pi+", "K+", "p", "l+",  "pi-", "K-", "barp", "l-", "Others"};
  25
  26
  27 AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager() : AliHistoListWrapper(), fHNameSuffix(""), fParticleSpecies(0){
  28   // standard ctor
  29
  30 }
  31
  32 AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager(const char * name, const char * title): AliHistoListWrapper(name,title), fHNameSuffix(""), fParticleSpecies(0)  {
  33   //named ctor
  34 };
  35
  36 AliAnalysisMultPbTrackHistoManager::AliAnalysisMultPbTrackHistoManager(const AliAnalysisMultPbTrackHistoManager& obj) : AliHistoListWrapper (obj) {
  37   // copy ctor
  38   AliError("Not Implemented");
  39 };
  40
  41 AliAnalysisMultPbTrackHistoManager::~AliAnalysisMultPbTrackHistoManager() {
  42   // dtor
  43
  44 }
  45
  46 TH2D * AliAnalysisMultPbTrackHistoManager::GetHistoElectronCutQA() {
  47   // Get a p vs dE/dx plot, to check the histos we are rejecting
  48   TString name = "histoElectronCutQA";
  49
  50   TH2D * h = (TH2D*) GetHisto(name);
  51
  52   if (!h) {
  53     AliInfo(Form("Booking %s", (name+fHNameSuffix).Data()));
  54     Bool_t oldStatus = TH1::AddDirectoryStatus();
  55     TH1::AddDirectory(kFALSE);
  56     h = new TH2D (name+fHNameSuffix, name+fHNameSuffix, 100, 0, 3, 50, 0, 200);
  57     fList->Add(h);
  58     TH1::AddDirectory(oldStatus);
  59
  60   }
  61
  62   return h;
  63 }
  64
  65
  66 TH3D * AliAnalysisMultPbTrackHistoManager::GetHistoPtEtaVz(Histo_t id, Int_t particle) {
  67   // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
  68
  69   TString name = kHistoPtEtaVzNames[id];
  70
  71   if(particle >= 0) name += kSpeciesName[particle];
  72
  73   TH3D * h = (TH3D*) GetHisto(name.Data());
  74   if (!h) {
  75     h = BookHistoPtEtaVz(name.Data(), Form("Pt Eta Vz distribution (%s)",kHistoPtEtaVzNames[id]));
  76   }
  77
  78   return h;
  79
  80 }
  81
  82 TH2D * AliAnalysisMultPbTrackHistoManager::GetHistoDCA(Histo_t id) {
  83   // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
  84
  85   TH2D * h = (TH2D*) GetHisto(kHistoDCANames[id]);
  86   if (!h) {
  87     h = BookHistoDCA(kHistoDCANames[id], Form("DCA vs pt distribution (%s)",kHistoDCANames[id]));
  88   }
  89
  90   return h;
  91
  92 }
  93
  94 TH2D * AliAnalysisMultPbTrackHistoManager::GetHistoV0vsNtracks(Histo_t id) {
  95   // Returns a histo of V0 vs ntracks
  96
  97
  98   TString name = TString(kHistoPrefix[id])+"_V0vsNtracks";
  99
 100   TH2D * h =  (TH2D*) GetHisto(name);
 101   if (!h) {
 102     name+=fHNameSuffix;
 103     Bool_t oldStatus = TH1::AddDirectoryStatus();
 104     TH1::AddDirectory(kFALSE);
 105
 106     AliInfo(Form("Booking histo %s",name.Data()));
 107
 108     h = new TH2D (name.Data(), Form("V0 vs Ntracks (%s)",kHistoPrefix[id]), 300,0,3000, 300, 0, 20000);
 109     h->SetXTitle("N_{Tracks}");
 110     h->SetYTitle("V0 Amplitude");
 111     TH1::AddDirectory(oldStatus);
 112     fList->Add(h);
 113
 114
 115   }
 116   return h;
 117
 118
 119 }
 120
 121
 122 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoMult(Histo_t id) {
 123   // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
 124
 125   TString name = kHistoPrefix[id];
 126   name += "Mult";
 127   TH1D * h = (TH1D*) GetHisto(name.Data());
 128   if (!h) {
 129     h = BookHistoMult(name.Data(), Form("Multiplicity distribution (%s)",kHistoPrefix[id]));
 130   }
 131
 132   return h;
 133
 134 }
 135
 136 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoPt (Histo_t id,
 137                                                        Float_t minEta, Float_t maxEta,
 138                                                        Float_t minVz, Float_t maxVz,
 139                                                        Bool_t scaleWidth) {
 140   // Returns a projection of the 3D histo pt/eta/vz.
 141   // WARNING: since that is a histo, the requested range will be discretized to the binning.
 142   // Always avoids under (over) flows
 143   // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
 144
 145   TH3D * h3D = GetHistoPtEtaVz(id);
 146
 147   // Get range in terms of bin numners.  If the float range is
 148   // less than -11111 take the range from the first to the last bin (i.e. no
 149   // under/over-flows)
 150
 151   Int_t min1 = minEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(minEta);
 152   Int_t min2  = minVz  < -11111 ? -1 : h3D ->GetZaxis()->FindBin(minVz) ;
 153
 154   Int_t max1 = maxEta  < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
 155   Int_t max2  = maxVz  < -11111 ? h3D->GetNbinsZ() : h3D ->GetZaxis()->FindBin(maxVz -0.00001);
 156   // Int_t max1 = maxEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
 157   // Int_t max2  = maxVz  < -11111 ? -1 : h3D ->GetZaxis()->FindBin(maxVz -0.00001);
 158
 159
 160   TString hname = h3D->GetName();
 161   hname = hname +  "_pt_" + long (min1)  + "_" + long(max1) + "_" + long (min2)  + "_" + long(max2);
 162
 163
 164   if (gROOT->FindObjectAny(hname.Data())){
 165     AliError(Form("An object called %s already exists,adding suffix",hname.Data()));
 166     hname += "_2";
 167   }
 168
 169   TH1D * h = h3D->ProjectionX(hname.Data(), min1, max1, min2, max2, "E");
 170   if(scaleWidth) h ->Scale(1.,"width");
 171
 172   return h;
 173
 174 }
 175
 176
 177 TH2D * AliAnalysisMultPbTrackHistoManager::GetHistoPtVz (Histo_t id,
 178                                                          Float_t minEta, Float_t maxEta,
 179                                                          Bool_t scaleWidth) {
 180   // Returns a projection of the 3D histo pt/eta/vz.
 181   // WARNING: since that is a histo, the requested range will be discretized to the binning.
 182   // Always avoids under (over) flows
 183   // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
 184
 185   // FIXME: what do I do here for the scaling?
 186
 187   TH3D * h3D = GetHistoPtEtaVz(id);
 188
 189   // Get range in terms of bin numners.  If the float range is
 190   // less than -11111 take the range from the first to the last bin (i.e. no
 191   // under/over-flows)
 192   Int_t min1 = minEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(minEta);
 193   Int_t max1 = maxEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
 194
 195
 196   TString hname = h3D->GetName();
 197   hname = hname +  "_ptvz_" + long (min1)  + "_" + long(max1);
 198
 199   if (gROOT->FindObjectAny(hname.Data())){
 200     AliError(Form("An object called %s already exists,adding suffix",hname.Data()));
 201     hname += "_2";
 202   }
 203
 204   h3D->GetYaxis()->SetRange(min1,max1);
 205
 206   TH2D * h =  (TH2D*) h3D->Project3D("zxe");
 207   if(scaleWidth) h ->Scale(1.,"width");
 208
 209   return h;
 210
 211 }
 212
 213
 214 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoVz (Histo_t id,
 215                                                        Float_t minPt, Float_t maxPt,
 216                                                        Float_t minEta, Float_t maxEta,
 217                                                        Bool_t scaleWidth) {
 218   // Returns a projection of the 3D histo pt/eta/vz.
 219   // WARNING: since that is a histo, the requested range will be discretized to the binning.
 220   // Always avoids under (over) flows
 221   // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
 222
 223   TH3D * h3D = GetHistoPtEtaVz(id);
 224
 225   // Get range in terms of bin numners.  If the float range is
 226   // less than -11111 take the range from the first to the last bin (i.e. no
 227   // under/over-flows)
 228   Int_t min1  = minPt  < -11111 ? 1 : h3D ->GetXaxis()->FindBin(minPt) ;
 229   Int_t min2  = minEta < -11111 ? 1 : h3D ->GetYaxis()->FindBin(minEta);
 230
 231   Int_t max1  = maxPt  < -11111 ? h3D->GetNbinsX() : h3D ->GetXaxis()->FindBin(maxPt -0.00001);
 232   Int_t max2  = maxEta < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
 233
 234
 235   TString hname = h3D->GetName();
 236   hname = hname +  "_Vz_" + long (min1)  + "_" + long(max1) + "_" + long (min2)  + "_" + long(max2);
 237
 238   if (gROOT->FindObjectAny(hname.Data())){
 239     AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
 240     hname+="_2";
 241   }
 242
 243   TH1D * h = h3D->ProjectionZ(hname.Data(), min1, max1, min2, max2, "E");
 244   if(scaleWidth) h ->Scale(1.,"width");
 245   return h;
 246
 247
 248 }
 249
 250 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoEta (Histo_t id,
 251                                                         Float_t minPt, Float_t maxPt,
 252                                                         Float_t minVz, Float_t maxVz,
 253                                                         Bool_t scaleWidth) {
 254   // Returns a projection of the 3D histo pt/eta/vz.
 255   // WARNING: since that is a histo, the requested range will be discretized to the binning.
 256   // Always avoids under (over) flows
 257   // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
 258
 259   TH3D * h3D = GetHistoPtEtaVz(id);
 260
 261   // Get range in terms of bin numners.  If the float range is
 262   // less than -11111 take the range from the first to the last bin (i.e. no
 263   // under/over-flows)
 264   Int_t min1 = minPt < -11111 ? 1 : h3D ->GetXaxis()->FindBin(minPt) ;
 265   Int_t min2 = minVz < -11111 ? 1 : h3D ->GetYaxis()->FindBin(minVz);
 266
 267   Int_t max1 = maxPt < -11111 ? h3D->GetNbinsX() : h3D ->GetXaxis()->FindBin(maxPt -0.00001);
 268   Int_t max2 = maxVz < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxVz-0.00001);
 269
 270   TString hname = h3D->GetName();
 271   hname = hname +  "_Eta_" + long (min1)  + "_" + long(max1) + "_" + long (min2)  + "_" + long(max2);
 272
 273   if (gROOT->FindObjectAny(hname.Data())){
 274     AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
 275     hname+="_2";
 276   }
 277
 278   TH1D * h = h3D->ProjectionY(hname.Data(), min1, max1, min2, max2, "E");
 279   if(scaleWidth) h ->Scale(1.,"width");
 280   return h;
 281 }
 282
 283 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoProcess(Histo_t id) {
 284
 285   // Returns histogram with particle specties
 286
 287   TString name = TString(kHistoPrefix[id])+"_Process";
 288
 289   TH1D * h =  (TH1D*) GetHisto(name);
 290   if (!h) {
 291     name+=fHNameSuffix;
 292     Bool_t oldStatus = TH1::AddDirectoryStatus();
 293     TH1::AddDirectory(kFALSE);
 294
 295     AliInfo(Form("Booking histo %s",name.Data()));
 296
 297     h = new TH1D (name.Data(), Form("Particle production process (%s)",kHistoPrefix[id]), kPNoProcess+1, -0.5, kPNoProcess+1-0.5);
 298     Int_t nbin = kPNoProcess+1;
 299     for(Int_t ibin = 0; ibin < nbin; ibin++){
 300       h->GetXaxis()->SetBinLabel(ibin+1,TMCProcessName[ibin]);
 301     }
 302     TH1::AddDirectory(oldStatus);
 303     fList->Add(h);
 304
 305
 306   }
 307   return h;
 308
 309
 310 }
 311
 312
 313 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoSpecies(Histo_t id) {
 314
 315   // Returns histogram with particle specties
 316
 317   TString name = TString(kHistoPrefix[id])+"_Species";
 318
 319   TH1D * h =  (TH1D*) GetHisto(name);
 320   if (!h) {
 321     name+=fHNameSuffix;
 322     Bool_t oldStatus = TH1::AddDirectoryStatus();
 323     TH1::AddDirectory(kFALSE);
 324
 325     AliInfo(Form("Booking histo %s",name.Data()));
 326
 327     h = new TH1D (name.Data(), Form("Particle species (%s)",kHistoPrefix[id]), kNPart+1, -0.5, kNPart+1-0.5);
 328     Int_t nbin = kNPart+1;
 329     for(Int_t ibin = 0; ibin < nbin; ibin++){
 330       h->GetXaxis()->SetBinLabel(ibin+1,kSpeciesName[ibin]);
 331     }
 332     TH1::AddDirectory(oldStatus);
 333     fList->Add(h);
 334
 335
 336   }
 337   return h;
 338
 339
 340 }
 341
 342
 343 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoMeanPt(Histo_t id) {
 344   // mean pt computed event by event
 345   TString name = TString(kHistoPrefix[id])+"_MeanPt";
 346
 347   TH1D * h =  (TH1D*) GetHisto(name);
 348   if (!h) {
 349     name+=fHNameSuffix;
 350     Bool_t oldStatus = TH1::AddDirectoryStatus();
 351     TH1::AddDirectory(kFALSE);
 352
 353     AliInfo(Form("Booking histo %s",name.Data()));
 354
 355     h = new TH1D (name.Data(), Form("Pt Event by Event(%s)",kHistoPrefix[id]), 100, 0, 1);
 356     h->Sumw2();
 357     TH1::AddDirectory(oldStatus);
 358     fList->Add(h);
 359
 360
 361   }
 362   return h;
 363
 364 }
 365
 366 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoPtEvent(Histo_t id) {
 367
 368   // Returns of dNdpt, used to compute mean pt event by event
 369
 370   TString name = TString(kHistoPrefix[id])+"_PtEvent";
 371
 372   TH1D * h =  (TH1D*) GetHisto(name);
 373   if (!h) {
 374     name+=fHNameSuffix;
 375     Bool_t oldStatus = TH1::AddDirectoryStatus();
 376     TH1::AddDirectory(kFALSE);
 377
 378     AliInfo(Form("Booking histo %s",name.Data()));
 379     const Int_t nptbins = 68;
 380     const Double_t binsPt[] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0};
 381
 382     h = new TH1D (name.Data(), Form("Pt Event by Event(%s)",kHistoPrefix[id]), nptbins,binsPt);
 383     h->SetBit(kKeepMaxMean,1);
 384     if(id == kHistoRecLowestMeanPt ) h->SetBit(kKeepMinMean,1);
 385     h->Sumw2();
 386     TH1::AddDirectory(oldStatus);
 387     fList->Add(h);
 388
 389
 390   }
 391   return h;
 392
 393
 394 }
 395
 396
 397 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoVzEvent(Histo_t id) {
 398
 399   // Returns histogram with Vz of the event
 400
 401   TString name = TString(kHistoPrefix[id])+"_VzEvent";
 402
 403   TH1D * h =  (TH1D*) GetHisto(name);
 404   if (!h) {
 405     name+=fHNameSuffix;
 406     Bool_t oldStatus = TH1::AddDirectoryStatus();
 407     TH1::AddDirectory(kFALSE);
 408
 409     AliInfo(Form("Booking histo %s",name.Data()));
 410
 411     h = new TH1D (name.Data(), Form("Vz of the event (%s)",kHistoPrefix[id]), 10, -10, 10);
 412     h->Sumw2();
 413     h->SetXTitle("V_{z}^{event} (cm)");
 414     TH1::AddDirectory(oldStatus);
 415     fList->Add(h);
 416
 417
 418   }
 419   return h;
 420
 421
 422 }
 423
 424
 425
 426 TH1I * AliAnalysisMultPbTrackHistoManager::GetHistoStats() {
 427   // Returns histogram with event statistiscs (processed events at each step)
 428
 429   TH1I * h =  (TH1I*) GetHisto("hStats");
 430   if (!h) h = BookHistoStats();
 431   return h;
 432
 433 }
 434
 435
 436
 437 TH1 * AliAnalysisMultPbTrackHistoManager::GetHisto(const char * name) {
 438   //Search list for histo
 439   // TODO: keep track of histo id rather than searching by name?
 440   return (TH1*) fList->FindObject(TString(name)+fHNameSuffix);
 441
 442 }
 443
 444 void AliAnalysisMultPbTrackHistoManager::ScaleHistos(Double_t nev, Option_t * option) {
 445   // Scales all histos in the list for nev
 446   // option can be used to pass further options to TH1::Scale
 447   TH1 * h = 0;
 448   TIter iter = fList->MakeIterator();
 449   while ((h = (TH1*) iter.Next())) {
 450     if (!h->InheritsFrom("TH1")) {
 451       AliFatal (Form("%s does not inherits from TH1, cannot scale",h->GetName()));
 452     }
 453     if (h->InheritsFrom("TH1I")) {
 454       AliInfo (Form("Not scaling integer histo %s",h->GetName()));
 455       continue;
 456     }
 457     AliInfo(Form("Scaling %s, nev %2.2f", h->GetName(), nev));
 458
 459     h->Scale(1./nev,option);
 460   }
 461
 462 }
 463
 464 TH3D * AliAnalysisMultPbTrackHistoManager::BookHistoPtEtaVz(const char * name, const char * title) {
 465   // Books a 3D histo of Pt/eta/vtx
 466   // TODO: make the binning settable, variable binning?
 467
 468   Bool_t oldStatus = TH1::AddDirectoryStatus();
 469   TH1::AddDirectory(kFALSE);
 470
 471   TString hname = name;
 472   hname+=fHNameSuffix;
 473
 474   AliInfo(Form("Booking %s",hname.Data()));
 475
 476   // Binning from Jacek task
 477   // const Int_t nptbins = 49;
 478   // const Double_t binsPt[] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0};//,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0};
 479   const Int_t nptbins = 68;
 480   const Double_t binsPt[] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0};
 481
 482   const Int_t netabins=20;
 483   Double_t binsEta[netabins+1];
 484   Float_t minEta = -1;
 485   Float_t maxEta =  1;
 486   Float_t etaStep = (maxEta-minEta)/netabins;
 487   for(Int_t ibin = 0; ibin < netabins; ibin++){
 488     binsEta[ibin]   = minEta + ibin*etaStep;
 489     binsEta[ibin+1] = minEta + ibin*etaStep + etaStep;
 490   }
 491
 492   const Int_t nvzbins=10;
 493   Double_t binsVz[nvzbins+1];
 494   Float_t minVz = -10;
 495   Float_t maxVz =  10;
 496   Float_t vzStep = (maxVz-minVz)/nvzbins;
 497   for(Int_t ibin = 0; ibin < nvzbins; ibin++){
 498     binsVz[ibin]   = minVz + ibin*vzStep;
 499     binsVz[ibin+1] = minVz + ibin*vzStep + vzStep;
 500   }
 501
 502
 503   TH3D * h = new TH3D (hname,title,
 504                        nptbins,  binsPt,
 505                        netabins, binsEta,
 506                        nvzbins,  binsVz
 507                        );
 508
 509   h->SetXTitle("p_{T} (GeV)");
 510   h->SetYTitle("#eta");
 511   h->SetZTitle("V_{z}^{tracks} (cm)");
 512   h->Sumw2();
 513
 514   TH1::AddDirectory(oldStatus);
 515   fList->Add(h);
 516   return h;
 517 }
 518
 519 TH2D * AliAnalysisMultPbTrackHistoManager::BookHistoDCA(const char * name, const char * title) {
 520   // Books a DCA histo
 521
 522   const Int_t nptbins = 68;
 523   const Double_t binsPt[] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,18.0,20.0,22.0,24.0,26.0,28.0,30.0,32.0,34.0,36.0,40.0,45.0,50.0};
 524   const Int_t ndcabins=100;
 525   Double_t binsDca[ndcabins+1];
 526   Float_t minDca = -3;
 527   Float_t maxDca =  3;
 528   //  const Double_t binsDCA[] = {-3,-2.8,-2.6,-2.4,-2.2,-2.0,-1.9,};
 529   Float_t dcaStep = (maxDca-minDca)/ndcabins;
 530   for(Int_t ibin = 0; ibin < ndcabins; ibin++){
 531     binsDca[ibin]   = minDca + ibin*dcaStep;
 532     binsDca[ibin+1] = minDca + ibin*dcaStep + dcaStep;
 533   }
 534
 535
 536
 537
 538   Bool_t oldStatus = TH1::AddDirectoryStatus();
 539   TH1::AddDirectory(kFALSE);
 540
 541   TString hname = name;
 542   hname+=fHNameSuffix;
 543
 544   AliInfo(Form("Booking %s",hname.Data()));
 545
 546
 547 #if defined WEIGHTED_DCA
 548   TH1D * h = new TH1D (hname,title, 200,0,200);
 549   h->SetXTitle("#Delta DCA");
 550 #elif defined TRANSVERSE_DCA
 551   TH2D * h = new TH2D (hname,title, ndcabins,binsDca,nptbins,binsPt);
 552   h->SetYTitle("p_{T} (GeV)");
 553   h->SetXTitle("d_{0} r#phi (cm)");
 554 #endif
 555   h->Sumw2();
 556
 557   fList->Add(h);
 558
 559   TH1::AddDirectory(oldStatus);
 560   return h;
 561 }
 562 TH1D * AliAnalysisMultPbTrackHistoManager::BookHistoMult(const char * name, const char * title) {
 563   // Books a multiplicity histo
 564
 565   Bool_t oldStatus = TH1::AddDirectoryStatus();
 566   TH1::AddDirectory(kFALSE);
 567
 568   TString hname = name;
 569   hname+=fHNameSuffix;
 570
 571   AliInfo(Form("Booking %s",hname.Data()));
 572
 573
 574   TH1D * h = new TH1D (hname,title, 600, 0,6000);
 575
 576   h->SetXTitle("N tracks");
 577   h->Sumw2();
 578
 579   fList->Add(h);
 580
 581   TH1::AddDirectory(oldStatus);
 582   return h;
 583 }
 584
 585 TH1I * AliAnalysisMultPbTrackHistoManager::BookHistoStats() {
 586   // Books histogram with event statistiscs (processed events at each step)
 587
 588   Bool_t oldStatus = TH1::AddDirectoryStatus();
 589   TH1::AddDirectory(kFALSE);
 590
 591   AliInfo(Form("Booking Stat histo"));
 592
 593   TH1I * h = new TH1I (TString("hStats")+fHNameSuffix, "Number of processed events", kNStatBins, -0.5, kNStatBins-0.5);
 594   for(Int_t istatbin = 0; istatbin < kNStatBins; istatbin++){
 595     h->GetXaxis()->SetBinLabel(istatbin+1,kStatStepNames[istatbin]);
 596   }
 597   TH1::AddDirectory(oldStatus);
 598   fList->Add(h);
 599   return h;
 600 }
 601
 602 Int_t AliAnalysisMultPbTrackHistoManager::GetLocalParticleID(AliMCParticle * part) {
 603   // returns the local code (Part_t)
 604
 605   Int_t pdgcode = part->PdgCode();
 606   switch(pdgcode) {
 607
 608   case 211:
 609     return kPartPiPlus;
 610     break;
 611   case -211:
 612     return kPartPiMinus;
 613     break;
 614   case 2212:
 615     return kPartP;
 616     break;
 617   case -2212:
 618     return kPartPBar;
 619     break;
 620   case 321:
 621     return kPartKPlus;
 622     break;
 623   case -321:
 624     return kPartKMinus;
 625     break;
 626   case -11:
 627     return kPartLMinus;
 628     break;
 629   case 11:
 630     return kPartLPlus;
 631     break;
 632   case -13:
 633     return kPartLMinus;
 634     break;
 635   case 13:
 636     return kPartLPlus;
 637     break;
 638   default:
 639     return kPartOther;
 640   }
 641 }
 642
 643
 644 void AliAnalysisMultPbTrackHistoManager::FillSpeciesMap(Int_t pdgCode) {
 645
 646   // Fills map of species vs ID
 647   if(!fParticleSpecies) fParticleSpecies = new TMap;
 648
 649   TObjString * key = new TObjString(Form("%s",TDatabasePDG::Instance()->GetParticle(pdgCode)->GetName()));
 650   TParameter<Long64_t>* count = dynamic_cast<TParameter<Long64_t>*> (fParticleSpecies->GetValue(key));
 651   if (!count)
 652   {
 653     count = new TParameter<Long64_t>(key->String().Data(), 0);
 654     fParticleSpecies->Add(key, count);
 655   }
 656   count->SetVal(count->GetVal() + 1);
 657
 658 }
 659
 660
 661 Long64_t AliAnalysisMultPbTrackHistoManager::Merge(TCollection* list)
 662 {
 663   // Merge a list of AliHistoListWrapper objects with this.
 664   // Returns the number of merged objects (including this).
 665
 666   // We have to make sure that all the list contain the same histos in
 667   // the same order. We thus also have to sort the list (sorting is
 668   // done by name in TList).
 669
 670   AliInfo("Merging");
 671   if (!list)
 672
 673   if (list->IsEmpty())
 674     return 1;
 675
 676   TIterator* iter = list->MakeIterator();
 677   TObject* obj;
 678
 679   // collections of all histograms
 680   TList collections;
 681
 682   Int_t count = 0;
 683
 684   while ((obj = iter->Next())) {
 685     cout << "1" << endl;
 686     //    AliHistoListWrapper* entry = dynamic_cast<AliHistoListWrapper*> (obj);
 687     AliAnalysisMultPbTrackHistoManager* entry = dynamic_cast<AliAnalysisMultPbTrackHistoManager*> (obj);
 688     cout << "2 " << endl;
 689     if (entry == 0)
 690       continue;
 691     cout << "3 " << entry->GetName() << endl;
 692
 693     // Merge map fParticleSpecies
 694     if(!entry->fParticleSpecies){
 695       AliInfo("Cannot get fParticleSpecies");
 696       continue;
 697     }
 698     TIterator* iter2 = entry->fParticleSpecies->MakeIterator();
 699     TObjString* obj2 = 0;
 700     cout << "4" << endl;
 701     while ((obj2 = dynamic_cast<TObjString*> (iter2->Next())))
 702     {
 703       cout << "5" << endl;
 704       TParameter<Long64_t>* param2 = static_cast<TParameter<Long64_t>*> (entry->fParticleSpecies->GetValue(obj2));
 705
 706       TParameter<Long64_t>* param1 = dynamic_cast<TParameter<Long64_t>*> (fParticleSpecies->GetValue(obj2));
 707       cout << " - (other)" << param2->GetName() << " " << param2->GetVal();
 708       if (param1)
 709       {
 710         cout << " (this) " << param1->GetName() << " " << param1->GetVal() << endl;
 711         param1->SetVal(param1->GetVal() + param2->GetVal());
 712       }
 713       else
 714       {
 715         cout << "" << endl;
 716         param1 = dynamic_cast<TParameter<Long64_t>*> (param2->Clone());
 717         fParticleSpecies->Add(new TObjString(obj2->String()), param1);
 718       }
 719     }
 720     delete iter2;
 721
 722   }
 723
 724   // merge histograms
 725   // We need a separate loop as this one could be restarted if needed.
 726   iter->Reset();
 727   while ((obj = iter->Next())) {
 728     AliAnalysisMultPbTrackHistoManager* entry = dynamic_cast<AliAnalysisMultPbTrackHistoManager*> (obj);
 729     if (entry == 0)
 730       continue;
 731
 732     Bool_t foundDiffinThisIterStep = kFALSE;
 733
 734     //    Printf("%d - %s",count, obj->GetName());
 735
 736     TList * hlist = entry->GetList();
 737
 738     // Check if all histos in this fList are also in the one from entry and viceversa
 739     // Use getters to automatically book non defined histos
 740
 741     Bool_t areListsDifferent=kTRUE;
 742     Int_t iloop = 0;
 743     Int_t max_loops = hlist->GetSize() + fList->GetSize(); // In the worst case all of the histos will be different...
 744     while(areListsDifferent) {
 745       if(iloop>max_loops) AliFatal("Infinite Loop?");
 746       iloop++;
 747       // sort
 748       hlist->Sort();
 749       fList->Sort();
 750       // loop over the largest
 751       TObject * hist =0;
 752       TIterator * iterlist = 0;
 753       TList * thislist  = 0; // the list over which I'm iterating
 754       TList * otherlist = 0; // the other
 755
 756       if (hlist->GetSize() >= fList->GetSize()) {
 757         thislist  = hlist;
 758         otherlist = fList;
 759       }
 760       else{
 761         thislist  = fList;
 762         otherlist = hlist;
 763       }
 764       iterlist = thislist->MakeIterator();
 765
 766       while ((hist= iterlist->Next())){
 767         if(!otherlist->FindObject(hist->GetName())){
 768           AliInfo(Form("Adding object %s",hist->GetName()));
 769           TH1 * hclone =  (TH1*) hist->Clone();
 770           if (!hclone->InheritsFrom("TH1")) AliFatal(Form("Found a %s. This class only supports objects inheriting from TH1",hclone->ClassName()));
 771           hclone->Reset();
 772           otherlist->Add(hclone);
 773           foundDiffinThisIterStep=kTRUE;
 774         }
 775       }
 776
 777       // re-sort before checking
 778       hlist->Sort();
 779       fList->Sort();
 780
 781       // check if everything is fine
 782       areListsDifferent=kFALSE;
 783       if (hlist->GetSize() == fList->GetSize()) {
 784         Int_t nhist =  fList->GetSize();
 785         for(Int_t ihist = 0; ihist < nhist; ihist++){
 786           if(strcmp(fList->At(ihist)->GetName(),hlist->At(ihist)->GetName())) areListsDifferent = kTRUE;
 787         }
 788       } else {
 789         areListsDifferent=kTRUE;
 790       }
 791     }
 792
 793     // last check: if something is not ok die loudly
 794     if (hlist->GetSize() != fList->GetSize()) {
 795       AliFatal("Mismatching size!");
 796     }
 797     Int_t nhist =  fList->GetSize();
 798     for(Int_t ihist = 0; ihist < nhist; ihist++){
 799       if(strcmp(fList->At(ihist)->GetName(),hlist->At(ihist)->GetName())){
 800         AliFatal(Form("Mismatching histos: %s -> %s", fList->At(ihist)->GetName(),hlist->At(ihist)->GetName()));
 801       } else {
 802         // Specific merging strategies, according to the bits...
 803         if (fList->At(ihist)->TestBits(kKeepMinMean|kKeepMaxMean)){
 804           TH1 * h1 = (TH1*) fList->At(ihist);
 805           TH1 * h2 = (TH1*) hlist->At(ihist);
 806           if (h1->GetEntries()>0 && h2->GetEntries()>0) {
 807             if (h1->TestBit(kKeepMinMean)) {
 808               AliInfo(Form("Keeping only the histo with the lowest mean [%s][%f][%f]",h1->GetName(), h1->GetMean(), h2->GetMean()));
 809               if(h1->GetMean() > h2->GetMean()) h1->Reset();
 810               else h2->Reset();
 811             }
 812             if (h1->TestBit(kKeepMaxMean)) {
 813               AliInfo(Form("Keeping only the histo with the highest mean [%s][%f][%f]",h1->GetName(), h1->GetMean(), h2->GetMean()));
 814               if(h2->GetMean() > h1->GetMean()) h1->Reset();
 815               else h2->Reset();
 816             }
 817           }
 818         }
 819       }
 820     }
 821
 822     if (foundDiffinThisIterStep){
 823       iter->Reset(); // We found a difference: previous lists could
 824                      // also be affected... We start from scratch
 825       collections.Clear();
 826       count = 0;
 827     }
 828     else {
 829
 830       collections.Add(hlist);
 831
 832       count++;
 833     }
 834   }
 835
 836   fList->Merge(&collections);
 837
 838   delete iter;
 839
 840   AliInfo("Merged");
 841
 842   return count+1;
 843 }