Use AliFMDMCTrackDensity for hits
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 26 Apr 2011 20:26:38 +0000 (20:26 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 26 Apr 2011 20:26:38 +0000 (20:26 +0000)
PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.cxx
PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.h

index cc2e866..a880a48 100644 (file)
@@ -43,15 +43,13 @@ ClassImp(AliFMDMCSharingFilter)
 //____________________________________________________________________
 AliFMDMCSharingFilter::AliFMDMCSharingFilter(const char* title)
   : AliFMDSharingFilter(title), 
+    fTrackDensity(title),
     fFMD1i(0),
     fFMD2i(0),
     fFMD2o(0),
     fFMD3i(0),
     fFMD3o(0),
-    fSumEta(0),
-    fOperComp(0),
-    fThetaVsNr(0), 
-    fOnlyPrimary(false)
+    fOperComp(0)
 {
   // 
   // Constructor 
@@ -79,15 +77,6 @@ AliFMDMCSharingFilter::AliFMDMCSharingFilter(const char* title)
   fFMD2o->SetDirectory(0);
   fFMD3i->SetDirectory(0);
   fFMD3o->SetDirectory(0);
-  fSumEta = new TH1D("mcSumEta", "MC INEL Truth", 200, -4, 6);
-  fSumEta->SetXTitle("#eta");
-  fSumEta->SetYTitle("dN_{ch}/d#eta");
-  fSumEta->SetDirectory(0);
-  fSumEta->Sumw2();
-  fSumEta->SetMarkerColor(kOrange+2);
-  fSumEta->SetMarkerStyle(22);
-  fSumEta->SetFillColor(0);
-  fSumEta->SetFillStyle(0);
 
   fOper     = new AliFMDFloatMap(0,0,0,0);
   fOperComp = new TH2I("operComp", "Operation vs # track refs", 
@@ -101,26 +90,18 @@ AliFMDMCSharingFilter::AliFMDMCSharingFilter(const char* title)
   fOperComp->GetXaxis()->SetBinLabel(kMergedWithOther, "Merged w/other");
   fOperComp->GetXaxis()->SetBinLabel(kMergedInto,      "Merged into");
   fOperComp->SetDirectory(0);
-  
-  fThetaVsNr = new TH2D("thetaVsNr", "#theta of track vs # track references",
-                       360, 0, 360, 20, -.5, 19.5);
-  fThetaVsNr->SetXTitle("#theta [degrees]");
-  fThetaVsNr->SetYTitle("# of track references");
-  fThetaVsNr->SetDirectory(0);
 }
 
 //____________________________________________________________________
 AliFMDMCSharingFilter::AliFMDMCSharingFilter(const AliFMDMCSharingFilter& o)
   : AliFMDSharingFilter(o), 
+    fTrackDensity(o.fTrackDensity),
     fFMD1i(o.fFMD1i),
     fFMD2i(o.fFMD2i),
     fFMD2o(o.fFMD2o),
     fFMD3i(o.fFMD3i),
     fFMD3o(o.fFMD3o),
-    fSumEta(o.fSumEta),
-    fOperComp(o.fOperComp),
-    fThetaVsNr(o.fThetaVsNr),
-    fOnlyPrimary(o.fOnlyPrimary)
+    fOperComp(o.fOperComp)
 {
   // 
   // Copy constructor 
@@ -136,12 +117,6 @@ AliFMDMCSharingFilter::~AliFMDMCSharingFilter()
   // 
   // Destructor
   //
-  if (fFMD1i)  delete fFMD1i;
-  if (fFMD2i)  delete fFMD2i;
-  if (fFMD2o)  delete fFMD2o;
-  if (fFMD3i)  delete fFMD3i;
-  if (fFMD3o)  delete fFMD3o;
-  if (fSumEta) delete fSumEta;
 }
 
 //____________________________________________________________________
@@ -158,39 +133,10 @@ AliFMDMCSharingFilter::operator=(const AliFMDMCSharingFilter& o)
   //    Reference to this 
   //
   AliFMDSharingFilter::operator=(o);
-  fOnlyPrimary = o.fOnlyPrimary;
+  fTrackDensity = o.fTrackDensity;
   return *this;
 }
 
-//____________________________________________________________________
-void
-AliFMDMCSharingFilter::StoreParticle(UShort_t   d, 
-                                    Char_t     r,
-                                    UShort_t   s, 
-                                    UShort_t   t, 
-                                    UShort_t   nR,
-                                    Double_t   theta,
-                                    AliESDFMD& output) const
-{
-  // 
-  // Store a particle hit in FMD<i>dr</i>[<i>s,t</i>] in @a output
-  // 
-  // Parameters:
-  //    d       Detector
-  //    r       Ring
-  //    s       Sector
-  //    t       Strip
-  //    nR      Number of references to this particle in this sector
-  //    output  Output ESD object
-  //
-  Double_t old = output.Multiplicity(d,r,s,t);
-  if (old == AliESDFMD::kInvalidMult) old = 0;
-  if (fOper) fOperComp->Fill(fOper->operator()(d,r,s,t), nR);
-  if (theta < 0) theta += 2*TMath::Pi();
-  theta *= 180. / TMath::Pi();
-  fThetaVsNr->Fill(theta, nR);
-  output.SetMultiplicity(d,r,s,t,old+1);
-}
 
 //____________________________________________________________________
 Bool_t
@@ -216,119 +162,9 @@ AliFMDMCSharingFilter::FilterMC(const AliESDFMD&  input,
   //
   output.Clear();
 
-  // Increase event count - stored in 
-  // underflow bin 
-  fSumEta->AddBinContent(0); 
-
-  // Copy eta values to output 
-  for (UShort_t ed = 1; ed <= 3; ed++) { 
-    UShort_t nq = (ed == 1 ? 1 : 2);
-    for (UShort_t eq = 0; eq < nq; eq++) {
-      Char_t   er = (eq == 0 ? 'I' : 'O');
-      UShort_t ns = (eq == 0 ?  20 :  40);
-      UShort_t nt = (eq == 0 ? 512 : 256);
-      for (UShort_t es = 0; es < ns; es++) 
-       for (UShort_t et = 0; et < nt; et++) 
-         output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et));
-    }
-  }
-  AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
-  Int_t nTracks   = stack->GetNtrack();//event.GetNumberOfTracks();
-  Int_t nPrim     = stack->GetNprimary();//event.GetNumberOfPrimary();
-  for (Int_t iTr = 0; iTr < nTracks; iTr++) { 
-    AliMCParticle* particle = 
-      static_cast<AliMCParticle*>(event.GetTrack(iTr));
-    
-    // Check the returned particle 
-    if (!particle) continue;
-    
-    // Check if this charged and a primary 
-    Bool_t isCharged = particle->Charge() != 0;
-    if (!isCharged) continue;
-    Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
-
-    // Pseudo rapidity and azimuthal angle 
-    Double_t eta = particle->Eta();
-    Double_t phi = particle->Phi();
-    
-    // Fill 'dn/deta' histogram 
-    if (isPrimary && iTr < nPrim) { 
-      // Avoid under count - used to store event count
-      if (eta >= fSumEta->GetXaxis()->GetXmin()) fSumEta->Fill(eta);
-      primary->Fill(eta, phi);
-    }
-
-    // Bail out if we're only processing primaries - perhaps we should
-    // track back to the original primary?
-    if (fOnlyPrimary && !isPrimary) continue;
-
-    Int_t    nTrRef  = particle->GetNumberOfTrackReferences();
-    Int_t    longest = -1;
-    Double_t angle   = 0;
-    UShort_t oD = 0, oS = 1024, oT = 1024;
-    Char_t   oR = '\0';
-    UShort_t nC = 0;
-    Double_t oTheta = 0;
-    for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) { 
-      AliTrackReference* ref = particle->GetTrackReference(iTrRef);
-      
-      // Check existence 
-      if (!ref) continue;
 
-      // Check that we hit an FMD element 
-      if (ref->DetectorId() != AliTrackReference::kFMD) 
-       continue;
-
-      // Count number of track refs in this sector 
-      nC++;
-
-      // Get the detector coordinates 
-      UShort_t d, s, t;
-      Char_t r;
-      AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
-      // If this is a new detector/ring, then reset the other one 
-      if (oD > 0 && oR != '\0' && oS != 1024 && 
-         (d != oD || r != oR || s != oS)) {
-       longest = -1;
-       angle   = 0;
-       StoreParticle(oD, oR, oS, oT, nC, oTheta, output);
-       nC = 0;
-       oD = 0;
-       oR = '\0';
-       oS = 1024;
-       oT = 1024;
-      }
+  fTrackDensity.Calculate(input, event, vz, output, primary);
 
-      oD = d;
-      oR = r;
-      oS = s;
-      oT = t;
-
-      // The longest passage is determined through the angle 
-      Double_t x    = ref->X();
-      Double_t y    = ref->Y();
-      Double_t z    = ref->Z()-vz;
-      Double_t rr   = TMath::Sqrt(x*x+y*y);
-      Double_t theta= TMath::ATan2(rr,z);
-      Double_t ang  = TMath::Abs(TMath::Pi()-theta);
-      if (ang > angle) {
-       longest = iTrRef;
-       angle   = ang;
-      }
-      oTheta = theta;
-    } // Loop over track references
-    if (longest < 0) continue;
-
-    // Get the reference corresponding to the longest path through the detector
-    AliTrackReference* ref = particle->GetTrackReference(longest);
-
-    // Get the detector coordinates 
-    UShort_t d, s, t;
-    Char_t r;
-    AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
-    
-    StoreParticle(d,r,s,t,nC,particle->Theta(),output);
-  } // Loop over tracks
   return kTRUE;
 }
 
