Updates to MC density calculator
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 5 Jan 2011 11:34:50 +0000 (11:34 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 5 Jan 2011 11:34:50 +0000 (11:34 +0000)
PWG2/FORWARD/analysis2/AliFMDMCDensityCalculator.cxx
PWG2/FORWARD/analysis2/AliFMDMCDensityCalculator.h

index 7810d52..26e41a9 100644 (file)
@@ -4,14 +4,31 @@
 #include "AliFMDStripIndex.h"
 #include "AliMCEvent.h"
 // #include "AliFMDAnaParameters.h"
+#include "AliESDFMD.h"
 #include "AliLog.h"
 #include <TH2D.h>
+#include <TProfile2D.h>
 
 ClassImp(AliFMDMCDensityCalculator)
 #if 0
 ; // For Emacs
 #endif 
 
+//____________________________________________________________________
+AliFMDMCDensityCalculator::~AliFMDMCDensityCalculator()
+{
+  if (fComps)  fComps->Clear();
+  if (fFMD1i)  delete fFMD1i;
+  if (fFMD2i)  delete fFMD2i;
+  if (fFMD2o)  delete fFMD2o;
+  if (fFMD3i)  delete fFMD3i;
+  if (fFMD3o)  delete fFMD3o;
+  if (fFMD1iC) delete fFMD1iC;
+  if (fFMD2iC) delete fFMD2iC;
+  if (fFMD2oC) delete fFMD2oC;
+  if (fFMD3iC) delete fFMD3iC;
+  if (fFMD3oC) delete fFMD3oC;
+}
 
 //____________________________________________________________________
 AliFMDMCDensityCalculator&
@@ -22,72 +39,160 @@ AliFMDMCDensityCalculator::operator=(const AliFMDMCDensityCalculator& o)
 }
 
     
+
 //____________________________________________________________________
 Bool_t
-AliFMDMCDensityCalculator::CalculateMC(const AliMCEvent&       event,
-                                      AliForwardUtil::Histos& hists,
-                                      Double_t                vz,
-                                      UShort_t                vtxbin)
+AliFMDMCDensityCalculator::CalculateMC(const AliESDFMD&        fmd,
+                                      AliForwardUtil::Histos& hists)
 {
-  Int_t nTracks = event.GetNumberOfTracks();
-  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;
+  for (UShort_t d=1; d<=3; d++) { 
+    UShort_t nr = (d == 1 ? 1 : 2);
+    for (UShort_t q=0; q<nr; q++) { 
+      Char_t      r = (q == 0 ? 'I' : 'O');
+      UShort_t    ns= (q == 0 ?  20 :  40);
+      UShort_t    nt= (q == 0 ? 512 : 256);
+      TH2D*       h = hists.Get(d,r);
 
-    Int_t nTrRef = particle->GetNumberOfTrackReferences();
-    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);
-
-      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 phi  = TMath::ATan2(y,x);
-      Double_t theta= TMath::ATan2(rr,z);
-      Double_t eta  = -TMath::Log(TMath::Tan(theta/2));
-
-      Float_t  c    = Correction(d,r,s,t,vtxbin,eta,false);
-      fCorrections->Fill(c);
-      
-      TH2D*    h    = hists.Get(d,r);
-      h->Fill(eta,phi, 1 * c);
+      for (UShort_t s=0; s<ns; s++) { 
+       for (UShort_t t=0; t<nt; t++) {
+         Float_t mult = fmd.Multiplicity(d,r,s,t);
+         
+         if (mult == 0 || mult > 20) continue;
+
+         Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
+         Float_t eta = fmd.Eta(d,r,s,t);
+
+         Float_t corr = Correction(d, r, s, t, eta, 0, false);
+         Float_t sig  = (corr <= 0 ? 0 : mult / corr);
+         h->Fill(eta,phi,sig);
+       }
+      }
     }
   }
   return kTRUE;
 }
 
 //____________________________________________________________________
