PWG2/FORWARD/analysis2/AliFMDMCTrackDensity.cxx

   1 #include "AliFMDMCTrackDensity.h"
   2 #include <AliESDFMD.h>
   3 #include <AliMCEvent.h>
   4 #include <AliTrackReference.h>
   5 #include <AliStack.h>
   6 #include <TMath.h>
   7 #include "AliFMDStripIndex.h"
   8 #include "AliFMDFloatMap.h"
   9 #include <AliLog.h>
  10 #include <TH2D.h>
  11 #include <TH1D.h>
  12 #include <TList.h>
  13 #include <TROOT.h>
  14 #include <iostream>
  15
  16 //____________________________________________________________________
  17 AliFMDMCTrackDensity::AliFMDMCTrackDensity()
  18   : TNamed(),
  19     fUseOnlyPrimary(false),
  20     fMaxConsequtiveStrips(3),
  21     fNr(0),
  22     fNt(0),
  23     fBinFlow(0),
  24     fEtaBinFlow(0),
  25     fPhiBinFlow(0),
  26     fDebug(false)
  27 {
  28   // Default constructor
  29 }
  30
  31 //____________________________________________________________________
  32 AliFMDMCTrackDensity::AliFMDMCTrackDensity(const char*)
  33   : TNamed("fmdMCTrackDensity","fmdMCTrackDensity"),
  34     fUseOnlyPrimary(false),
  35     fMaxConsequtiveStrips(3),
  36     fNr(0),
  37     fNt(0),
  38     fBinFlow(0),
  39     fEtaBinFlow(0),
  40     fPhiBinFlow(0),
  41     fDebug(false)
  42 {
  43   // Normal constructor constructor
  44 }
  45
  46 //____________________________________________________________________
  47 AliFMDMCTrackDensity::AliFMDMCTrackDensity(const AliFMDMCTrackDensity& o)
  48   : TNamed(o),
  49     fUseOnlyPrimary(o.fUseOnlyPrimary),
  50     fMaxConsequtiveStrips(o.fMaxConsequtiveStrips),
  51     fNr(o.fNr),
  52     fNt(o.fNt),
  53     fBinFlow(o.fBinFlow),
  54     fEtaBinFlow(o.fEtaBinFlow),
  55     fPhiBinFlow(o.fPhiBinFlow),
  56     fDebug(o.fDebug)
  57 {
  58   // Normal constructor constructor
  59 }
  60
  61 //____________________________________________________________________
  62 AliFMDMCTrackDensity&
  63 AliFMDMCTrackDensity::operator=(const AliFMDMCTrackDensity& o)
  64 {
  65   // Assignment operator
  66   TNamed::operator=(o);
  67   fUseOnlyPrimary       = o.fUseOnlyPrimary;
  68   fMaxConsequtiveStrips = o.fMaxConsequtiveStrips;
  69   fNr                   = o.fNr;
  70   fNt                   = o.fNt;
  71   fBinFlow              = o.fBinFlow;
  72   fEtaBinFlow           = o.fEtaBinFlow;
  73   fPhiBinFlow           = o.fPhiBinFlow;
  74   fDebug                = o.fDebug;
  75   return *this;
  76 }
  77
  78 //____________________________________________________________________
  79 void
  80 AliFMDMCTrackDensity::DefineOutput(TList* l)
  81 {
  82   fNr = new TH1D("clusterRefs", "# track references per cluster",
  83                  21, -.5, 20.5);
  84   fNr->SetXTitle("# of track references/cluster");
  85   fNr->SetDirectory(0);
  86   l->Add(fNr);
  87
  88   fNt = new TH1D("clusterSize", "cluster length in strips", 21, -.5, 20.5);
  89   fNt->SetXTitle("Cluster size [strips]");
  90   fNt->SetDirectory(0);
  91   l->Add(fNt);
  92
  93   fBinFlow = new TH2D("binFlow", "#eta and #varphi bin flow",
  94                       200, -5, 5, 40, -180, 180);
  95   fBinFlow->SetXTitle("#Delta#eta");
  96   fBinFlow->SetYTitle("#Delta#varphi");
  97   fBinFlow->SetDirectory(0);
  98   l->Add(fBinFlow);
  99
 100   fEtaBinFlow = new TH2D("binFlowEta", "#eta bin flow vs #eta",
 101                          200, -4, 6, 200, -5, 5);
 102   fEtaBinFlow->SetXTitle("#eta of strip");
 103   fEtaBinFlow->SetYTitle("#Delta#eta");
 104   fEtaBinFlow->SetDirectory(0);
 105   l->Add(fEtaBinFlow);
 106
 107   fPhiBinFlow = new TH2D("binFlowPhi", "#phi bin flow vs #phi",
 108                          40, 0, 360, 40, -180, 180);
 109   fPhiBinFlow->SetXTitle("#varphi of strip");
 110   fPhiBinFlow->SetYTitle("#Delta#varphi");
 111   fPhiBinFlow->SetDirectory(0);
 112   l->Add(fPhiBinFlow);
 113 }
 114
 115 //____________________________________________________________________
 116 void
 117 AliFMDMCTrackDensity::StoreParticle(AliMCParticle*       particle,
 118                                     const AliMCParticle* mother,
 119                                     Int_t                longest,
 120                                     Double_t             vz,
 121                                     UShort_t             nC,
 122                                     UShort_t             nT,
 123                                     AliESDFMD&           output) const
 124 {
 125   // Store a particle.
 126   if (longest < 0) return;
 127
 128   AliTrackReference* ref = particle->GetTrackReference(longest);
 129   if (!ref) return;
 130
 131   // Get the detector coordinates
 132   UShort_t d, s, t;
 133   Char_t r;
 134   AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
 135
 136   // Check if we have value already
 137   Double_t old = output.Multiplicity(d,r,s,t);
 138
 139   // If invalid, force it valid
 140   if (old == AliESDFMD::kInvalidMult) old = 0;
 141
 142   // Increment count
 143   output.SetMultiplicity(d,r,s,t,old+1);
 144
 145   // Get track-reference stuff
 146   Double_t x      = ref->X();
 147   Double_t y      = ref->Y();
 148   Double_t z      = ref->Z()-vz;
 149   Double_t rr     = TMath::Sqrt(x*x+y*y);
 150   Double_t thetaR = TMath::ATan2(rr,z);
 151   Double_t phiR   = TMath::ATan2(y,x);
 152   Double_t etaR   = -TMath::Log(TMath::Tan(thetaR/2));
 153
 154   // Correct angle and convert to degrees
 155   if (thetaR < 0) thetaR += 2*TMath::Pi();
 156   thetaR *= 180. / TMath::Pi();
 157   if (phiR < 0) phiR += 2*TMath::Pi();
 158   phiR *= 180. / TMath::Pi();
 159
 160   // Fill histograms
 161   fNr->Fill(nC);
 162   fNt->Fill(nT);
 163
 164   const AliMCParticle* mp = (mother ? mother : particle);
 165   Double_t dEta = mp->Eta() - etaR;
 166   Double_t dPhi = mp->Phi() * 180 / TMath::Pi() - phiR;
 167   if (dPhi >  180) dPhi -= 360;
 168   if (dPhi < -180) dPhi += 360;
 169   fBinFlow->Fill(dEta, dPhi);
 170   fEtaBinFlow->Fill(etaR, dEta);
 171   fPhiBinFlow->Fill(phiR, dPhi);
 172
 173   if (fDebug)
 174     Info("StoreParticle", "Store @ FMD%d%c[%2d,%3d] nStrips=%3d, "
 175          "dEta=%7.4f, dPhi=%d",
 176          d, r, s, t, nT, dEta, int(dPhi+.5));
 177 }
 178
 179 //____________________________________________________________________
 180 Double_t
 181 AliFMDMCTrackDensity::GetTrackRefTheta(const AliTrackReference* ref,
 182                                        Double_t vz) const
 183 {
 184   // Get the incidient angle of the track reference.
 185   Double_t x    = ref->X();
 186   Double_t y    = ref->Y();
 187   Double_t z    = ref->Z()-vz;
 188   Double_t rr   = TMath::Sqrt(x*x+y*y);
 189   Double_t theta= TMath::ATan2(rr,z);
 190   Double_t ang  = TMath::Abs(TMath::Pi()-theta);
 191   return ang;
 192 }
 193
 194 //____________________________________________________________________
 195 const AliMCParticle*
 196 AliFMDMCTrackDensity::GetMother(Int_t     iTr,
 197                                 const AliMCEvent& event) const
 198 {
 199   //
 200   // Track down primary mother
 201   //
 202   Int_t i  = iTr;
 203   do {
 204     const AliMCParticle* p = static_cast<AliMCParticle*>(event.GetTrack(i));
 205     if (const_cast<AliMCEvent&>(event).Stack()->IsPhysicalPrimary(i)) return p;
 206
 207     i = p->GetMother();
 208   } while (i > 0);
 209
 210   return 0;
 211 }
 212
 213 //____________________________________________________________________
 214 Bool_t
 215 AliFMDMCTrackDensity::Calculate(const AliESDFMD&  input,
 216                                 const AliMCEvent& event,
 217                                 Double_t          vz,
 218                                 AliESDFMD&        output,
 219                                 TH2D*             primary)
 220 {
 221   //
 222   // Filter the input kinematics and track references, using
 223   // some of the ESD information
 224   //
 225   // Parameters:
 226   //    input   Input ESD event
 227   //    event   Input MC event
 228   //    vz      Vertex position
 229   //    output  Output ESD-like object
 230   //    primary Per-event histogram of primaries
 231   //
 232   // Return:
 233   //    True on succes, false otherwise
 234   //
 235   output.Clear();
 236
 237   // Copy eta values to output
 238   for (UShort_t ed = 1; ed <= 3; ed++) {
 239     UShort_t nq = (ed == 1 ? 1 : 2);
 240     for (UShort_t eq = 0; eq < nq; eq++) {
 241       Char_t   er = (eq == 0 ? 'I' : 'O');
 242       UShort_t ns = (eq == 0 ?  20 :  40);
 243       UShort_t nt = (eq == 0 ? 512 : 256);
 244       for (UShort_t es = 0; es < ns; es++)
 245         for (UShort_t et = 0; et < nt; et++)
 246           output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et));
 247     }
 248   }
 249   AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
 250   Int_t nTracks   = stack->GetNtrack();//event.GetNumberOfTracks();
 251   Int_t nPrim     = stack->GetNprimary();//event.GetNumberOfPrimary();
 252   for (Int_t iTr = 0; iTr < nTracks; iTr++) {
 253     AliMCParticle* particle =
 254       static_cast<AliMCParticle*>(event.GetTrack(iTr));
 255
 256     // Check the returned particle
 257     if (!particle) continue;
 258
 259     // Check if this charged and a primary
 260     Bool_t isCharged = particle->Charge() != 0;
 261     if (!isCharged) continue;
 262     Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
 263
 264     // Pseudo rapidity and azimuthal angle
 265     Double_t eta = particle->Eta();
 266     Double_t phi = particle->Phi();
 267
 268     // Fill 'dn/deta' histogram
 269     if (isPrimary && iTr < nPrim) {
 270       if (primary) primary->Fill(eta, phi);
 271     }
 272
 273     // Bail out if we're only processing primaries - perhaps we should
 274     // track back to the original primary?
 275     if (fUseOnlyPrimary && !isPrimary) continue;
 276
 277     Int_t    nTrRef   = particle->GetNumberOfTrackReferences();
 278     Int_t    longest  = -1;
 279     Double_t angle    = 0;
 280     UShort_t oD       = 0, oS = 1024, oT = 1024, ooT=1024;
 281     Char_t   oR       = '\0';
 282     UShort_t nC       = 0;
 283     UShort_t nT       = 0;
 284     // Double_t oTheta= 0;
 285     for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
 286       AliTrackReference* ref = particle->GetTrackReference(iTrRef);
 287
 288       // Check existence
 289       if (!ref) continue;
 290
 291       // Check that we hit an FMD element
 292       if (ref->DetectorId() != AliTrackReference::kFMD)
 293         continue;
 294
 295       // Get the detector coordinates
 296       UShort_t d, s, t;
 297       Char_t r;
 298       AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
 299
 300       // Calculate length of last and second to last strips.
 301
 302       // IF base of cluster not set, set it here.
 303       if (ooT == 1024) ooT = t;
 304
 305       // Calculate distance of previous reference to base of cluster
 306       nT = TMath::Abs(t - ooT) + 1;
 307
 308       // Count number of track refs in this sector
 309       nC++;
 310
 311       Bool_t used = false;
 312       // If this is a new detector/ring, then reset the other one
 313       // Check if we have a valid old detectorm ring, and sector
 314       if (oD >  0 && oR != '\0' && oS != 1024) {
 315         // New detector, new ring, or new sector
 316         if (d != oD   || r != oR   || s != oS) {
 317           if (fDebug) Info("Process", "New because new sector");
 318           used = true;
 319         }
 320       }
 321       if (nT > fMaxConsequtiveStrips) {
 322         if (fDebug) Info("Process", "New because too long: %d (%d,%d,%d)",
 323                          nT, t, oT, ooT);
 324         used = true;
 325       }
 326       if (used) {
 327         if (fDebug)
 328           Info("Process", "I=%3d L=%3d D=%d (was %d), R=%c (was %c), "
 329                "S=%2d (was %2d) t=%3d (was %3d) nT=%3d/%4d",
 330                iTr, longest, d, oD, r, oR, s, oS, t, oT, nT,
 331                fMaxConsequtiveStrips);
 332         Int_t nnT   = TMath::Abs(oT - ooT) + 1;
 333         const AliMCParticle* mother = GetMother(iTr, event);
 334         StoreParticle(particle, mother, longest, vz, nC, nnT, output);
 335         longest = -1;
 336         angle   = 0;
 337         nC  = 1;    // Reset track-ref counter - we have this->1
 338         nT  = 0;    // Reset cluster size
 339         oD  = 0;    // Reset detector
 340         oR  = '\0'; // Reset ring
 341         oS  = 1024; // Reset sector
 342         oT  = 1024; // Reset old strip
 343         ooT = t;    // Reset base
 344       }
 345       else if (fDebug) {
 346         if (ooT == t)
 347           Info("Process", "New cluster starting at FMD%d%c[%2d,%3d]",
 348                d, r, s, t);
 349         else
 350           Info("Process", "Adding to cluster starting at FMD%d%c[%2d,%3d], "
 351                "length=%3d (now in %3d, previous %3d)",
 352                d, r, s, ooT, nT, t, oT);
 353       }
 354
 355       oD = d;
 356       oR = r;
 357       oS = s;
 358       oT = t;
 359       nT = TMath::Abs(t-ooT)+1;
 360
 361       // The longest passage is determined through the angle
 362       Double_t ang  = GetTrackRefTheta(ref, vz);
 363       if (ang > angle) {
 364         longest = iTrRef;
 365         angle   = ang;
 366       }
 367       // oTheta = ang;
 368     } // Loop over track references
 369     if (longest < 0) continue; // Nothing found
 370
 371     // Get the reference corresponding to the longest path through the detector
 372     if (fDebug)
 373       Info("Process", "I=%3d L=%3d nT=%3d (out of %3d)",
 374            iTr, longest, nT, fMaxConsequtiveStrips);
 375     const AliMCParticle* mother = GetMother(iTr, event);
 376     StoreParticle(particle, mother, longest, vz, nC, nT, output);
 377   } // Loop over tracks
 378   return kTRUE;
 379 }
 380 //____________________________________________________________________
 381 void
 382 AliFMDMCTrackDensity::Print(Option_t* /*option*/) const
 383 {
 384   char ind[gROOT->GetDirLevel()+1];
 385   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 386   ind[gROOT->GetDirLevel()] = '\0';
 387   std::cout << ind << ClassName() << ": " << GetName() << '\n'
 388             << std::boolalpha
 389             << ind << " Only primary tracks:    " << fUseOnlyPrimary << '\n'
 390             << ind << " Max cluster size:       " << fMaxConsequtiveStrips
 391             << std::noboolalpha << std::endl;
 392
 393 }
 394
 395 //____________________________________________________________________
 396 //
 397 // EOF
 398 //