PWGLF/FORWARD/analysis2/AliCentralMCCorrectionsTask.cxx

   1 //
   2 // Calculate the corrections in the central regions
   3 //
   4 // Inputs:
   5 //   - AliESDEvent
   6 //
   7 // Outputs:
   8 //   - AliAODCentralMult
   9 //
  10 // Histograms
  11 //
  12 // Corrections used
  13 //
  14 #include "AliCentralMCCorrectionsTask.h"
  15 #include "AliCentralCorrectionManager.h"
  16 #include "AliTriggerAnalysis.h"
  17 #include "AliPhysicsSelection.h"
  18 #include "AliLog.h"
  19 #include "AliHeader.h"
  20 #include "AliGenEventHeader.h"
  21 #include "AliESDEvent.h"
  22 #include "AliAODHandler.h"
  23 #include "AliMultiplicity.h"
  24 #include "AliInputEventHandler.h"
  25 #include "AliStack.h"
  26 #include "AliMCEvent.h"
  27 #include "AliAODForwardMult.h"
  28 #include "AliCentralCorrSecondaryMap.h"
  29 #include "AliCentralCorrAcceptance.h"
  30 #include "AliForwardUtil.h"
  31 #include "AliMultiplicity.h"
  32 #include <TH1.h>
  33 #include <TH2D.h>
  34 #include <TDirectory.h>
  35 #include <TList.h>
  36 #include <TROOT.h>
  37 #include <iostream>
  38
  39 //====================================================================
  40 namespace {
  41   const char* GetEventName(Bool_t tr, Bool_t vtx)
  42   {
  43     return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
  44   }
  45 }
  46
  47 //====================================================================
  48 AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask()
  49   : AliBaseMCCorrectionsTask(),
  50     fTrackDensity(),
  51     fSecCorr(0),
  52     fAccCorr(0),
  53     fNPhiBins(20),
  54     fEffectiveCorr(true),
  55     fEtaCut(1.9),
  56     fCorrCut(0.6)
  57 {
  58   //
  59   // Constructor
  60   //
  61   // Parameters:
  62   //    name Name of task
  63   //
  64   DGUARD(fDebug, 3,"Default CTOR of AliCentralMCCorrectionsTask");
  65 }
  66
  67 //____________________________________________________________________
  68 AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const char* name)
  69   : AliBaseMCCorrectionsTask(name, &(AliCentralCorrectionManager::Instance())),
  70     fTrackDensity("trackDensity"),
  71     fSecCorr(0),
  72     fAccCorr(0),
  73     fNPhiBins(20),
  74     fEffectiveCorr(true),
  75     fEtaCut(1.9),
  76     fCorrCut(0.6)
  77 {
  78   //
  79   // Constructor
  80   //
  81   // Parameters:
  82   //    name Name of task
  83   //
  84   DGUARD(fDebug, 3,"Named CTOR of AliCentralMCCorrectionsTask: %s",name);
  85 }
  86
  87 //____________________________________________________________________
  88 AliBaseMCCorrectionsTask::VtxBin*
  89 AliCentralMCCorrectionsTask::CreateVtxBin(Double_t low, Double_t high)
  90 {
  91   return new AliCentralMCCorrectionsTask::VtxBin(low,high,fEtaAxis,fNPhiBins);
  92 }
  93
  94 //____________________________________________________________________
  95 Bool_t
  96 AliCentralMCCorrectionsTask::ProcessESD(const AliESDEvent& esd,
  97                                         const AliMCEvent& mc,
  98                                         AliBaseMCCorrectionsTask::VtxBin& bin,
  99                                         Double_t          vz)
 100 {
 101   AliCentralMCCorrectionsTask::VtxBin& vb =
 102     static_cast<AliCentralMCCorrectionsTask::VtxBin&>(bin);
 103
 104   // Now process our input data and store in own ESD object
 105   fTrackDensity.Calculate(mc, vz, *vb.fHits, bin.fPrimary);
 106
 107   // Get the ESD object
 108   const AliMultiplicity* spdmult = esd.GetMultiplicity();
 109
 110   // Count number of tracklets per bin
 111   for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
 112     vb.fClusters->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
 113   //...and then the unused clusters in layer 1
 114   for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
 115     Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
 116     vb.fClusters->Fill(eta, spdmult->GetPhiSingle(j));
 117   }
 118
 119   // Count events
 120   bin.fCounts->Fill(0.5);
 121
 122   return true;
 123 }
 124 //____________________________________________________________________
 125 void
 126 AliCentralMCCorrectionsTask::CreateCorrections(TList* results)
 127 {
 128   fSecCorr = new AliCentralCorrSecondaryMap;
 129   fSecCorr->SetVertexAxis(fVtxAxis);
 130
 131   fAccCorr = new AliCentralCorrAcceptance;
 132   fAccCorr->SetVertexAxis(fVtxAxis);
 133
 134   results->Add(fSecCorr);
 135   results->Add(fAccCorr);
 136 }
 137
 138 //____________________________________________________________________
 139 Bool_t
 140 AliCentralMCCorrectionsTask::FinalizeVtxBin(AliBaseMCCorrectionsTask::VtxBin*
 141                                             bin,  UShort_t iVz)
 142 {
 143
 144   AliCentralMCCorrectionsTask::VtxBin* vb =
 145     static_cast<AliCentralMCCorrectionsTask::VtxBin*>(bin);
 146   vb->Terminate(fList, fResults, iVz, fEffectiveCorr,
 147                 fEtaCut, fCorrCut, fSecCorr,fAccCorr);
 148   return true;
 149 }
 150
 151
 152 //____________________________________________________________________
 153 void
 154 AliCentralMCCorrectionsTask::Print(Option_t* option) const
 155 {
 156   AliBaseMCCorrectionsTask::Print(option);
 157   std::cout << "  # of phi bins:    " << fNPhiBins << "\n"
 158             << "  Effective corr.:  " << fEffectiveCorr << "\n"
 159             << "  Eta cut-off:      " << fEtaCut << "\n"
 160             << "  Acceptance cut:   " << fCorrCut
 161             << std::endl;
 162   gROOT->IncreaseDirLevel();
 163   fTrackDensity.Print(option);
 164   gROOT->DecreaseDirLevel();
 165 }
 166
 167 //====================================================================
 168 AliCentralMCCorrectionsTask::VtxBin::VtxBin()
 169   : AliBaseMCCorrectionsTask::VtxBin(),
 170     fHits(0),
 171     fClusters(0)
 172 {
 173 }
 174 //____________________________________________________________________
 175 AliCentralMCCorrectionsTask::VtxBin::VtxBin(Double_t     low,
 176                                             Double_t     high,
 177                                             const TAxis& axis,
 178                                             UShort_t     nPhi)
 179   : AliBaseMCCorrectionsTask::VtxBin(low, high, axis, nPhi),
 180     fHits(0),
 181     fClusters(0)
 182 {
 183   fHits = static_cast<TH2D*>(fPrimary->Clone("hits"));
 184   fHits->SetTitle("Hits");
 185   fHits->SetDirectory(0);
 186
 187   fClusters = static_cast<TH2D*>(fPrimary->Clone("clusters"));
 188   fClusters->SetTitle("Clusters");
 189   fClusters->SetDirectory(0);
 190 }
 191
 192
 193 //____________________________________________________________________
 194 TList*
 195 AliCentralMCCorrectionsTask::VtxBin::CreateOutputObjects(TList* l)
 196 {
 197   TList* d = AliBaseMCCorrectionsTask::VtxBin::CreateOutputObjects(l);
 198
 199   d->Add(fHits);
 200   d->Add(fClusters);
 201
 202   return d;
 203 }
 204 //____________________________________________________________________
 205 void
 206 AliCentralMCCorrectionsTask::VtxBin::Terminate(const TList* input,
 207                                                TList* output,
 208                                                UShort_t iVz,
 209                                                Bool_t effectiveCorr,
 210                                                Double_t etaCut,
 211                                                Double_t accCut,
 212                                                AliCentralCorrSecondaryMap* map,
 213                                                AliCentralCorrAcceptance* acorr)
 214 {
 215   TList* out = new TList;
 216   out->SetName(GetName());
 217   out->SetOwner();
 218   output->Add(out);
 219
 220   TList* l = static_cast<TList*>(input->FindObject(GetName()));
 221   if (!l) {
 222     AliError(Form("List %s not found in %s", GetName(), input->GetName()));
 223     return;
 224   }
 225
 226   // Get the sums
 227   TH2D*   hits  = static_cast<TH2D*>(l->FindObject("hits"));
 228   TH2D*   clus  = static_cast<TH2D*>(l->FindObject("clusters"));
 229   TH2D*   prim  = static_cast<TH2D*>(l->FindObject("primary"));
 230   if (!hits || !prim) {
 231     AliError(Form("Missing histograms: %p, %p", hits, prim));
 232     return;
 233   }
 234
 235   // Clone cluster and hit map
 236   TH2D* secMapEff = static_cast<TH2D*>(clus->Clone("secMapEff"));
 237   TH2D* secMapHit = static_cast<TH2D*>(hits->Clone("secMapHit"));
 238   secMapEff->SetTitle("2^{nd} map from clusters");
 239   secMapEff->SetDirectory(0);
 240   secMapHit->SetTitle("2^{nd} map from MC hits");
 241   secMapHit->SetDirectory(0);
 242
 243   // Divide cluster and hit map with number of primaries
 244   secMapEff->Divide(prim);
 245   secMapHit->Divide(prim);
 246
 247   // Create acceptance histograms
 248   TH1D* accEff = new TH1D("accEff",
 249                           "Acceptance correction for SPD (from 2^{nd} map)" ,
 250                           fPrimary->GetXaxis()->GetNbins(),
 251                           fPrimary->GetXaxis()->GetXmin(),
 252                           fPrimary->GetXaxis()->GetXmax());
 253   TH1D* accHit = static_cast<TH1D*>(accEff->Clone("accHit"));
 254   accHit->SetTitle("Acceptance correction for SPD (from clusters)");
 255
 256   // Diagnostics histogra,
 257   TH2*  dia    = static_cast<TH2D*>(clus->Clone("diagnostics"));
 258   dia->SetTitle("Scaled cluster density");
 259
 260   // Total number of channels along phi and # of eta bins
 261   Int_t nTotal = secMapHit->GetNbinsY();
 262   Int_t nEta   = secMapHit->GetNbinsX();
 263
 264   for(Int_t xx = 1; xx <= nEta; xx++) {
 265     Double_t eta = secMapHit->GetXaxis()->GetBinCenter(xx);
 266     Bool_t   ins = TMath::Abs(eta) <= etaCut;
 267     Double_t mm  = 0;
 268     if (ins) {
 269       // Find the maximum cluster signal in this phi range
 270       for (Int_t yy = 1; yy <= nTotal; yy++) {
 271         Double_t c = clus->GetBinContent(xx,yy);
 272         mm         = TMath::Max(mm, c);
 273       }
 274     }
 275     // Count number of phi bins with enough clusters or high enough
 276     // correction.
 277     Int_t nOKEff    = 0;
 278     Int_t nOKHit    = 0;
 279     for(Int_t yy = 1; yy <=nTotal; yy++) {
 280       if (!ins) { // Not inside Eta cut
 281         secMapEff->SetBinContent(xx,yy,0.);
 282         secMapEff->SetBinError(xx,yy,0.);
 283         secMapHit->SetBinContent(xx,yy,0.);
 284         secMapHit->SetBinError(xx,yy,0.);
 285         dia->SetBinContent(xx,yy,0);
 286         continue;
 287       }
 288
 289       // Check if the background correction is large enough, or zero map
 290       if(secMapEff->GetBinContent(xx,yy) > 0.9) {
 291         // acc->Fill(h->GetXaxis()->GetBinCenter(xx));
 292         nOKEff++;
 293       }
 294       else {
 295         secMapEff->SetBinContent(xx,yy,0.);
 296         secMapEff->SetBinError(xx,yy,0.);
 297       }
 298
 299       // Check if the number of cluster is large enough, or zero map
 300       Double_t c = clus->GetBinContent(xx,yy);
 301       Double_t s = (mm < 1e-6) ? 0 : c / mm;
 302       dia->SetBinContent(xx,yy,s);
 303       if (s >= accCut) {
 304         nOKHit++;
 305       }
 306       else {
 307         secMapHit->SetBinContent(xx,yy,0);
 308         secMapHit->SetBinError(xx,yy,0);
 309       }
 310     }
 311
 312     // Calculate acceptance as ratio of bins with enough clusters and
 313     // total number of phi bins.
 314     Double_t accXX = float(nOKHit) / nTotal;
 315     if (accXX < 0.2) accXX = 0;
 316     accHit->SetBinContent(xx, accXX);
 317
 318     // Calculate acceptance as ratio of bins with large enough
 319     // correction and total number of phi bins.
 320     accXX = float(nOKEff) / nTotal;
 321     if (accXX < 0.2) accXX = 0;
 322     accEff->SetBinContent(xx, accXX);
 323   }
 324
 325   TH2D* secMap    = (effectiveCorr ? secMapEff : secMapHit);
 326   TH2D* secMapAlt = (effectiveCorr ? secMapHit : secMapEff);
 327   TH1D* acc       = (effectiveCorr ? accEff    : accHit);
 328   TH1D* accAlt    = (effectiveCorr ? accHit    : accEff);
 329   out->Add(secMap->Clone("secMap"));
 330   out->Add(secMapAlt->Clone());
 331   out->Add(acc->Clone("acc"));
 332   out->Add(accAlt->Clone());
 333   out->Add(dia->Clone());
 334
 335   map->SetCorrection(iVz, secMap);
 336   acorr->SetCorrection(iVz, acc);
 337 }
 338
 339 //
 340 // EOF
 341 //