Added sub-tasks to do hit counting based on track-refs.
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 26 Apr 2011 20:17:56 +0000 (20:17 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 26 Apr 2011 20:17:56 +0000 (20:17 +0000)
This is used in the MC analysis and correction tasks to
ensure consistency.

PWG2/FORWARD/analysis2/AliFMDMCTrackDensity.cxx [new file with mode: 0644]
PWG2/FORWARD/analysis2/AliFMDMCTrackDensity.h [new file with mode: 0644]
PWG2/FORWARD/analysis2/AliSPDMCTrackDensity.cxx [new file with mode: 0644]
PWG2/FORWARD/analysis2/AliSPDMCTrackDensity.h [new file with mode: 0644]

diff --git a/PWG2/FORWARD/analysis2/AliFMDMCTrackDensity.cxx b/PWG2/FORWARD/analysis2/AliFMDMCTrackDensity.cxx
new file mode 100644 (file)
index 0000000..66b5f54
--- /dev/null
@@ -0,0 +1,398 @@
+#include "AliFMDMCTrackDensity.h"
+#include <AliESDFMD.h>
+#include <AliMCEvent.h>
+#include <AliTrackReference.h>
+#include <AliStack.h>
+#include <TMath.h>
+#include "AliFMDStripIndex.h"
+#include "AliFMDFloatMap.h"
+#include <AliLog.h>
+#include <TH2D.h>
+#include <TH1D.h>
+#include <TList.h>
+#include <TROOT.h>
+#include <iostream>
+
+//____________________________________________________________________
+AliFMDMCTrackDensity::AliFMDMCTrackDensity()
+  : TNamed(), 
+    fUseOnlyPrimary(false), 
+    fMaxConsequtiveStrips(3), 
+    fNr(0), 
+    fNt(0),
+    fBinFlow(0), 
+    fEtaBinFlow(0),
+    fPhiBinFlow(0),
+    fDebug(false)
+{
+  // Default constructor 
+}
+
+//____________________________________________________________________
+AliFMDMCTrackDensity::AliFMDMCTrackDensity(const char*)
+  : TNamed("fmdMCTrackDensity","fmdMCTrackDensity"), 
+    fUseOnlyPrimary(false), 
+    fMaxConsequtiveStrips(3), 
+    fNr(0), 
+    fNt(0),
+    fBinFlow(0), 
+    fEtaBinFlow(0),
+    fPhiBinFlow(0),
+    fDebug(false)
+{
+  // Normal constructor constructor 
+}
+
+//____________________________________________________________________
+AliFMDMCTrackDensity::AliFMDMCTrackDensity(const AliFMDMCTrackDensity& o)
+  : TNamed(o),
+    fUseOnlyPrimary(o.fUseOnlyPrimary), 
+    fMaxConsequtiveStrips(o.fMaxConsequtiveStrips), 
+    fNr(o.fNr), 
+    fNt(o.fNt),
+    fBinFlow(o.fBinFlow), 
+    fEtaBinFlow(o.fEtaBinFlow),
+    fPhiBinFlow(o.fPhiBinFlow),
+    fDebug(o.fDebug)
+{
+  // Normal constructor constructor 
+}
+
+//____________________________________________________________________
+AliFMDMCTrackDensity&
+AliFMDMCTrackDensity::operator=(const AliFMDMCTrackDensity& o)
+{
+  // Assignment operator 
+  TNamed::operator=(o);
+  fUseOnlyPrimary       = o.fUseOnlyPrimary;
+  fMaxConsequtiveStrips = o.fMaxConsequtiveStrips;
+  fNr                   = o.fNr;
+  fNt                   = o.fNt;
+  fBinFlow              = o.fBinFlow;
+  fEtaBinFlow           = o.fEtaBinFlow;
+  fPhiBinFlow           = o.fPhiBinFlow;
+  fDebug                = o.fDebug;
+  return *this;
+}
+
+//____________________________________________________________________
+void
+AliFMDMCTrackDensity::DefineOutput(TList* l)
+{
+  fNr = new TH1D("clusterRefs", "# track references per cluster",
+                21, -.5, 20.5);
+  fNr->SetXTitle("# of track references/cluster");
+  fNr->SetDirectory(0);
+  l->Add(fNr);
+
+  fNt = new TH1D("clusterSize", "cluster length in strips", 21, -.5, 20.5);
+  fNt->SetXTitle("Cluster size [strips]");
+  fNt->SetDirectory(0);
+  l->Add(fNt);
+
+  fBinFlow = new TH2D("binFlow", "#eta and #varphi bin flow", 
+                     200, -5, 5, 40, -180, 180);
+  fBinFlow->SetXTitle("#Delta#eta");
+  fBinFlow->SetYTitle("#Delta#varphi");
+  fBinFlow->SetDirectory(0);
+  l->Add(fBinFlow);
+
+  fEtaBinFlow = new TH2D("binFlowEta", "#eta bin flow vs #eta", 
+                        200, -4, 6, 200, -5, 5);
+  fEtaBinFlow->SetXTitle("#eta of strip");
+  fEtaBinFlow->SetYTitle("#Delta#eta");
+  fEtaBinFlow->SetDirectory(0);
+  l->Add(fEtaBinFlow);
+
+  fPhiBinFlow = new TH2D("binFlowPhi", "#phi bin flow vs #phi", 
+                        40, 0, 360, 40, -180, 180);
+  fPhiBinFlow->SetXTitle("#varphi of strip");
+  fPhiBinFlow->SetYTitle("#Delta#varphi");
+  fPhiBinFlow->SetDirectory(0);
+  l->Add(fPhiBinFlow);
+}
+
+//____________________________________________________________________
+void
+AliFMDMCTrackDensity::StoreParticle(AliMCParticle*       particle, 
+                                   const AliMCParticle* mother, 
+                                   Int_t                longest,
+                                   Double_t             vz,
+                                   UShort_t             nC, 
+                                   UShort_t             nT,
+                                   AliESDFMD&           output) const
+{
+  // Store a particle. 
+  if (longest < 0) return;
+
+  AliTrackReference* ref = particle->GetTrackReference(longest);
+  if (!ref) return;
+    
+  // Get the detector coordinates 
+  UShort_t d, s, t;
+  Char_t r;
+  AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
+
+  // Check if we have value already 
+  Double_t old = output.Multiplicity(d,r,s,t);
+
+  // If invalid, force it valid 
+  if (old == AliESDFMD::kInvalidMult) old = 0;
+
+  // Increment count 
+  output.SetMultiplicity(d,r,s,t,old+1);
+
+  // Get track-reference stuff 
+  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 thetaR = TMath::ATan2(rr,z);
+  Double_t phiR   = TMath::ATan2(y,x);
+  Double_t etaR   = -TMath::Log(TMath::Tan(thetaR/2));
+  
+  // Correct angle and convert to degrees 
+  if (thetaR < 0) thetaR += 2*TMath::Pi();
+  thetaR *= 180. / TMath::Pi();
+  if (phiR < 0) phiR += 2*TMath::Pi();
+  phiR *= 180. / TMath::Pi();
+
+  // Fill histograms 
+  fNr->Fill(nC);
+  fNt->Fill(nT);
+
+  const AliMCParticle* mp = (mother ? mother : particle);
+  Double_t dEta = mp->Eta() - etaR;
+  Double_t dPhi = mp->Phi() * 180 / TMath::Pi() - phiR;
+  if (dPhi >  180) dPhi -= 360;
+  if (dPhi < -180) dPhi += 360;
+  fBinFlow->Fill(dEta, dPhi);
+  fEtaBinFlow->Fill(etaR, dEta);
+  fPhiBinFlow->Fill(phiR, dPhi);
+
+  if (fDebug) 
+    Info("StoreParticle", "Store @ FMD%d%c[%2d,%3d] nStrips=%3d, "
+        "dEta=%7.4f, dPhi=%d",
+        d, r, s, t, nT, dEta, int(dPhi+.5));
+}
+
+//____________________________________________________________________
+Double_t
+AliFMDMCTrackDensity::GetTrackRefTheta(const AliTrackReference* ref,
+                                      Double_t vz) const
+{
+  // Get the incidient angle of the track reference. 
+  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);
+  return ang;
+}
+                                   
+//____________________________________________________________________
+const AliMCParticle*
+AliFMDMCTrackDensity::GetMother(Int_t     iTr,
+                               const AliMCEvent& event) const
+{
+  // 
+  // Track down primary mother 
+  // 
+  Int_t i  = iTr;
+  do { 
+    const AliMCParticle* p = static_cast<AliMCParticle*>(event.GetTrack(i));
+    if (const_cast<AliMCEvent&>(event).Stack()->IsPhysicalPrimary(i)) return p;
+    
+    i = p->GetMother();
+  } while (i > 0);
+
+  return 0;
+}  
+
+//____________________________________________________________________
+Bool_t
+AliFMDMCTrackDensity::Calculate(const AliESDFMD&  input, 
+                               const AliMCEvent& event,
+                               Double_t          vz,
+                               AliESDFMD&        output, 
+                               TH2D*             primary)
+{
+  // 
+  // Filter the input kinematics and track references, using 
+  // some of the ESD information
+  // 
+  // Parameters:
+  //    input   Input ESD event
+  //    event   Input MC event
+  //    vz      Vertex position 
+  //    output  Output ESD-like object
+  //    primary Per-event histogram of primaries 
+  //
+  // Return:
+  //    True on succes, false otherwise 
+  //
+  output.Clear();
+
+  // 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) {
+      if (primary) primary->Fill(eta, phi);
+    }
+
+    // Bail out if we're only processing primaries - perhaps we should
+    // track back to the original primary?
+    if (fUseOnlyPrimary && !isPrimary) continue;
+
+    Int_t    nTrRef   = particle->GetNumberOfTrackReferences();
+    Int_t    longest  = -1;
+    Double_t angle    = 0;
+    UShort_t oD       = 0, oS = 1024, oT = 1024, ooT=1024;
+    Char_t   oR       = '\0';
+    UShort_t nC       = 0;
+    UShort_t nT       = 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;
+
+      // Get the detector coordinates 
+      UShort_t d, s, t;
+      Char_t r;
+      AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
+
+      // Calculate length of last and second to last strips. 
+
+      // IF base of cluster not set, set it here. 
+      if (ooT == 1024) ooT = t;
+
+      // Calculate distance of previous reference to base of cluster 
+      nT = TMath::Abs(t - ooT) + 1;
+
+      // Count number of track refs in this sector 
+      nC++;
+
+      Bool_t used = false;
+      // If this is a new detector/ring, then reset the other one 
+      // Check if we have a valid old detectorm ring, and sector 
+      if (oD >  0 && oR != '\0' && oS != 1024) {
+       // New detector, new ring, or new sector 
+       if (d != oD   || r != oR   || s != oS) {
+         if (fDebug) Info("Process", "New because new sector");
+         used = true;
+       }
+      }
+      if (nT > fMaxConsequtiveStrips) {
+       if (fDebug) Info("Process", "New because too long: %d (%d,%d,%d)", 
+                        nT, t, oT, ooT);
+       used = true;
+      }
+      if (used) {
+       if (fDebug) 
+         Info("Process", "I=%3d L=%3d D=%d (was %d), R=%c (was %c), "
+              "S=%2d (was %2d) t=%3d (was %3d) nT=%3d/%4d",
+              iTr, longest, d, oD, r, oR, s, oS, t, oT, nT, 
+              fMaxConsequtiveStrips);
+       Int_t nnT   = TMath::Abs(oT - ooT) + 1;
+       const AliMCParticle* mother = GetMother(iTr, event);
+       StoreParticle(particle, mother, longest, vz, nC, nnT, output);
+       longest = -1;
+       angle   = 0;
+       nC  = 1;    // Reset track-ref counter - we have this->1
+       nT  = 0;    // Reset cluster size 
+       oD  = 0;    // Reset detector
+       oR  = '\0'; // Reset ring 
+       oS  = 1024; // Reset sector 
+       oT  = 1024; // Reset old strip
+       ooT = t;    // Reset base 
+      }
+      else if (fDebug) {
+       if (ooT == t) 
+         Info("Process", "New cluster starting at FMD%d%c[%2d,%3d]", 
+              d, r, s, t);
+       else 
+         Info("Process", "Adding to cluster starting at FMD%d%c[%2d,%3d], "
+              "length=%3d (now in %3d, previous %3d)", 
+              d, r, s, ooT, nT, t, oT);
+      }
+       
+      oD = d;
+      oR = r;
+      oS = s;
+      oT = t;
+      nT = TMath::Abs(t-ooT)+1;
+
+      // The longest passage is determined through the angle 
+      Double_t ang  = GetTrackRefTheta(ref, vz);
+      if (ang > angle) {
+       longest = iTrRef;
+       angle   = ang;
+      }
+      oTheta = ang;
+    } // Loop over track references
+    if (longest < 0) continue; // Nothing found
+
+    // Get the reference corresponding to the longest path through the detector
+    if (fDebug) 
+      Info("Process", "I=%3d L=%3d nT=%3d (out of %3d)", 
+          iTr, longest, nT, fMaxConsequtiveStrips);
+    const AliMCParticle* mother = GetMother(iTr, event);
+    StoreParticle(particle, mother, longest, vz, nC, nT, output);
+  } // Loop over tracks
+  return kTRUE;
+}
+//____________________________________________________________________
+void
+AliFMDMCTrackDensity::Print(Option_t* /*option*/) const 
+{
+  char ind[gROOT->GetDirLevel()+1];
+  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
+  ind[gROOT->GetDirLevel()] = '\0';
+  std::cout << ind << ClassName() << ": " << GetName() << '\n'
+           << std::boolalpha 
+           << ind << " Only primary tracks:    " << fUseOnlyPrimary << '\n'
+           << ind << " Max cluster size:       " << fMaxConsequtiveStrips
+           << std::noboolalpha << std::endl;
+  
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/PWG2/FORWARD/analysis2/AliFMDMCTrackDensity.h b/PWG2/FORWARD/analysis2/AliFMDMCTrackDensity.h
new file mode 100644 (file)
index 0000000..dfd894a
--- /dev/null
@@ -0,0 +1,176 @@
+#ifndef ALIFMDMCTRACKDENSITY_MC
+#define ALIFMDMCTRACKDENSITY_MC
+#include "AliForwardUtil.h"
+#include <TNamed.h>
+class TList;
+class TH1D;
+class TH2D;
+class AliMCEvent;
+class AliESDFMD;
+class AliMCParticle;
+class AliTrackReference;
+class AliStack;
+
+/**
+ * A class to calculate the particle density from track references.
+ * This code is used both in AliForwardMCCorrectionsTask and
+ * AliFMDMCDensity calculator. 
+ * 
+ * @par Input: 
+ *    - AliESDFMD object  - from reconstruction
+ *    - Kinematics
+ *    - Track-References
+ *
+ * @par Output: 
+ *    - AliESDFMD object  - content is # of track references/strip
+ *
+ * @par Corrections used: 
+ *    - None
+ *
+ * @par Histograms: 
+ *    - Incident angle vs number of track references
+ *    - Incident angle vs number of strips/cluster
+ *
+ * @ingroup pwg2_forward_algo
+ * @ingroup pwg2_forward_mc
+ * @ingroup pwg2_forward_aod
+ */
+class AliFMDMCTrackDensity : public TNamed
+{
+public:
+  /** 
+   * Default constructor.  Do not use - for ROOT I/O system use only 
+   */
+  AliFMDMCTrackDensity();
+  /** 
+   * Normal constructor 
+   * 
+   * @param name Not used
+   */
+  AliFMDMCTrackDensity(const char* name);
+  /** 
+   * Copy constructor 
+   * 
+   * @param o Object to copy from 
+   */
+  AliFMDMCTrackDensity(const AliFMDMCTrackDensity& o);
+  /** 
+   * Assignment operator
+   * 
+   * @param o Object to assign from 
+   * 
+   * @return Reference to this 
+   */
+  AliFMDMCTrackDensity& operator=(const AliFMDMCTrackDensity& o);
+  /** 
+   * Destructor. 
+   */
+  virtual ~AliFMDMCTrackDensity() {}
+
+  /** 
+   * Set maximum number of strips per 'cluster' 
+   * 
+   * @param n  Maximum number of strips per 'cluster' 
+   */
+  void SetMaxConsequtiveStrips(UShort_t n) { fMaxConsequtiveStrips = n; }
+  /** 
+   * Set whether to only consider primaries 
+   * 
+   * @param use If true, consider only primaries
+   */
+  void SetUseOnlyPrimary(Bool_t use) { fUseOnlyPrimary = use; }
+  
+  /** 
+   * Set whether to print debug messages.  Please note this will
+   * produce a lot of output. 
+   * 
+   * @param debug Whether to enable debug messages or not 
+   */
+  void SetDebug(Bool_t debug=true) { fDebug = debug; }
+  /** 
+   * Loops over all the particles in the passed event.  If @a primary
+   * is not null, then that histogram is filled with the primary
+   * particle information - irrespective of whether the particle
+   * actually hits the FMD or not.  For each track (primary or
+   * secondary, unless only primary information is requested - see
+   * SetUseOnlyPrimary) loop over all track references to that
+   * particle and check if they come from the FMD.  In that case,
+   * figure out which strip(s) to assign the track to, and fill the @a
+   * hits structure.
+   * 
+   * @param esd      FMD ESD structure 
+   * @param event    MC event 
+   * @param vz       IP z-coordinate
+   * @param output   Output of FMD hits
+   * @param primary  Primary information, if available. 
+   * 
+   * @return true 
+   */
+  Bool_t Calculate(const AliESDFMD&    esd, 
+                  const AliMCEvent&   event, 
+                  Double_t            vz,
+                  AliESDFMD&          output,
+                  TH2D*               primary);
+  /** 
+   * Define ouputs 
+   * 
+   * @param list List to add outputs to
+   */
+  void DefineOutput(TList* list);
+  
+  void Print(Option_t* option="") const;
+protected:
+  /** 
+   * Store a particle hit in FMD<i>dr</i>[<i>s,t</i>] in @a output
+   * 
+   * 
+   * @param particle  Particle to store
+   * @param mother    Ultimate mother of particle 
+   * @param longest   Longest track reference
+   * @param vz        Z coordinate of IP
+   * @param nC        Total number of track-references in this sector  
+   * @param nT               Number of distint strips hit in this sector
+   * @param output    Output structure 
+   */  
+  void StoreParticle(AliMCParticle* particle, 
+                    const AliMCParticle* mother,
+                    Int_t          longest,
+                    Double_t       vz,
+                    UShort_t       nC, 
+                    UShort_t       nT,
+                    AliESDFMD&     output) const;
+  /** 
+   * Get incident angle of this track reference
+   * 
+   * @param ref Track reference
+   * @param vz  Z coordinate of the IP
+   * 
+   * @return incident angle (in radians)
+   */
+  Double_t GetTrackRefTheta(const AliTrackReference* ref,
+                           Double_t vz) const;
+  /** 
+   * Get ultimate mother of a track 
+   * 
+   * @param iTr   Track number 
+   * @param stack Stack 
+   * 
+   * @return Pointer to mother or null 
+   */
+  const AliMCParticle* GetMother(Int_t iTr, const AliMCEvent& event) const;
+  Bool_t   fUseOnlyPrimary;       // Only use primaries 
+  UShort_t fMaxConsequtiveStrips; // Max 'cluster' size
+  TH1D*    fNr;                   // Number of track-refs per cluster
+  TH1D*    fNt;                   // Size of cluster in strips 
+  TH2D*    fBinFlow;              // eta,phi bin flow 
+  TH2D*    fEtaBinFlow;           // dEta vs eta of strip
+  TH2D*    fPhiBinFlow;           // dPhi vs phi of strip
+  Bool_t   fDebug;                // Debug flag
+
+  ClassDef(AliFMDMCTrackDensity,1); // Calculate track-ref density
+};
+
+#endif
+// Local Variables:
+//  mode: C++ 
+// End:
diff --git a/PWG2/FORWARD/analysis2/AliSPDMCTrackDensity.cxx b/PWG2/FORWARD/analysis2/AliSPDMCTrackDensity.cxx
new file mode 100644 (file)
index 0000000..d49f00a
--- /dev/null
@@ -0,0 +1,272 @@
+#include "AliSPDMCTrackDensity.h"
+#include <AliMCEvent.h>
+#include <AliTrackReference.h>
+#include <AliStack.h>
+#include <TMath.h>
+#include <AliLog.h>
+#include <TH3D.h>
+#include <TH2D.h>
+#include <TH1D.h>
+#include <TList.h>
+#include <TROOT.h>
+#include <iostream>
+
+//____________________________________________________________________
+AliSPDMCTrackDensity::AliSPDMCTrackDensity()
+  : TNamed(), 
+    fUseOnlyPrimary(false), 
+    fMinR(3.5), 
+    fMaxR(4.5),
+    fMinZ(-14.1), 
+    fMaxZ(+14.1),
+    fRZ(0), 
+    fXYZ(0),
+    fNRefs(0),
+    fBinFlow(0)
+{
+  // Default constructor 
+}
+
+//____________________________________________________________________
+AliSPDMCTrackDensity::AliSPDMCTrackDensity(const char*)
+  : TNamed("spdMCTrackDensity","spdMCTrackDensity"), 
+    fUseOnlyPrimary(false), 
+    fMinR(3.5), 
+    fMaxR(4.5),
+    fMinZ(-14.1), 
+    fMaxZ(+14.1),
+    fRZ(0), 
+    fXYZ(0),
+    fNRefs(0),
+    fBinFlow(0)
+{
+  // Normal constructor constructor 
+}
+
+//____________________________________________________________________
+AliSPDMCTrackDensity::AliSPDMCTrackDensity(const AliSPDMCTrackDensity& o)
+  : TNamed(o),
+    fUseOnlyPrimary(o.fUseOnlyPrimary), 
+    fMinR(o.fMinR), 
+    fMaxR(o.fMaxR),
+    fMinZ(o.fMinZ), 
+    fMaxZ(o.fMaxZ),
+    fRZ(o.fRZ), 
+    fXYZ(o.fXYZ),
+    fNRefs(o.fNRefs),
+    fBinFlow(o.fBinFlow)
+{
+  // Normal constructor constructor 
+}
+
+//____________________________________________________________________
+AliSPDMCTrackDensity&
+AliSPDMCTrackDensity::operator=(const AliSPDMCTrackDensity& o)
+{
+  // Assignment operator 
+  TNamed::operator=(o);
+  fUseOnlyPrimary       = o.fUseOnlyPrimary;
+  fMinR             = o.fMinR;
+  fMaxR             = o.fMaxR;
+  fMinZ             = o.fMinZ;
+  fMaxZ             = o.fMaxZ;
+  fRZ               = o.fRZ;
+  fXYZ              = o.fXYZ;
+  fNRefs            = o.fNRefs;
+  fBinFlow          = o.fBinFlow;
+  return *this;
+}
+
+//____________________________________________________________________
+void
+AliSPDMCTrackDensity::DefineOutput(TList* l)
+{
+  fRZ = new TH2D("rz", "(r,z) of used track references", 
+                50, 0, 5, 30, -15, 15);
+  fRZ->SetXTitle("z [cm]");
+  fRZ->SetYTitle("r [cm]");
+  fRZ->SetDirectory(0);
+  l->Add(fRZ);
+
+  fXYZ = new TH3D("xyz", "(x,y,z) of used track references", 
+                 100, -5., 5., 100, -5., 5., 30, -15., 15.);
+  fXYZ->SetXTitle("x [cm]");
+  fXYZ->SetYTitle("y [cm]");
+  fXYZ->SetZTitle("z [cm]");
+  fXYZ->SetDirectory(0);
+  l->Add(fXYZ);
+
+  fNRefs = new TH1D("nrefs", "Number of references used per track", 
+                   11, -.5, 10.5); 
+  fNRefs->SetXTitle("# of references per track");
+  fNRefs->SetDirectory(0);
+  l->Add(fNRefs);
+
+  fBinFlow = new TH2D("binFlow", "#eta and #varphi bin flow", 
+                     120, -3, 3, 40, -180, 180);
+  fBinFlow->SetXTitle("#Delta#eta");
+  fBinFlow->SetYTitle("#Delta#varphi");
+  fBinFlow->SetDirectory(0);
+  l->Add(fBinFlow);
+}
+
+//____________________________________________________________________
+void
+AliSPDMCTrackDensity::StoreParticle(AliMCParticle* particle, 
+                                   const AliMCParticle* mother, 
+                                   Int_t          refNo,
+                                   Double_t       vz, 
+                                   TH2D&     output) const
+{
+  // Store a particle. 
+  if (refNo < 0) return;
+
+  AliTrackReference* ref = particle->GetTrackReference(refNo);
+  if (!ref) return;
+    
+  Double_t r = ref->R();
+  Double_t x = ref->X();
+  Double_t y = ref->Y();
+  Double_t z = ref->Z();
+
+  Double_t zr = z-vz;
+  Double_t th = TMath::ATan2(r,zr);
+  if (th < 0) th += 2*TMath::Pi();
+  Double_t et = -TMath::Log(TMath::Tan(th/2));
+  Double_t ph = TMath::ATan2(y,x);
+  if (ph < 0) ph += 2*TMath::Pi();
+  output.Fill(et,ph);
+
+  const AliMCParticle* mp = (mother ? mother : particle);
+  Double_t dEta = mp->Eta() - et;
+  Double_t dPhi = (mp->Phi() - ph) * 180 / TMath::Pi();
+  if (dPhi >  180) dPhi -= 360;
+  if (dPhi < -180) dPhi += 360;
+  fBinFlow->Fill(dEta, dPhi);
+}
+
+
+//____________________________________________________________________
+const AliMCParticle*
+AliSPDMCTrackDensity::GetMother(Int_t     iTr,
+                               const AliMCEvent& event) const
+{
+  // 
+  // Track down primary mother 
+  // 
+  Int_t i  = iTr;
+  do { 
+    const AliMCParticle* p = static_cast<AliMCParticle*>(event.GetTrack(i));
+    if (const_cast<AliMCEvent&>(event).Stack()->IsPhysicalPrimary(i)) return p;
+    
+    i = p->GetMother();
+  } while (i > 0);
+
+  return 0;
+}  
+
+//____________________________________________________________________
+Bool_t
+AliSPDMCTrackDensity::Calculate(const AliMCEvent& event, 
+                               Double_t          vz,
+                               TH2D&             output, 
+                               TH2D*             primary)
+{
+  // 
+  // Filter the input kinematics and track references, using 
+  // some of the ESD information
+  // 
+  // Parameters:
+  //    input   Input ESD event
+  //    event   Input MC event
+  //    vz      Vertex position 
+  //    output  Output ESD-like object
+  //    primary Per-event histogram of primaries 
+  //
+  // Return:
+  //    True on succes, false otherwise 
+  //
+
+  AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
+  Int_t nTracks   = stack->GetNtrack();
+  Int_t nPrim     = stack->GetNtrack();
+  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);
+
+    // Fill 'dn/deta' histogram 
+    if (isPrimary && iTr < nPrim) {
+      if (primary) {
+       Double_t eta = particle->Eta();
+       Double_t phi = particle->Phi();
+       primary->Fill(eta, phi);
+      }
+    }
+
+    // Bail out if we're only processing primaries - perhaps we should
+    // track back to the original primary?
+    if (fUseOnlyPrimary && !isPrimary) continue;
+
+    Int_t nTrRef  = particle->GetNumberOfTrackReferences();
+    Int_t nRef    = 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::kITS) 
+       continue;
+      
+      // Get radius and z where the track reference was made 
+      Double_t r = ref->R();
+      Double_t x = ref->X();
+      Double_t y = ref->Y();
+      Double_t z = ref->Z();
+      if (r > fMaxR || r < fMinR) continue;
+      if (z > fMaxZ || z < fMinZ) continue;
+
+      fRZ->Fill(z,r);
+      fXYZ->Fill(x, y, z);
+      nRef++;
+      // Only fill first reference 
+      if (nRef == 1) { 
+       const AliMCParticle* mother = GetMother(iTrRef, event);
+       StoreParticle(particle, mother, iTrRef, vz, output);
+      }
+    }
+    fNRefs->Fill(nRef);
+  }
+  return kTRUE;
+}
+//____________________________________________________________________
+void
+AliSPDMCTrackDensity::Print(Option_t* /*option*/) const 
+{
+  char ind[gROOT->GetDirLevel()+1];
+  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
+  ind[gROOT->GetDirLevel()] = '\0';
+  std::cout << ind << ClassName() << ": " << GetName() << '\n'
+           << std::boolalpha 
+           << ind << " Only primary tracks:    " << fUseOnlyPrimary << '\n'
+           << ind << " R range:                [" << fMinR << ',' << fMaxR
+           << "]\n"
+           << ind << " Z range:                [" << fMinZ << ',' << fMaxZ
+           << "]" << std::noboolalpha << std::endl;
+  
+}
+
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/PWG2/FORWARD/analysis2/AliSPDMCTrackDensity.h b/PWG2/FORWARD/analysis2/AliSPDMCTrackDensity.h
new file mode 100644 (file)
index 0000000..d983f59
--- /dev/null
@@ -0,0 +1,154 @@
+#ifndef ALISPDMCTRACKDENSITY_MC
+#define ALISPDMCTRACKDENSITY_MC
+#include <TNamed.h>
+class TList;
+class AliMCEvent;
+class AliMultiplicity;
+class AliMCParticle;
+class AliTrackReference;
+class TH3D;
+class TH2D;
+class TH1D;
+
+/**
+ * A class to calculate the particle density from track references.
+ * This code is used both in AliForwardMCCorrectionsTask and
+ * AliSPDMCDensity calculator. 
+ * 
+ * @par Input: 
+ *    - AliMultiplicity object  - from reconstruction
+ *    - Kinematics
+ *    - Track-References
+ *
+ * @par Output: 
+ *    - AliESDSPD object  - content is # of track references/strip
+ *
+ * @par Corrections used: 
+ *    - None
+ *
+ * @par Histograms: 
+ *    - Incident angle vs number of track references
+ *    - Incident angle vs number of strips/cluster
+ *
+ * @ingroup pwg2_forward_algo
+ * @ingroup pwg2_forward_mc
+ * @ingroup pwg2_forward_aod
+ */
+class AliSPDMCTrackDensity : public TNamed
+{
+public:
+  /** 
+   * Default constructor.  Do not use - for ROOT I/O system use only 
+   */
+  AliSPDMCTrackDensity();
+  /** 
+   * Normal constructor 
+   * 
+   * @param name Not used
+   */
+  AliSPDMCTrackDensity(const char* name);
+  /** 
+   * Copy constructor 
+   * 
+   * @param o Object to copy from 
+   */
+  AliSPDMCTrackDensity(const AliSPDMCTrackDensity& o);
+  /** 
+   * Assignment operator
+   * 
+   * @param o Object to assign from 
+   * 
+   * @return Reference to this 
+   */
+  AliSPDMCTrackDensity& operator=(const AliSPDMCTrackDensity& o);
+  /** 
+   * Destructor. 
+   */
+  virtual ~AliSPDMCTrackDensity() {}
+
+  /** 
+   * Set whether to only consider primaries 
+   * 
+   * @param use If true, consider only primaries
+   */
+  void SetUseOnlyPrimary(Bool_t use) { fUseOnlyPrimary = use; }
+  /** 
+   * Loops over all the particles in the passed event.  If @a primary
+   * is not null, then that histogram is filled with the primary
+   * particle information - irrespective of whether the particle
+   * actually hits the SPD or not.  For each track (primary or
+   * secondary, unless only primary information is requested - see
+   * SetUseOnlyPrimary) loop over all track references to that
+   * particle and check if they come from the SPD.  In that case,
+   * figure out which @f$(\eta,\varphi)@f$-bin to assign the track to,
+   * and fill the @a output histogram
+   * 
+   * @param event    MC event 
+   * @param vz       IP z--coordinate
+   * @param output   Output of SPD hits
+   * @param primary  Primary information, if available. 
+   * 
+   * @return true 
+   */
+  Bool_t Calculate(const AliMCEvent&   event, 
+                  Double_t            vz,
+                  TH2D&               output,
+                  TH2D*               primary);
+  /** 
+   * Define ouputs 
+   * 
+   * @param list List to add outputs to
+   */
+  void DefineOutput(TList* list);
+  
+  void Print(Option_t* option="") const;
+protected:
+  /** 
+   * Store a particle hit in FMD<i>dr</i>[<i>s,t</i>] in @a output
+   * 
+   * 
+   * @param particle  Particle to store
+   * @param mother    Ultimate mother 
+   * @param output    Output structure 
+   */  
+  void StoreParticle(AliMCParticle*       particle, 
+                    const AliMCParticle* mother,
+                    Int_t                refNo,
+                    Double_t             vz, 
+                    TH2D&                output) const;
+  /** 
+   * Get ultimate mother of a track 
+   * 
+   * @param iTr   Track number 
+   * @param stack Stack 
+   * 
+   * @return Pointer to mother or null 
+   */
+  const AliMCParticle* GetMother(Int_t iTr, const AliMCEvent& event) const;
+  /** 
+   * Get incident angle of this track reference
+   * 
+   * @param ref Track reference
+   * @param vz  Z coordinate of the IP
+   * 
+   * @return incident angle (in radians)
+   */
+  Double_t GetTrackRefTheta(const AliTrackReference* ref,
+                           Double_t vz) const;
+  Bool_t   fUseOnlyPrimary;       // Only use primaries 
+  Double_t fMinR;             // Min radius 
+  Double_t fMaxR;             // Max radius 
+  Double_t fMinZ;             // Min z
+  Double_t fMaxZ;             // Max z
+  TH2D*    fRZ;               // Location in (r,z)
+  TH3D*    fXYZ;              // Location in (x,y,z)
+  TH1D*    fNRefs;            // Refs per track 
+  TH2D*    fBinFlow;          // eta,phi bin flow 
+
+  ClassDef(AliSPDMCTrackDensity,1); // Calculate track-ref density
+};
+
+#endif
+// Local Variables:
+//  mode: C++ 
+// End: