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 <AliFMDDigit.h>
  17 #include <AliFMDRecPoint.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 //____________________________________________________________________
  32 /** @class DrawHitsRecs
  33     @brief Draw hit energy loss versus rec point mult
  34     @code
  35     Root> .L Compile.C
  36     Root> Compile("DrawHitsRecs.C")
  37     Root> DrawHitsRecs c
  38     Root> c.Run();
  39     @endcode
  40     @ingroup FMD_script
  41  */
  42 class DrawHitsRecs : public AliFMDInput
  43 {
  44 private:
  45   TH2D* fHitEvsAdc;
  46   TH2D* fHitEvsRecM;  // Histogram
  47   TH2D* fHitEvsRecE;  // Histogram
  48   TH1D* fDiffE;       // Histogram
  49   TH2D* fHitsVsRecM;  // Histogram
  50   AliFMDEdepMap  fMap;
  51   AliFMDFloatMap fEta;
  52   AliFMDFloatMap fPhi;
  53   AliFMDFloatMap fMult;
  54   Bool_t fPrimary;
  55 public:
  56   //__________________________________________________________________
  57   TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max)
  58   {
  59     TArrayF bins(n+1);
  60     Float_t dp   = n / TMath::Log10(max / min);
  61     Float_t pmin = TMath::Log10(min);
  62     bins[0]      = min;
  63     for (Int_t i = 1; i < n+1; i++) {
  64       Float_t p = pmin + i / dp;
  65       bins[i]   = TMath::Power(10, p);
  66     }
  67     return bins;
  68   }
  69   //__________________________________________________________________
  70   DrawHitsRecs(Bool_t primary=kFALSE,
  71                Int_t n=900, Double_t emin=1e-3, Double_t emax=10,
  72                Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5)
  73   {
  74     fPrimary = primary;
  75     AddLoad(kRecPoints);
  76     AddLoad(kHits);
  77     AddLoad(kDigits);
  78     if (fPrimary) AddLoad(kKinematics);
  79     TArrayF eloss(MakeLogScale(n, emin, emax));
  80     TArrayF mults(m+1);
  81     mults[0] = mmin;
  82     for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;
  83
  84     fHitEvsAdc  = new TH2D("hitEvsAdc", "#Delta E_{sim} vs. ADC",
  85                            n, emin, emax, 1025, -.5, 1024.5);
  86     fHitEvsAdc->SetXTitle("#Delta E_{sim} [MeV]");
  87     fHitEvsAdc->SetYTitle("ADC");
  88
  89     fHitEvsRecM = new TH2D("hitEvsRecM", "#Delta E_{sim} vs. M_{rec}",
  90                            eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray);
  91     fHitEvsRecM->SetXTitle("#Delta E_{sim} [MeV]");
  92     fHitEvsRecM->SetYTitle("M_{rec}");
  93
  94     fHitEvsRecE = new TH2D("hitEvsRecE", "#Delta E_{sim} vs. #Delta E_{rec}",
  95                             n, emin, emax, n, emin, emax);
  96     fHitEvsRecE->SetXTitle("#Delta E_{sim} [MeV]");
  97     fHitEvsRecE->SetYTitle("#Delta E_{rec} [MeV]");
  98
  99
 100     fDiffE = new TH1D("diffE",
 101                       "#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}",
 102                       1100, -1, 1.1);
 103     fDiffE->SetXTitle("#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}");
 104
 105     Double_t omin = -.5;
 106     Double_t omax = 7.5;
 107     Int_t    o    = 8;
 108     fHitsVsRecM = new TH2D("hitsVsStrip", "# of Hits vs. M_{rec}",
 109                            o, omin, omax, m, mmin, mmax);
 110     fHitsVsRecM->SetXTitle("# of Hits");
 111     fHitsVsRecM->SetYTitle("M_{rec}");
 112   }
 113   //__________________________________________________________________
 114   /** Begining of event
 115       @param ev Event number
 116       @return @c false on error */
 117   Bool_t Begin(Int_t ev)
 118   {
 119     fMap.Reset();
 120     return AliFMDInput::Begin(ev);
 121   }
 122   //__________________________________________________________________
 123   Bool_t ProcessHit(AliFMDHit* hit, TParticle*)
 124   {
 125     // Cache the energy loss
 126     if (!hit) return kFALSE;
 127     UShort_t det = hit->Detector();
 128     Char_t   rng = hit->Ring();
 129     UShort_t sec = hit->Sector();
 130     UShort_t str = hit->Strip();
 131     if (str > 511) {
 132       AliWarning(Form("Bad strip number %d in hit", str));
 133       return kTRUE;
 134     }
 135     if (fPrimary) {
 136       TParticle* kine = fStack->Particle(hit->Track());
 137       if (!kine) return kTRUE;
 138       if (!kine->IsPrimary()) return kTRUE;
 139     }
 140
 141     fMap(det, rng, sec, str).fEdep += hit->Edep();
 142     fMap(det, rng, sec, str).fN++;
 143     return kTRUE;
 144   }
 145
 146   //__________________________________________________________________
 147   Bool_t ProcessDigit(AliFMDDigit* digit)
 148   {
 149     if (!digit) return kTRUE;
 150
 151     UShort_t det = digit->Detector();
 152     Char_t   rng = digit->Ring();
 153     UShort_t sec = digit->Sector();
 154     UShort_t str = digit->Strip();
 155     if (str > 511) {
 156       AliWarning(Form("Bad strip number %d in digit", str));
 157       return kTRUE;
 158     }
 159     Double_t edep = fMap(det, rng, sec, str).fEdep;
 160     if (edep > 0) fHitEvsAdc->Fill(edep, digit->Counts());
 161
 162     return kTRUE;
 163   }
 164
 165   //__________________________________________________________________
 166   Bool_t ProcessRecPoint(AliFMDRecPoint* single)
 167   {
 168     if (!single) return kTRUE;
 169     UShort_t det = single->Detector();
 170     Char_t   rng = single->Ring();
 171     UShort_t sec = single->Sector();
 172     UShort_t str = single->Strip();
 173     if (str > 511) {
 174       AliWarning(Form("Bad strip number %d in single", str));
 175       return kTRUE;
 176     }
 177     Double_t edep = fMap(det, rng, sec, str).fEdep;
 178     Int_t    nhit = fMap(det, rng, sec, str).fN;
 179     if (edep > 0) {
 180       fHitEvsRecM->Fill(edep, single->Particles());
 181       fHitEvsRecE->Fill(edep, single->Edep());
 182       fDiffE->Fill((single->Edep() - edep) / edep);
 183     }
 184     if (nhit > 0) fHitsVsRecM->Fill(nhit, single->Particles());
 185     return kTRUE;
 186   }
 187   //__________________________________________________________________
 188   Bool_t Finish()
 189   {
 190     gStyle->SetPalette(1);
 191     gStyle->SetOptTitle(0);
 192     gStyle->SetCanvasColor(0);
 193     gStyle->SetCanvasBorderSize(0);
 194     gStyle->SetPadColor(0);
 195     gStyle->SetPadBorderSize(0);
 196
 197     new TCanvas("c0", fHitEvsAdc->GetTitle());
 198     fHitEvsAdc->SetStats(kFALSE);
 199     fHitEvsAdc->Draw("COLZ");
 200
 201     new TCanvas("c1", fHitEvsRecM->GetTitle());
 202     fHitEvsRecM->SetStats(kFALSE);
 203     fHitEvsRecM->Draw("COLZ");
 204
 205     new TCanvas("c2", fHitEvsRecE->GetTitle());
 206     fHitEvsRecE->SetStats(kFALSE);
 207     fHitEvsRecE->Draw("COLZ");
 208
 209     new TCanvas("c3", fDiffE->GetTitle());
 210     fDiffE->Draw();
 211
 212     new TCanvas("c4", fHitsVsRecM->GetTitle());
 213     fHitsVsRecM->SetStats(kFALSE);
 214     fHitsVsRecM->Draw("COLZ");
 215
 216     return kTRUE;
 217   }
 218
 219   ClassDef(DrawHitsRecs,0);
 220 };
 221
 222 //____________________________________________________________________
 223 //
 224 // EOF
 225 //