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   rng,
  69                     UShort_t sec,
  70                     UShort_t str,
  71                     Float_t  /* eta */,
  72                     Float_t  mult)
  73   {
  74     // Cache the energy loss
  75     if (mult > 0)
  76       Info("ProcessESD", "FMD%d%c[%2d,%3d]=%f", det, rng, sec, str, mult);
  77     if (mult/fCorr > 0.001) fMult->Fill(mult/fCorr);
  78     return kTRUE;
  79   }
  80   //__________________________________________________________________
  81   TF1* FitPeak(Int_t n, TH1D* hist, Double_t min, Double_t& max)
  82   {
  83     if (TMath::Abs(max-min) < .25) return 0;
  84     std::cout << "Fit peack in range " << min << " to " << max << std::endl;
  85     TF1* l = new TF1(Form("l%02d", n), "landau", min, max);
  86     hist->GetXaxis()->SetRangeUser(0, 4);
  87     hist->Fit(l, "0Q", "", min, max);
  88     Double_t mpv   = l->GetParameter(1);
  89     Double_t empv  = l->GetParError(1);
  90     Double_t sigma = l->GetParameter(2);
  91     l->SetRange(mpv-empv, mpv+3*sigma);
  92     hist->Fit(l, "EMQ0", "", mpv-3*empv, mpv+3*sigma);
  93     std::cout << "Peak # " << n << " [" << min << "," << max << "]\n"
  94               << " MPV: " << l->GetParameter(1)
  95               << " +/- "  << l->GetParError(1)
  96               << " Var: " << l->GetParameter(2)
  97               << " +/- "  << l->GetParError(2)
  98               << " Chi^2/NDF: " << l->GetChisquare() / l->GetNDF()
  99               << std::endl;
 100     mpv   = l->GetParameter(1);
 101     sigma = l->GetParameter(2);
 102     min   = mpv - sigma * 2; // (n==1 ? 3 : 2);
 103     max   = mpv + sigma * 3;
 104     // l->SetRange(min, max);
 105     l->Draw("same");
 106     return  l;
 107   }
 108   //__________________________________________________________________
 109   void MaxInRange(TH1D* hist, Double_t min, Double_t& mean, Double_t& var)
 110   {
 111     hist->GetXaxis()->SetRangeUser(min, 4);
 112     mean = hist->GetMean();
 113     var  = hist->GetRMS();
 114   }
 115
 116   //__________________________________________________________________
 117   Bool_t Finish()
 118   {
 119     Info("Finish", "Will draw results");
 120     gStyle->SetPalette(1);
 121     gStyle->SetOptTitle(0);
 122     gStyle->SetOptFit(1111111);
 123     gStyle->SetCanvasColor(0);
 124     gStyle->SetCanvasBorderSize(0);
 125     gStyle->SetPadColor(0);
 126     gStyle->SetPadBorderSize(0);
 127
 128     if (fMult->GetEntries() <= 0) return kFALSE;
 129
 130     TCanvas* c = new TCanvas("c", "C");
 131     c->cd();
 132     // c->SetLogy();
 133     fMult->GetXaxis()->SetRangeUser(0,4);
 134     fMult->Scale(1. / fMult->GetEntries());
 135     fMult->SetStats(kFALSE);
 136     fMult->SetFillColor(2);
 137     fMult->SetFillStyle(3001);
 138     fMult->Draw("hist e");
 139
 140     return kTRUE;
 141
 142     Double_t mean, rms;
 143     MaxInRange(fMult, 0.2, mean, rms);
 144     Double_t x1   = mean-rms/2; // .75; // .8;  // .65 / fCorr;
 145     Double_t x2   = mean+rms/2; // 1.3; // 1.7; // fCorr;
 146     TF1*     l1   = FitPeak(1, fMult, x1, x2);
 147     x2            = TMath::Max(mean+rms/2, x2);
 148     MaxInRange(fMult, x2, mean, rms);
 149     Double_t x3   = mean + rms;
 150     TF1*     l2   = FitPeak(2, fMult, x2, x3);
 151     TF1*     f    = 0;
 152     Double_t diff = 0;
 153     if (l2) {
 154       diff          = l2->GetParameter(1)-l1->GetParameter(1);
 155       f             = new TF1("user", "landau(0)+landau(3)", x1, x3);
 156       f->SetParNames("A_{1}", "Mpv_{1}", "#sigma_{1}",
 157                      "A_{2}", "Mpv_{2}", "#sigma_{2}",
 158                      "A_{3}", "Mpv_{3}", "#sigma_{3}");
 159       f->SetParameters(l1->GetParameter(0),
 160                        l1->GetParameter(1),
 161                        l1->GetParameter(2),
 162                        l2 ? l2->GetParameter(0) : 0,
 163                        l2 ? l2->GetParameter(1) : 0,
 164                        l2 ? l2->GetParameter(2) : 0,
 165                        l2->GetParameter(0)/10,
 166                        l2->GetParameter(1) + diff,
 167                        l2->GetParameter(2));
 168     }
 169     else {
 170       x3 = x2;
 171       f  = new TF1("usr", "landau", x1, x3);
 172     }
 173
 174     std::cout << "Range: " << x1 << "-" << x3 << std::endl;
 175
 176     fMult->GetXaxis()->SetRangeUser(0, 4);
 177     fMult->Fit(f, "0QR", "", x1, x3);
 178     fMult->Fit(f, "ME0R", "E1", x1, x3);
 179     fMult->DrawClone("same hist");
 180
 181     l1->SetLineColor(3);
 182     l1->SetLineWidth(2);
 183     l1->SetRange(0, 4);
 184     l1->Draw("same");
 185     if (l2) {
 186       l2->SetLineColor(4);
 187       l2->SetLineWidth(2);
 188       l2->SetRange(0, 4);
 189       l2->Draw("same");
 190     }
 191     f->SetLineWidth(2);
 192     f->SetRange(0, 4);
 193     f->Draw("same");
 194
 195     TLegend* l = new TLegend(0.6, 0.6, .89, .89);
 196     l->AddEntry(l1, "1 particle Landau", "l");
 197     if (l2) l->AddEntry(l2, "2 particle Landau", "l");
 198     l->AddEntry(f,  "1+2 particle Landau", "l");
 199     l->SetFillColor(0);
 200     l->Draw("same");
 201
 202
 203 #if 0
 204     c = new TCanvas("c2", "Landaus");
 205     // c->SetLogy();
 206     fMult->DrawClone("axis");
 207     f->Draw("same");
 208     Double_t* p1 = f->GetParameters();
 209     Double_t* p2 = &(p1[3]);
 210     TF1* ll1 = new TF1("ll1", "landau", 0, 4);
 211     ll1->SetParameters(p1);
 212     ll1->SetLineColor(3);
 213     ll1->Draw("same");
 214     TF1* ll2 = new TF1("ll2", "landau", 0, 4);
 215     ll2->SetParameters(p2);
 216     ll2->SetLineColor(4);
 217     ll2->Draw("same");
 218
 219     Double_t y1  = fMult->GetMinimum() * 1.1;
 220     Double_t y2  = fMult->GetMaximum() * .9;
 221     Double_t xc1 = p1[1]-3*p1[2];
 222     Double_t xc2 = p2[1]-2*p2[2];
 223     Double_t xc3 = p2[1]-2*p2[2]+diff;
 224     TLine* c1 = new TLine(xc1, y1, xc1, y2);
 225     c1->Draw("same");
 226     TLine* c2 = new TLine(xc2, y1, xc2, y2);
 227     c2->Draw("same");
 228     TLine* c3 = new TLine(xc3, y1, xc3, y2);
 229     c3->Draw("same");
 230
 231     l = new TLegend(0.6, 0.6, .89, .89);
 232     l->AddEntry(ll1, "1 particle Landau", "l");
 233     l->AddEntry(ll2, "2 particle Landau", "l");
 234     l->AddEntry(f,  "1+2 particle Landau", "l");
 235     l->SetFillColor(0);
 236     l->Draw("same");
 237 #endif
 238
 239
 240     return kTRUE;
 241   }
 242
 243   ClassDef(DrawESD,0);
 244
 245 };
 246
 247 //____________________________________________________________________
 248 //
 249 // EOF
 250 //