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