FMD/scripts/DrawHitsRecs.C

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