@@ -387,6 +223,7 @@ AliFMDMCSharingFilter::DefineOutput(TList* dir)
   AliFMDSharingFilter::DefineOutput(dir);
   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
   TList* cd = new TList;
+  cd->SetOwner();
   cd->SetName("esd_mc_comparion");
   d->Add(cd);
   cd->Add(fFMD1i);
@@ -394,9 +231,8 @@ AliFMDMCSharingFilter::DefineOutput(TList* dir)
   cd->Add(fFMD2o);
   cd->Add(fFMD3i);
   cd->Add(fFMD3o);
-  dir->Add(fSumEta);
   cd->Add(fOperComp);
-  cd->Add(fThetaVsNr);
+  fTrackDensity.DefineOutput(d);
 }
 
 //____________________________________________________________________
@@ -411,13 +247,6 @@ AliFMDMCSharingFilter::ScaleHistograms(const TList* dir, Int_t nEvents)
   //    nEvents Number of events 
   //
   AliFMDSharingFilter::ScaleHistograms(dir, nEvents);
-  TH1D* sumEta = static_cast<TH1D*>(dir->FindObject("mcSumEta"));
-  if (!sumEta) { 
-    AliWarning(Form("No mcSumEta histogram found in output list"));
-    return;
-  }
-  Double_t n = nEvents; // sumEta->GetBinContent(0);
-  sumEta->Scale(1. / n, "width");
 }
 
 //____________________________________________________________________
