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