FMD/scripts/DrawHits.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 <TH2D.h>
  15 #include <AliFMDHit.h>
  16 #include <AliFMDInput.h>
  17 #include <iostream>
  18 #include <TStyle.h>
  19 #include <TArrayF.h>
  20 #include <TParticle.h>
  21
  22 /** @class DrawHits
  23     @brief Draw hit energy loss
  24     @code
  25     Root> .L Compile.C
  26     Root> Compile("DrawHits.C")
  27     Root> DrawHits c
  28     Root> c.Run();
  29     @endcode
  30     @ingroup FMD_script
  31  */
  32 class DrawHits : public AliFMDInput
  33 {
  34 private:
  35   TH2D* fElossVsPMQ; // Histogram
  36   Int_t fPdg;
  37 public:
  38   //__________________________________________________________________
  39   TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max)
  40   {
  41     TArrayF bins(n+1);
  42     bins[0]      = min;
  43     if (n <= 20) {
  44       for (Int_t i = 1; i < n+1; i++) bins[i] = bins[i-1] + (max-min)/n;
  45       return bins;
  46     }
  47     Float_t dp   = n / TMath::Log10(max / min);
  48     Float_t pmin = TMath::Log10(min);
  49     for (Int_t i = 1; i < n+1; i++) {
  50       Float_t p = pmin + i / dp;
  51       bins[i]   = TMath::Power(10, p);
  52     }
  53     return bins;
  54   }
  55   //__________________________________________________________________
  56   DrawHits(Int_t m=1000, Double_t emin=1, Double_t emax=1000,
  57            Int_t n=900, Double_t tmin=1e-2, Double_t tmax=1e3,
  58            Int_t pdg=211)
  59     : fPdg(pdg)
  60   {
  61     AddLoad(kKinematics);
  62     AddLoad(kHits);
  63     TArrayF tkine(MakeLogScale(n, tmin, tmax));
  64     TArrayF eloss(MakeLogScale(m, emin, emax));
  65     TString name(Form("elossVsP%s", (fPdg == 0 ? "MQ" : "")));
  66     TString title(Form("#Delta E/#Delta x vs. p%s",
  67                        (fPdg == 0 ? "/(mq^{2})" : "")));
  68     fElossVsPMQ = new TH2D(name.Data(), title.Data(),
  69                            tkine.fN-1, tkine.fArray,
  70                            eloss.fN-1, eloss.fArray);
  71     fElossVsPMQ->SetXTitle(Form("p%s", (fPdg == 0 ? "/(mq^{2})" : " [GeV]")));
  72     fElossVsPMQ->SetYTitle("#Delta E/#Delta x [MeV/cm]");
  73   }
  74   //__________________________________________________________________
  75   Bool_t ProcessHit(AliFMDHit* hit, TParticle* p)
  76   {
  77     if (!hit) {
  78       std::cout << "No hit" << std::endl;
  79       return kFALSE;
  80     }
  81
  82     if (!p) {
  83       std::cout << "No track" << std::endl;
  84       return kFALSE;
  85     }
  86     if (!p->IsPrimary()) return kTRUE;
  87     if (hit->IsStop()) return kTRUE;
  88     Float_t x = hit->P();
  89     Float_t y = hit->Edep()/hit->Length();
  90     if (fPdg != 0) {
  91       if (TMath::Abs(hit->Pdg()) != fPdg) return kTRUE;
  92     }
  93     else {
  94       Float_t q = hit->Q() / 3.;
  95       if (hit->M() == 0 || q == 0) return kTRUE;
  96       x /= hit->M();
  97       y /= q * q;
  98     }
  99     fElossVsPMQ->Fill(x, y);
 100     return kTRUE;
 101   }
 102   //__________________________________________________________________
 103   Bool_t Finish()
 104   {
 105     gStyle->SetPalette(1);
 106     gStyle->SetOptTitle(0);
 107     gStyle->SetCanvasColor(0);
 108     gStyle->SetCanvasBorderSize(0);
 109     gStyle->SetPadColor(0);
 110     gStyle->SetPadBorderSize(0);
 111     fElossVsPMQ->SetStats(kFALSE);
 112     fElossVsPMQ->Draw("COLZ box");
 113     return kTRUE;
 114   }
 115   ClassDef(DrawHits,0);
 116 };
 117
 118 //____________________________________________________________________
 119 //
 120 // EOF
 121 //