PWG2/FORWARD/analysis2/AliSPDMCTrackDensity.cxx

   1 #include "AliSPDMCTrackDensity.h"
   2 #include <AliMCEvent.h>
   3 #include <AliTrackReference.h>
   4 #include <AliStack.h>
   5 #include <TMath.h>
   6 #include <AliLog.h>
   7 #include <TH3D.h>
   8 #include <TH2D.h>
   9 #include <TH1D.h>
  10 #include <TList.h>
  11 #include <TROOT.h>
  12 #include <iostream>
  13
  14 //____________________________________________________________________
  15 AliSPDMCTrackDensity::AliSPDMCTrackDensity()
  16   : TNamed(),
  17     fUseOnlyPrimary(false),
  18     fMinR(3.5),
  19     fMaxR(4.5),
  20     fMinZ(-15), // -14.1),
  21     fMaxZ(+15), // +14.1),
  22     fRZ(0),
  23     fXYZ(0),
  24     fNRefs(0),
  25     fBinFlow(0)
  26 {
  27   // Default constructor
  28 }
  29
  30 //____________________________________________________________________
  31 AliSPDMCTrackDensity::AliSPDMCTrackDensity(const char*)
  32   : TNamed("spdMCTrackDensity","spdMCTrackDensity"),
  33     fUseOnlyPrimary(false),
  34     fMinR(3.5),
  35     fMaxR(4.5),
  36     fMinZ(-14.1),
  37     fMaxZ(+14.1),
  38     fRZ(0),
  39     fXYZ(0),
  40     fNRefs(0),
  41     fBinFlow(0)
  42 {
  43   // Normal constructor constructor
  44 }
  45
  46 //____________________________________________________________________
  47 AliSPDMCTrackDensity::AliSPDMCTrackDensity(const AliSPDMCTrackDensity& o)
  48   : TNamed(o),
  49     fUseOnlyPrimary(o.fUseOnlyPrimary),
  50     fMinR(o.fMinR),
  51     fMaxR(o.fMaxR),
  52     fMinZ(o.fMinZ),
  53     fMaxZ(o.fMaxZ),
  54     fRZ(o.fRZ),
  55     fXYZ(o.fXYZ),
  56     fNRefs(o.fNRefs),
  57     fBinFlow(o.fBinFlow)
  58 {
  59   // Normal constructor constructor
  60 }
  61
  62 //____________________________________________________________________
  63 AliSPDMCTrackDensity&
  64 AliSPDMCTrackDensity::operator=(const AliSPDMCTrackDensity& o)
  65 {
  66   // Assignment operator
  67   TNamed::operator=(o);
  68   fUseOnlyPrimary       = o.fUseOnlyPrimary;
  69   fMinR             = o.fMinR;
  70   fMaxR             = o.fMaxR;
  71   fMinZ             = o.fMinZ;
  72   fMaxZ             = o.fMaxZ;
  73   fRZ               = o.fRZ;
  74   fXYZ              = o.fXYZ;
  75   fNRefs            = o.fNRefs;
  76   fBinFlow          = o.fBinFlow;
  77   return *this;
  78 }
  79
  80 //____________________________________________________________________
  81 void
  82 AliSPDMCTrackDensity::DefineOutput(TList* l)
  83 {
  84   fRZ = new TH2D("rz", "(r,z) of used track references",
  85                  30, -15, 15, 50, 0, 5);
  86   fRZ->SetXTitle("z [cm]");
  87   fRZ->SetYTitle("r [cm]");
  88   fRZ->SetDirectory(0);
  89   l->Add(fRZ);
  90
  91   fXYZ = new TH3D("xyz", "(x,y,z) of used track references",
  92                   100, -5., 5., 100, -5., 5., 30, -15., 15.);
  93   fXYZ->SetXTitle("x [cm]");
  94   fXYZ->SetYTitle("y [cm]");
  95   fXYZ->SetZTitle("z [cm]");
  96   fXYZ->SetDirectory(0);
  97   l->Add(fXYZ);
  98
  99   fNRefs = new TH1D("nrefs", "Number of references used per track",
 100                     11, -.5, 10.5);
 101   fNRefs->SetXTitle("# of references per track");
 102   fNRefs->SetDirectory(0);
 103   l->Add(fNRefs);
 104
 105   fBinFlow = new TH2D("binFlow", "#eta and #varphi bin flow",
 106                       120, -3, 3, 40, -180, 180);
 107   fBinFlow->SetXTitle("#Delta#eta");
 108   fBinFlow->SetYTitle("#Delta#varphi");
 109   fBinFlow->SetDirectory(0);
 110   l->Add(fBinFlow);
 111 }
 112
 113 //____________________________________________________________________
 114 void
 115 AliSPDMCTrackDensity::StoreParticle(AliMCParticle* particle,
 116                                     const AliMCParticle* mother,
 117                                     Int_t          refNo,
 118                                     Double_t       vz,
 119                                     TH2D&     output) const
 120 {
 121   // Store a particle.
 122   if (refNo < 0) return;
 123
 124   AliTrackReference* ref = particle->GetTrackReference(refNo);
 125   if (!ref) return;
 126
 127   Double_t r = ref->R();
 128   Double_t x = ref->X();
 129   Double_t y = ref->Y();
 130   Double_t z = ref->Z();
 131
 132   Double_t zr = z-vz;
 133   Double_t th = TMath::ATan2(r,zr);
 134   if (th < 0) th += 2*TMath::Pi();
 135   Double_t et = -TMath::Log(TMath::Tan(th/2));
 136   Double_t ph = TMath::ATan2(y,x);
 137   if (ph < 0) ph += 2*TMath::Pi();
 138   output.Fill(et,ph);
 139
 140   const AliMCParticle* mp = (mother ? mother : particle);
 141   Double_t dEta = mp->Eta() - et;
 142   Double_t dPhi = (mp->Phi() - ph) * 180 / TMath::Pi();
 143   if (dPhi >  180) dPhi -= 360;
 144   if (dPhi < -180) dPhi += 360;
 145   fBinFlow->Fill(dEta, dPhi);
 146 }
 147
 148
 149 //____________________________________________________________________
 150 const AliMCParticle*
 151 AliSPDMCTrackDensity::GetMother(Int_t     iTr,
 152                                 const AliMCEvent& event) const
 153 {
 154   //
 155   // Track down primary mother
 156   //
 157   Int_t i  = iTr;
 158   do {
 159     const AliMCParticle* p = static_cast<AliMCParticle*>(event.GetTrack(i));
 160     if (const_cast<AliMCEvent&>(event).Stack()->IsPhysicalPrimary(i)) return p;
 161
 162     i = p->GetMother();
 163   } while (i > 0);
 164
 165   return 0;
 166 }
 167
 168 //____________________________________________________________________
 169 Bool_t
 170 AliSPDMCTrackDensity::Calculate(const AliMCEvent& event,
 171                                 Double_t          vz,
 172                                 TH2D&             output,
 173                                 TH2D*             primary)
 174 {
 175   //
 176   // Filter the input kinematics and track references, using
 177   // some of the ESD information
 178   //
 179   // Parameters:
 180   //    input   Input ESD event
 181   //    event   Input MC event
 182   //    vz      Vertex position
 183   //    output  Output ESD-like object
 184   //    primary Per-event histogram of primaries
 185   //
 186   // Return:
 187   //    True on succes, false otherwise
 188   //
 189
 190   AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
 191   Int_t nTracks   = stack->GetNtrack();
 192   Int_t nPrim     = stack->GetNtrack();
 193   for (Int_t iTr = 0; iTr < nTracks; iTr++) {
 194     AliMCParticle* particle =
 195       static_cast<AliMCParticle*>(event.GetTrack(iTr));
 196
 197     // Check the returned particle
 198     if (!particle) continue;
 199
 200     // Check if this charged and a primary
 201     Bool_t isCharged = particle->Charge() != 0;
 202     if (!isCharged) continue;
 203
 204     Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
 205
 206     // Fill 'dn/deta' histogram
 207     if (isPrimary && iTr < nPrim) {
 208       if (primary) {
 209         Double_t eta = particle->Eta();
 210         Double_t phi = particle->Phi();
 211         primary->Fill(eta, phi);
 212       }
 213     }
 214
 215     // Bail out if we're only processing primaries - perhaps we should
 216     // track back to the original primary?
 217     if (fUseOnlyPrimary && !isPrimary) continue;
 218
 219     Int_t nTrRef  = particle->GetNumberOfTrackReferences();
 220     Int_t nRef    = 0;
 221     for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
 222       AliTrackReference* ref = particle->GetTrackReference(iTrRef);
 223
 224       // Check existence
 225       if (!ref) continue;
 226
 227       // Check that we hit an FMD element
 228       if (ref->DetectorId() != AliTrackReference::kITS)
 229         continue;
 230
 231       // Get radius and z where the track reference was made
 232       Double_t r = ref->R();
 233       Double_t x = ref->X();
 234       Double_t y = ref->Y();
 235       Double_t z = ref->Z();
 236       if (r > fMaxR || r < fMinR) continue;
 237       if (z > fMaxZ || z < fMinZ) continue;
 238
 239       fRZ->Fill(z,r);
 240       fXYZ->Fill(x, y, z);
 241       nRef++;
 242       // Only fill first reference
 243       if (nRef == 1) {
 244         const AliMCParticle* mother = GetMother(iTr, event);
 245         StoreParticle(particle, mother, iTrRef, vz, output);
 246       }
 247     }
 248     fNRefs->Fill(nRef);
 249   }
 250   return kTRUE;
 251 }
 252 //____________________________________________________________________
 253 void
 254 AliSPDMCTrackDensity::Print(Option_t* /*option*/) const
 255 {
 256   char ind[gROOT->GetDirLevel()+1];
 257   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 258   ind[gROOT->GetDirLevel()] = '\0';
 259   std::cout << ind << ClassName() << ": " << GetName() << '\n'
 260             << std::boolalpha
 261             << ind << " Only primary tracks:    " << fUseOnlyPrimary << '\n'
 262             << ind << " R range:                [" << fMinR << ',' << fMaxR
 263             << "]\n"
 264             << ind << " Z range:                [" << fMinZ << ',' << fMaxZ
 265             << "]" << std::noboolalpha << std::endl;
 266
 267 }
 268
 269 //____________________________________________________________________
 270 //
 271 // EOF
 272 //