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