PWG0/dNdEta/AlidNdEtaCorrection.cxx

   1 /* $Id$ */
   2
   3 #include "AlidNdEtaCorrection.h"
   4
   5 #include <TCanvas.h>
   6 #include <TH3F.h>
   7 #include <TH1D.h>
   8
   9 //____________________________________________________________________
  10 ClassImp(AlidNdEtaCorrection)
  11
  12 //____________________________________________________________________
  13 AlidNdEtaCorrection::AlidNdEtaCorrection()
  14   : TNamed(),
  15   fTrack2ParticleCorrection(0),
  16   fVertexRecoCorrection(0),
  17   fTriggerBiasCorrectionMBToINEL(0),
  18   fTriggerBiasCorrectionMBToNSD(0),
  19   fTriggerBiasCorrectionMBToND(0)
  20 {
  21   // default constructor
  22 }
  23
  24 //____________________________________________________________________
  25 AlidNdEtaCorrection::AlidNdEtaCorrection(const Char_t* name, const Char_t* title)
  26   : TNamed(name, title),
  27   fTrack2ParticleCorrection(0),
  28   fVertexRecoCorrection(0),
  29   fTriggerBiasCorrectionMBToINEL(0),
  30   fTriggerBiasCorrectionMBToNSD(0),
  31   fTriggerBiasCorrectionMBToND(0)
  32 {
  33   // constructor
  34   //
  35
  36   Float_t binLimitsPt[] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 5.0, 10.0, 100.0};
  37
  38   TString matrixName;
  39   matrixName.Form("%s_nTrackToNPart", name);
  40
  41   fTrack2ParticleCorrection = new AliCorrectionMatrix3D(matrixName, matrixName, 40, -20, 20, 20, -2, 2, 15, binLimitsPt);
  42
  43   Float_t binLimitsN[]   = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5,
  44                             10.5, 12.5, 14.5, 16.5, 18.5, 20.5, 25.5, 30.5, 40.5, 50.5, 100.5, 300.5};
  45   Float_t binLimitsVtx[] = {-20,-15,-10,-6,-3,0,3,6,10,15,20};
  46
  47   matrixName.Form("%s_vtxReco", name);
  48   fVertexRecoCorrection        = new AliCorrectionMatrix2D(matrixName, matrixName, 10,binLimitsVtx ,22,binLimitsN);
  49
  50   matrixName.Form("%s_triggerBias_MBToINEL", name);
  51   fTriggerBiasCorrectionMBToINEL       = new AliCorrectionMatrix2D(matrixName, matrixName, 10,binLimitsVtx ,22,binLimitsN);
  52   matrixName.Form("%s_triggerBias_MBToNSD", name);
  53   fTriggerBiasCorrectionMBToNSD        = new AliCorrectionMatrix2D(matrixName, matrixName, 10,binLimitsVtx ,22,binLimitsN);
  54   matrixName.Form("%s_triggerBias_MBToND", name);
  55   fTriggerBiasCorrectionMBToND         = new AliCorrectionMatrix2D(matrixName, matrixName, 10,binLimitsVtx ,22,binLimitsN);
  56
  57   fTrack2ParticleCorrection      ->SetAxisTitles("vtx z [cm]", "#eta", "p_{T} [GeV/c]");
  58   fVertexRecoCorrection          ->SetAxisTitles("vtx z [cm]", "Ntracks");
  59   fTriggerBiasCorrectionMBToINEL ->SetAxisTitles("vtx z [cm]", "Ntracks");
  60   fTriggerBiasCorrectionMBToNSD  ->SetAxisTitles("vtx z [cm]", "Ntracks");
  61   fTriggerBiasCorrectionMBToND   ->SetAxisTitles("vtx z [cm]", "Ntracks");
  62 }
  63
  64 //____________________________________________________________________
  65 AlidNdEtaCorrection::~AlidNdEtaCorrection()
  66 {
  67   // destructor
  68
  69   if (fTrack2ParticleCorrection) {
  70     delete fTrack2ParticleCorrection;
  71     fTrack2ParticleCorrection = 0;
  72   }
  73
  74   if (fVertexRecoCorrection) {
  75     delete fVertexRecoCorrection;
  76     fVertexRecoCorrection = 0;
  77   }
  78
  79   if (fTriggerBiasCorrectionMBToINEL) {
  80     delete fTriggerBiasCorrectionMBToINEL;
  81     fTriggerBiasCorrectionMBToINEL = 0;
  82   }
  83
  84   if (fTriggerBiasCorrectionMBToNSD) {
  85     delete fTriggerBiasCorrectionMBToNSD;
  86     fTriggerBiasCorrectionMBToNSD = 0;
  87   }
  88
  89   if (fTriggerBiasCorrectionMBToND) {
  90     delete fTriggerBiasCorrectionMBToND;
  91     fTriggerBiasCorrectionMBToND = 0;
  92   }
  93 }
  94
  95 //____________________________________________________________________
  96 void
  97 AlidNdEtaCorrection::Finish() {
  98   //
  99   // finish method
 100   //
 101   // divide the histograms in the AliCorrectionMatrix2D objects to get the corrections
 102
 103   fTrack2ParticleCorrection->Divide();
 104
 105   TH3F* hist = fTrack2ParticleCorrection->GetCorrectionHistogram();
 106   Int_t emptyBins = 0;
 107   for (Int_t x=hist->GetXaxis()->FindBin(-10); x<=hist->GetXaxis()->FindBin(10); ++x)
 108     for (Int_t y=hist->GetYaxis()->FindBin(-0.8); y<=hist->GetYaxis()->FindBin(0.8); ++y)
 109       for (Int_t z=hist->GetZaxis()->FindBin(0.3); z<=hist->GetZaxis()->FindBin(9.9); ++z)
 110         if (hist->GetBinContent(x, y, z) == 0)
 111         {
 112           printf("Empty bin in fTrack2ParticleCorrection at vtx = %f, eta = %f, pt = %f\n", hist->GetXaxis()->GetBinCenter(x), hist->GetYaxis()->GetBinCenter(y), hist->GetZaxis()->GetBinCenter(z));
 113           ++emptyBins;
 114         }
 115
 116   printf("INFO: In the central region fTrack2ParticleCorrection has %d empty bins\n", emptyBins);
 117
 118   fVertexRecoCorrection->Divide();
 119   fTriggerBiasCorrectionMBToINEL->Divide();
 120   fTriggerBiasCorrectionMBToNSD->Divide();
 121   fTriggerBiasCorrectionMBToND->Divide();
 122 }
 123
 124 //____________________________________________________________________
 125 Long64_t
 126 AlidNdEtaCorrection::Merge(TCollection* list) {
 127   // Merge a list of dNdEtaCorrection objects with this (needed for
 128   // PROOF).
 129   // Returns the number of merged objects (including this).
 130
 131   if (!list)
 132     return 0;
 133
 134   if (list->IsEmpty())
 135     return 1;
 136
 137   TIterator* iter = list->MakeIterator();
 138   TObject* obj;
 139
 140   // collections of measured and generated histograms
 141   TList* collectionNtrackToNparticle    = new TList;
 142   TList* collectionVertexReco           = new TList;
 143   TList* collectionTriggerBiasMBToINEL  = new TList;
 144   TList* collectionTriggerBiasMBToNSD   = new TList;
 145   TList* collectionTriggerBiasMBToND    = new TList;
 146
 147   Int_t count = 0;
 148   while ((obj = iter->Next())) {
 149
 150     AlidNdEtaCorrection* entry = dynamic_cast<AlidNdEtaCorrection*> (obj);
 151     if (entry == 0)
 152       continue;
 153
 154     collectionNtrackToNparticle  ->Add(entry->GetTrack2ParticleCorrection());
 155     collectionVertexReco         ->Add(entry->GetVertexRecoCorrection());
 156     collectionTriggerBiasMBToINEL->Add(entry->GetTriggerBiasCorrection("INEL"));
 157     collectionTriggerBiasMBToNSD ->Add(entry->GetTriggerBiasCorrection("NSD"));
 158     collectionTriggerBiasMBToND  ->Add(entry->GetTriggerBiasCorrection("ND"));
 159
 160     count++;
 161
 162     //fNEvents += entry->fNEvents;
 163     //fNTriggeredEvents += entry->fNTriggeredEvents;
 164   }
 165   fTrack2ParticleCorrection      ->Merge(collectionNtrackToNparticle);
 166   fVertexRecoCorrection          ->Merge(collectionVertexReco);
 167   fTriggerBiasCorrectionMBToINEL ->Merge(collectionTriggerBiasMBToINEL);
 168   fTriggerBiasCorrectionMBToNSD  ->Merge(collectionTriggerBiasMBToNSD);
 169   fTriggerBiasCorrectionMBToND   ->Merge(collectionTriggerBiasMBToND);
 170
 171   delete collectionNtrackToNparticle;
 172   delete collectionVertexReco;
 173   delete collectionTriggerBiasMBToINEL;
 174   delete collectionTriggerBiasMBToNSD;
 175   delete collectionTriggerBiasMBToND;
 176
 177   return count+1;
 178 }
 179
 180
 181
 182 //____________________________________________________________________
 183 Bool_t
 184 AlidNdEtaCorrection::LoadHistograms(const Char_t* fileName, const Char_t* dir) {
 185   //
 186   // loads the histograms
 187   //
 188
 189   fTrack2ParticleCorrection      ->LoadHistograms(fileName, dir);
 190   fVertexRecoCorrection          ->LoadHistograms(fileName, dir);
 191   fTriggerBiasCorrectionMBToINEL ->LoadHistograms(fileName, dir);
 192   fTriggerBiasCorrectionMBToNSD  ->LoadHistograms(fileName, dir);
 193   fTriggerBiasCorrectionMBToND   ->LoadHistograms(fileName, dir);
 194
 195
 196
 197   return kTRUE;
 198 }
 199
 200 //____________________________________________________________________
 201 void
 202 AlidNdEtaCorrection::SaveHistograms() {
 203   //
 204   // save the histograms
 205   //
 206
 207   gDirectory->mkdir(fName.Data());
 208   gDirectory->cd(fName.Data());
 209
 210   fTrack2ParticleCorrection     ->SaveHistograms();
 211   fVertexRecoCorrection         ->SaveHistograms();
 212   fTriggerBiasCorrectionMBToINEL->SaveHistograms();
 213   fTriggerBiasCorrectionMBToNSD ->SaveHistograms();
 214   fTriggerBiasCorrectionMBToND  ->SaveHistograms();
 215
 216   gDirectory->cd("../");
 217 }
 218
 219 //____________________________________________________________________
 220 void AlidNdEtaCorrection::DrawHistograms()
 221 {
 222   //
 223   // call the draw histogram method of the two AliCorrectionMatrix2D objects
 224
 225   fTrack2ParticleCorrection     ->DrawHistograms();
 226   fVertexRecoCorrection         ->DrawHistograms();
 227   fTriggerBiasCorrectionMBToINEL->DrawHistograms();
 228   fTriggerBiasCorrectionMBToNSD ->DrawHistograms();
 229   fTriggerBiasCorrectionMBToND  ->DrawHistograms();
 230
 231 }
 232
 233 //____________________________________________________________________
 234 void AlidNdEtaCorrection::FillEventWithTrigger(Float_t vtx, Float_t n)
 235 {
 236   // fill events with trigger.
 237   // used to calculate vertex reco and trigger bias corrections
 238
 239   fVertexRecoCorrection->FillGene(vtx, n);
 240   fTriggerBiasCorrectionMBToINEL ->FillMeas(vtx, n);
 241   fTriggerBiasCorrectionMBToNSD  ->FillMeas(vtx, n);
 242   fTriggerBiasCorrectionMBToND   ->FillMeas(vtx, n);
 243 }
 244
 245 //____________________________________________________________________
 246 void AlidNdEtaCorrection::FillEventAll(Float_t vtx, Float_t n, Char_t* opt)
 247 {
 248   // fill all events
 249   // used to calculate trigger bias corrections
 250
 251   if (strcmp(opt,"INEL")==0)
 252     fTriggerBiasCorrectionMBToINEL->FillGene(vtx, n);
 253   else if (strcmp(opt,"NSD")==0)
 254     fTriggerBiasCorrectionMBToNSD->FillGene(vtx, n);
 255   else if (strcmp(opt,"ND")==0)
 256     fTriggerBiasCorrectionMBToND->FillGene(vtx, n);
 257   else
 258     AliDebug(AliLog::kWarning, Form(" event type %s unknown (use INEL, NSD or ND)",opt));
 259 }
 260
 261 //____________________________________________________________________
 262 AliCorrectionMatrix2D* AlidNdEtaCorrection::GetTriggerBiasCorrection(Char_t* opt)
 263 {
 264   // returning the trigger bias correction matrix
 265   // option can be used to specify to which collision process (INEL, NSD or ND)
 266   if (strcmp(opt,"INEL")==0)
 267     return fTriggerBiasCorrectionMBToINEL;
 268   else if (strcmp(opt,"NSD")==0)
 269     return fTriggerBiasCorrectionMBToNSD;
 270   else if (strcmp(opt,"ND")==0)
 271     return fTriggerBiasCorrectionMBToND;
 272   else
 273     {
 274       AliDebug(AliLog::kWarning, Form(" %s is unknown (use INEL, NSD or ND). returning INEL ",opt));
 275       return fTriggerBiasCorrectionMBToINEL;
 276     }
 277 }
 278
 279 //____________________________________________________________________
 280 Float_t AlidNdEtaCorrection::GetTriggerBiasCorrection(Float_t vtx, Float_t n, Char_t* opt)
 281 {
 282   // returning the trigger bias correction matrix
 283   // 3rd option can be used to specify to which collision process (INEL, NSD or ND)
 284
 285   if (strcmp(opt,"INEL")==0)
 286     return fTriggerBiasCorrectionMBToINEL->GetCorrection(vtx, n);
 287   else if (strcmp(opt,"NSD")==0)
 288     return fTriggerBiasCorrectionMBToNSD->GetCorrection(vtx, n);
 289   else if (strcmp(opt,"ND")==0)
 290     return fTriggerBiasCorrectionMBToND->GetCorrection(vtx, n);
 291   else {
 292     AliDebug(AliLog::kWarning, Form(" %s unknown (use INEL, NSD or ND). returning corr for INEL",opt));
 293     return fTriggerBiasCorrectionMBToINEL->GetCorrection(vtx, n);
 294   }
 295 }
 296
 297
 298 //____________________________________________________________________
 299 Float_t AlidNdEtaCorrection::GetMeasuredFraction(Float_t ptCutOff, Float_t eta, Bool_t debug)
 300 {
 301   // calculates the fraction of particles measured (some are missed due to the pt cut off)
 302   // uses the generated particle histogram from fTrack2ParticleCorrection
 303
 304   const TH3F* generated = fTrack2ParticleCorrection->GetGeneratedHistogram();
 305
 306   // find eta borders, if eta is negative assume -0.8 ... 0.8
 307   Int_t etaBegin = 0;
 308   Int_t etaEnd = 0;
 309   if (eta < -99)
 310   {
 311     etaBegin = generated->GetYaxis()->FindBin(-0.8);
 312     etaEnd = generated->GetYaxis()->FindBin(0.8);
 313   }
 314   else
 315   {
 316     etaBegin = generated->GetYaxis()->FindBin(eta);
 317     etaEnd = etaBegin;
 318   }
 319
 320   Int_t vertexBegin = generated->GetXaxis()->FindBin(-10);
 321   Int_t vertexEnd = generated->GetXaxis()->FindBin(10);
 322
 323   TH1D* ptProj = dynamic_cast<TH1D*> (generated->ProjectionZ(Form("%s_pt", generated->GetName()), vertexBegin, vertexEnd, etaBegin, etaEnd));
 324   ptProj->GetXaxis()->SetTitle(generated->GetZaxis()->GetTitle());
 325
 326   Int_t ptBin = ptProj->FindBin(ptCutOff);
 327   Float_t abovePtCut = ptProj->Integral(ptBin, ptProj->GetNbinsX());
 328   Float_t all = ptProj->Integral();
 329
 330   if (all == 0)
 331     return -1;
 332
 333   Float_t fraction = abovePtCut / all;
 334
 335   if (debug)
 336   {
 337     new TCanvas;
 338     ptProj->Draw();
 339   }
 340   else
 341     delete ptProj;
 342
 343   if (debug)
 344     printf("AlidNdEtaCorrection::GetMeasuredFraction: pt cut off = %f, eta = %f, => fraction = %f\n", ptCutOff, eta, fraction);
 345
 346   return fraction;
 347 }
 348
 349 void AlidNdEtaCorrection::ReduceInformation()
 350 {
 351   // this function deletes the measured and generated histograms from the corrections to reduce the amount of data
 352   // in memory
 353
 354   // these are needed for GetMeasuredFraction(): fTrack2ParticleCorrection->ReduceInformation();
 355   fVertexRecoCorrection          ->ReduceInformation();
 356   fTriggerBiasCorrectionMBToINEL ->ReduceInformation();
 357   fTriggerBiasCorrectionMBToNSD  ->ReduceInformation();
 358   fTriggerBiasCorrectionMBToND   ->ReduceInformation();
 359 }
 360