PWGLF/FORWARD/analysis2/AliFMDMCTrackELoss.cxx

   1 #include "AliFMDMCTrackELoss.h"
   2 #include "AliESDFMD.h"
   3 #include "AliTrackReference.h"
   4 #include <TMath.h>
   5 #include "AliFMDStripIndex.h"
   6 #include "AliFMDEncodedEdx.h"
   7 #include "AliLog.h"
   8 #include "AliMCParticle.h"
   9 #include <TH2D.h>
  10 #include <TH1D.h>
  11 #include <TList.h>
  12 #include <TROOT.h>
  13 #include <TTree.h>
  14 #include <TLorentzVector.h>
  15 #include <TParticle.h>
  16 #include <TClonesArray.h>
  17 #include <iostream>
  18
  19 //====================================================================
  20 void
  21 AliFMDMCTrackELoss::State::Clear(Bool_t alsoCount)
  22 {
  23   angle       = 0;
  24   oldDetector = 0;
  25   oldRing     = '\0';
  26   oldSector   = 1024;
  27   oldStrip    = 1024;
  28   startStrip  = 1024;
  29   nRefs       = 0;
  30   nStrips     = 0;
  31   longest     = 0x0;
  32   de          = 0;
  33   dx          = 0;
  34   if (alsoCount) count = 0;
  35 }
  36
  37 //____________________________________________________________________
  38 AliFMDMCTrackELoss::State&
  39 AliFMDMCTrackELoss::State::operator=(const State& o)
  40 {
  41   if (&o == this) return *this;
  42   angle          = o.angle;
  43   oldDetector    = o.oldDetector;
  44   oldRing        = o.oldRing;
  45   oldSector      = o.oldSector;
  46   oldStrip       = o.oldStrip;
  47   startStrip     = o.startStrip;
  48   nRefs          = o.nRefs;
  49   nStrips        = o.nStrips;
  50   count          = o.count;
  51   longest        = o.longest;
  52   de             = o.de;
  53   dx             = o.dx;
  54   return *this;
  55 }
  56
  57 //====================================================================
  58 AliFMDMCTrackELoss::Hit::Hit()
  59   : fGamma(0), fBeta(0), fEta(0), fDe(0), fDx(0), fDetId(0), fPdg(0),
  60     fPrimary(true), fDetector(0), fRing('\0'),
  61     fSector(0xFFFF), fStrip(0xFFFF)
  62 {}
  63 //____________________________________________________________________
  64 void
  65 AliFMDMCTrackELoss::Hit::Decode() const
  66 {
  67   if (fDetector > 0) return;
  68   AliFMDStripIndex::Unpack(fDetId, fDetector, fRing, fSector, fStrip);
  69 }
  70 //____________________________________________________________________
  71 UInt_t
  72 AliFMDMCTrackELoss::Hit::AbsPdg() const
  73 {
  74   return TMath::Abs(fPdg);
  75 }
  76
  77 //====================================================================
  78 AliFMDMCTrackELoss::AliFMDMCTrackELoss()
  79   : AliBaseMCTrackDensity(),
  80     fState(),
  81     fEvent(),
  82     fMaxConsequtiveStrips(3),
  83     fUseTree(false),
  84     fHits(0),
  85     fTree(0),
  86     fNr(0),
  87     fNt(0),
  88     fNc(0),
  89     fNcr(0),
  90     fBetaGammadEdx(0),
  91     fBetaGammaEta(0),
  92     fDEdxEta(0),
  93     fPrimaries(0),
  94     fSecondaries(0),
  95     fAll(0),
  96     fEta(0)
  97 {
  98   // Default constructor
  99 }
 100
 101 //____________________________________________________________________
 102 AliFMDMCTrackELoss::AliFMDMCTrackELoss(const char*)
 103   : AliBaseMCTrackDensity("fmdMCTrackELoss"),
 104     fState(),
 105     fEvent(),
 106     fMaxConsequtiveStrips(3),
 107     fUseTree(false),
 108     fHits(0),
 109     fTree(0),
 110     fNr(0),
 111     fNt(0),
 112     fNc(0),
 113     fNcr(0),
 114     fBetaGammadEdx(0),
 115     fBetaGammaEta(0),
 116     fDEdxEta(0),
 117     fPrimaries(0),
 118     fSecondaries(0),
 119     fAll(0),
 120     fEta(0)
 121 {
 122   // Normal constructor constructor
 123   SetTitle("mcTrackELoss");
 124   fHits = new TClonesArray("AliFMDMCTrackELoss::Hit");
 125 }
 126
 127 #if 0
 128 //____________________________________________________________________
 129 AliFMDMCTrackELoss::AliFMDMCTrackELoss(const AliFMDMCTrackELoss& o)
 130   : AliBaseMCTrackDensity(o),
 131     fState(o.fState),
 132     fMaxConsequtiveStrips(o.fMaxConsequtiveStrips),
 133     fNr(o.fNr),
 134     fNt(o.fNt),
 135     fNc(o.fNc),
 136     fNcr(o.fNcr),
 137     fOutput(o.fOutput)
 138 {
 139   // Normal constructor constructor
 140 }
 141
 142 //____________________________________________________________________
 143 AliFMDMCTrackELoss&
 144 AliFMDMCTrackELoss::operator=(const AliFMDMCTrackELoss& o)
 145 {
 146   // Assignment operator
 147   if (&o == this) return *this;
 148   AliBaseMCTrackDensity::operator=(o);
 149   fMaxConsequtiveStrips = o.fMaxConsequtiveStrips;
 150   fUseTree             = o.fUseTree;
 151   fNr                   = o.fNr;
 152   fNt                   = o.fNt;
 153   fNc                   = o.fNc;
 154   fNcr                  = o.fNcr;
 155   fState                = o.fState;
 156
 157   if (o.fTree)
 158     fTree               = static_cast<TTree*>(o.fTree->Clone());
 159
 160   return *this;
 161 }
 162 #endif
 163
 164 //____________________________________________________________________
 165 void
 166 AliFMDMCTrackELoss::CreateOutputObjects(TList* l)
 167 {
 168   AliBaseMCTrackDensity::CreateOutputObjects(l);
 169   TList* ll = static_cast<TList*>(l->FindObject(GetTitle()));
 170   if (!ll) ll = l;
 171
 172   fNr = new TH1D("clusterRefs", "# track references per cluster",
 173                  21, -.5, 20.5);
 174   fNr->SetXTitle("# of track references/cluster");
 175   fNr->SetFillColor(kRed+1);
 176   fNr->SetFillStyle(3001);
 177   fNr->SetDirectory(0);
 178   ll->Add(fNr);
 179
 180   fNt = new TH1D("clusterSize", "cluster length in strips", 21, -.5, 20.5);
 181   fNt->SetXTitle("Cluster size [strips]");
 182   fNt->SetFillColor(kBlue+1);
 183   fNt->SetFillStyle(3001);
 184   fNt->SetDirectory(0);
 185   ll->Add(fNt);
 186
 187   fNc = new TH1D("nClusters", "# clusters per track", 21, -.5, 20.5);
 188   fNc->SetXTitle("# clusters");
 189   fNc->SetFillColor(kGreen+1);
 190   fNc->SetFillStyle(3001);
 191   fNc->SetDirectory(0);
 192   ll->Add(fNc);
 193
 194   fNcr = new TH2D("clusterVsRefs", "# clusters vs. # refs",
 195                   21, -.5, 20.5, 21, -.5, 20.5);
 196   fNcr->SetXTitle("# References");
 197   fNcr->SetYTitle("# Clusters");
 198   fNcr->SetOption("COLZ");
 199   fNcr->SetDirectory(0);
 200   ll->Add(fNcr);
 201
 202
 203   TArrayD edepArray;
 204   TArrayD bgArray(601);
 205   AliFMDEncodedEdx::GetdEAxis().FillBinArray(edepArray);
 206   AliForwardUtil::MakeLogScale(600, -2, 5, bgArray);
 207   for (Int_t i = 0; i < edepArray.GetSize(); i++)
 208     edepArray[i] *= 10;
 209
 210   fBetaGammadEdx = new TH2D("betaGammadEdx", "Energy loss",
 211                             bgArray.GetSize()-1, bgArray.GetArray(),
 212                             edepArray.GetSize()-1, edepArray.GetArray());
 213   fBetaGammadEdx->SetDirectory(0);
 214   fBetaGammadEdx->SetXTitle("#it{#beta#gamma}");
 215   fBetaGammadEdx->SetYTitle("d#it{#Delta}/d#it{x}");
 216   fBetaGammadEdx->Sumw2();
 217   ll->Add(fBetaGammadEdx);
 218
 219   fBetaGammaEta = new TH2D("betaGammaEta", "#beta#gamma",
 220                            200, -4, 6, bgArray.GetSize()-1, bgArray.GetArray());
 221   fBetaGammaEta->SetXTitle("#eta");
 222   fBetaGammaEta->SetYTitle("#it{#beta#gamma}");
 223   fBetaGammaEta->Sumw2();
 224   ll->Add(fBetaGammaEta);
 225
 226   fDEdxEta = new TH2D("dEdxEta", "d#it{#Delta}/d#it{x}",
 227                       200, -4, 6, edepArray.GetSize()-1, edepArray.GetArray());
 228   fDEdxEta->SetXTitle("#eta");
 229   fDEdxEta->SetYTitle("d#it{#Delta}/d#it{x}");
 230   fDEdxEta->Sumw2();
 231   ll->Add(fDEdxEta);
 232
 233   if (!fUseTree) return;
 234
 235   fTree  = new TTree("tree", "Tree of hits");
 236   fTree->Branch("event", &fEvent, "ipZ/D:cent");
 237   fTree->Branch("hits", &fHits);
 238 }
 239
 240 //____________________________________________________________________
 241 void
 242 AliFMDMCTrackELoss::Clear(Option_t*)
 243 {
 244   fPrimaries.Reset(0);
 245   fSecondaries.Reset(0);
 246   fAll.Reset(0);
 247   fEta.Reset(AliESDFMD::kInvalidEta);
 248   fHits->Clear();
 249 }
 250
 251 //____________________________________________________________________
 252 Int_t
 253 AliFMDMCTrackELoss::GetDetectorId() const
 254 {
 255   return AliTrackReference::kFMD;
 256 }
 257
 258 //____________________________________________________________________
 259 void
 260 AliFMDMCTrackELoss::BeginTrackRefs()
 261 {
 262   fState.Clear(true);
 263 }
 264
 265 //____________________________________________________________________
 266 void
 267 AliFMDMCTrackELoss::EndTrackRefs(Int_t nRefs)
 268 {
 269   fNc->Fill(fState.count);
 270   fNcr->Fill(nRefs, fState.count);
 271   fState.Clear(true);
 272 }
 273
 274 //____________________________________________________________________
 275 AliTrackReference*
 276 AliFMDMCTrackELoss::ProcessRef(AliMCParticle*       particle,
 277                                const AliMCParticle* mother,
 278                                AliTrackReference*   ref)
 279 {
 280   // Get the detector coordinates
 281   UShort_t d, s, t;
 282   Char_t r;
 283   AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
 284   Double_t edep, length, dEdep, dLength;
 285   AliFMDEncodedEdx::Decode((ref->UserId() >> 19), edep, length, dEdep, dLength);
 286
 287
 288
 289
 290   // Calculate distance of previous reference to base of cluster
 291   UShort_t nT = TMath::Abs(t - fState.startStrip) + 1;
 292
 293   // Now check if we should flush to output
 294   Bool_t used = false;
 295
 296   // If this is a new detector/ring, then reset the other one
 297   // Check if we have a valid old detectorm ring, and sector
 298   if (fState.oldDetector >  0 &&
 299       fState.oldRing     != '\0' &&
 300       fState.oldSector   != 1024) {
 301     // New detector, new ring, or new sector
 302     if (d != fState.oldDetector   ||
 303         r != fState.oldRing       ||
 304         s != fState.oldSector) {
 305       if (fDebug) Info("Process", "New because new sector");
 306       used = true;
 307     }
 308     else if (nT > fMaxConsequtiveStrips) {
 309       if (fDebug) Info("Process", "New because too long: %d (%d,%d,%d)",
 310                        fState.nStrips, t, fState.oldStrip, fState.startStrip);
 311       used = true;
 312     }
 313   }
 314   if (used) {
 315     if (fDebug)
 316       Info("Process", "I=%p L=%p D=%d (was %d), R=%c (was %c), "
 317            "S=%2d (was %2d) t=%3d (was %3d) nT=%3d/%4d",
 318            ref, fState.longest,
 319            d, fState.oldDetector,
 320            r, fState.oldRing,
 321            s, fState.oldSector,
 322            t, fState.oldStrip,
 323            fState.nStrips, fMaxConsequtiveStrips);
 324     // Int_t nnT   = TMath::Abs(fState.oldStrip - fState.startStrip) + 1;
 325     StoreParticle(particle, mother, fState.longest);
 326     fState.Clear(false);
 327   }
 328
 329   // If base of cluster not set, set it here.
 330   if (fState.startStrip == 1024) fState.startStrip = t;
 331
 332   // Calculate distance of previous reference to base of cluster
 333   fState.nStrips = TMath::Abs(t - fState.startStrip) + 1;
 334
 335   // Count number of track refs in this sector
 336   fState.nRefs++;
 337
 338   fState.oldDetector = d;
 339   fState.oldRing     = r;
 340   fState.oldSector   = s;
 341   fState.oldStrip    = t;
 342   fState.de          += edep;
 343   fState.dx          += length;
 344
 345   // Debug output
 346   if (fDebug) {
 347     if (t == fState.startStrip)
 348       Info("Process", "New cluster starting at FMD%d%c[%2d,%3d]",
 349            d, r, s, t);
 350     else
 351       Info("Process", "Adding to cluster starting at FMD%d%c[%2d,%3d], "
 352            "length=%3d (now in %3d, previous %3d)",
 353            d, r, s, fState.startStrip, fState.nStrips, t, fState.oldStrip);
 354   }
 355
 356   // The longest passage is determined through the angle
 357   Double_t ang  = GetTrackRefTheta(ref);
 358   if (ang > fState.angle) {
 359     fState.longest = ref;
 360     fState.angle   = ang;
 361   }
 362
 363   return fState.longest;
 364 }
 365
 366 //____________________________________________________________________
 367 Double_t
 368 AliFMDMCTrackELoss::StoreParticle(AliMCParticle*       particle,
 369                                   const AliMCParticle* mother,
 370                                   AliTrackReference*   ref) const
 371 {
 372   Double_t w =
 373     AliBaseMCTrackDensity::StoreParticle(particle, mother, ref);
 374   if (w <= 0) return w;
 375
 376   // Get the detector coordinates
 377   UShort_t d, s, t;
 378   Char_t r;
 379   AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
 380
 381   Double_t de   = fState.de;
 382   Double_t dx   = fState.dx;
 383   Bool_t   prim = particle == mother;
 384   fAll(d,r,s,t) += de;
 385   if (prim) fPrimaries(d,r,s,t)   += de;
 386   else      fSecondaries(d,r,s,t) += de;
 387
 388   // Fill histograms
 389   fNr->Fill(fState.nRefs);
 390   fNt->Fill(fState.nStrips);
 391
 392   fState.count++;
 393
 394   if (de <= 0 || dx <= 0) return w;
 395
 396   TLorentzVector v;
 397   particle->Particle()->Momentum(v);
 398   if (v.E() <= 0) return w;
 399   if (v.Beta() < 0 || v.Beta() > 1) return w;
 400
 401   Double_t eta   = fEta(d,r,s,t);
 402   Double_t beta  = v.Beta();
 403   Double_t gamma = v.Gamma();
 404   Double_t dEdx  = de/dx;
 405
 406
 407   fBetaGammadEdx->Fill(beta*gamma, dEdx);
 408   fBetaGammaEta->Fill(eta, beta*gamma);
 409   fDEdxEta->Fill(eta, dEdx);
 410
 411   Int_t nHits = fHits->GetEntries();
 412   Hit* hit = new((*fHits)[nHits]) Hit;
 413   hit->fPrimary   = prim;
 414   hit->fEta       = eta;
 415   hit->fDe        = de;
 416   hit->fDx        = dx;
 417   hit->fDetId     = ref->UserId() & AliFMDStripIndex::kIdMask;
 418   hit->fBeta      = beta;
 419   hit->fGamma     = gamma;
 420   hit->fPdg       = particle->PdgCode();
 421
 422   return w;
 423 }
 424
 425 //____________________________________________________________________
 426 Bool_t
 427 AliFMDMCTrackELoss::Calculate(const AliESDFMD&  input,
 428                               const AliMCEvent& event,
 429                               Double_t          vz,
 430                               Double_t          cent)
 431 {
 432   //
 433   // Filter the input kinematics and track references, using
 434   // some of the ESD information
 435   //
 436   // Parameters:
 437   //    input   Input ESD event
 438   //    event   Input MC event
 439   //    vz      Vertex position
 440   //    output  Output ESD-like object
 441   //    primary Per-event histogram of primaries
 442   //
 443   // Return:
 444   //    True on succes, false otherwise
 445   //
 446   Clear();
 447   fEvent.fCent = cent;
 448   fEvent.fIpZ  = vz;
 449
 450   // Copy eta values to output
 451   for (UShort_t ed = 1; ed <= 3; ed++) {
 452     UShort_t nq = (ed == 1 ? 1 : 2);
 453     for (UShort_t eq = 0; eq < nq; eq++) {
 454       Char_t   er = (eq == 0 ? 'I' : 'O');
 455       UShort_t ns = (eq == 0 ?  20 :  40);
 456       UShort_t nt = (eq == 0 ? 512 : 256);
 457       for (UShort_t es = 0; es < ns; es++)
 458         for (UShort_t et = 0; et < nt; et++)
 459           fEta(ed, er, es, et) = input.Eta(ed, er, es, et);
 460     }
 461   }
 462
 463   Bool_t ret = ProcessTracks(event, vz, 0);
 464
 465   if (fTree) fTree->Fill();
 466
 467   return ret;
 468 }
 469
 470 #define PFV(N,VALUE)                                    \
 471   do {                                                  \
 472     AliForwardUtil::PrintName(N);                       \
 473     std::cout << (VALUE) << std::endl; } while(false)
 474 //____________________________________________________________________
 475 void
 476 AliFMDMCTrackELoss::Print(Option_t* option) const
 477 {
 478   AliBaseMCTrackDensity::Print(option);
 479   gROOT->IncreaseDirLevel();
 480   PFV("Max cluster size", fMaxConsequtiveStrips);
 481   gROOT->DecreaseDirLevel();
 482 }
 483
 484 //____________________________________________________________________
 485 //
 486 // EOF
 487 //