PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraBothHistoManager.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 //-----------------------------------------------------------------
  18 //         AliSpectraBothHistoManager class
  19 //-----------------------------------------------------------------
  20
  21 #include "TChain.h"
  22 #include "TTree.h"
  23 #include "TLegend.h"
  24 #include "TH1.h"
  25 #include "TH1F.h"
  26 #include "TH2F.h"
  27 #include "TCanvas.h"
  28 #include "AliAnalysisTask.h"
  29 #include "AliAnalysisManager.h"
  30 #include "AliAODTrack.h"
  31 #include "AliAODMCParticle.h"
  32 #include "AliAODEvent.h"
  33 #include "AliAODInputHandler.h"
  34 #include "AliAnalysisTaskESDfilter.h"
  35 #include "AliAnalysisDataContainer.h"
  36 #include "AliSpectraBothHistoManager.h"
  37 #include <iostream>
  38
  39 using namespace std;
  40
  41 ClassImp(AliSpectraBothHistoManager)
  42
  43
  44 using namespace AliSpectraNameSpaceBoth;
  45 #include "HistogramNamesBoth.cxx" // generate this automatically running createNames.py
  46
  47
  48 //const char* kParticleSpecies[] =
  49  // {
  50    // "PionPlus",
  51   //  "KaonPlus",
  52   //  "ProtonPlus",
  53   //  "PionMinus",
  54   //  "KaonMinus",
  55   //  "ProtonMinus",
  56  // };
  57
  58
  59 AliSpectraBothHistoManager::AliSpectraBothHistoManager(const char *name,Int_t nrebin,Bool_t pidqa): TNamed(name, "AOD Spectra Histo Manager"), fOutputList(0), fNRebin(0)
  60 {
  61   // ctor
  62   fNRebin=nrebin;
  63   fOutputList = new TList;
  64   fOutputList->SetOwner();
  65   Bool_t oldStatus = TH1::AddDirectoryStatus();
  66   TH1::AddDirectory(kFALSE);
  67   for (Int_t ihist  = 0; ihist < kNHist ; ihist++)
  68     {
  69       if (ihist <= kNPtGenHist) BookPtGenHistogram(kHistNameBoth[ihist]);  // PT histos
  70       if (ihist > kNPtGenHist && ihist <= kNPtGenAllChHist) BookPtGenAllChHistogram(kHistNameBoth[ihist]);  // PT histos
  71       if (ihist > kNPtGenAllChHist && ihist <= kNPtRecHist) BookPtRecHistogram(kHistNameBoth[ihist]);  // PT histos
  72       if (ihist > kNPtRecHist && ihist <= kNPtRecAllChHist) BookPtRecAllChHistogram(kHistNameBoth[ihist]);  // PT histos
  73       if (ihist > kNPtRecAllChHist && ihist <= kNHistPID && pidqa) BookPIDHistogram(kHistNameBoth[ihist]);  // PID histos
  74       if (ihist > kNHistPID && ihist <= kNHistNSig && pidqa) BookNSigHistogram(kHistNameBoth[ihist]);  // NSigmaSep histos
  75       if(ihist==kHistGenMulvsRawMul) BookGenMulvsRawMulHistogram(kHistNameBoth[ihist]);
  76     }
  77
  78   TH1::AddDirectory(oldStatus);
  79 }
  80 //_______________________________________________________________________________________
  81 AliSpectraBothHistoManager::~AliSpectraBothHistoManager()
  82 {
  83         if (fOutputList)
  84         {
  85                 fOutputList->Delete();
  86                 delete fOutputList;
  87         }
  88 }
  89 //_______________________________________________________
  90
  91 TH2F* AliSpectraBothHistoManager::BookPtGenHistogram(const char * name)
  92 {
  93   // Return a pt histogram with predefined binning, set the ID and add it to the output list
  94   AliInfo(Form("Booking pt gen histogram %s", name));
  95
  96   //standard histo
  97   const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,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.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0};
  98   Int_t nbinsTempl=52;
  99
 100   TH2F * hist = new TH2F(name,Form("P_{T} distribution (%s)", name),nbinsTempl,templBins,2,-0.5,1.5);//need to be at least 1 becuase the generated are filled with (pt,IsPhysPrim)
 101   hist->GetXaxis()->SetTitle("generated P_{T} (GeV / c)");
 102   hist->GetYaxis()->SetTitle("IsPhysicalPrimary()");
 103   hist->SetMarkerStyle(kFullCircle);
 104   hist->Sumw2();
 105   fOutputList->Add(hist);
 106
 107   return hist;
 108 }
 109
 110 //_______________________________________________________
 111
 112 TH2F* AliSpectraBothHistoManager::BookPtGenAllChHistogram(const char * name)
 113 {
 114   // Return a pt histogram with predefined binning, set the ID and add it to the output list
 115   AliInfo(Form("Booking pt gen histogram - no PID %s", name));
 116
 117   //standard histo
 118   const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,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.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,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,7.5,8,8.5,9,9.5,10};
 119   Int_t nbinsTempl=62;
 120
 121   TH2F * hist = new TH2F(name,Form("P_{T} distribution (All Ch) (%s)", name),nbinsTempl,templBins,2,-0.5,1.5);//need to be at least 1 becuase the generated are filled with (pt,IsPhysPrim)
 122   hist->GetXaxis()->SetTitle("generated P_{T} (GeV / c)");
 123   hist->GetYaxis()->SetTitle("IsPhysicalPrimary()");
 124   hist->SetMarkerStyle(kFullCircle);
 125   hist->Sumw2();
 126   fOutputList->Add(hist);
 127
 128   return hist;
 129 }
 130
 131
 132 //_______________________________________________________
 133 TH2F* AliSpectraBothHistoManager::BookPtRecHistogram(const char * name)
 134 {
 135   // Return a pt histogram with predefined binning, set the ID and add it to the output list
 136   AliInfo(Form("Booking pt rec histogram %s,  rebin:%d", name, fNRebin));
 137
 138   //standard histo
 139   const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,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.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0};
 140   Int_t nbinsTempl=52;
 141
 142   TH2F * hist = new TH2F(name,Form("reconstructed P_{T} distribution (%s)", name),nbinsTempl,templBins,3000,-3,3);//need to be at least 1 becuase the generated are filled with (pt,IsPhysPrim)
 143   if(fNRebin!=0)hist->RebinY(fNRebin);
 144   hist->GetXaxis()->SetTitle("P_{T} (GeV / c)");
 145   hist->GetYaxis()->SetTitle("DCA xy");
 146   hist->SetMarkerStyle(kFullCircle);
 147   hist->Sumw2();
 148   fOutputList->Add(hist);
 149
 150   return hist;
 151 }
 152
 153 //_______________________________________________________
 154 TH2F* AliSpectraBothHistoManager::BookPtRecAllChHistogram(const char * name)
 155 {
 156   // Return a pt histogram with predefined binning, set the ID and add it to the output list
 157   AliInfo(Form("Booking pt rec histogram %s,  rebin:%d", name, fNRebin));
 158
 159   //standard histo
 160   const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,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.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,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,7.5,8,8.5,9,9.5,10};
 161   Int_t nbinsTempl=62;
 162
 163   TH2F * hist = new TH2F(name,Form("reconstructed P_{T} distribution (All Ch)  (%s)", name),nbinsTempl,templBins,3000,-3,3);//need to be at least 1 becuase the generated are filled with (pt,IsPhysPrim)
 164   if(fNRebin!=0)hist->RebinY(fNRebin);
 165   hist->GetXaxis()->SetTitle("P_{T} (GeV / c)");
 166   hist->GetYaxis()->SetTitle("DCA xy");
 167   hist->SetMarkerStyle(kFullCircle);
 168   hist->Sumw2();
 169   fOutputList->Add(hist);
 170
 171   return hist;
 172 }
 173
 174 //_____________________________________________________________________________
 175
 176 TH2F* AliSpectraBothHistoManager::BookPIDHistogram(const char * name)
 177 {
 178   // Return a pt histogram with predefined binning, set the ID and add it to the output list
 179   AliInfo(Form("Booking PID histogram %s, rebin:%d", name, fNRebin));
 180 TString tmp(name);
 181   TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 50, 0, 2.5, 200, -1000, 1000);
 182   if(fNRebin!=0){
 183     hist->RebinY(fNRebin);
 184   //  hist->RebinX(fNRebin);
 185   }
 186  if(tmp.Contains("Pt"))
 187         hist->GetXaxis()->SetTitle("P_{t} (GeV / c)");
 188   else
 189         hist->GetXaxis()->SetTitle("P (GeV / c)");
 190
 191   hist->GetYaxis()->SetTitle("PID signal");
 192   //  hist->Sumw2();
 193   fOutputList->Add(hist);
 194
 195   return hist;
 196 }
 197
 198 //_____________________________________________________________________________
 199
 200 TH2F* AliSpectraBothHistoManager::BookNSigHistogram(const char * name)
 201 {
 202   // Return a pt histogram with predefined binning, set the ID and add it to the output list
 203   AliInfo(Form("Booking NSigma histogram %s, rebin:%d", name, fNRebin));
 204   Int_t nbins=200;
 205   Float_t miny=-20;
 206   TString tmp(name);
 207   if(tmp.Contains("TPCTOF"))
 208   {
 209         nbins=100;
 210         miny=0.0;
 211   }
 212   TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 50, 0, 2.5,nbins,miny, 20);
 213   if(fNRebin!=0){
 214     hist->RebinY(fNRebin);
 215     //hist->RebinX(fNRebin);
 216   }
 217   if(tmp.Contains("Pt"))
 218         hist->GetXaxis()->SetTitle("P_{t} (GeV / c)");
 219   else
 220         hist->GetXaxis()->SetTitle("P (GeV / c)");
 221   //hist->GetYaxis()->SetTitle("TPC");
 222   //hist->Sumw2();
 223   fOutputList->Add(hist);
 224
 225   return hist;
 226 }
 227
 228 //_____________________________________________________________________________
 229
 230  TH2F*   AliSpectraBothHistoManager::BookGenMulvsRawMulHistogram(const char * name)
 231 {
 232   // Return a pt histogram with predefined binning, set the ID and add it to the output list
 233   AliInfo(Form("Booking  Gen vs Raw multilicity histogram %s", name));
 234
 235   TH2F * hist = new TH2F(name, Form("Gen vs Raw multilicity (%s)", name), 100, -0.5, 99.5, 100, -0.5, 99.5);
 236   hist->GetXaxis()->SetTitle("gen mul ");
 237   hist->GetYaxis()->SetTitle("raw mul ");
 238   //  hist->Sumw2();
 239   fOutputList->Add(hist);
 240
 241   return hist;
 242 }
 243 //_____________________________________________________________________________
 244
 245 TH1F* AliSpectraBothHistoManager::GetPtHistogram1D(const char * name,Double_t minDCA,Double_t maxDCA)
 246 {
 247   //   //return the projection of the TH2 (pt,DCA) in the DCA bin range [firstDCAbin,lastDCAbin]
 248   //   //if minDCA=-1 && maxDCA=-1 the projection is done using the full DCA range
 249   TH2F *hist=(TH2F*)fOutputList->FindObject(name);
 250   TH1F *outhist=0x0;
 251   AliDebug(2,Form("--- Projecting %s on Xaxis[%f,%f]:",name,minDCA,maxDCA));
 252   if(minDCA==-1 && maxDCA==-1){
 253     outhist=(TH1F*)hist->ProjectionX("_px",0,-1,"e");
 254     AliDebug(2,"Full Range");
 255   }else {
 256     Int_t firstbin=hist->GetYaxis()->FindBin(minDCA);
 257     Int_t lastbin=hist->GetYaxis()->FindBin(maxDCA);
 258     AliDebug(2,Form("firstbin: %d lastbin: %d",firstbin,lastbin));
 259     outhist=(TH1F*)hist->ProjectionX("_px",firstbin,lastbin,"e");
 260   }
 261   AliDebug(2,Form("Entries outhist: %.0f   Entries hist: %.0f",outhist->GetEntries(),hist->GetEntries()));
 262   return outhist;
 263 }
 264
 265 //_____________________________________________________________________________
 266
 267 TH1F* AliSpectraBothHistoManager::GetDCAHistogram1D(const char * name,Double_t minPt,Double_t maxPt)
 268 {
 269   //   //return the projection of the TH2 (pt,DCA) in the DCA bin range [firstDCAbin,lastDCAbin]
 270   //   //if minPt=-1 && maxPt=-1 the projection is done using the full DCA range
 271   TH2F *hist=(TH2F*)fOutputList->FindObject(name);
 272   TH1F *outhist=0x0;
 273   AliDebug(2,Form("--- Projecting %s on Yaxis[%f,%f]:",name,minPt,maxPt));
 274   if(minPt==-1 && maxPt==-1){
 275     outhist=(TH1F*)hist->ProjectionY("_py",0,-1,"e");
 276     AliDebug(2,"Full Range");
 277   }else {
 278     Int_t firstbin=hist->GetXaxis()->FindBin(minPt);
 279     Int_t lastbin=hist->GetXaxis()->FindBin(maxPt);
 280     AliDebug(2,Form("firstbin: %d lastbin: %d",firstbin,lastbin));
 281     outhist=(TH1F*)hist->ProjectionY("_py",firstbin,lastbin,"e");
 282     AliDebug(2,Form("GetDCAHistogram1D(%s) BinRange:%d  %d  Pt Range: %f %f",hist->GetName(),firstbin,lastbin,hist->GetXaxis()->GetBinLowEdge(firstbin),hist->GetXaxis()->GetBinLowEdge(firstbin)+hist->GetXaxis()->GetBinWidth(lastbin)));
 283   }
 284   AliDebug(2,Form("Entries outhist: %.0f   Entries hist: %.0f",outhist->GetEntries(),hist->GetEntries()));
 285   return outhist;
 286 }
 287
 288 //_____________________________________________________________________________
 289
 290 Long64_t AliSpectraBothHistoManager::Merge(TCollection* list)
 291 {
 292   // Merge a list of AliSpectraBothHistoManager objects with this.
 293   // Returns the number of merged objects (including this).
 294
 295   //  AliInfo("Merging");
 296
 297   if (!list)
 298     return 0;
 299
 300   if (list->IsEmpty())
 301     return 1;
 302
 303   TIterator* iter = list->MakeIterator();
 304   TObject* obj;
 305
 306   // collections of all histograms
 307   TList collections;
 308
 309   Int_t count = 0;
 310
 311   while ((obj = iter->Next())) {
 312     AliSpectraBothHistoManager* entry = dynamic_cast<AliSpectraBothHistoManager*> (obj);
 313     if (entry == 0)
 314       continue;
 315
 316     TList * hlist = entry->GetOutputList();
 317     collections.Add(hlist);
 318     count++;
 319   }
 320
 321   fOutputList->Merge(&collections);
 322
 323   delete iter;
 324
 325   return count+1;
 326 }
 327
 328 //________________________________________________________________________________
 329 TH1* AliSpectraBothHistoManager::GetHistogram1D(UInt_t histoType, UInt_t particleType, UInt_t charge) {
 330   // GetHistogram using particle ID and histogram type
 331   Int_t baseId = -1;
 332
 333   if (particleType == kSpUndefined) {
 334     AliError ("Trying to get histo for undefined particle");
 335     return 0;
 336   }
 337
 338   switch(histoType) {
 339   case kHistPtGenTruePrimary:
 340     baseId = kHistPtGenTruePrimaryPionPlus;
 341     break;
 342   case kHistPtRecSigma:
 343     baseId = kHistPtRecSigmaPionPlus;
 344     break;
 345   case kHistPtRecTrue:
 346     baseId = kHistPtRecTruePionPlus;
 347     break;
 348   case kHistPtRecTruePrimary:
 349     baseId = kHistPtRecTruePrimaryPionPlus;
 350     break;
 351   case kHistPtRecPrimary:
 352     baseId = kHistPtRecPrimaryPionPlus;
 353     break;
 354   case kHistPtRecSigmaPrimary:
 355     baseId = kHistPtRecSigmaPrimaryPionPlus;
 356     break;
 357   case kHistPtRecSigmaSecondaryMaterial:
 358     baseId = kHistPtRecSigmaSecondaryMaterialPionPlus;
 359     break;
 360   case kHistPtRecSigmaSecondaryWeakDecay:
 361     baseId = kHistPtRecSigmaSecondaryWeakDecayPionPlus;
 362     break;
 363   case kHistNSigTPC:
 364     baseId = kHistNSigPionTPC;
 365     break;
 366   case kHistNSigTOF:
 367     baseId = kHistNSigPionTOF;
 368     break;
 369   case kHistNSigTPCTOF:
 370     baseId = kHistNSigPionTPCTOF;
 371     break;
 372   default:
 373     baseId = -1;
 374   }
 375
 376   if (baseId < 0)
 377     AliFatal(Form("Wrong histogram type %d", histoType));
 378
 379   //cout << "T[" << histoType << "] ID["<< baseId <<"] P["<<particleType<<"] C[" << charge
 380   //     << " --> ["<< baseId + particleType + 3*(charge) <<"] = " ;
 381
 382   baseId = baseId + particleType + 3*(charge);
 383
 384   //cout <<  GetHistogram(baseId)->GetName() << endl;
 385
 386   return GetHistogram(baseId);
 387 }
 388 //____________________________________________________________________________________________________
 389 TH2* AliSpectraBothHistoManager::GetHistogram2D(UInt_t histoType, UInt_t particleType, UInt_t charge){
 390   // returns histo based on ids, casting it to TH2*
 391   return (TH2*) GetHistogram1D(histoType,particleType,charge);
 392
 393
 394 }