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