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 physics Selection", "After centrality selection", "With Vertex" };
  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"};
  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 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoDCA(Histo_t id) {
  61   // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
  62
  63   TH1D * h = (TH1D*) GetHisto(kHistoDCANames[id]);
  64   if (!h) {
  65     h = BookHistoDCA(kHistoDCANames[id], Form("Pt Eta Vz distribution (%s)",kHistoDCANames[id]));
  66   }
  67
  68   return h;
  69
  70 }
  71
  72 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoMult(Histo_t id) {
  73   // Returns a 3D histo of Pt/eta/vtx. It it does not exist, books it.
  74
  75   TString name = kHistoPrefix[id];
  76   name += "Mult";
  77   TH1D * h = (TH1D*) GetHisto(name.Data());
  78   if (!h) {
  79     h = BookHistoMult(name.Data(), Form("Multiplicity distribution (%s)",kHistoPrefix[id]));
  80   }
  81
  82   return h;
  83
  84 }
  85
  86 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoPt (Histo_t id,
  87                                                        Float_t minEta, Float_t maxEta,
  88                                                        Float_t minVz, Float_t maxVz,
  89                                                        Bool_t scaleWidth) {
  90   // Returns a projection of the 3D histo pt/eta/vz.
  91   // WARNING: since that is a histo, the requested range will be discretized to the binning.
  92   // Always avoids under (over) flows
  93   // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
  94
  95   TH3D * h3D = GetHistoPtEtaVz(id);
  96
  97   // Get range in terms of bin numners.  If the float range is
  98   // less than -11111 take the range from the first to the last bin (i.e. no
  99   // under/over-flows)
 100   // FIXME: UNDERFLOWS
 101   Int_t min1 = minEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(minEta);
 102   Int_t min2  = minVz  < -11111 ? -1 : h3D ->GetZaxis()->FindBin(minVz) ;
 103
 104   // Int_t max1 = maxEta  < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
 105   // Int_t max2  = maxVz  < -11111 ? h3D->GetNbinsZ() : h3D ->GetZaxis()->FindBin(maxVz -0.00001);
 106   Int_t max1 = maxEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
 107   Int_t max2  = maxVz  < -11111 ? -1 : h3D ->GetZaxis()->FindBin(maxVz -0.00001);
 108
 109
 110   TString hname = h3D->GetName();
 111   hname = hname +  "_pt_" + long (min1)  + "_" + long(max1) + "_" + long (min2)  + "_" + long(max2);
 112
 113
 114   if (gROOT->FindObjectAny(hname.Data())){
 115     AliError(Form("An object called %s already exists,adding suffix",hname.Data()));
 116     hname += "_2";
 117   }
 118
 119   TH1D * h = h3D->ProjectionX(hname.Data(), min1, max1, min2, max2, "E");
 120   if(scaleWidth) h ->Scale(1.,"width");
 121
 122   return h;
 123
 124 }
 125
 126
 127 TH2D * AliAnalysisMultPbTrackHistoManager::GetHistoPtVz (Histo_t id,
 128                                                          Float_t minEta, Float_t maxEta,
 129                                                          Bool_t scaleWidth) {
 130   // Returns a projection of the 3D histo pt/eta/vz.
 131   // WARNING: since that is a histo, the requested range will be discretized to the binning.
 132   // Always avoids under (over) flows
 133   // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
 134
 135   // FIXME: what do I do here for the scaling?
 136
 137   TH3D * h3D = GetHistoPtEtaVz(id);
 138
 139   // Get range in terms of bin numners.  If the float range is
 140   // less than -11111 take the range from the first to the last bin (i.e. no
 141   // under/over-flows)
 142   Int_t min1 = minEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(minEta);
 143   Int_t max1 = maxEta  < -11111 ? -1 : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
 144
 145
 146   TString hname = h3D->GetName();
 147   hname = hname +  "_ptvz_" + long (min1)  + "_" + long(max1);
 148
 149   if (gROOT->FindObjectAny(hname.Data())){
 150     AliError(Form("An object called %s already exists,adding suffix",hname.Data()));
 151     hname += "_2";
 152   }
 153
 154   h3D->GetYaxis()->SetRange(min1,max1);
 155
 156   TH2D * h =  (TH2D*) h3D->Project3D("zxe");
 157   if(scaleWidth) h ->Scale(1.,"width");
 158
 159   return h;
 160
 161 }
 162
 163
 164 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoVz (Histo_t id,
 165                                                        Float_t minPt, Float_t maxPt,
 166                                                        Float_t minEta, Float_t maxEta,
 167                                                        Bool_t scaleWidth) {
 168   // Returns a projection of the 3D histo pt/eta/vz.
 169   // WARNING: since that is a histo, the requested range will be discretized to the binning.
 170   // Always avoids under (over) flows
 171   // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
 172
 173   TH3D * h3D = GetHistoPtEtaVz(id);
 174
 175   // Get range in terms of bin numners.  If the float range is
 176   // less than -11111 take the range from the first to the last bin (i.e. no
 177   // under/over-flows)
 178   Int_t min1  = minPt  < -11111 ? 1 : h3D ->GetXaxis()->FindBin(minPt) ;
 179   Int_t min2  = minEta < -11111 ? 1 : h3D ->GetYaxis()->FindBin(minEta);
 180
 181   Int_t max1  = maxPt  < -11111 ? h3D->GetNbinsX() : h3D ->GetXaxis()->FindBin(maxPt -0.00001);
 182   Int_t max2  = maxEta < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxEta-0.00001);
 183
 184
 185   TString hname = h3D->GetName();
 186   hname = hname +  "_Vz_" + long (min1)  + "_" + long(max1) + "_" + long (min2)  + "_" + long(max2);
 187
 188   if (gROOT->FindObjectAny(hname.Data())){
 189     AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
 190     hname+="_2";
 191   }
 192
 193   TH1D * h = h3D->ProjectionZ(hname.Data(), min1, max1, min2, max2, "E");
 194   if(scaleWidth) h ->Scale(1.,"width");
 195   return h;
 196
 197
 198 }
 199
 200 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoEta (Histo_t id,
 201                                                         Float_t minPt, Float_t maxPt,
 202                                                         Float_t minVz, Float_t maxVz,
 203                                                         Bool_t scaleWidth) {
 204   // Returns a projection of the 3D histo pt/eta/vz.
 205   // WARNING: since that is a histo, the requested range will be discretized to the binning.
 206   // Always avoids under (over) flows
 207   // If scaleWidth = kTRUE, the projection is scaled for the bin width (default)
 208
 209   TH3D * h3D = GetHistoPtEtaVz(id);
 210
 211   // Get range in terms of bin numners.  If the float range is
 212   // less than -11111 take the range from the first to the last bin (i.e. no
 213   // under/over-flows)
 214   Int_t min1 = minPt < -11111 ? 1 : h3D ->GetXaxis()->FindBin(minPt) ;
 215   Int_t min2 = minVz < -11111 ? 1 : h3D ->GetYaxis()->FindBin(minVz);
 216
 217   Int_t max1 = maxPt < -11111 ? h3D->GetNbinsX() : h3D ->GetXaxis()->FindBin(maxPt -0.00001);
 218   Int_t max2 = maxVz < -11111 ? h3D->GetNbinsY() : h3D ->GetYaxis()->FindBin(maxVz-0.00001);
 219
 220   TString hname = h3D->GetName();
 221   hname = hname +  "_Eta_" + long (min1)  + "_" + long(max1) + "_" + long (min2)  + "_" + long(max2);
 222
 223   if (gROOT->FindObjectAny(hname.Data())){
 224     AliError(Form("An object called %s already exists, adding suffix",hname.Data()));
 225     hname+="_2";
 226   }
 227
 228   TH1D * h = h3D->ProjectionY(hname.Data(), min1, max1, min2, max2, "E");
 229   if(scaleWidth) h ->Scale(1.,"width");
 230   return h;
 231 }
 232
 233 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoProcess(Histo_t id) {
 234
 235   // Returns histogram with particle specties
 236
 237   TString name = TString(kHistoPrefix[id])+"_Process";
 238
 239   TH1D * h =  (TH1D*) GetHisto(name);
 240   if (!h) {
 241     name+=fHNameSuffix;
 242     Bool_t oldStatus = TH1::AddDirectoryStatus();
 243     TH1::AddDirectory(kFALSE);
 244
 245     AliInfo(Form("Booking histo %s",name.Data()));
 246
 247     h = new TH1D (name.Data(), Form("Particle production process (%s)",kHistoPrefix[id]), kPNoProcess+1, -0.5, kPNoProcess+1-0.5);
 248     Int_t nbin = kPNoProcess+1;
 249     for(Int_t ibin = 0; ibin < nbin; ibin++){
 250       h->GetXaxis()->SetBinLabel(ibin+1,TMCProcessName[ibin]);
 251     }
 252     TH1::AddDirectory(oldStatus);
 253     fList->Add(h);
 254
 255
 256   }
 257   return h;
 258
 259
 260 }
 261
 262
 263 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoSpecies(Histo_t id) {
 264
 265   // Returns histogram with particle specties
 266
 267   TString name = TString(kHistoPrefix[id])+"_Species";
 268
 269   TH1D * h =  (TH1D*) GetHisto(name);
 270   if (!h) {
 271     name+=fHNameSuffix;
 272     Bool_t oldStatus = TH1::AddDirectoryStatus();
 273     TH1::AddDirectory(kFALSE);
 274
 275     AliInfo(Form("Booking histo %s",name.Data()));
 276
 277     h = new TH1D (name.Data(), Form("Particle species (%s)",kHistoPrefix[id]), kNPart+1, -0.5, kNPart+1-0.5);
 278     Int_t nbin = kNPart+1;
 279     for(Int_t ibin = 0; ibin < nbin; ibin++){
 280       h->GetXaxis()->SetBinLabel(ibin+1,kSpeciesName[ibin]);
 281     }
 282     TH1::AddDirectory(oldStatus);
 283     fList->Add(h);
 284
 285
 286   }
 287   return h;
 288
 289
 290 }
 291
 292 TH1D * AliAnalysisMultPbTrackHistoManager::GetHistoVzEvent(Histo_t id) {
 293
 294   // Returns histogram with Vz of the event
 295
 296   TString name = TString(kHistoPrefix[id])+"_VzEvent";
 297
 298   TH1D * h =  (TH1D*) GetHisto(name);
 299   if (!h) {
 300     name+=fHNameSuffix;
 301     Bool_t oldStatus = TH1::AddDirectoryStatus();
 302     TH1::AddDirectory(kFALSE);
 303
 304     AliInfo(Form("Booking histo %s",name.Data()));
 305
 306     h = new TH1D (name.Data(), Form("Vz of the event (%s)",kHistoPrefix[id]), 10, -10, 10);
 307     h->Sumw2();
 308     h->SetXTitle("V_{z}^{event} (cm)");
 309     TH1::AddDirectory(oldStatus);
 310     fList->Add(h);
 311
 312
 313   }
 314   return h;
 315
 316
 317 }
 318
 319
 320
 321 TH1I * AliAnalysisMultPbTrackHistoManager::GetHistoStats() {
 322   // Returns histogram with event statistiscs (processed events at each step)
 323
 324   TH1I * h =  (TH1I*) GetHisto("hStats");
 325   if (!h) h = BookHistoStats();
 326   return h;
 327
 328 }
 329
 330
 331
 332 TH1 * AliAnalysisMultPbTrackHistoManager::GetHisto(const char * name) {
 333   //Search list for histo
 334   // TODO: keep track of histo id rather than searching by name?
 335   return (TH1*) fList->FindObject(TString(name)+fHNameSuffix);
 336
 337 }
 338
 339 void AliAnalysisMultPbTrackHistoManager::ScaleHistos(Double_t nev, Option_t * option) {
 340   // Scales all histos in the list for nev
 341   // option can be used to pass further options to TH1::Scale
 342   TH1 * h = 0;
 343   TIter iter = fList->MakeIterator();
 344   while ((h = (TH1*) iter.Next())) {
 345     if (!h->InheritsFrom("TH1")) {
 346       AliFatal (Form("%s does not inherits from TH1, cannot scale",h->GetName()));
 347     }
 348     AliInfo(Form("Scaling %s, nev %2.2f", h->GetName(), nev));
 349
 350     h->Scale(1./nev,option);
 351   }
 352
 353 }
 354
 355 TH3D * AliAnalysisMultPbTrackHistoManager::BookHistoPtEtaVz(const char * name, const char * title) {
 356   // Books a 3D histo of Pt/eta/vtx
 357   // TODO: make the binning settable, variable binning?
 358
 359   Bool_t oldStatus = TH1::AddDirectoryStatus();
 360   TH1::AddDirectory(kFALSE);
 361
 362   TString hname = name;
 363   hname+=fHNameSuffix;
 364
 365   AliInfo(Form("Booking %s",hname.Data()));
 366
 367   // Binning from Jacek task
 368   // const Int_t nptbins = 49;
 369   // 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};
 370   const Int_t nptbins = 68;
 371   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};
 372
 373   const Int_t netabins=20;
 374   Double_t binsEta[netabins+1];
 375   Float_t minEta = -1;
 376   Float_t maxEta =  1;
 377   Float_t etaStep = (maxEta-minEta)/netabins;
 378   for(Int_t ibin = 0; ibin < netabins; ibin++){
 379     binsEta[ibin]   = minEta + ibin*etaStep;
 380     binsEta[ibin+1] = minEta + ibin*etaStep + etaStep;
 381   }
 382
 383   const Int_t nvzbins=10;
 384   Double_t binsVz[nvzbins+1];
 385   Float_t minVz = -10;
 386   Float_t maxVz =  10;
 387   Float_t vzStep = (maxVz-minVz)/nvzbins;
 388   for(Int_t ibin = 0; ibin < nvzbins; ibin++){
 389     binsVz[ibin]   = minVz + ibin*vzStep;
 390     binsVz[ibin+1] = minVz + ibin*vzStep + vzStep;
 391   }
 392
 393
 394   TH3D * h = new TH3D (hname,title,
 395                        nptbins,  binsPt,
 396                        netabins, binsEta,
 397                        nvzbins,  binsVz
 398                        );
 399
 400   h->SetXTitle("p_{T} (GeV)");
 401   h->SetYTitle("#eta");
 402   h->SetZTitle("V_{z}^{tracks} (cm)");
 403   h->Sumw2();
 404
 405   fList->Add(h);
 406
 407   TH1::AddDirectory(oldStatus);
 408   return h;
 409 }
 410
 411 TH1D * AliAnalysisMultPbTrackHistoManager::BookHistoDCA(const char * name, const char * title) {
 412   // Books a DCA histo
 413
 414   Bool_t oldStatus = TH1::AddDirectoryStatus();
 415   TH1::AddDirectory(kFALSE);
 416
 417   TString hname = name;
 418   hname+=fHNameSuffix;
 419
 420   AliInfo(Form("Booking %s",hname.Data()));
 421
 422
 423   TH1D * h = new TH1D (hname,title, 200,0,200);
 424
 425   h->SetXTitle("#Delta DCA");
 426   h->Sumw2();
 427
 428   fList->Add(h);
 429
 430   TH1::AddDirectory(oldStatus);
 431   return h;
 432 }
 433 TH1D * AliAnalysisMultPbTrackHistoManager::BookHistoMult(const char * name, const char * title) {
 434   // Books a multiplicity histo
 435
 436   Bool_t oldStatus = TH1::AddDirectoryStatus();
 437   TH1::AddDirectory(kFALSE);
 438
 439   TString hname = name;
 440   hname+=fHNameSuffix;
 441
 442   AliInfo(Form("Booking %s",hname.Data()));
 443
 444
 445   TH1D * h = new TH1D (hname,title, 600, 0,6000);
 446
 447   h->SetXTitle("N tracks");
 448   h->Sumw2();
 449
 450   fList->Add(h);
 451
 452   TH1::AddDirectory(oldStatus);
 453   return h;
 454 }
 455
 456 TH1I * AliAnalysisMultPbTrackHistoManager::BookHistoStats() {
 457   // Books histogram with event statistiscs (processed events at each step)
 458
 459   Bool_t oldStatus = TH1::AddDirectoryStatus();
 460   TH1::AddDirectory(kFALSE);
 461
 462   AliInfo(Form("Booking Stat histo"));
 463
 464   TH1I * h = new TH1I (TString("hStats")+fHNameSuffix, "Number of processed events", kNStatBins, -0.5, kNStatBins-0.5);
 465   for(Int_t istatbin = 0; istatbin < kNStatBins; istatbin++){
 466     h->GetXaxis()->SetBinLabel(istatbin+1,kStatStepNames[istatbin]);
 467   }
 468   TH1::AddDirectory(oldStatus);
 469   fList->Add(h);
 470   return h;
 471 }
 472
 473 Int_t AliAnalysisMultPbTrackHistoManager::GetLocalParticleID(AliMCParticle * part) {
 474   // returns the local code (Part_t)
 475
 476   Int_t pdgcode = part->PdgCode();
 477   switch(pdgcode) {
 478
 479   case 211:
 480     return kPartPiPlus;
 481     break;
 482   case -211:
 483     return kPartPiMinus;
 484     break;
 485   case 2212:
 486     return kPartP;
 487     break;
 488   case -2212:
 489     return kPartPBar;
 490     break;
 491   case 321:
 492     return kPartKPlus;
 493     break;
 494   case -321:
 495     return kPartKMinus;
 496     break;
 497   case -11:
 498     return kPartLMinus;
 499     break;
 500   case 11:
 501     return kPartLPlus;
 502     break;
 503   case -13:
 504     return kPartLMinus;
 505     break;
 506   case 13:
 507     return kPartLPlus;
 508     break;
 509   default:
 510     return kPartOther;
 511   }
 512 }
 513
 514
 515