FMD/scripts/DrawHitsRecs.C

   1 //____________________________________________________________________
   2 //
   3 // $Id$
   4 //
   5 // Script that contains a class to draw eloss from hits, versus ADC
   6 // counts from digits, using the 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 <AliFMDMultStrip.h>
  17 #include <AliFMDMultRegion.h>
  18 #include <AliFMDDigit.h>
  19 #include <AliFMDInput.h>
  20 #include <AliFMDEdepMap.h>
  21 #include <AliFMDFloatMap.h>
  22 #include <iostream>
  23 #include <TStyle.h>
  24 #include <TArrayF.h>
  25 #include <TCanvas.h>
  26 #include <TParticle.h>
  27 #include <TClonesArray.h>
  28 #include <TTree.h>
  29 #include <AliStack.h>
  30 #include <AliLog.h>
  31
  32 //____________________________________________________________________
  33 class DrawHitsRecs : public AliFMDInputHits
  34 {
  35 private:
  36   TH2D* fElossVsMult; // Histogram
  37   TH2D* fHitsVsStrip;  // Histogram
  38   TH2D* fHitsVsRegion;  // Histogram
  39   AliFMDEdepMap fMap;
  40   AliFMDFloatMap fEta;
  41   AliFMDFloatMap fPhi;
  42   AliFMDFloatMap fMult;
  43   Bool_t fPrimary;
  44 public:
  45   //__________________________________________________________________
  46   TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max)
  47   {
  48     TArrayF bins(n+1);
  49     Float_t dp   = n / TMath::Log10(max / min);
  50     Float_t pmin = TMath::Log10(min);
  51     bins[0]      = min;
  52     for (Int_t i = 1; i < n+1; i++) {
  53       Float_t p = pmin + i / dp;
  54       bins[i]   = TMath::Power(10, p);
  55     }
  56     return bins;
  57   }
  58   //__________________________________________________________________
  59   DrawHitsRecs(Bool_t primary=kFALSE,
  60                Int_t n=900, Double_t emin=1e-3, Double_t emax=10,
  61                Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5)
  62   {
  63     fPrimary = primary;
  64     AddLoad(kRecPoints);
  65     if (fPrimary) AddLoad(kKinematics);
  66     TArrayF eloss(MakeLogScale(n, emin, emax));
  67     TArrayF mults(m+1);
  68     mults[0] = mmin;
  69     for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;
  70     fElossVsMult = new TH2D("elossVsMult",
  71                             "#Delta E vs. Multiplicity (single)",
  72                            eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray);
  73     fElossVsMult->SetXTitle("#Delta E/#Delta x [MeV/cm]");
  74     fElossVsMult->SetYTitle("Strip Multiplicity");
  75
  76     Double_t omin = -.5;
  77     Double_t omax = 7.5;
  78     Int_t    o    = 8;
  79     fHitsVsStrip = new TH2D("hitsVsStrip",
  80                             "# of Hits vs. Multiplicity (strip)",
  81                            o, omin, omax, m, mmin, mmax);
  82     fHitsVsStrip->SetXTitle("# of Hits");
  83     fHitsVsStrip->SetYTitle("Strip Multiplicity");
  84   }
  85   //__________________________________________________________________
  86   Bool_t Begin(Int_t ev)
  87   {
  88     fMap.Reset();
  89     return AliFMDInputHits::Begin(ev);
  90   }
  91   //__________________________________________________________________
  92   Bool_t ProcessHit(AliFMDHit* hit, TParticle*)
  93   {
  94     // Cache the energy loss
  95     if (!hit) return kFALSE;
  96     UShort_t det = hit->Detector();
  97     Char_t   rng = hit->Ring();
  98     UShort_t sec = hit->Sector();
  99     UShort_t str = hit->Strip();
 100     if (str > 511) {
 101       AliWarning(Form("Bad strip number %d in hit", str));
 102       return kTRUE;
 103     }
 104     if (fPrimary) {
 105       TParticle* kine = fStack->Particle(hit->Track());
 106       if (!kine) return kTRUE;
 107       if (!kine->IsPrimary()) return kTRUE;
 108     }
 109
 110     fMap(det, rng, sec, str).fEdep += hit->Edep();
 111     fMap(det, rng, sec, str).fN++;
 112     return kTRUE;
 113   }
 114   //__________________________________________________________________
 115   Bool_t ProcessRecPoint(AliFMDRecPoint* single)
 116   {
 117     if (!single) continue;
 118     UShort_t det = single->Detector();
 119     Char_t   rng = single->Ring();
 120     UShort_t sec = single->Sector();
 121     UShort_t str = single->Strip();
 122     if (str > 511) {
 123       AliWarning(Form("Bad strip number %d in single", str));
 124       continue;
 125     }
 126     if (fMap(det, rng, sec, str).fEdep > 0)
 127       fElossVsMult->Fill(fMap(det, rng, sec, str).fEdep,
 128                          single->Particles());
 129     if (fMap(det, rng, sec, str).fN > 0)
 130       fHitsVsStrip->Fill(fMap(det, rng, sec, str).fN,
 131                          single->Particles());
 132     return kTRUE;
 133   }
 134   //__________________________________________________________________
 135   Bool_t Finish()
 136   {
 137     gStyle->SetPalette(1);
 138     gStyle->SetOptTitle(0);
 139     gStyle->SetCanvasColor(0);
 140     gStyle->SetCanvasBorderSize(0);
 141     gStyle->SetPadColor(0);
 142     gStyle->SetPadBorderSize(0);
 143
 144     new TCanvas("c1", "Energy loss vs. Strip Multiplicity");
 145     fElossVsMult->SetStats(kFALSE);
 146     fElossVsMult->Draw("COLZ");
 147
 148     new TCanvas("c2", "# of Hits vs. Strip Multiplicity");
 149     fHitsVsStrip->SetStats(kFALSE);
 150     fHitsVsStrip->Draw("COLZ");
 151
 152     return kTRUE;
 153   }
 154
 155   ClassDef(DrawHitsRecs,0);
 156 };
 157
 158 //____________________________________________________________________
 159 //
 160 // EOF
 161 //