FMD/scripts/DrawESD.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 <TH1D.h>
  15 #include <AliFMDHit.h>
  16 #include <AliFMDDigit.h>
  17 #include <AliFMDInput.h>
  18 #include <AliFMDUShortMap.h>
  19 #include <AliFMDFloatMap.h>
  20 #include <AliFMDRecPoint.h>
  21 #include <AliESDFMD.h>
  22 #include <AliLog.h>
  23 #include <iostream>
  24 #include <TStyle.h>
  25 #include <TArrayF.h>
  26 #include <TCanvas.h>
  27 #include <TMath.h>
  28 #include <TF1.h>
  29 #include <TSpectrum.h>
  30 #include <TLegend.h>
  31 #include <TLine.h>
  32
  33 /** @class DrawESD
  34     @brief Draw digit ADC versus Rec point mult
  35     @code
  36     Root> .L Compile.C
  37     Root> Compile("DrawESD.C")
  38     Root> DrawESD c
  39     Root> c.Run();
  40     @endcode
  41     @ingroup FMD_script
  42  */
  43 class DrawESD : public AliFMDInput
  44 {
  45 private:
  46   TH1D* fMult; // Histogram
  47   const Double_t fCorr;
  48 public:
  49   //__________________________________________________________________
  50   DrawESD(Int_t n=1000, Double_t mmin=-0.5, Double_t mmax=20.5)
  51     : fCorr(1) // 0.68377 / 1.1)
  52   {
  53     AddLoad(kESD);
  54     fMult = new TH1D("mult", " Multiplicity (strip)", n, mmin, mmax);
  55     fMult->Sumw2();
  56     fMult->SetXTitle("Strip Multiplicity");
  57   }
  58   //__________________________________________________________________
  59   /** Begining of event
  60       @param ev Event number
  61       @return @c false on error */
  62   Bool_t Begin(Int_t ev)
  63   {
  64     return AliFMDInput::Begin(ev);
  65   }
  66   //__________________________________________________________________
  67   Bool_t ProcessESD(UShort_t /* det */,
  68                     Char_t   /* ring */,
  69                     UShort_t /* sec */,
  70                     UShort_t /* strip */,
  71                     Float_t  /* eta */,
  72                     Float_t  mult)
  73   {
  74     // Cache the energy loss
  75     if (mult/fCorr > 0.01) fMult->Fill(mult/fCorr);
  76     return kTRUE;
  77   }
  78   //__________________________________________________________________
  79   TF1* FitPeak(Int_t n, TH1D* hist, Double_t min, Double_t& max)
  80   {
  81     TF1* l = new TF1(Form("l%02d", n), "landau", min, max);
  82     hist->GetXaxis()->SetRangeUser(0, 4);
  83     hist->Fit(l, "0Q", "", min, max);
  84     Double_t mpv   = l->GetParameter(1);
  85     Double_t empv  = l->GetParError(1);
  86     Double_t sigma = l->GetParameter(2);
  87     l->SetRange(mpv-empv, mpv+3*sigma);
  88     hist->Fit(l, "EMQ0", "", mpv-3*empv, mpv+3*sigma);
  89     std::cout << "Peak # " << n << " [" << min << "," << max << "]\n"
  90               << " MPV: " << l->GetParameter(1)
  91               << " +/- "  << l->GetParError(1)
  92               << " Var: " << l->GetParameter(2)
  93               << " +/- "  << l->GetParError(2)
  94               << " Chi^2/NDF: " << l->GetChisquare() / l->GetNDF()
  95               << std::endl;
  96     mpv   = l->GetParameter(1);
  97     sigma = l->GetParameter(2);
  98     min   = mpv - sigma * 2; // (n==1 ? 3 : 2);
  99     max   = mpv + sigma * 3;
 100     // l->SetRange(min, max);
 101     l->Draw("same");
 102     return  l;
 103   }
 104   //__________________________________________________________________
 105   void MaxInRange(TH1D* hist, Double_t min, Double_t& mean, Double_t& var)
 106   {
 107     hist->GetXaxis()->SetRangeUser(min, 4);
 108     mean = hist->GetMean();
 109     var  = hist->GetRMS();
 110   }
 111
 112   //__________________________________________________________________
 113   Bool_t Finish()
 114   {
 115     gStyle->SetPalette(1);
 116     gStyle->SetOptTitle(0);
 117     gStyle->SetOptFit(1111111);
 118     gStyle->SetCanvasColor(0);
 119     gStyle->SetCanvasBorderSize(0);
 120     gStyle->SetPadColor(0);
 121     gStyle->SetPadBorderSize(0);
 122
 123     TCanvas* c = new TCanvas("c", "C");
 124     c->cd();
 125     c->SetLogy();
 126     fMult->GetXaxis()->SetRangeUser(0,4);
 127     fMult->Scale(1. / fMult->GetEntries());
 128     fMult->SetStats(kFALSE);
 129     fMult->SetFillColor(2);
 130     fMult->SetFillStyle(3001);
 131     fMult->Draw("hist e");
 132
 133     Double_t mean, rms;
 134     MaxInRange(fMult, 0.2, mean, rms);
 135     Double_t x1   = mean-rms/2; // .75; // .8;  // .65 / fCorr;
 136     Double_t x2   = mean+rms/2; // 1.3; // 1.7; // fCorr;
 137     TF1*     l1   = FitPeak(1, fMult, x1, x2);
 138     x2            = TMath::Max(mean+rms/2, x2);
 139     MaxInRange(fMult, x2, mean, rms);
 140     Double_t x3   = mean + rms;
 141     TF1*     l2   = FitPeak(2, fMult, x2, x3);
 142     Double_t diff = l2->GetParameter(1)-l1->GetParameter(1);
 143     TF1*     f    = new TF1("user", "landau(0)+landau(3)", x1, x3);
 144
 145     fMult->GetXaxis()->SetRangeUser(0, 4);
 146     f->SetParNames("A_{1}", "Mpv_{1}", "#sigma_{1}",
 147                    "A_{2}", "Mpv_{2}", "#sigma_{2}",
 148                    "A_{3}", "Mpv_{3}", "#sigma_{3}");
 149     f->SetParameters(l1->GetParameter(0),
 150                      l1->GetParameter(1),
 151                      l1->GetParameter(2),
 152                      l2->GetParameter(0),
 153                      l2->GetParameter(1),
 154                      l2->GetParameter(2),
 155                      l2->GetParameter(0)/10,
 156                      l2->GetParameter(1) + diff,
 157                      l2->GetParameter(2));
 158     fMult->Fit(f, "0Q", "", x1, x3);
 159     fMult->Fit(f, "ME0", "E1", x1, x3);
 160     fMult->DrawClone("same hist");
 161
 162     l1->SetLineColor(3);
 163     l1->SetLineWidth(2);
 164     l1->SetRange(0, 4);
 165     l1->Draw("same");
 166     l2->SetLineColor(4);
 167     l2->SetLineWidth(2);
 168     l2->SetRange(0, 4);
 169     l2->Draw("same");
 170     f->SetLineWidth(2);
 171     f->SetRange(0, 4);
 172     f->Draw("same");
 173
 174     TLegend* l = new TLegend(0.6, 0.6, .89, .89);
 175     l->AddEntry(l1, "1 particle Landau", "l");
 176     l->AddEntry(l2, "2 particle Landau", "l");
 177     l->AddEntry(f,  "1+2 particle Landau", "l");
 178     l->SetFillColor(0);
 179     l->Draw("same");
 180
 181
 182 #if 1
 183     c = new TCanvas("c2", "Landaus");
 184     c->SetLogy();
 185     fMult->DrawClone("axis");
 186     f->Draw("same");
 187     Double_t* p1 = f->GetParameters();
 188     Double_t* p2 = &(p1[3]);
 189     TF1* ll1 = new TF1("ll1", "landau", 0, 4);
 190     ll1->SetParameters(p1);
 191     ll1->SetLineColor(3);
 192     ll1->Draw("same");
 193     TF1* ll2 = new TF1("ll2", "landau", 0, 4);
 194     ll2->SetParameters(p2);
 195     ll2->SetLineColor(4);
 196     ll2->Draw("same");
 197
 198     Double_t y1  = fMult->GetMinimum() * 1.1;
 199     Double_t y2  = fMult->GetMaximum() * .9;
 200     Double_t xc1 = p1[1]-3*p1[2];
 201     Double_t xc2 = p2[1]-2*p2[2];
 202     Double_t xc3 = p2[1]-2*p2[2]+diff;
 203     TLine* c1 = new TLine(xc1, y1, xc1, y2);
 204     c1->Draw("same");
 205     TLine* c2 = new TLine(xc2, y1, xc2, y2);
 206     c2->Draw("same");
 207     TLine* c3 = new TLine(xc3, y1, xc3, y2);
 208     c3->Draw("same");
 209
 210     l = new TLegend(0.6, 0.6, .89, .89);
 211     l->AddEntry(ll1, "1 particle Landau", "l");
 212     l->AddEntry(ll2, "2 particle Landau", "l");
 213     l->AddEntry(f,  "1+2 particle Landau", "l");
 214     l->SetFillColor(0);
 215     l->Draw("same");
 216 #endif
 217
 218
 219     return kTRUE;
 220   }
 221
 222   ClassDef(DrawESD,0);
 223
 224 };
 225
 226 //____________________________________________________________________
 227 //
 228 // EOF
 229 //