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 #include <TF1.h>
  31
  32 //____________________________________________________________________
  33 /** @class DrawHitsRecs
  34     @brief Draw hit energy loss versus rec point mult
  35     @code
  36     Root> .L Compile.C
  37     Root> Compile("DrawHitsRecs.C")
  38     Root> DrawHitsRecs c
  39     Root> c.Run();
  40     @endcode
  41     @ingroup FMD_script
  42  */
  43 class DrawHitsRecs : public AliFMDInput
  44 {
  45 private:
  46   TH2D* fHitEvsAdc;
  47   TH2D* fHitEvsRecM;  // Histogram
  48   TH2D* fHitEvsRecE;  // Histogram
  49   TH1D* fDiffE;       // Histogram
  50   TH2D* fHitsVsRecM;  // Histogram
  51   TH2D* fDiffM;       // Histogram
  52   TH1*  fHitEloss;
  53   TH1*  fRecEloss;
  54   AliFMDEdepMap  fMap;
  55   AliFMDFloatMap fEta;
  56   AliFMDFloatMap fPhi;
  57   AliFMDFloatMap fMult;
  58   Bool_t fPrimary;
  59 public:
  60   //__________________________________________________________________
  61   DrawHitsRecs(Bool_t primary=kFALSE,
  62                Int_t n=900, Double_t emin=1e-3, Double_t emax=10,
  63                Int_t m=21, Double_t mmin=-0.5, Double_t mmax=20.5)
  64   {
  65     fPrimary = primary;
  66     AddLoad(kRecPoints);
  67     AddLoad(kHits);
  68     AddLoad(kDigits);
  69     if (fPrimary) AddLoad(kKinematics);
  70     TArrayF eloss(MakeLogScale(n, emin, emax));
  71     TArrayF mults(m+1);
  72     mults[0] = mmin;
  73     for (Int_t i = 1; i < m+1; i++) mults[i] = mults[i-1] + (mmax-mmin)/m;
  74
  75     fHitEvsAdc  = new TH2D("hitEvsAdc", "#Delta E_{sim} vs. ADC",
  76                            n, emin, emax, 1025, -.5, 1024.5);
  77     fHitEvsAdc->SetXTitle("#Delta E_{sim} [MeV]");
  78     fHitEvsAdc->SetYTitle("ADC");
  79     fHitEvsAdc->Sumw2();
  80
  81     fHitEvsRecM = new TH2D("hitEvsRecM", "#Delta E_{sim} vs. M_{rec}",
  82                            eloss.fN-1, eloss.fArray, mults.fN-1, mults.fArray);
  83     fHitEvsRecM->SetXTitle("#Delta E_{sim} [MeV]");
  84     fHitEvsRecM->SetYTitle("M_{rec}");
  85     fHitEvsRecM->Sumw2();
  86
  87     fHitEvsRecE = new TH2D("hitEvsRecE", "#Delta E_{sim} vs. #Delta E_{rec}",
  88                             n, emin, emax, n, emin, emax);
  89     fHitEvsRecE->SetXTitle("#Delta E_{sim} [MeV]");
  90     fHitEvsRecE->SetYTitle("#Delta E_{rec} [MeV]");
  91     fHitEvsRecE->Sumw2();
  92
  93
  94     fDiffE = new TH1D("diffE",
  95                       "#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}",
  96                       1100, -1, 1.1);
  97     fDiffE->SetXTitle("#frac{#Delta E_{sim}-#Delta E_{rec}}{#Delta E_{sim}}");
  98     fDiffE->Sumw2();
  99
 100     Double_t omin = mmin; // -.5;
 101     Double_t omax = mmax; // 7.5;
 102     Int_t    o    = m;    // 8;
 103     fHitsVsRecM = new TH2D("hitsVsStrip", "# of Hits vs. M_{rec}",
 104                            o, omin, omax, m, mmin, mmax);
 105     fHitsVsRecM->SetXTitle("# of Hits");
 106     fHitsVsRecM->SetYTitle("M_{rec}");
 107     fHitsVsRecM->Sumw2();
 108
 109     fDiffM = new TH2D("diffM", "M_{sim} - M_{rec}",
 110                       41, -20.5, 20.5, 70, 1.5, 5);
 111     // 36, -TMath::Pi(),TMath::Pi());
 112     fDiffM->SetXTitle("M_{sim} - M_{rec}");
 113     fDiffM->SetYTitle("|#eta|");
 114     fDiffM->Sumw2();
 115     // fDiffM->SetYTitle("Detector");
 116
 117     fHitEloss = new TH1D("hitEloss","#frac{#Delta E_{sim}}{#Delta x} (MeV/cm)",
 118                          100, 0, 10);
 119     fHitEloss->SetFillColor(2);
 120     fHitEloss->SetFillStyle(3001);
 121     fHitEloss->Sumw2();
 122
 123     fRecEloss = new TH1D("recEloss","#frac{#Delta E_{rec}}{#Delta x} (MeV/cm)",
 124                          100, 0, 10);
 125     fRecEloss->SetFillColor(4);
 126     fRecEloss->SetFillStyle(3001);
 127     fRecEloss->Sumw2();
 128   }
 129   //__________________________________________________________________
 130   /** Begining of event
 131       @param ev Event number
 132       @return @c false on error */
 133   Bool_t Begin(Int_t ev)
 134   {
 135     fMap.Reset();
 136     return AliFMDInput::Begin(ev);
 137   }
 138   //__________________________________________________________________
 139   Bool_t ProcessHit(AliFMDHit* hit, TParticle*)
 140   {
 141     // Cache the energy loss
 142     if (!hit) return kFALSE;
 143     UShort_t det = hit->Detector();
 144     Char_t   rng = hit->Ring();
 145     UShort_t sec = hit->Sector();
 146     UShort_t str = hit->Strip();
 147     if (str > 511) {
 148       AliWarning(Form("Bad strip number %d in hit", str));
 149       return kTRUE;
 150     }
 151     if (fPrimary) {
 152       TParticle* kine = fStack->Particle(hit->Track());
 153       if (!kine) return kTRUE;
 154       if (!kine->IsPrimary()) return kTRUE;
 155     }
 156
 157     if (hit->Edep()/hit->Length() > 0.1)
 158       fHitEloss->Fill(hit->Edep() / hit->Length());
 159     fMap(det, rng, sec, str).fEdep += hit->Edep();
 160     fMap(det, rng, sec, str).fN++;
 161     return kTRUE;
 162   }
 163
 164   //__________________________________________________________________
 165   Bool_t ProcessDigit(AliFMDDigit* digit)
 166   {
 167     if (!digit) return kTRUE;
 168
 169     UShort_t det = digit->Detector();
 170     Char_t   rng = digit->Ring();
 171     UShort_t sec = digit->Sector();
 172     UShort_t str = digit->Strip();
 173     if (str > 511) {
 174       AliWarning(Form("Bad strip number %d in digit", str));
 175       return kTRUE;
 176     }
 177     Double_t edep = fMap(det, rng, sec, str).fEdep;
 178     if (edep > 0) fHitEvsAdc->Fill(edep, digit->Counts());
 179
 180     return kTRUE;
 181   }
 182
 183   //__________________________________________________________________
 184   Bool_t ProcessRecPoint(AliFMDRecPoint* single)
 185   {
 186     if (!single) return kTRUE;
 187     UShort_t det = single->Detector();
 188     Char_t   rng = single->Ring();
 189     UShort_t sec = single->Sector();
 190     UShort_t str = single->Strip();
 191     if (str > 511) {
 192       AliWarning(Form("Bad strip number %d in single", str));
 193       return kTRUE;
 194     }
 195     Double_t edep = fMap(det, rng, sec, str).fEdep;
 196     Int_t    nhit = fMap(det, rng, sec, str).fN;
 197     if (edep > 0) {
 198       fHitEvsRecM->Fill(edep, single->Particles());
 199       fHitEvsRecE->Fill(edep, single->Edep());
 200       fDiffE->Fill((single->Edep() - edep) / edep);
 201     }
 202     if (nhit > 0) fHitsVsRecM->Fill(nhit, single->Particles());
 203     fDiffM->Fill(nhit - single->Particles(), TMath::Abs(single->Eta()));
 204     if (single->Edep()/.03 > 0.1) fRecEloss->Fill(single->Edep() / 0.0300);
 205     return kTRUE;
 206   }
 207   //__________________________________________________________________
 208   Bool_t Finish()
 209   {
 210     gStyle->SetPalette(1);
 211     gStyle->SetOptTitle(0);
 212     gStyle->SetCanvasColor(0);
 213     gStyle->SetCanvasBorderSize(0);
 214     gStyle->SetPadColor(0);
 215     gStyle->SetPadBorderSize(0);
 216     TCanvas* c = 0;
 217
 218     c = new TCanvas("c0", fHitEvsAdc->GetTitle());
 219     fHitEvsAdc->SetStats(kFALSE);
 220     fHitEvsAdc->Draw("COLZ");
 221
 222     c = new TCanvas("c1", fHitEvsRecM->GetTitle());
 223     fHitEvsRecM->SetStats(kFALSE);
 224     fHitEvsRecM->Draw("COLZ");
 225
 226     c = new TCanvas("c2", fHitEvsRecE->GetTitle());
 227     fHitEvsRecE->SetStats(kFALSE);
 228     fHitEvsRecE->Draw("COLZ");
 229
 230     c = new TCanvas("c3", fDiffE->GetTitle());
 231     c->SetLogz();
 232     fDiffE->Draw();
 233
 234     c = new TCanvas("c4", fHitsVsRecM->GetTitle());
 235     c->SetLogz();
 236     fHitsVsRecM->SetStats(kFALSE);
 237     fHitsVsRecM->Draw("COLZ");
 238
 239     c = new TCanvas("c5", fDiffM->GetTitle());
 240     fDiffM->SetFillColor(2);
 241     fDiffM->SetFillStyle(3001);
 242     c->SetLogz();
 243     fDiffM->Draw("colz");
 244
 245     c = new TCanvas("c6", "Hit Eloss, Reco Eloss");
 246     fRecEloss->Scale(1./fRecEloss->GetEntries());
 247     fRecEloss->Draw("hist e");
 248     fRecEloss->Fit("landau", "", "SAME", 2, 4);
 249     TF1* recResp = new TF1(*fRecEloss->GetFunction("landau"));
 250     fHitEloss->Scale(1./fHitEloss->GetEntries());
 251     fHitEloss->Draw("same hist e");
 252     fHitEloss->Fit("landau", "", "SAME", 2, 10);
 253     TF1* hitResp = new TF1(*fHitEloss->GetFunction("landau"));
 254     std::cout << "From hits: MPV="
 255               << hitResp->GetParameter(1) << "+/-"
 256               << hitResp->GetParError(1) << "  Width="
 257               << hitResp->GetParameter(2) << "+/-"
 258               << hitResp->GetParError(2) << "\nFrom recs: MPV="
 259               << recResp->GetParameter(1) << "+/-"
 260               << recResp->GetParError(1) << "  Width="
 261               << recResp->GetParameter(2) << "+/-"
 262               << recResp->GetParError(2) << "\nRatio: MPV(hit/rec)="
 263               << hitResp->GetParameter(1) / recResp->GetParameter(1)
 264               << std::endl;
 265     c->SetLogy();
 266
 267     return kTRUE;
 268   }
 269
 270   ClassDef(DrawHitsRecs,0);
 271 };
 272
 273 //____________________________________________________________________
 274 //
 275 // EOF
 276 //