FMD/scripts/PoissonHit.C

   1 //____________________________________________________________________
   2 //
   3 // $Id$
   4 //
   5 // Script that contains a class to draw hits, using the
   6 // AliFMDInputHits class in the util library.
   7 //
   8 // It draws the energy loss versus the p/(mq^2).  It can be overlayed
   9 // with the Bethe-Bloc curve to show how the simulation behaves
  10 // relative to the expected.
  11 //
  12 // Use the script `Compile.C' to compile this class using ACLic.
  13 //
  14 #include "Poisson.C"
  15 #include <TMath.h>
  16 #include <TCanvas.h>
  17 #include <AliFMDHit.h>
  18 #include <AliFMDGeometry.h>
  19
  20 /** @class PoissonHit
  21     @brief Make a poisson reconstruction and compare to simulated hits
  22     @code
  23     Root> .L Compile.C
  24     Root> Compile("Poisson.C")
  25     Root> Compile("PoissonHit.C")
  26     Root> PoissonHit c
  27     Root> c.Run();
  28     @endcode
  29     @ingroup FMD_script
  30  */
  31 class PoissonHit : public Poisson
  32 {
  33 protected:
  34   TH2D*  fHits; // Histogram
  35   TH2D*  fDiff; // Histogram
  36 public:
  37   /** Constructor
  38       @param threshold Threshold
  39       @param nEta      # of @f$ \eta@f$ bins
  40       @param minEta    minimum @f$ \eta@f$
  41       @param maxEta    maximum @f$ \eta@f$
  42       @param nPhi      # of @f$ \eta@f$ bins
  43       @param minPhi    minimum @f$ \varphi@f$
  44       @param maxPhi    maximum @f$ \varphi@f$ */
  45   PoissonHit(Double_t threshold=.3,
  46              Int_t nEta=120, Float_t minEta=-6, Float_t maxEta=6,
  47              Int_t nPhi=4,   Float_t minPhi=0,  Float_t maxPhi=2*TMath::Pi())
  48     : Poisson(threshold, nEta, minEta, maxEta, nPhi, minPhi, maxPhi)
  49   {
  50     AddLoad(kHits);
  51     AddLoad(kGeometry);
  52     fHits  = new TH2D(*fEmpty);
  53     fHits->SetName("hits");
  54     fHits->SetTitle("# of hits");
  55     fDiff  = new TH2D(*fEmpty);
  56     fDiff->SetName("diff");
  57     fDiff->SetTitle("Difference between poisson and hits");
  58     fHits->SetXTitle("#eta");
  59     fHits->SetYTitle("#phi");
  60     fHits->SetZTitle("N");
  61     fDiff->SetXTitle("#eta");
  62     fDiff->SetYTitle("#phi");
  63     fDiff->SetZTitle("#frac{N_{hit}-N_{poisson}}{N_{hit}}");
  64   }
  65   /** Initialize the analyser. Opens the output file.
  66       @return @c true on success. */
  67   virtual Bool_t Init()
  68   {
  69     if (!Poisson::Init()) return kFALSE;
  70     AliFMDGeometry::Instance()->Init();
  71     AliFMDGeometry::Instance()->InitTransformations();
  72     return kTRUE;
  73   }
  74   /** Get the @f$ \eta@f$ and @f$\varphi@f$ corresponding to the
  75       spatial coordinates @f$ \mathbf{v} = (x,y,z)@f$
  76       @param x    X coordinate
  77       @param y    Y coordinate
  78       @param z    Z coordinate
  79       @param eta  Psuedo rapidity @f$ \eta@f$
  80       @param phi  Azimuthal angle @f$\varphi@f$
  81   */
  82   void PhysicalCoordinates(Double_t x, Double_t y, Double_t z,
  83                            Double_t& eta, Double_t& phi)
  84   {
  85     Double_t r, theta;
  86     phi   =  TMath::ATan2(y, x);
  87     r     =  TMath::Sqrt(y * y + x * x);
  88     theta =  TMath::ATan2(r, z);
  89     eta   = -TMath::Log(TMath::Tan(theta / 2));
  90     if (phi < 0) phi += 2 * TMath::Pi();
  91   }
  92   /** Process one hit.  Increment bin corresponding to strip.
  93       @param hit Hit.
  94       @return @c true on success */
  95   virtual Bool_t ProcessHit(AliFMDHit* hit, TParticle*)
  96   {
  97     Double_t x, y, z;
  98 #if 0
  99     AliFMDGeometry* geom = AliFMDGeometry::Instance();
 100     geom->Detector2XYZ(hit->Detector(),hit->Ring(),hit->Sector(),
 101                        hit->Strip(),x,y,z);
 102 #else
 103     x = hit->X();
 104     y = hit->Y();
 105     z = hit->Z();
 106 #endif
 107     Double_t eta, phi;
 108     PhysicalCoordinates(x, y, z, eta, phi);
 109     fHits->Fill(eta, phi);
 110     return kTRUE;
 111   }
 112   /** Begining of event
 113       @param event Event number
 114       @return @c false on error */
 115   virtual Bool_t Begin(Int_t event)
 116   {
 117     if (!Poisson::Begin(event)) return kFALSE;
 118     fHits->Clear();
 119     fDiff->Clear();
 120     return kTRUE;
 121   }
 122   /** Let the poisson code do it's job, and then compare to the
 123       numberr of hits.
 124       @return @c true  */
 125   virtual Bool_t End()
 126   {
 127     if (!Poisson::End()) return kFALSE;
 128     fDiff->Add(fMult,fHits,-1.,1.);
 129     fDiff->Divide(fHits);
 130     if (!gROOT->IsBatch()) {
 131       gStyle->SetPalette(1);
 132       TCanvas* c1 = new TCanvas("hits", "Hit multiplicity");
 133       c1->SetFillColor(0);
 134       fHits->Draw("colz");
 135       TCanvas* c2 = new TCanvas("diff", "Difference between Hit and poisson");
 136       c2->SetFillColor(0);
 137       fDiff->Draw("colz");
 138       TCanvas* c3 = new TCanvas("empty", "# of Empty strips");
 139       c3->SetFillColor(0);
 140       fEmpty->Draw("colz");
 141       TCanvas* c4 = new TCanvas("total", "Total # of strips");
 142       c4->SetFillColor(0);
 143       fTotal->Draw("colz");
 144     }
 145     return kTRUE;
 146   }
 147
 148   ClassDef(PoissonHit,0);
 149 };
 150
 151 //____________________________________________________________________
 152 //
 153 // EOF
 154 //