PWGLF/FORWARD/analysis2/AliFMDMCTrackDensity.cxx

   1 #include "AliFMDMCTrackDensity.h"
   2 #include "AliESDFMD.h"
   3 #include "AliTrackReference.h"
   4 #include <TMath.h>
   5 #include "AliFMDStripIndex.h"
   6 #include "AliLog.h"
   7 #include <TH2D.h>
   8 #include <TH1D.h>
   9 #include <TList.h>
  10 #include <TROOT.h>
  11 #include <iostream>
  12
  13 //____________________________________________________________________
  14 void
  15 AliFMDMCTrackDensity::State::Clear(Bool_t alsoCount)
  16 {
  17   angle       = 0;
  18   oldDetector = 0;
  19   oldRing     = '\0';
  20   oldSector   = 1024;
  21   oldStrip    = 1024;
  22   startStrip  = 1024;
  23   nRefs       = 0;
  24   nStrips     = 0;
  25   longest     = 0x0;
  26   if (alsoCount) count = 0;
  27 }
  28
  29 //____________________________________________________________________
  30 AliFMDMCTrackDensity::State&
  31 AliFMDMCTrackDensity::State::operator=(const State& o)
  32 {
  33   if (&o == this) return *this;
  34   angle          = o.angle;
  35   oldDetector    = o.oldDetector;
  36   oldRing        = o.oldRing;
  37   oldSector      = o.oldSector;
  38   oldStrip       = o.oldStrip;
  39   startStrip     = o.startStrip;
  40   nRefs          = o.nRefs;
  41   nStrips        = o.nStrips;
  42   count          = o.count;
  43   longest        = o.longest;
  44   return *this;
  45 }
  46
  47 //____________________________________________________________________
  48 AliFMDMCTrackDensity::AliFMDMCTrackDensity()
  49   : AliBaseMCTrackDensity(),
  50     fState(),
  51     fMaxConsequtiveStrips(3),
  52     fNr(0),
  53     fNt(0),
  54     fNc(0),
  55     fNcr(0),
  56     fOutput(0)
  57 {
  58   // Default constructor
  59 }
  60
  61 //____________________________________________________________________
  62 AliFMDMCTrackDensity::AliFMDMCTrackDensity(const char*)
  63   : AliBaseMCTrackDensity("fmdMCTrackDensity"),
  64     fState(),
  65     fMaxConsequtiveStrips(3),
  66     fNr(0),
  67     fNt(0),
  68     fNc(0),
  69     fNcr(0),
  70     fOutput(0)
  71 {
  72   // Normal constructor constructor
  73 }
  74
  75 //____________________________________________________________________
  76 AliFMDMCTrackDensity::AliFMDMCTrackDensity(const AliFMDMCTrackDensity& o)
  77   : AliBaseMCTrackDensity(o),
  78     fState(o.fState),
  79     fMaxConsequtiveStrips(o.fMaxConsequtiveStrips),
  80     fNr(o.fNr),
  81     fNt(o.fNt),
  82     fNc(o.fNc),
  83     fNcr(o.fNcr),
  84     fOutput(o.fOutput)
  85 {
  86   // Normal constructor constructor
  87 }
  88
  89 //____________________________________________________________________
  90 AliFMDMCTrackDensity&
  91 AliFMDMCTrackDensity::operator=(const AliFMDMCTrackDensity& o)
  92 {
  93   // Assignment operator
  94   if (&o == this) return *this;
  95   AliBaseMCTrackDensity::operator=(o);
  96   fMaxConsequtiveStrips = o.fMaxConsequtiveStrips;
  97   fNr                   = o.fNr;
  98   fNt                   = o.fNt;
  99   fNc                   = o.fNc;
 100   fNcr                  = o.fNcr;
 101   fState                = o.fState;
 102   fOutput               = o.fOutput;
 103
 104   return *this;
 105 }
 106
 107 //____________________________________________________________________
 108 void
 109 AliFMDMCTrackDensity::DefineOutput(TList* l)
 110 {
 111   AliBaseMCTrackDensity::DefineOutput(l);
 112   TList* ll = static_cast<TList*>(l->FindObject(GetTitle()));
 113   if (!ll) ll = l;
 114
 115   fNr = new TH1D("clusterRefs", "# track references per cluster",
 116                  21, -.5, 20.5);
 117   fNr->SetXTitle("# of track references/cluster");
 118   fNr->SetFillColor(kRed+1);
 119   fNr->SetFillStyle(3001);
 120   fNr->SetDirectory(0);
 121   ll->Add(fNr);
 122
 123   fNt = new TH1D("clusterSize", "cluster length in strips", 21, -.5, 20.5);
 124   fNt->SetXTitle("Cluster size [strips]");
 125   fNt->SetFillColor(kBlue+1);
 126   fNt->SetFillStyle(3001);
 127   fNt->SetDirectory(0);
 128   ll->Add(fNt);
 129
 130   fNc = new TH1D("nClusters", "# clusters per track", 21, -.5, 20.5);
 131   fNc->SetXTitle("# clusters");
 132   fNc->SetFillColor(kGreen+1);
 133   fNc->SetFillStyle(3001);
 134   fNc->SetDirectory(0);
 135   ll->Add(fNc);
 136
 137   fNcr = new TH2D("clusterVsRefs", "# clusters vs. # refs",
 138                   21, -.5, 20.5, 21, -.5, 20.5);
 139   fNcr->SetXTitle("# References");
 140   fNcr->SetYTitle("# Clusters");
 141   fNcr->SetOption("COLZ");
 142   fNcr->SetDirectory(0);
 143   ll->Add(fNcr);
 144
 145
 146 }
 147 //____________________________________________________________________
 148 Int_t
 149 AliFMDMCTrackDensity::GetDetectorId() const
 150 {
 151   return AliTrackReference::kFMD;
 152 }
 153
 154 //____________________________________________________________________
 155 void
 156 AliFMDMCTrackDensity::BeginTrackRefs()
 157 {
 158   fState.Clear(true);
 159 }
 160
 161 //____________________________________________________________________
 162 void
 163 AliFMDMCTrackDensity::EndTrackRefs(Int_t nRefs)
 164 {
 165   fNc->Fill(fState.count);
 166   fNcr->Fill(nRefs, fState.count);
 167   fState.Clear(true);
 168 }
 169
 170 //____________________________________________________________________
 171 AliTrackReference*
 172 AliFMDMCTrackDensity::ProcessRef(AliMCParticle*       particle,
 173                                  const AliMCParticle* mother,
 174                                  AliTrackReference*   ref)
 175 {
 176   // Get the detector coordinates
 177   UShort_t d, s, t;
 178   Char_t r;
 179   AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
 180
 181   // Calculate distance of previous reference to base of cluster
 182   UShort_t nT = TMath::Abs(t - fState.startStrip) + 1;
 183
 184   // Now check if we should flush to output
 185   Bool_t used = false;
 186
 187   // If this is a new detector/ring, then reset the other one
 188   // Check if we have a valid old detectorm ring, and sector
 189   if (fState.oldDetector >  0 &&
 190       fState.oldRing     != '\0' &&
 191       fState.oldSector   != 1024) {
 192     // New detector, new ring, or new sector
 193     if (d != fState.oldDetector   ||
 194         r != fState.oldRing       ||
 195         s != fState.oldSector) {
 196       if (fDebug) Info("Process", "New because new sector");
 197       used = true;
 198     }
 199     else if (nT > fMaxConsequtiveStrips) {
 200       if (fDebug) Info("Process", "New because too long: %d (%d,%d,%d)",
 201                        fState.nStrips, t, fState.oldStrip, fState.startStrip);
 202       used = true;
 203     }
 204   }
 205   if (used) {
 206     if (fDebug)
 207       Info("Process", "I=%p L=%p D=%d (was %d), R=%c (was %c), "
 208            "S=%2d (was %2d) t=%3d (was %3d) nT=%3d/%4d",
 209            ref, fState.longest,
 210            d, fState.oldDetector,
 211            r, fState.oldRing,
 212            s, fState.oldSector,
 213            t, fState.oldStrip,
 214            fState.nStrips, fMaxConsequtiveStrips);
 215     // Int_t nnT   = TMath::Abs(fState.oldStrip - fState.startStrip) + 1;
 216     StoreParticle(particle, mother, fState.longest);
 217     fState.Clear(false);
 218   }
 219
 220   // If base of cluster not set, set it here.
 221   if (fState.startStrip == 1024) fState.startStrip = t;
 222
 223   // Calculate distance of previous reference to base of cluster
 224   fState.nStrips = TMath::Abs(t - fState.startStrip) + 1;
 225
 226   // Count number of track refs in this sector
 227   fState.nRefs++;
 228
 229   fState.oldDetector = d;
 230   fState.oldRing     = r;
 231   fState.oldSector   = s;
 232   fState.oldStrip    = t;
 233
 234   // Debug output
 235   if (fDebug) {
 236     if (t == fState.startStrip)
 237       Info("Process", "New cluster starting at FMD%d%c[%2d,%3d]",
 238            d, r, s, t);
 239     else
 240       Info("Process", "Adding to cluster starting at FMD%d%c[%2d,%3d], "
 241            "length=%3d (now in %3d, previous %3d)",
 242            d, r, s, fState.startStrip, fState.nStrips, t, fState.oldStrip);
 243   }
 244
 245   // The longest passage is determined through the angle
 246   Double_t ang  = GetTrackRefTheta(ref);
 247   if (ang > fState.angle) {
 248     fState.longest = ref;
 249     fState.angle   = ang;
 250   }
 251
 252   return fState.longest;
 253 }
 254
 255 //____________________________________________________________________
 256 Double_t
 257 AliFMDMCTrackDensity::StoreParticle(AliMCParticle*       particle,
 258                                     const AliMCParticle* mother,
 259                                     AliTrackReference*   ref) const
 260 {
 261   Double_t w =
 262     AliBaseMCTrackDensity::StoreParticle(particle, mother, ref);
 263   if (w <= 0) return w;
 264
 265   // Get the detector coordinates
 266   UShort_t d, s, t;
 267   Char_t r;
 268   AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
 269
 270   // Check if we have value already
 271   Double_t old = fOutput->Multiplicity(d,r,s,t);
 272
 273   // If invalid, force it valid
 274   if (old == AliESDFMD::kInvalidMult) old = 0;
 275
 276   // Increment count
 277   fOutput->SetMultiplicity(d,r,s,t,old+w);
 278
 279   // Fill histograms
 280   fNr->Fill(fState.nRefs);
 281   fNt->Fill(fState.nStrips);
 282
 283   fState.count++;
 284
 285   return w;
 286 }
 287
 288 //____________________________________________________________________
 289 Bool_t
 290 AliFMDMCTrackDensity::Calculate(const AliESDFMD&  input,
 291                                 const AliMCEvent& event,
 292                                 Double_t          vz,
 293                                 AliESDFMD&        output,
 294                                 TH2D*             primary)
 295 {
 296   //
 297   // Filter the input kinematics and track references, using
 298   // some of the ESD information
 299   //
 300   // Parameters:
 301   //    input   Input ESD event
 302   //    event   Input MC event
 303   //    vz      Vertex position
 304   //    output  Output ESD-like object
 305   //    primary Per-event histogram of primaries
 306   //
 307   // Return:
 308   //    True on succes, false otherwise
 309   //
 310   output.Clear();
 311   fOutput = &output;
 312
 313   // Copy eta values to output
 314   for (UShort_t ed = 1; ed <= 3; ed++) {
 315     UShort_t nq = (ed == 1 ? 1 : 2);
 316     for (UShort_t eq = 0; eq < nq; eq++) {
 317       Char_t   er = (eq == 0 ? 'I' : 'O');
 318       UShort_t ns = (eq == 0 ?  20 :  40);
 319       UShort_t nt = (eq == 0 ? 512 : 256);
 320       for (UShort_t es = 0; es < ns; es++)
 321         for (UShort_t et = 0; et < nt; et++)
 322           output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et));
 323     }
 324   }
 325
 326   return ProcessTracks(event, vz, primary);
 327 }
 328 //____________________________________________________________________
 329 void
 330 AliFMDMCTrackDensity::Print(Option_t* option) const
 331 {
 332   AliBaseMCTrackDensity::Print(option);
 333   char ind[gROOT->GetDirLevel()+1];
 334   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 335   ind[gROOT->GetDirLevel()] = '\0';
 336   std::cout << ind << " Max cluster size:       " << fMaxConsequtiveStrips
 337             << std::endl;
 338
 339 }
 340
 341 //____________________________________________________________________
 342 //
 343 // EOF
 344 //