PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraAODHistoManager.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 //         AliSpectraAODHistoManager 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 "AliSpectraAODHistoManager.h"
  37
  38 #include <iostream>
  39
  40 using namespace std;
  41
  42 ClassImp(AliSpectraAODHistoManager)
  43
  44
  45 using namespace AliSpectraNameSpace;
  46 #include "HistogramNames.cxx" // generate this automatically running createNames.py
  47
  48 const char * kParticleSpecies[] =
  49   {
  50     "PionPlus",
  51     "KaonPlus",
  52     "ProtonPlus",
  53     "PionMinus",
  54     "KaonMinus",
  55     "ProtonMinus",
  56   };
  57
  58
  59 AliSpectraAODHistoManager::AliSpectraAODHistoManager(const char *name,Int_t nrebin): 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(kHistName[ihist]);  // PT histos
  70       if (ihist > kNPtGenHist && ihist <= kNPtGenAllChHist) BookPtGenAllChHistogram(kHistName[ihist]);  // PT histos
  71       if (ihist > kNPtGenAllChHist && ihist <= kNPtRecHist) BookPtRecHistogram(kHistName[ihist]);  // PT histos
  72       if (ihist > kNPtRecHist && ihist <= kNPtRecAllChHist) BookPtRecAllChHistogram(kHistName[ihist]);  // PT histos
  73       if (ihist > kNPtRecAllChHist && ihist <= kNHistPID) BookPIDHistogram(kHistName[ihist]);  // PID histos
  74       if (ihist > kNHistPID && ihist <= kNHistNSig) BookNSigHistogram(kHistName[ihist]);  // NSigmaSep histos
  75     }
  76
  77   TH1::AddDirectory(oldStatus);
  78 }
  79
  80 //_______________________________________________________
  81
  82 TH2F* AliSpectraAODHistoManager::BookPtGenHistogram(const char * name)
  83 {
  84   // Return a pt histogram with predefined binning, set the ID and add it to the output list
  85   AliInfo(Form("Booking pt gen histogram %s", name));
  86
  87   //standard histo
  88   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};
  89   Int_t nbinsTempl=52;
  90
  91   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)
  92   hist->GetXaxis()->SetTitle("generated P_{T} (GeV / c)");
  93   hist->GetYaxis()->SetTitle("IsPhysicalPrimary()");
  94   hist->SetMarkerStyle(kFullCircle);
  95   hist->Sumw2();
  96   fOutputList->Add(hist);
  97
  98   return hist;
  99 }
 100
 101 //_______________________________________________________
 102
 103 TH2F* AliSpectraAODHistoManager::BookPtGenAllChHistogram(const char * name)
 104 {
 105   // Return a pt histogram with predefined binning, set the ID and add it to the output list
 106   AliInfo(Form("Booking pt gen histogram - no PID %s", name));
 107
 108   //standard histo
 109   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};
 110   Int_t nbinsTempl=62;
 111
 112   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)
 113   hist->GetXaxis()->SetTitle("generated P_{T} (GeV / c)");
 114   hist->GetYaxis()->SetTitle("IsPhysicalPrimary()");
 115   hist->SetMarkerStyle(kFullCircle);
 116   hist->Sumw2();
 117   fOutputList->Add(hist);
 118
 119   return hist;
 120 }
 121
 122
 123 //_______________________________________________________
 124 TH2F* AliSpectraAODHistoManager::BookPtRecHistogram(const char * name)
 125 {
 126   // Return a pt histogram with predefined binning, set the ID and add it to the output list
 127   AliInfo(Form("Booking pt rec histogram %s,  rebin:%d", name, fNRebin));
 128
 129   //standard histo
 130   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};
 131   Int_t nbinsTempl=52;
 132
 133   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)
 134   if(fNRebin!=0)hist->RebinY(fNRebin);
 135   hist->GetXaxis()->SetTitle("P_{T} (GeV / c)");
 136   hist->GetYaxis()->SetTitle("DCA xy");
 137   hist->SetMarkerStyle(kFullCircle);
 138   hist->Sumw2();
 139   fOutputList->Add(hist);
 140
 141   return hist;
 142 }
 143
 144 //_______________________________________________________
 145 TH2F* AliSpectraAODHistoManager::BookPtRecAllChHistogram(const char * name)
 146 {
 147   // Return a pt histogram with predefined binning, set the ID and add it to the output list
 148   AliInfo(Form("Booking pt rec histogram %s,  rebin:%d", name, fNRebin));
 149
 150   //standard histo
 151   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};
 152   Int_t nbinsTempl=62;
 153
 154   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)
 155   if(fNRebin!=0)hist->RebinY(fNRebin);
 156   hist->GetXaxis()->SetTitle("P_{T} (GeV / c)");
 157   hist->GetYaxis()->SetTitle("DCA xy");
 158   hist->SetMarkerStyle(kFullCircle);
 159   hist->Sumw2();
 160   fOutputList->Add(hist);
 161
 162   return hist;
 163 }
 164
 165 //_____________________________________________________________________________
 166
 167 TH2F* AliSpectraAODHistoManager::BookPIDHistogram(const char * name)
 168 {
 169   // Return a pt histogram with predefined binning, set the ID and add it to the output list
 170   AliInfo(Form("Booking PID histogram %s, rebin:%d", name, fNRebin));
 171
 172   TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 1000, 0, 2.5, 2000, -1000, 1000);
 173   if(fNRebin!=0){
 174     hist->RebinY(fNRebin);
 175     hist->RebinX(fNRebin);
 176   }
 177   hist->GetXaxis()->SetTitle("(GeV / c)");
 178   hist->GetYaxis()->SetTitle("PID signal");
 179   //  hist->Sumw2();
 180   fOutputList->Add(hist);
 181
 182   return hist;
 183 }
 184
 185 //_____________________________________________________________________________
 186
 187 TH2F* AliSpectraAODHistoManager::BookNSigHistogram(const char * name)
 188 {
 189   // Return a pt histogram with predefined binning, set the ID and add it to the output list
 190   AliInfo(Form("Booking NSigma histogram %s, rebin:%d", name, fNRebin));
 191
 192   TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 1000, 0, 2.5, 2000,-40, 40);
 193   if(fNRebin!=0){
 194     hist->RebinY(fNRebin);
 195     hist->RebinX(fNRebin);
 196   }
 197   hist->GetXaxis()->SetTitle("P (GeV / c)");
 198   hist->GetYaxis()->SetTitle("TPC");
 199   //hist->Sumw2();
 200   fOutputList->Add(hist);
 201
 202   return hist;
 203 }
 204
 205 //_____________________________________________________________________________
 206
 207 TH1F* AliSpectraAODHistoManager::GetPtHistogram1D(const char * name,Double_t minDCA,Double_t maxDCA)
 208 {
 209   //   //return the projection of the TH2 (pt,DCA) in the DCA bin range [firstDCAbin,lastDCAbin]
 210   //   //if minDCA=-1 && maxDCA=-1 the projection is done using the full DCA range
 211   TH2F *hist=(TH2F*)fOutputList->FindObject(name);
 212   TH1F *outhist=0x0;
 213   Printf("--- Projecting %s on Xaxis[%f,%f]:",name,minDCA,maxDCA);
 214   if(minDCA==-1 && maxDCA==-1){
 215     outhist=(TH1F*)hist->ProjectionX("_px",0,-1,"e");
 216     Printf("Full Range");
 217   }else {
 218     Int_t firstbin=hist->GetYaxis()->FindBin(minDCA);
 219     Int_t lastbin=hist->GetYaxis()->FindBin(maxDCA);
 220     Printf("firstbin: %d lastbin: %d",firstbin,lastbin);
 221     outhist=(TH1F*)hist->ProjectionX("_px",firstbin,lastbin,"e");
 222   }
 223   Printf("Entries outhist: %.0f   Entries hist: %.0f",outhist->GetEntries(),hist->GetEntries());
 224   return outhist;
 225 }
 226
 227 //_____________________________________________________________________________
 228
 229 TH1F* AliSpectraAODHistoManager::GetDCAHistogram1D(const char * name,Double_t minPt,Double_t maxPt)
 230 {
 231   //   //return the projection of the TH2 (pt,DCA) in the DCA bin range [firstDCAbin,lastDCAbin]
 232   //   //if minPt=-1 && maxPt=-1 the projection is done using the full DCA range
 233   TH2F *hist=(TH2F*)fOutputList->FindObject(name);
 234   TH1F *outhist=0x0;
 235   Printf("--- Projecting %s on Yaxis[%f,%f]:",name,minPt,maxPt);
 236   if(minPt==-1 && maxPt==-1){
 237     outhist=(TH1F*)hist->ProjectionY("_py",0,-1,"e");
 238     Printf("Full Range");
 239   }else {
 240     Int_t firstbin=hist->GetXaxis()->FindBin(minPt);
 241     Int_t lastbin=hist->GetXaxis()->FindBin(maxPt);
 242     Printf("firstbin: %d lastbin: %d",firstbin,lastbin);
 243     outhist=(TH1F*)hist->ProjectionY("_py",firstbin,lastbin,"e");
 244     Printf("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));
 245   }
 246   Printf("Entries outhist: %.0f   Entries hist: %.0f",outhist->GetEntries(),hist->GetEntries());
 247   return outhist;
 248 }
 249
 250 //_____________________________________________________________________________
 251
 252 Long64_t AliSpectraAODHistoManager::Merge(TCollection* list)
 253 {
 254   // Merge a list of AliSpectraAODHistoManager objects with this.
 255   // Returns the number of merged objects (including this).
 256
 257   //  AliInfo("Merging");
 258
 259   if (!list)
 260     return 0;
 261
 262   if (list->IsEmpty())
 263     return 1;
 264
 265   TIterator* iter = list->MakeIterator();
 266   TObject* obj;
 267
 268   // collections of all histograms
 269   TList collections;
 270
 271   Int_t count = 0;
 272
 273   while ((obj = iter->Next())) {
 274     AliSpectraAODHistoManager* entry = dynamic_cast<AliSpectraAODHistoManager*> (obj);
 275     if (entry == 0)
 276       continue;
 277
 278     TList * hlist = entry->GetOutputList();
 279     collections.Add(hlist);
 280     count++;
 281   }
 282
 283   fOutputList->Merge(&collections);
 284
 285   delete iter;
 286
 287   return count+1;
 288 }
 289
 290
 291 TH1* AliSpectraAODHistoManager::GetHistogram1D(UInt_t histoType, UInt_t particleType, UInt_t charge) {
 292   // GetHistogram using particle ID and histogram type
 293   Int_t baseId = -1;
 294
 295   if (particleType == kSpUndefined) {
 296     AliError ("Trying to get histo for undefined particle");
 297     return 0;
 298   }
 299
 300   switch(histoType) {
 301   case kHistPtGenTruePrimary:
 302     baseId = kHistPtGenTruePrimaryPionPlus;
 303     break;
 304   case kHistPtRecSigma:
 305     baseId = kHistPtRecSigmaPionPlus;
 306     break;
 307   case kHistPtRecTrue:
 308     baseId = kHistPtRecTruePionPlus;
 309     break;
 310   case kHistPtRecTruePrimary:
 311     baseId = kHistPtRecTruePrimaryPionPlus;
 312     break;
 313   case kHistPtRecPrimary:
 314     baseId = kHistPtRecPrimaryPionPlus;
 315     break;
 316   case kHistPtRecSigmaPrimary:
 317     baseId = kHistPtRecSigmaPrimaryPionPlus;
 318     break;
 319   case kHistPtRecSigmaSecondaryMaterial:
 320     baseId = kHistPtRecSigmaSecondaryMaterialPionPlus;
 321     break;
 322   case kHistPtRecSigmaSecondaryWeakDecay:
 323     baseId = kHistPtRecSigmaSecondaryWeakDecayPionPlus;
 324     break;
 325   case kHistNSigTPC:
 326     baseId = kHistNSigPionTPC;
 327     break;
 328   case kHistNSigTOF:
 329     baseId = kHistNSigPionTOF;
 330     break;
 331   case kHistNSigTPCTOF:
 332     baseId = kHistNSigPionTPCTOF;
 333     break;
 334   default:
 335     baseId = -1;
 336   }
 337
 338   if (baseId < 0)
 339     AliFatal(Form("Wrong histogram type %d", histoType));
 340
 341   //cout << "T[" << histoType << "] ID["<< baseId <<"] P["<<particleType<<"] C[" << charge
 342   //     << " --> ["<< baseId + particleType + 3*(charge) <<"] = " ;
 343
 344   baseId = baseId + particleType + 3*(charge);
 345
 346   //cout <<  GetHistogram(baseId)->GetName() << endl;
 347
 348   return GetHistogram(baseId);
 349 }
 350
 351 TH2* AliSpectraAODHistoManager::GetHistogram2D(UInt_t histoType, UInt_t particleType, UInt_t charge){
 352   // returns histo based on ids, casting it to TH2*
 353   return (TH2*) GetHistogram1D(histoType,particleType,charge);
 354
 355
 356 }