+void
+AliFMDMCDensityCalculator::Init(const TAxis& eAxis)
+{
+  fFMD1i  = Make(1,'I',eAxis);
+  fFMD2i  = Make(2,'I',eAxis);
+  fFMD2o  = Make(2,'O',eAxis);
+  fFMD3i  = Make(3,'I',eAxis);
+  fFMD3o  = Make(3,'O',eAxis); 
+  fFMD1iC = Make(1,'I');
+  fFMD2iC = Make(2,'I');
+  fFMD2oC = Make(2,'O');
+  fFMD3iC = Make(3,'I');
+  fFMD3oC = Make(3,'O');
+  fComps->Add(fFMD1i);
+  fComps->Add(fFMD2i);
+  fComps->Add(fFMD2o);
+  fComps->Add(fFMD3i);
+  fComps->Add(fFMD3o);
+  fComps->Add(fFMD1iC);
+  fComps->Add(fFMD2iC);
+  fComps->Add(fFMD2oC);
+  fComps->Add(fFMD3iC);
+  fComps->Add(fFMD3oC);
+}
+
+//____________________________________________________________________
+TProfile2D*
+AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r, 
+                               const TAxis& axis) const
+{
+  TProfile2D* ret = new TProfile2D(Form("FMD%d%c_esd_vs_mc", d, r),
+                                  Form("ESD/MC signal for FMD%d%c", d, r),
+                                  axis.GetNbins(), 
+                                  axis.GetXmin(),
+                                  axis.GetXmax(), 
+                                  (r == 'I' || r == 'i') ? 20 : 40,
+                                  0, 2*TMath::Pi());
+  ret->GetXaxis()->SetTitle("#eta");
+  ret->GetYaxis()->SetTitle("#varphi [degrees]");
+  ret->GetZaxis()->SetTitle("#LT incusive density ESD/MC#GT");
+  ret->SetDirectory(0);
+  return ret;
+}
+//____________________________________________________________________
+TH2D*
+AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r) const
+{
+  TH2D* ret = new TH2D(Form("FMD%d%c_corr_mc_esd", d, r),
+                      Form("ESD-MC correlation for FMD%d%c", d, r),
+                      200, 0, 20, 200, 0, 20);
+  ret->GetXaxis()->SetTitle("m_{incl} (MC)");
+  ret->GetYaxis()->SetTitle("#Delta/#Delta_{mp} (ESD)");
+  ret->GetZaxis()->SetTitle("Correlation");
+  ret->SetDirectory(0);
+  return ret;
+}
+//____________________________________________________________________
+void
+AliFMDMCDensityCalculator::Fill(UShort_t d, Char_t r, TH2* esd, TH2* mc)
+{
+  if (!esd || !mc) return;
+  TProfile2D* p = 0;
+  TH2D*       h = 0;
+  switch (d) { 
+  case 1:  
+    p = fFMD1i;                                   
+    h = fFMD1iC;
+    break;
+  case 2:  
+    p = (r == 'I' || r == 'i' ? fFMD2i  : fFMD2o); 
+    h = (r == 'I' || r == 'i' ? fFMD2iC : fFMD2oC); 
+    break;
+  case 3:  
+    p = (r == 'I' || r == 'i' ? fFMD3i  : fFMD3o); 
+    h = (r == 'I' || r == 'i' ? fFMD3iC : fFMD3oC); 
+    break;
+  }
+  if (!p) return;
+
+  for (Int_t iEta = 1; iEta <= esd->GetNbinsX(); iEta++) { 
+    Double_t eta = esd->GetXaxis()->GetBinCenter(iEta);
+    for (Int_t iPhi = 1; iPhi <= esd->GetNbinsY(); iPhi++) { 
+      Double_t phi  = esd->GetYaxis()->GetBinCenter(iPhi);
+      Double_t mEsd = esd->GetBinContent(iEta,iPhi);
+      Double_t mMc  = mc->GetBinContent(iEta,iPhi);
+      
+      p->Fill(eta, phi, (mMc > 0 ? mEsd / mMc : 0));
+      h->Fill(mMc,mEsd);
+    }
+  }
+}
+
+//____________________________________________________________________
 Bool_t
