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