@@ -434,9 +263,10 @@ AliFMDMCSharingFilter::Print(Option_t* option) const
   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
   ind[gROOT->GetDirLevel()] = '\0';
   AliFMDSharingFilter::Print(option);
-  std::cout << std::boolalpha 
-           << ind << " Only primary tracks:    " << fOnlyPrimary 
-           << std::noboolalpha << std::endl;
+  gROOT->IncreaseDirLevel();
+  fTrackDensity.Print(option);
+  gROOT->DecreaseDirLevel();
+
 }
 
 //____________________________________________________________________
index 09aa4e0..6e14e17 100644 (file)
@@ -14,7 +14,7 @@
  * @ingroup pwg2_forward_aod
  */
 #include "AliFMDSharingFilter.h"
-class AliMCEvent;
+#include "AliFMDMCTrackDensity.h"
 
 /**
  * Class to do the sharing correction for MC data.
@@ -53,15 +53,13 @@ public:
    */
   AliFMDMCSharingFilter()
   : AliFMDSharingFilter(), 
+    fTrackDensity(),
     fFMD1i(0),
     fFMD2i(0),
     fFMD2o(0),
     fFMD3i(0),
     fFMD3o(0), 
-    fSumEta(0), 
-    fOperComp(0), 
-    fThetaVsNr(0), 
-    fOnlyPrimary(false)
+    fOperComp(0)
   {}
   /** 
    * Constructor 
@@ -83,13 +81,20 @@ public:
    * @return Reference to this 
    */
   AliFMDMCSharingFilter& operator=(const AliFMDMCSharingFilter& o);
-  
+
+  /** 
+   * Return the track density calculator 
+   * 
+   * @return Track density calculator 
+   */
+  const AliFMDMCTrackDensity& GetTrackDensity() const { return fTrackDensity; }
   /** 
-   * If set, then only process primary tracks 
+   * Return the track density calculator 
    * 
-   * @param use 
+   * @return Track density calculator 
    */
-  void SetOnlyPrimary(Bool_t use) { fOnlyPrimary = use; }
+  AliFMDMCTrackDensity& GetTrackDensity() { return fTrackDensity; }
+
   /** 
    * Filter the input kinematics and track references, using 
    * some of the ESD information
@@ -138,26 +143,13 @@ public:
    */
   void Print(Option_t* option="") const;
 protected:
-  /** 
-   * Store a particle hit in FMD<i>dr</i>[<i>s,t</i>] in @a output
-   * 
-   * @param d       Detector
-   * @param r       Ring
-   * @param s       Sector
-   * @param t       Strip
-   * @param output  Output ESD object
-   */
-  void StoreParticle(UShort_t d, Char_t r, UShort_t s, UShort_t t, 
-                    UShort_t nr, Double_t theta, AliESDFMD& output) const;
-  TH2D* fFMD1i;  // ESD-MC correlation 
-  TH2D* fFMD2i;  // ESD-MC correlation 
-  TH2D* fFMD2o;  // ESD-MC correlation 
-  TH2D* fFMD3i;  // ESD-MC correlation 
-  TH2D* fFMD3o;  // ESD-MC correlation 
-  TH1D* fSumEta; // MC dN/deta 
-  TH2I* fOperComp; // Operation vs # trackrefs
-  TH2D* fThetaVsNr; // Theta vs # trackrefs
-  Bool_t fOnlyPrimary; // Only process primary tracks 
+  AliFMDMCTrackDensity fTrackDensity;
+  TH2D* fFMD1i;      // ESD-MC correlation 
+  TH2D* fFMD2i;      // ESD-MC correlation 
+  TH2D* fFMD2o;      // ESD-MC correlation 
+  TH2D* fFMD3i;      // ESD-MC correlation 
+  TH2D* fFMD3o;      // ESD-MC correlation 
+  TH2I* fOperComp;   // Operation vs # trackrefs
   ClassDef(AliFMDMCSharingFilter,1); //
 };