-AliFMDMCDensityCalculator::Calculate(const AliESDFMD&,
-                                    AliForwardUtil::Histos&,
-                                    UShort_t,
-                                    Bool_t)
+AliFMDMCDensityCalculator::CompareResults(AliForwardUtil::Histos& esd,
+                                         AliForwardUtil::Histos& mc)
 {
-  AliWarning("Method Calculate disabled for this class. If you need this, "
-            "make an AliFMDDensityCalculator object instead");
-  return kFALSE;
+  Fill(1, 'I', esd.Get(1,'I'), mc.Get(1,'I'));
+  Fill(2, 'I', esd.Get(2,'I'), mc.Get(2,'I'));
+  Fill(2, 'O', esd.Get(2,'O'), mc.Get(2,'O'));
+  Fill(3, 'I', esd.Get(3,'I'), mc.Get(3,'I'));
+  Fill(3, 'O', esd.Get(3,'O'), mc.Get(3,'O'));
+
+  return kTRUE;
 }
 
 //____________________________________________________________________
+void
+AliFMDMCDensityCalculator::DefineOutput(TList* dir)
+{
+  AliFMDDensityCalculator::DefineOutput(dir);
+  TList* d = static_cast<TList*>(dir->FindObject(GetName()));
+
+  fComps = new TList;
+  fComps->SetName("esd_mc_comparison");
+  d->Add(fComps);
+
+}
+//____________________________________________________________________
 //
 // EOF
 //
index 6f39f9b..55d9794 100644 (file)
@@ -1,11 +1,12 @@
-#ifndef ALIROOT_PWG2_FORWARD_ANALYSIS2_ALIFMDMCDENSITYCALCULATOR_H
-#define ALIROOT_PWG2_FORWARD_ANALYSIS2_ALIFMDMCDENSITYCALCULATOR_H
+#ifndef ALIFMDMCDENSITYCALCULATOR_H
+#define ALIFMDMCDENSITYCALCULATOR_H
 #include "AliFMDDensityCalculator.h"
 #include <TList.h>
 #include "AliForwardUtil.h"
 class AliMCEvent;
+class TH2;
 class TH2D;
