PWGLF/FORWARD/analysis2/AliBaseMCTrackDensity.cxx

   1 #include "AliBaseMCTrackDensity.h"
   2 #include <AliMCEvent.h>
   3 #include <AliTrackReference.h>
   4 #include <AliStack.h>
   5 #include <TMath.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 #include "AliCollisionGeometry.h"
  13 #include "AliGenEventHeader.h"
  14 #include "AliForwardUtil.h"
  15 #include <TF1.h>
  16 #include <TGraph.h>
  17
  18 //____________________________________________________________________
  19 AliBaseMCTrackDensity::AliBaseMCTrackDensity()
  20   : TNamed(),
  21     fUseOnlyPrimary(false),
  22     fUseFlowWeights(false),
  23     fBinFlow(0),
  24     fEtaBinFlow(0),
  25     fPhiBinFlow(0),
  26     fNRefs(0),
  27     fWeights(),
  28     fVz(0),
  29     fB(0),
  30     fPhiR(0),
  31     fDebug(false)
  32 {
  33   // Default constructor
  34   DGUARD(fDebug, 3,"Default CTOR of AliBasMCTrackDensity");
  35 }
  36
  37 //____________________________________________________________________
  38 AliBaseMCTrackDensity::AliBaseMCTrackDensity(const char* name)
  39   : TNamed(name,"mcTrackDensity"),
  40     fUseOnlyPrimary(false),
  41     fUseFlowWeights(false),
  42     fBinFlow(0),
  43     fEtaBinFlow(0),
  44     fPhiBinFlow(0),
  45     fNRefs(0),
  46     fWeights(),
  47     fVz(0),
  48     fB(0),
  49     fPhiR(0),
  50     fDebug(false)
  51 {
  52   // Normal constructor constructor
  53   DGUARD(fDebug, 3,"Named CTOR of AliBasMCTrackDensity: %s", name);
  54 }
  55
  56 //____________________________________________________________________
  57 AliBaseMCTrackDensity::AliBaseMCTrackDensity(const AliBaseMCTrackDensity& o)
  58   : TNamed(o),
  59     fUseOnlyPrimary(o.fUseOnlyPrimary),
  60     fUseFlowWeights(o.fUseFlowWeights),
  61     fBinFlow(o.fBinFlow),
  62     fEtaBinFlow(o.fEtaBinFlow),
  63     fPhiBinFlow(o.fPhiBinFlow),
  64     fNRefs(o.fNRefs),
  65     fWeights(o.fWeights),
  66     fVz(o.fVz),
  67     fB(o.fB),
  68     fPhiR(o.fPhiR),
  69     fDebug(o.fDebug)
  70 {
  71   // Normal constructor constructor
  72   DGUARD(fDebug, 3,"Copy CTOR of AliBasMCTrackDensity");
  73 }
  74
  75 //____________________________________________________________________
  76 AliBaseMCTrackDensity&
  77 AliBaseMCTrackDensity::operator=(const AliBaseMCTrackDensity& o)
  78 {
  79   // Assignment operator
  80   DGUARD(fDebug,3,"MC track density assignmetn");
  81   if (&o == this) return *this;
  82   TNamed::operator=(o);
  83   fUseOnlyPrimary       = o.fUseOnlyPrimary;
  84   fBinFlow              = o.fBinFlow;
  85   fEtaBinFlow           = o.fEtaBinFlow;
  86   fPhiBinFlow           = o.fPhiBinFlow;
  87   fNRefs                = o.fNRefs;
  88   fDebug                = o.fDebug;
  89   fUseFlowWeights       = o.fUseFlowWeights;
  90   fWeights              = o.fWeights;
  91   fVz                   = o.fVz;
  92   fB                    = o.fB;
  93   fPhiR                 = o.fPhiR;
  94   return *this;
  95 }
  96
  97 //____________________________________________________________________
  98 void
  99 AliBaseMCTrackDensity::CreateOutputObjects(TList* l)
 100 {
 101   DGUARD(fDebug,1,"MC track defines output");
 102   TList* ll = new TList;
 103   ll->SetName(GetTitle());
 104   ll->SetOwner();
 105   l->Add(ll);
 106
 107   fBinFlow = new TH2D("binFlow", "#eta and #varphi bin flow",
 108                       200, -5, 5, 40, -180, 180);
 109   fBinFlow->SetXTitle("#Delta#eta");
 110   fBinFlow->SetYTitle("#Delta#varphi");
 111   fBinFlow->SetOption("colz");
 112   fBinFlow->SetDirectory(0);
 113   ll->Add(fBinFlow);
 114
 115   fEtaBinFlow = new TH2D("binFlowEta", "#eta bin flow vs #eta",
 116                          200, -4, 6, 200, -5, 5);
 117   fEtaBinFlow->SetXTitle("#eta");
 118   fEtaBinFlow->SetYTitle("#Delta#eta");
 119   fEtaBinFlow->SetOption("colz");
 120   fEtaBinFlow->SetDirectory(0);
 121   ll->Add(fEtaBinFlow);
 122
 123   fPhiBinFlow = new TH2D("binFlowPhi", "#phi bin flow vs #phi",
 124                          40, 0, 360, 40, -180, 180);
 125   fPhiBinFlow->SetXTitle("#varphi");
 126   fPhiBinFlow->SetYTitle("#Delta#varphi");
 127   fPhiBinFlow->SetOption("colz");
 128   fPhiBinFlow->SetDirectory(0);
 129   ll->Add(fPhiBinFlow);
 130
 131   fNRefs = new TH1D("nRefs", "# references per track", 21, -.5, 20.5);
 132   fNRefs->SetXTitle("# references");
 133   fNRefs->SetFillColor(kMagenta+1);
 134   fNRefs->SetFillStyle(3001);
 135   fNRefs->SetDirectory(0);
 136   ll->Add(fNRefs);
 137
 138   fWeights.Init(ll);
 139 }
 140
 141
 142 //____________________________________________________________________
 143 Double_t
 144 AliBaseMCTrackDensity::StoreParticle(AliMCParticle*       particle,
 145                                      const AliMCParticle* mother,
 146                                      AliTrackReference*   ref) const
 147 {
 148   DGUARD(fDebug,3,"MC track density store particle");
 149   // Store a particle.
 150   if (!ref) return 0;
 151
 152   Double_t weight = 1;
 153   if (fUseFlowWeights) {
 154     Double_t phi = (mother ? mother->Phi() : particle->Phi());
 155     Double_t eta = (mother ? mother->Eta() : particle->Eta());
 156     Double_t pt  = (mother ? mother->Pt() : particle->Pt());
 157     Int_t    id  = (mother ? mother->PdgCode() : 2212);
 158     weight       = CalculateWeight(eta, pt, phi, id);
 159   }
 160
 161   // Get track-reference stuff
 162   Double_t x      = ref->X();
 163   Double_t y      = ref->Y();
 164   Double_t z      = ref->Z()-fVz;
 165   Double_t rr     = TMath::Sqrt(x*x+y*y);
 166   Double_t thetaR = TMath::ATan2(rr,z);
 167   Double_t phiR   = TMath::ATan2(y,x);
 168   Double_t etaR   = -TMath::Log(TMath::Tan(thetaR/2));
 169
 170   // Correct angle and convert to degrees
 171   if (thetaR < 0) thetaR += 2*TMath::Pi();
 172   thetaR *= 180. / TMath::Pi();
 173   if (phiR < 0) phiR += 2*TMath::Pi();
 174   phiR *= 180. / TMath::Pi();
 175
 176   const AliMCParticle* mp = (mother ? mother : particle);
 177   Double_t dEta = mp->Eta() - etaR;
 178   Double_t dPhi = mp->Phi() * 180 / TMath::Pi() - phiR;
 179   if (dPhi >  180) dPhi -= 360;
 180   if (dPhi < -180) dPhi += 360;
 181   fBinFlow->Fill(dEta, dPhi);
 182   fEtaBinFlow->Fill(etaR, dEta);
 183   fPhiBinFlow->Fill(phiR, dPhi);
 184
 185   return weight;
 186 }
 187
 188 //____________________________________________________________________
 189 Double_t
 190 AliBaseMCTrackDensity::GetTrackRefTheta(const AliTrackReference* ref) const
 191 {
 192   // Get the incidient angle of the track reference.
 193   Double_t x    = ref->X();
 194   Double_t y    = ref->Y();
 195   Double_t z    = ref->Z()-fVz;
 196   Double_t rr   = TMath::Sqrt(x*x+y*y);
 197   Double_t theta= TMath::ATan2(rr,z);
 198   Double_t ang  = TMath::Abs(TMath::Pi()-theta);
 199   return ang;
 200 }
 201
 202 //____________________________________________________________________
 203 const AliMCParticle*
 204 AliBaseMCTrackDensity::GetMother(Int_t     iTr,
 205                                 const AliMCEvent& event) const
 206 {
 207   //
 208   // Track down primary mother
 209   //
 210   Int_t i  = iTr;
 211   do {
 212     const AliMCParticle* p = static_cast<AliMCParticle*>(event.GetTrack(i));
 213     if (const_cast<AliMCEvent&>(event).Stack()->IsPhysicalPrimary(i)) return p;
 214
 215     i = p->GetMother();
 216   } while (i > 0);
 217
 218   return 0;
 219 }
 220
 221 //____________________________________________________________________
 222 Bool_t
 223 AliBaseMCTrackDensity::GetCollisionParameters(const AliMCEvent& event)
 224 {
 225   DGUARD(fDebug,3,"MC track density get collision parameters");
 226   AliCollisionGeometry* hd =
 227     dynamic_cast<AliCollisionGeometry*>(event.GenEventHeader());
 228   fPhiR = (hd ? hd->ReactionPlaneAngle() : 0.);
 229   fB    = (hd ? hd->ImpactParameter() : -1 );
 230   return hd != 0;
 231 }
 232
 233 //____________________________________________________________________
 234 Bool_t
 235 AliBaseMCTrackDensity::ProcessTrack(AliMCParticle* particle,
 236                                     const AliMCParticle* mother)
 237 {
 238   // Check the returned particle
 239   DGUARD(fDebug,3,"MC track density Process a track");
 240   if (!particle) return false;
 241
 242   Int_t              nTrRef = particle->GetNumberOfTrackReferences();
 243   AliTrackReference* store  = 0;
 244
 245   BeginTrackRefs();
 246
 247   // Double_t oTheta= 0;
 248   Int_t nRefs = 0;
 249   for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
 250     AliTrackReference* ref = particle->GetTrackReference(iTrRef);
 251
 252     // Check existence
 253     if (!ref) continue;
 254
 255     // Check that we hit an Base element
 256     if (ref->DetectorId() != GetDetectorId()) continue;
 257     if (!CheckTrackRef(ref)) continue;
 258
 259     nRefs++;
 260
 261     AliTrackReference* test = ProcessRef(particle, mother, ref);
 262     if (test) store = test;
 263
 264   } // Loop over track references
 265   if (!store) return true; // Nothing found
 266
 267   StoreParticle(particle, mother, store);
 268
 269   fNRefs->Fill(nRefs);
 270
 271   EndTrackRefs(nRefs);
 272
 273   return true;
 274 }
 275
 276
 277 //____________________________________________________________________
 278 Bool_t
 279 AliBaseMCTrackDensity::ProcessTracks(const AliMCEvent& event,
 280                                      Double_t          vz,
 281                                      TH2D*             primary)
 282 {
 283   //
 284   // Filter the input kinematics and track references, using
 285   // some of the ESD information
 286   //
 287   // Parameters:
 288   //    input   Input ESD event
 289   //    event   Input MC event
 290   //    vz      Vertex position
 291   //    output  Output ESD-like object
 292   //    primary Per-event histogram of primaries
 293   //
 294   // Return:
 295   //    True on succes, false otherwise
 296   //
 297   DGUARD(fDebug,3,"MC track density Process a tracks");
 298   fVz = vz;
 299   GetCollisionParameters(event);
 300
 301   AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
 302   Int_t nTracks   = stack->GetNtrack();//event.GetNumberOfTracks();
 303   Int_t nPrim     = stack->GetNprimary();//event.GetNumberOfPrimary();
 304   for (Int_t iTr = 0; iTr < nTracks; iTr++) {
 305     AliMCParticle* particle =
 306       static_cast<AliMCParticle*>(event.GetTrack(iTr));
 307
 308     // Check if this charged and a primary
 309     if (!particle->Charge() != 0) continue;
 310
 311     Bool_t isPrimary = stack->IsPhysicalPrimary(iTr) && iTr < nPrim;
 312
 313     // Fill 'dn/deta' histogram
 314     if (isPrimary && primary)
 315       primary->Fill(particle->Eta(), particle->Phi());
 316
 317     // Bail out if we're only processing primaries - perhaps we should
 318     // track back to the original primary?
 319     if (fUseOnlyPrimary && !isPrimary) continue;
 320
 321     const AliMCParticle* mother = GetMother(iTr, event);
 322     ProcessTrack(particle, mother);
 323
 324   } // Loop over tracks
 325   return kTRUE;
 326 }
 327
 328
 329 //____________________________________________________________________
 330 Double_t
 331 AliBaseMCTrackDensity::CalculateWeight(Double_t eta, Double_t pt,
 332                                        Double_t phi, Int_t id) const
 333 {
 334   return fWeights.CalcWeight(eta, pt, phi, id, fPhiR, fB);
 335 }
 336 //____________________________________________________________________
 337 void
 338 AliBaseMCTrackDensity::Print(Option_t* /*option*/) const
 339 {
 340   char ind[gROOT->GetDirLevel()+1];
 341   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 342   ind[gROOT->GetDirLevel()] = '\0';
 343   std::cout << ind << ClassName() << ": " << GetName() << '\n'
 344             << std::boolalpha
 345             << ind << " Only primary tracks:    " << fUseOnlyPrimary << '\n'
 346             << ind << " Use flow after burner:  " << fUseFlowWeights
 347             << std::noboolalpha << std::endl;
 348
 349 }
 350
 351 //____________________________________________________________________
 352 //
 353 // EOF
 354 //