PWGLF/FORWARD/analysis2/AliFMDMCHitEnergyFitter.cxx

   1 #include "AliFMDMCHitEnergyFitter.h"
   2 #include "AliMCEvent.h"
   3 #include "AliESDEvent.h"
   4 #include "AliStack.h"
   5 #include "AliFMDHit.h"
   6 #include "AliESDFMD.h"
   7 #include "AliMCAuxHandler.h"
   8 #include "AliGenEventHeader.h"
   9 #include "AliLog.h"
  10 #include "AliFMDEncodedEdx.h"
  11 #include <TH1.h>
  12 #include <TH2.h>
  13 #include <TNtuple.h>
  14 #include <TMath.h>
  15
  16
  17 //____________________________________________________________________
  18 AliFMDMCHitEnergyFitter::AliFMDMCHitEnergyFitter()
  19   : AliFMDEnergyFitter(),
  20     fSumPrimary(0),
  21     fSumSecondary(0),
  22     fIp(),
  23     fNTrack(0),
  24     fNPrimary(0),
  25     fTuple(0)
  26 {
  27 }
  28
  29 //____________________________________________________________________
  30 AliFMDMCHitEnergyFitter::AliFMDMCHitEnergyFitter(const char* title,
  31                                                  Bool_t useTuple)
  32   : AliFMDEnergyFitter(title),
  33     fSumPrimary(0),
  34     fSumSecondary(0),
  35     fIp(),
  36     fNTrack(0),
  37     fNPrimary(0),
  38     fTuple(0)
  39 {
  40   // Some defaults
  41   fUseIncreasingBins = true;
  42   fDoFits            = true;
  43   fDoMakeObject      = false;
  44   fResidualMethod    = kNoResiduals;
  45
  46   if (!useTuple) return;
  47
  48   fTuple = new TNtuple("hits", "Hits",
  49                        "primary:eta:phi:edep:length:"
  50                        "detector:ring:sector:strip:ipz");
  51
  52 }
  53 //____________________________________________________________________
  54 AliFMDMCHitEnergyFitter::~AliFMDMCHitEnergyFitter()
  55 {}
  56
  57 //____________________________________________________________________
  58 void
  59 AliFMDMCHitEnergyFitter::CreateOutputObjects(TList* dir)
  60 {
  61   DGUARD(fDebug, 2, "Create output objects w/MC hits (%p)", dir);
  62   AliFMDEnergyFitter::CreateOutputObjects(dir);
  63   // TList* d = static_cast<TList*>(l->FindObject(GetName()));
  64 }
  65 //____________________________________________________________________
  66 Bool_t
  67 AliFMDMCHitEnergyFitter::PreEvent(const AliMCEvent& mcInput)
  68 {
  69   DGUARD(fDebug, 5, "Reset for event w/MC hits (%p)", &mcInput);
  70   fSumPrimary.Reset(0);
  71   fSumSecondary.Reset(0);
  72
  73   AliMCEvent&        mc        = const_cast<AliMCEvent&>(mcInput);
  74   AliHeader*         header    = mc.Header();
  75   AliStack*          stack     = mc.Stack();
  76   AliGenEventHeader* genHeader = header->GenEventHeader();
  77
  78   genHeader->PrimaryVertex(fIp);
  79   fNTrack   = stack->GetNtrack();
  80   fNPrimary = stack->GetNprimary();
  81
  82   return true;
  83 }
  84
  85 //____________________________________________________________________
  86 Bool_t
  87 AliFMDMCHitEnergyFitter::Event(const AliESDEvent& esdInput,
  88                                const AliMCEvent&  mcInput,
  89                                AliMCAuxHandler&   handler)
  90 {
  91   DGUARD(fDebug, 5, "Process an event for MC hits (%p,%p,%p)",
  92          &esdInput, &mcInput, &handler);
  93   PreEvent(mcInput);
  94
  95   Int_t  nEntries = handler.GetNEntry();
  96   DMSG(fDebug,5, "We got a total of %d particles", nEntries);
  97
  98   for (Int_t i = 0; i < nEntries; i++) {
  99     TClonesArray* hits = handler.GetEntryArray(i);
 100     if (!hits) continue;
 101
 102     AccumulateHits(mcInput, *hits);
 103   }
 104   return PostEvent(esdInput);
 105 }
 106
 107 //____________________________________________________________________
 108 Bool_t
 109 AliFMDMCHitEnergyFitter::AccumulateHits(const AliMCEvent& mcInput,
 110                                         const TClonesArray& hits)
 111 {
 112   DGUARD(fDebug, 5, "Accumulate MC hit energy loss (%p,%p,%d)",
 113          &mcInput, &hits, hits.GetEntries());
 114
 115   AliStack*  stack = const_cast<AliMCEvent&>(mcInput).Stack();
 116   Float_t    cache[10];
 117   Int_t      nHit  = hits.GetEntries();
 118   Int_t      oldTr = -1;
 119   UShort_t   oldD  = 0;
 120   Char_t     oldR  = '\0';
 121   for (Int_t j = 0; j < nHit; j++) {
 122
 123     // Check the hit
 124     AliFMDHit* hit  = static_cast<AliFMDHit*>(hits.At(j));
 125     if (!hit) continue;
 126
 127     // Zero fill array
 128     if (fTuple) for (Int_t k = 0; k < 10; k++) cache[k] = 0;
 129
 130     // Get the track number
 131     Int_t iTr = hit->Track();
 132     if (iTr < 0 || iTr >= fNTrack)
 133       AliWarningF("Track # %d seems out of bounds [0,%d]",
 134                   iTr, fNTrack-1);
 135
 136     // Get the track
 137     AliMCParticle* particle =
 138       static_cast<AliMCParticle*>(mcInput.GetTrack(iTr));
 139
 140     // Particle type 0: unknown, 1: primary, 2: secondary
 141     // Int_t flag = 0; - not used at the moment
 142
 143     // Check if this charged and a primary
 144     if (particle) {
 145       // Only charged tracks
 146       if (!particle->Charge() != 0) continue;
 147     }
 148     Bool_t   isPrimary = stack->IsPhysicalPrimary(iTr) && iTr < fNPrimary;
 149     Double_t eloss     = hit->Edep();
 150     Double_t length    = hit->Length();
 151     // Double_t eta  = particle->Eta();
 152     Double_t x         = hit->X() - fIp[0];
 153     Double_t y         = hit->Y() - fIp[1];
 154     Double_t z         = hit->Z() - fIp[2];
 155     Double_t r         = TMath::Sqrt(x*x + y*y);
 156     Double_t phi       = TMath::ATan2(y, x);
 157     Double_t theta     = TMath::ATan2(r, z);
 158     Double_t eta       = -TMath::Log(TMath::Tan(theta/2));
 159     UShort_t det       = hit->Detector();
 160     Char_t   rng       = hit->Ring();
 161     UShort_t sec       = hit->Sector();
 162     UShort_t str       = hit->Strip();
 163
 164     Bool_t isDistinct = (iTr != oldTr) || (det != oldD) || (rng != oldR);
 165     if (isDistinct) {
 166       RingHistos* h = static_cast<RingHistos*>(GetRingHistos(det, rng));
 167       h->fKind->Fill(eta, isPrimary ? .5 : 1.5);
 168     }
 169     oldTr = iTr;
 170     oldD  = det;
 171     oldR  = rng;
 172
 173     AliFMDFloatMap& sum     =  (isPrimary ? fSumPrimary : fSumSecondary);
 174     // Float_t         old     =  sum(det, rng, sec, str);
 175     sum(det, rng, sec, str) += eloss;
 176 #if 0
 177     Info("", "FMD%d%c[%02d,%03d]-%s: Adding %10.7f to %10.7f -> %10.7f",
 178          det, rng, sec, str, (isPrimary ? "1st" : "2nd"),
 179          eloss, old, sum(det, rng, sec, str));
 180 #endif
 181     if (!fTuple) continue;
 182
 183     // Leaves primary:eta:phi:edep:length:detector:ring:sector:strip:ipz
 184     cache[0] = isPrimary ? 1 : 0;
 185     cache[1] = eta;
 186     cache[2] = phi;
 187     cache[3] = eloss;
 188     cache[4] = length;
 189     cache[5] = hit->Detector();
 190     cache[6] = Int_t(hit->Ring());
 191     cache[7] = hit->Sector();
 192     cache[8] = hit->Strip();
 193     cache[9] = fIp[2];
 194     fTuple->Fill(cache);
 195   }
 196   return true;
 197 }
 198
 199 //____________________________________________________________________
 200 Bool_t AliFMDMCHitEnergyFitter::PostEvent(const AliESDEvent& input)
 201 {
 202   DGUARD(fDebug, 5, "Convert MC hit energy loss (%p)", &input);
 203
 204   AliESDFMD*  esdFMD    = input.GetFMDData();
 205   for (UShort_t d = 1; d <= 3; d++) {
 206     UShort_t nQ = (d == 1 ? 1 : 2);
 207     for (UShort_t q = 0; q < nQ; q++) {
 208       UShort_t    nS = (q == 0 ?  20 :  40);
 209       UShort_t    nT = (q == 0 ? 512 : 256);
 210       Char_t      r  = (q == 0 ? 'I' : 'O');
 211       RingHistos* h  =
 212         static_cast<AliFMDMCHitEnergyFitter::RingHistos*>(GetRingHistos(d,r));
 213       if (!h) continue;
 214
 215       for (UShort_t s = 0; s < nS; s++) {
 216         for (UShort_t t = 0; t < nT; t++) {
 217           Float_t  totPrim = fSumPrimary(d, r, s, t);
 218           Float_t  totSec  = fSumSecondary(d, r, s, t);
 219           Float_t  totAll  = totPrim + totSec;
 220           Double_t esdEta  = esdFMD->Eta(d, r, s, t);
 221           if (totAll  > 0) h->FillMC(0,  esdEta, totAll);
 222           if (totPrim > 0) h->FillMC(1,  esdEta, totPrim);
 223           if (totSec  > 0) h->FillMC(2,  esdEta, totSec);
 224         }
 225       }
 226     }
 227   }
 228   return true;
 229 }
 230
 231 //____________________________________________________________________
 232 AliFMDEnergyFitter::RingHistos*
 233 AliFMDMCHitEnergyFitter::CreateRingHistos(UShort_t d, Char_t r) const
 234 {
 235   // DGUARD(fDebug, 1, "Create Ring cache of MC hit energy loss (FMD%d%c)",d,r);
 236   return new AliFMDMCHitEnergyFitter::RingHistos(d,r);
 237 }
 238
 239 //====================================================================
 240 AliFMDMCHitEnergyFitter::RingHistos::RingHistos()
 241   : AliFMDEnergyFitter::RingHistos(),
 242     fPrimary(0),
 243     fSecondary(0),
 244     fKind(0)
 245 {}
 246
 247 //____________________________________________________________________
 248 AliFMDMCHitEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
 249   : AliFMDEnergyFitter::RingHistos(d,r),
 250     fPrimary(0),
 251     fSecondary(0),
 252     fKind(0)
 253 {}
 254 //____________________________________________________________________
 255 TArrayD
 256 AliFMDMCHitEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t,
 257                                                         Double_t,
 258                                                         Double_t) const
 259 {
 260   // Make an increasing axis for ELoss distributions.
 261   //
 262   // We use the service function of AliFMDEncodedEdx to do this
 263   TArrayD ret;
 264   const AliFMDEncodedEdx::Spec& s = AliFMDEncodedEdx::GetdEAxis();
 265   s.FillBinArray(ret);
 266
 267   return ret;
 268 }
 269 //____________________________________________________________________
 270 void
 271 AliFMDMCHitEnergyFitter::RingHistos::SetupForData(const TAxis& eAxis,
 272                                                   const TAxis& /*cAxis*/,
 273                                                   Double_t maxDE,
 274                                                   Int_t    nDEbins,
 275                                                   Bool_t   useIncrBin)
 276 {
 277   // AliFMDEnergyFitter::RingHistos::SetupForData(eAxis, cAxis, maxDE, nDEbins,
 278   // useIncrBin);
 279
 280   fHist      = Make("eloss", "#sum#Delta_{true} of all",
 281                     eAxis, maxDE, nDEbins, useIncrBin);
 282   fPrimary   = Make("primary", "#sum#Delta_{true} of primaries",
 283                     eAxis, maxDE, nDEbins, useIncrBin);
 284   fSecondary = Make("secondary","#sum#Delta_{true} of secondaries",
 285                     eAxis, maxDE, nDEbins, useIncrBin);
 286   fHist->SetXTitle("#Delta_{true}");
 287   fSecondary->SetXTitle("#Delta_{true}");
 288   fPrimary->SetXTitle("#Delta_{true}");
 289
 290   fKind       = 0;
 291   if (eAxis.GetXbins()->GetArray())
 292     fKind = new TH2I("kind", "Particle types",
 293                      eAxis.GetNbins(), eAxis.GetXbins()->GetArray(),
 294                      2, 0, 2);
 295   else
 296     fKind = new TH2I("kind", "Particle types",
 297                      eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(),
 298                      2, 0, 2);
 299   fKind->SetXTitle("#eta");
 300   fKind->GetYaxis()->SetBinLabel(1, "Primary");
 301   fKind->GetYaxis()->SetBinLabel(2, "Secondary");
 302   fKind->SetDirectory(0);
 303   fKind->Sumw2();
 304
 305   fList->Add(fHist);
 306   fList->Add(fPrimary);
 307   fList->Add(fSecondary);
 308   fList->Add(fKind);
 309 }
 310
 311 //____________________________________________________________________
 312 void
 313 AliFMDMCHitEnergyFitter::RingHistos::FillMC(UShort_t flag, Double_t eta,
 314                                             Double_t mult)
 315 {
 316   switch (flag) {
 317     // case 0: AliFMDEnergyFitter::RingHistos::Fill(false, eta, 0, mult); break;
 318   case 0: fHist->Fill(eta, mult); break;
 319   case 1: fPrimary->Fill(eta, mult); break;
 320   case 2: fSecondary->Fill(eta, mult); break;
 321   }
 322 }
 323
 324 //____________________________________________________________________
 325 TObjArray*
 326 AliFMDMCHitEnergyFitter::RingHistos::Fit(TList*           dir,
 327                                          Double_t         lowCut,
 328                                          UShort_t         nParticles,
 329                                          UShort_t         minEntries,
 330                                          UShort_t         minusBins,
 331                                          Double_t         relErrorCut,
 332                                          Double_t         chi2nuCut,
 333                                          Double_t         minWeight,
 334                                          Double_t         regCut,
 335                                          EResidualMethod  residuals) const
 336 {
 337   TObjArray* all  = FitSlices(dir, "eloss", lowCut, nParticles, minEntries,
 338                               minusBins, relErrorCut, chi2nuCut, minWeight,
 339                               regCut, residuals, false, 0);
 340   TObjArray* prim = FitSlices(dir, "primary", lowCut, nParticles, minEntries,
 341                               minusBins, relErrorCut, chi2nuCut, minWeight,
 342                               regCut, residuals, false, 0);
 343   TObjArray* sec  = FitSlices(dir, "secondary", lowCut, nParticles,
 344                               minEntries, minusBins, relErrorCut, chi2nuCut,
 345                               minWeight, regCut, residuals, false, 0);
 346   if (!all || !prim || !sec) {
 347     AliWarningF("Failed to get results for all(%p), primary(%p), and/or "
 348                 "secondary(%p)", all, prim, sec);
 349     return 0;
 350   }
 351   // Already added to the sub-folders
 352   // dir->Add(all);
 353   // dir->Add(prim);
 354   // dir->Add(sec);
 355
 356   Int_t      nPrim = prim->GetEntriesFast();
 357   TObjArray* out   = new TObjArray;
 358   for (Int_t i = 0; i < nPrim; i++) {
 359     TH1* h = static_cast<TH1*>(prim->At(i));
 360     if (!h) continue;
 361
 362     TAxis* xAxis  = h->GetXaxis();
 363     TH2*   hh     = 0;
 364     if (xAxis->GetXbins()->GetArray())
 365       hh = new TH2D(h->GetName(), h->GetTitle(),
 366                     xAxis->GetNbins(), xAxis->GetXbins()->GetArray(),
 367                     3, 0, 3);
 368     else
 369       hh = new TH2D(h->GetName(), h->GetTitle(),
 370                     xAxis->GetNbins(), xAxis->GetXmin(), xAxis->GetXmax(),
 371                     3, 0, 3);
 372     hh->GetYaxis()->SetBinLabel(1, "Primary");
 373     hh->GetYaxis()->SetBinLabel(2, "Secondary");
 374     hh->GetYaxis()->SetBinLabel(3, "All");
 375     hh->GetXaxis()->SetTitle("#eta");
 376     out->Add(hh);
 377   }
 378   for (Int_t i = 0; i < nPrim; i++) {
 379     TH2* h = static_cast<TH2*>(out->At(i));
 380     if (!h) continue;
 381
 382     TH1* hp = static_cast<TH1*>(prim->At(i));
 383     TH1* hs = static_cast<TH1*>(sec->At(i));
 384     TH1* ha = static_cast<TH1*>(all->At(i));
 385     TH1* hh[] = { hp, hs, ha };
 386     for (Int_t j = 0; j < 3; j++) {
 387       TH1* ph = hh[j];
 388       if (!ph) continue;
 389
 390       for (Int_t k = 1; k <= ph->GetNbinsX(); k++) {
 391         Double_t c = ph->GetBinContent(k);
 392         Double_t e = ph->GetBinError(k);
 393         h->SetBinContent(k, j+1, c);
 394         h->SetBinError(k, j+1, e);
 395       }
 396     }
 397   }
 398   TList* l = GetOutputList(dir);
 399   if (!l) return 0;
 400
 401   out->SetName("compartive");
 402   out->SetOwner();
 403   l->Add(out);
 404   return out;
 405 }