-class TH1D;
+class TProfile2D;
 
 /** 
  * This class calculates the inclusive charged particle density
@@ -20,7 +21,8 @@ class TH1D;
  * @par Corrections used: 
  *   - None
  *
- * @ingroup pwg2_forward 
+ * @ingroup pwg2_forward_algo
+ * @ingroup pwg2_forward_mc
  */
 class AliFMDMCDensityCalculator : public AliFMDDensityCalculator
 {
@@ -28,14 +30,38 @@ public:
   /** 
    * Constructor 
    */
-  AliFMDMCDensityCalculator() : AliFMDDensityCalculator() {}
+  AliFMDMCDensityCalculator() 
+    : AliFMDDensityCalculator(),
+      fFMD1i(0), 
+      fFMD2i(0),
+      fFMD2o(0),
+      fFMD3i(0),
+      fFMD3o(0),
+      fFMD1iC(0), 
+      fFMD2iC(0),
+      fFMD2oC(0),
+      fFMD3iC(0),
+      fFMD3oC(0),
+      fComps(0)
+  {}
   /** 
    * Constructor 
    * 
    * @param name Name of object
    */
   AliFMDMCDensityCalculator(const char* name) 
-   : AliFMDDensityCalculator(name)
+   : AliFMDDensityCalculator(name),
+      fFMD1i(0), 
+      fFMD2i(0),
+      fFMD2o(0),
+      fFMD3i(0),
+      fFMD3o(0),
+      fFMD1iC(0), 
+      fFMD2iC(0),
+      fFMD2oC(0),
+      fFMD3iC(0),
+      fFMD3oC(0),
+      fComps(0)
   {}
   /** 
    * Copy constructor 
@@ -43,12 +69,23 @@ public:
    * @param o Object to copy from 
    */
   AliFMDMCDensityCalculator(const AliFMDMCDensityCalculator& o)
-   : AliFMDDensityCalculator(o) 
+   : AliFMDDensityCalculator(o) ,
+      fFMD1i(o.fFMD1i), 
+      fFMD2i(o.fFMD2i),
+      fFMD2o(o.fFMD2o),
+      fFMD3i(o.fFMD3i),
+      fFMD3o(o.fFMD3o),
+      fFMD1iC(o.fFMD1iC), 
+      fFMD2iC(o.fFMD2iC),
+      fFMD2oC(o.fFMD2oC),
+      fFMD3iC(o.fFMD3iC),
+      fFMD3oC(o.fFMD3oC),
+      fComps(0)
   {}
   /** 
    * Destructor 
    */
-  virtual ~AliFMDMCDensityCalculator() {}
+  virtual ~AliFMDMCDensityCalculator();
   /** 
    * Assignement operator
    * 
@@ -58,18 +95,11 @@ public:
    */
   AliFMDMCDensityCalculator& operator=(const AliFMDMCDensityCalculator& o);
   /** 
-   * Do the calculations 
-   * 
-   * @param fmd      AliESDFMD object (possibly) corrected for sharing
-   * @param hists    Histogram cache
-   * @param vtxBin   Vertex bin 
-   * @param lowFlux  Low flux flag. 
+   * Initialize this object 
    * 
-   * @return true on successs 
+   * @param etaAxis Eta axis to use 
    */
-  virtual Bool_t Calculate(const AliESDFMD& fmd, 
-                          AliForwardUtil::Histos& hists, 
-                          UShort_t vtxBin, Bool_t lowFlux);
+  void Init(const TAxis& etaAxis);
   /** 
    * Calculate the charged particle density from the MC track references. 
    * 
@@ -80,11 +110,68 @@ public:
    * 
    * @return true on success
    */
-  virtual Bool_t CalculateMC(const AliMCEvent& event, 
-                            AliForwardUtil::Histos& hists, 
-                            Double_t vz,
-                            UShort_t vtxBin);                       
+  virtual Bool_t CalculateMC(const AliESDFMD&        fmd,
+                            AliForwardUtil::Histos& hists);
+                            
+  /** 
+   * Compare the result of analysing the ESD for 
+   * the inclusive charged particle density to analysing 
+   * MC truth 
+   * 
+   * @param esd 
+   * @param mc 
+   * 
+   * @return 
+   */
+  virtual Bool_t CompareResults(AliForwardUtil::Histos& esd, 
+                               AliForwardUtil::Histos& mc);
+  /** 
+   * Output diagnostic histograms to directory 
+   * 
+   * @param dir List to write in
+   */  
+  void DefineOutput(TList* dir);
 protected:
+  /** 
+   * MAke comparison profiles
+   * 
+   * @param d     Detector 
+   * @param r     Ring 
+   * @param axis  Eta axis 
+   * 
+   * @return Newly allocated profile object
+   */
+  TProfile2D* Make(UShort_t d, Char_t r, const TAxis& axis) const;
+  /** 
+   * MAke comparison profiles
+   * 
+   * @param d     Detector 
+   * @param r     Ring 
+   * 
+   * @return Newly allocated profile object
+   */
+  TH2D* Make(UShort_t d, Char_t r) const;
+  /** 
+   * Fill comparison profiles
+   * 
+   * @param d    Detector 
+   * @param r    Ring 
+   * @param esd  ESD histogram
+   * @param mc   MC histogram
+   */
+  void Fill(UShort_t d, Char_t r, TH2* esd, TH2* mc);
+
+  TProfile2D* fFMD1i; // Comparison
+  TProfile2D* fFMD2i; // Comparison
+  TProfile2D* fFMD2o; // Comparison
+  TProfile2D* fFMD3i; // Comparison
+  TProfile2D* fFMD3o; // Comparison
+  TH2D*       fFMD1iC; // Correlation in FMD1i
+  TH2D*       fFMD2iC; // Correlation in FMD2i
+  TH2D*       fFMD2oC; // Correlation in FMD2o
+  TH2D*       fFMD3iC; // Correlation in FMD3i
+  TH2D*       fFMD3oC; // Correlation in FMD3o
+  TList*      fComps; // List of comparisons 
 
   ClassDef(AliFMDMCDensityCalculator,1); // Calculate Nch density 
 };