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 #include <iostream>
  38
  39 using namespace std;
  40
  41 ClassImp(AliSpectraAODHistoManager)
  42
  43
  44 using namespace AliSpectraNameSpace;
  45
  46 AliSpectraAODHistoManager::AliSpectraAODHistoManager(const char *name): TNamed(name, "AOD Spectra Histo Manager"), fOutputList(0)
  47 {
  48   // ctor
  49    fOutputList = new TList;
  50    fOutputList->SetOwner();
  51    Bool_t oldStatus = TH1::AddDirectoryStatus();
  52    TH1::AddDirectory(kFALSE);
  53    for (Int_t ihist  = 0; ihist < kNHist ; ihist++)
  54    {
  55       if (ihist <= kNPtGenHist) BookPtGenHistogram(kHistName[ihist]);  // PT histos
  56       if (ihist > kNPtGenHist && ihist <= kNPtRecHist) BookPtRecHistogram(kHistName[ihist]);  // PT histos
  57       if (ihist > kNPtRecHist && ihist <= kNHistPID) BookPIDHistogram(kHistName[ihist]);  // PID histos
  58       if (ihist > kNHistPID && ihist <= kNHistNSig) BookNSigHistogram(kHistName[ihist]);  // NSigmaSep histos
  59       if (ihist > kNHistNSig) BookqVecHistogram(kHistName[ihist]);  // qDistr histos
  60    }
  61    //adding quickly o plot to check the centrality distr in two different ways
  62    TH2F *hist=new TH2F("CentCheck","CentCheck",1000,0,100,1000,0,100);
  63    hist->SetXTitle("fAOD->GetCentrality()->GetCentralityPercentile(V0M)");
  64    hist->SetYTitle("fAOD->GetHeader()->GetCentralityP()->GetCentralityPercentileUnchecked(V0M)");
  65    hist->Sumw2();
  66    fOutputList->Add(hist);
  67
  68    TH1::AddDirectory(oldStatus);
  69
  70 }
  71
  72 //_______________________________________________________
  73
  74 TH2F* AliSpectraAODHistoManager::BookPtGenHistogram(const char * name)
  75 {
  76    // Return a pt histogram with predefined binning, set the ID and add it to the output list
  77    AliInfo(Form("Booking pt histogram %s", name));
  78
  79    //standard histo
  80    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};
  81    Int_t nbinsTempl=52;
  82
  83    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)
  84    hist->GetXaxis()->SetTitle("generated P_{T} (GeV / c)");
  85    hist->GetYaxis()->SetTitle("DCA xy");
  86    hist->SetMarkerStyle(kFullCircle);
  87    hist->Sumw2();
  88    fOutputList->Add(hist);
  89
  90    return hist;
  91 }
  92
  93
  94 //_______________________________________________________
  95 TH2F* AliSpectraAODHistoManager::BookPtRecHistogram(const char * name)
  96 {
  97    // Return a pt histogram with predefined binning, set the ID and add it to the output list
  98    AliInfo(Form("Booking pt histogram %s", name));
  99
 100    //standard histo
 101    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};
 102    Int_t nbinsTempl=52;
 103
 104    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)
 105    hist->GetXaxis()->SetTitle("P_{T} (GeV / c)");
 106    hist->GetYaxis()->SetTitle("DCA xy");
 107    hist->SetMarkerStyle(kFullCircle);
 108    hist->Sumw2();
 109    fOutputList->Add(hist);
 110
 111    return hist;
 112 }
 113
 114 //_____________________________________________________________________________
 115
 116 TH2F* AliSpectraAODHistoManager::BookPIDHistogram(const char * name)
 117 {
 118    // Return a pt histogram with predefined binning, set the ID and add it to the output list
 119    AliInfo(Form("Booking PID histogram %s", name));
 120
 121    TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 200, 0, 2.5, 2000, -1000, 1000);
 122    hist->GetXaxis()->SetTitle("(GeV / c)");
 123    hist->GetYaxis()->SetTitle("PID signal");
 124 //  hist->Sumw2();
 125    fOutputList->Add(hist);
 126
 127    return hist;
 128 }
 129
 130 //_____________________________________________________________________________
 131
 132 TH2F* AliSpectraAODHistoManager::BookNSigHistogram(const char * name)
 133 {
 134   // Return a pt histogram with predefined binning, set the ID and add it to the output list
 135   AliInfo(Form("Booking NSigma histogram %s", name));
 136
 137   TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 200, 0, 2.5, 2000,-40, 40);
 138   hist->GetXaxis()->SetTitle("P (GeV / c)");
 139   hist->GetYaxis()->SetTitle("TPC");
 140   //hist->Sumw2();
 141   fOutputList->Add(hist);
 142
 143   return hist;
 144 }
 145
 146 //_____________________________________________________________________________
 147
 148 TH2F* AliSpectraAODHistoManager::BookqVecHistogram(const char * name)
 149 {
 150   // Return a pt histogram with predefined binning, set the ID and add it to the output list
 151   AliInfo(Form("Booking q Vector histogram %s", name));
 152
 153   TH2F * hist = new TH2F(name, Form("q Vector distribution vs Centrality (%s)", name), 100, 0, 10, 100, 0, 100);
 154   hist->GetXaxis()->SetTitle("q vector");
 155   hist->GetYaxis()->SetTitle("Centrality (V0)");
 156   //  hist->Sumw2();
 157   fOutputList->Add(hist);
 158
 159   return hist;
 160 }
 161
 162
 163 //_____________________________________________________________________________
 164
 165 TH1F* AliSpectraAODHistoManager::GetPtHistogram1D(const char * name,Double_t minDCA,Double_t maxDCA)
 166 {
 167   //   //return the projection of the TH2 (pt,DCA) in the DCA bin range [firstDCAbin,lastDCAbin]
 168   //   //if minDCA=-1 && maxDCA=-1 the projection is done using the full DCA range
 169   TH2F *hist=(TH2F*)fOutputList->FindObject(name);
 170   TH1F *outhist=0x0;
 171   Printf("--- Projecting %s on Xaxis[%f,%f]:",name,minDCA,maxDCA);
 172   if(minDCA==-1 && maxDCA==-1){
 173     outhist=(TH1F*)hist->ProjectionX("_px",0,-1,"e");
 174     Printf("Full Range");
 175   }else {
 176     Int_t firstbin=hist->GetYaxis()->FindBin(minDCA);
 177     Int_t lastbin=hist->GetYaxis()->FindBin(maxDCA);
 178     Printf("firstbin: %d lastbin: %d",firstbin,lastbin);
 179     outhist=(TH1F*)hist->ProjectionX("_px",firstbin,lastbin,"e");
 180   }
 181   Printf("Entries outhist: %.0f   Entries hist: %.0f",outhist->GetEntries(),hist->GetEntries());
 182   return outhist;
 183 }
 184
 185 //_____________________________________________________________________________
 186
 187 TH1F* AliSpectraAODHistoManager::GetDCAHistogram1D(const char * name,Double_t minPt,Double_t maxPt)
 188 {
 189   //   //return the projection of the TH2 (pt,DCA) in the DCA bin range [firstDCAbin,lastDCAbin]
 190   //   //if minPt=-1 && maxPt=-1 the projection is done using the full DCA range
 191   TH2F *hist=(TH2F*)fOutputList->FindObject(name);
 192   TH1F *outhist=0x0;
 193   Printf("--- Projecting %s on Yaxis[%f,%f]:",name,minPt,maxPt);
 194   if(minPt==-1 && maxPt==-1){
 195     outhist=(TH1F*)hist->ProjectionY("_py",0,-1,"e");
 196     Printf("Full Range");
 197   }else {
 198     Int_t firstbin=hist->GetXaxis()->FindBin(minPt);
 199     Int_t lastbin=hist->GetXaxis()->FindBin(maxPt);
 200     Printf("firstbin: %d lastbin: %d",firstbin,lastbin);
 201     outhist=(TH1F*)hist->ProjectionY("_py",firstbin,lastbin,"e");
 202     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));
 203   }
 204   Printf("Entries outhist: %.0f   Entries hist: %.0f",outhist->GetEntries(),hist->GetEntries());
 205   return outhist;
 206 }
 207
 208 //_____________________________________________________________________________
 209
 210 Long64_t AliSpectraAODHistoManager::Merge(TCollection* list)
 211 {
 212   // Merge a list of AliSpectraAODHistoManager objects with this.
 213   // Returns the number of merged objects (including this).
 214
 215   //  AliInfo("Merging");
 216
 217   if (!list)
 218     return 0;
 219
 220   if (list->IsEmpty())
 221     return 1;
 222
 223   TIterator* iter = list->MakeIterator();
 224   TObject* obj;
 225
 226   // collections of all histograms
 227   TList collections;
 228
 229   Int_t count = 0;
 230
 231   while ((obj = iter->Next())) {
 232     AliSpectraAODHistoManager* entry = dynamic_cast<AliSpectraAODHistoManager*> (obj);
 233     if (entry == 0)
 234       continue;
 235
 236     TList * hlist = entry->GetOutputList();
 237     collections.Add(hlist);
 238     count++;
 239   }
 240
 241   fOutputList->Merge(&collections);
 242
 243   delete iter;
 244
 245   return count+1;
 246 }
 247
 248
 249 TH1* AliSpectraAODHistoManager::GetHistogram1D(UInt_t histoType, UInt_t particleType, UInt_t charge) {
 250   // GetHistogram using particle ID and histogram type
 251   Int_t baseId = -1;
 252
 253   switch(histoType) {
 254   case kHistPtGenTruePrimary:
 255     baseId = kHistPtGenTruePrimaryPionPlus;
 256     break;
 257   case kHistPtRecSigma:
 258     baseId = kHistPtRecSigmaPionPlus;
 259     break;
 260   case kHistPtRecTruePrimary:
 261     baseId = kHistPtRecTruePrimaryPionPlus;
 262     break;
 263   case kHistPtRecSigmaPrimary:
 264     baseId = kHistPtRecSigmaPrimaryPionPlus;
 265     break;
 266   case kHistPtRecSigmaSecondaryMaterial:
 267     baseId = kHistPtRecSigmaSecondaryMaterialPionPlus;
 268     break;
 269   case kHistPtRecSigmaSecondaryWeakDecay:
 270     baseId = kHistPtRecSigmaSecondaryWeakDecayPionPlus;
 271     break;
 272   default:
 273     baseId = -1;
 274   }
 275
 276   if (baseId < 0)
 277     AliFatal("Wrong histogram type");
 278
 279   baseId = baseId + particleType + 3*(charge);
 280   //  cout << "ID " << baseId << endl;
 281
 282   return GetHistogram(baseId);
 283 }
 284
 285 TH2* AliSpectraAODHistoManager::GetHistogram2D(UInt_t histoType, UInt_t particleType, UInt_t charge){
 286   // returns histo based on ids, casting it to TH2*
 287   return (TH2*) GetHistogram1D(histoType,particleType,charge);
 288
 289
 290 }