Fixed some warning messages on request from FCA
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 1 Nov 2010 10:48:41 +0000 (10:48 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 1 Nov 2010 10:48:41 +0000 (10:48 +0000)
FMD/AliFMDBaseDigit.cxx
FMD/AliFMDBaseDigitizer.cxx
FMD/AliFMDBaseDigitizer.h
FMD/AliFMDHitDigitizer.cxx
FMD/AliFMDInput.cxx
FMD/AliFMDInput.h
FMD/scripts/DrawESD.C
FMD/scripts/DrawHits.C

index 6f1356c..712e5b5 100644 (file)
@@ -203,8 +203,9 @@ AliFMDBaseDigit::AddTrack(Int_t track)
   else if (fTracks[1] == -1) fTracks[1] = track;
   else if (fTracks[2] == -1) fTracks[2] = track;
   else 
-    AliWarning(Form("While adding track label to %s for %s: "
-                   "All 3 track labels used, can't add reference to track %d",
+    AliFMDDebug(1, ("While adding track label to %s for %s: "
+                   "All 3 track labels used, can't add "
+                   "reference to track %d",
                    ClassName(), GetName(), track));
 }
 
index 1970eb3..35f7edf 100644 (file)
@@ -231,7 +231,8 @@ AliFMDBaseDigitizer::AliFMDBaseDigitizer()
          AliFMDMap::kMaxSectors, 
          AliFMDMap::kMaxStrips),
     fShapingTime(6),
-    fStoreTrackRefs(kTRUE)
+    fStoreTrackRefs(kTRUE), 
+    fIgnoredLabels(0)
 {
   AliFMDDebug(1, ("Constructed"));
   // Default ctor - don't use it
@@ -244,7 +245,8 @@ AliFMDBaseDigitizer::AliFMDBaseDigitizer(AliRunDigitizer* manager)
     fRunLoader(0),
     fEdep(0),        // nDet==0 means 51200 slots
     fShapingTime(6),
-    fStoreTrackRefs(kTRUE)
+    fStoreTrackRefs(kTRUE), 
+    fIgnoredLabels(0)
 {
   // Normal CTOR
   AliFMDDebug(1, ("Constructed"));
@@ -259,7 +261,8 @@ AliFMDBaseDigitizer::AliFMDBaseDigitizer(const Char_t* name,
     fRunLoader(0),
     fEdep(0),        // nDet==0 means 51200 slots
     fShapingTime(6),
-    fStoreTrackRefs(kTRUE)
+    fStoreTrackRefs(kTRUE), 
+    fIgnoredLabels(0)
 {
   // Normal CTOR
   AliFMDDebug(1, (" Constructed"));
@@ -374,6 +377,7 @@ AliFMDBaseDigitizer::DigitizeHits() const
   // the digits array (AliFMD::fDigits)
   //
   AliFMDDebug(5, ("Will now digitize all the summed signals"));
+  fIgnoredLabels = 0;
   AliFMDGeometry* geometry = AliFMDGeometry::Instance();
   
   TArrayI counts(4);
@@ -453,6 +457,10 @@ AliFMDBaseDigitizer::DigitizeHits() const
       } // Sector 
     } // Ring 
   } // Detector 
+  if (fIgnoredLabels > 0) 
+    AliWarning(Form("%d track labels could not be associated with digits "
+                   "due to limited storage facilities in AliDigit", 
+                   fIgnoredLabels));
 }
 
 //____________________________________________________________________
@@ -568,6 +576,7 @@ AliFMDBaseDigitizer::AddDigit(UShort_t        detector,
   fFMD->AddDigitByFields(detector, ring, sector, strip, 
                         count1, count2, count3, count4, 
                         ntot, fStoreTrackRefs ? refs.fArray : 0);
+  if (fStoreTrackRefs && ntot > 3) fIgnoredLabels += ntot - 3;
 }
 
 //____________________________________________________________________
index 351e554..58da125 100644 (file)
@@ -258,6 +258,7 @@ protected:
   AliFMDEdepMap   fEdep;             // Cache of Energy from hits 
   Float_t         fShapingTime;      // Shaping profile parameter
   Bool_t          fStoreTrackRefs;   // Wether to store track references
+  mutable Int_t   fIgnoredLabels;    //! Number of labels not assigned 
   
   /** Copy CTOR 
       @param o object to copy from  */
@@ -267,7 +268,8 @@ protected:
       fRunLoader(0),
       fEdep(o.fEdep),
       fShapingTime(o.fShapingTime),
-      fStoreTrackRefs(o.fStoreTrackRefs)
+      fStoreTrackRefs(o.fStoreTrackRefs), 
+      fIgnoredLabels(o.fIgnoredLabels)
   {}
   /** Assignment operator
       @return Reference to this object */
@@ -278,9 +280,11 @@ protected:
     fEdep           = o.fEdep;
     fShapingTime    = o.fShapingTime;
     fStoreTrackRefs = o.fStoreTrackRefs;
+    fIgnoredLabels  = o.fIgnoredLabels;
     return *this; 
   }
-  ClassDef(AliFMDBaseDigitizer,3) // Base class for FMD digitizers
+
+  ClassDef(AliFMDBaseDigitizer,4) // Base class for FMD digitizers
 };
 
 
index 5e400b1..5c46bc0 100644 (file)
@@ -505,6 +505,7 @@ AliFMDHitDigitizer::AddDigit(UShort_t        detector,
   fFMD->AddSDigitByFields(detector, ring, sector, strip, edep,
                          count1, count2, count3, count4, 
                          ntotal, nprim, fStoreTrackRefs ? refs.fArray : 0);
+  if (fStoreTrackRefs && nprim > 3) fIgnoredLabels += nprim - 3;
 }
 
 //____________________________________________________________________
index cc7c89b..1f31bbe 100644 (file)
@@ -105,7 +105,8 @@ AliFMDInput::AliFMDInput()
     fTreeMask(0), 
     fRawFile(""),
     fIsInit(kFALSE),
-    fEventCount(0)
+    fEventCount(0), 
+    fNEvents(-1)
 {
 
   // Constructor of an FMD input object.  Specify what data to read in
@@ -149,7 +150,8 @@ AliFMDInput::AliFMDInput(const char* gAliceFile)
     fTreeMask(0), 
     fRawFile(""),
     fIsInit(kFALSE),
-    fEventCount(0)
+    fEventCount(0),
+    fNEvents(-1)
 {
   
   // Constructor of an FMD input object.  Specify what data to read in
@@ -351,7 +353,6 @@ AliFMDInput::Begin(Int_t event)
 
   // Get the event 
   if (fLoader && fLoader->GetEvent(event)) return kFALSE;
-  AliInfo(Form("Now in event %8d/%8d", event, NEvents()));
 
   // Possibly load global kinematics information 
   if (TESTBIT(fTreeMask, kKinematics)) {
@@ -852,19 +853,25 @@ AliFMDInput::End()
 
 //____________________________________________________________________
 Bool_t
-AliFMDInput::Run()
+AliFMDInput::Run(UInt_t maxEvents)
 {
   // Run over all events and files references in galice.root 
 
   Bool_t retval;
   if (!(retval = Init())) return retval;
 
-  Int_t nEvents = NEvents();
-  for (Int_t event = 0; nEvents < 0 || event < nEvents; event++) {
+  fNEvents = NEvents();
+  if (fNEvents < 0)       fNEvents = maxEvents;
+  else if (maxEvents > 0) fNEvents = TMath::Min(fNEvents,Int_t(maxEvents));
+
+  Int_t event = 0;
+  for (; fNEvents < 0 || event < fNEvents; event++) {
+    printf("\rEvent %8d/%8d ...", event, fNEvents);
     if (!(retval = Begin(event))) break;
     if (!(retval = Event())) break;
     if (!(retval = End())) break;
   }
+  printf("Looped over %8d events\n", event+1);
   if (!retval) return retval;
   retval = Finish();
   return retval;
index 70f6cae..214879c 100644 (file)
@@ -158,7 +158,7 @@ public:
   virtual Bool_t Finish() { return kTRUE; }
   /** Run a full job. 
       @return @c false on error  */
-  virtual Bool_t Run();
+  virtual Bool_t Run(UInt_t maxEvents=0);
 
   /** Loop over all hits, and call ProcessHit with that hit, and
       optionally the corresponding kinematics track. 
@@ -318,7 +318,8 @@ protected:
       fTreeMask(0),
       fRawFile(""),      
       fIsInit(kFALSE),
-      fEventCount(0)
+      fEventCount(0),
+      fNEvents(-1)
   {}
   /** Assignement operator 
       @return  REference to this */
@@ -366,6 +367,7 @@ protected:
   TString          fRawFile;    // Raw input file
   Bool_t           fIsInit;     // Have we been initialized 
   Int_t            fEventCount; // Event counter 
+  Int_t            fNEvents;    // The maximum number of events
   ClassDef(AliFMDInput,0)  //Hits for detector FMD
 };
 
index 94924fc..ed0c488 100644 (file)
 #include <TF1.h>
 #include <TSpectrum.h>
 #include <TLegend.h>
+#include <TLatex.h>
 #include <TLine.h>
+#include <TPolyMarker.h>
+#include <TSpectrum.h>
+#include <TList.h>
+#include <TGraph.h>
+
+Double_t landau(Double_t* xp, Double_t* pp)
+{
+  Double_t x     = xp[0];
+  Double_t A     = pp[0];
+  Double_t mpv   = pp[1];
+  Double_t w     = pp[2];
+
+  Double_t v     = mpv; //  + w * 0.22278;
+  return A * TMath::Landau(x, v, w); 
+}
+
+Double_t foldLandau(Double_t* xp, Double_t* pp)
+{
+  Double_t x     = xp[0];
+  Double_t A     = pp[0];
+  Double_t mpv   = pp[1];
+  Double_t w     = pp[2];
+  Double_t B     = pp[3];
+  Double_t sigma = pp[4];
+  
+  return A * (TMath::Landau(x,mpv,w) + B * TMath::Gaus(x, mpv, sigma));
+}
 
 /** @class DrawESD
     @brief Draw digit ADC versus Rec point mult
@@ -45,15 +73,23 @@ class DrawESD : public AliFMDInput
 private:
   TH1D* fMult; // Histogram 
   const Double_t fCorr;
+  TList fCleanup;
 public:
   //__________________________________________________________________
-  DrawESD(Int_t n=1000, Double_t mmin=-0.5, Double_t mmax=20.5) 
+  DrawESD(Int_t n=2000, Double_t mmin=-0.5, Double_t mmax=15.5) 
     : fCorr(1) // 0.68377 / 1.1)
   { 
     AddLoad(kESD);
-    fMult = new TH1D("mult", " Multiplicity (strip)", n, mmin, mmax);
+    fMult = new TH1D("mult", "#DeltaE/#DeltaE_{MIP)", n, mmin, mmax);
     fMult->Sumw2();
-    fMult->SetXTitle("Strip Multiplicity");
+    fMult->SetFillColor(kRed+1);
+    fMult->SetFillStyle(3001);
+    fMult->SetXTitle("#DeltaE/#DeltaE_{MIP}");
+    fCleanup.Add(fMult);
+  }
+  ~DrawESD() 
+  {
+    fCleanup.Delete();
   }
   //__________________________________________________________________
   /** Begining of event
@@ -64,24 +100,40 @@ public:
     return AliFMDInput::Begin(ev);
   }
   //__________________________________________________________________
-  Bool_t ProcessESD(UShort_t det, 
-                   Char_t   rng, 
-                   UShort_t sec, 
-                   UShort_t str, 
-                   Float_t  /* eta */, 
+  Bool_t ProcessESD(UShort_t /* det */, 
+                   Char_t   /* rng */, 
+                   UShort_t /* sec */, 
+                   UShort_t /* str */, 
+                   Float_t  eta, 
                    Float_t  mult)
   {
     // Cache the energy loss 
-    if (mult > 0) 
-      Info("ProcessESD", "FMD%d%c[%2d,%3d]=%f", det, rng, sec, str, mult);
-    if (mult/fCorr > 0.001) fMult->Fill(mult/fCorr);
+    // if (mult > 0) 
+    //   Info("ProcessESD", "FMD%d%c[%2d,%3d]=%f", det, rng, sec, str, mult);
+    if (mult <= 0 || mult == AliESDFMD::kInvalidMult) return kTRUE;
+    Double_t x = mult;
+    if (!fESD->IsAngleCorrected()) {
+      Double_t theta = 2 * TMath::ATan(TMath::Exp(-eta));
+      Double_t corr  = TMath::Abs(TMath::Cos(theta));
+      Double_t cmult = corr * mult;
+      x = cmult;
+    }
+    if (x > 0.001) fMult->Fill(x);
+    return kTRUE;
+  }
+  //__________________________________________________________________
+  Bool_t ProcessESDs()
+  {
+    if (!AliFMDInput::ProcessESDs()) return kFALSE;
+    // Info("ProcessESDs", "ESD is %sangle corrected", 
+    //      fESD->IsAngleCorrected() ? "" : "not ");
     return kTRUE;
   }
   //__________________________________________________________________
   TF1* FitPeak(Int_t n, TH1D* hist, Double_t min, Double_t& max)
   {
     if (TMath::Abs(max-min) < .25) return 0;
-    std::cout << "Fit peack in range " << min << " to " << max << std::endl;
+    std::cout << "Fit peak in range " << min << " to " << max << std::endl;
     TF1* l = new TF1(Form("l%02d", n), "landau", min, max);
     hist->GetXaxis()->SetRangeUser(0, 4);
     hist->Fit(l, "0Q", "", min, max);
@@ -114,6 +166,47 @@ public:
   }
 
   //__________________________________________________________________
+  const char* PrettyFloat(float x)
+  {
+    if (x == 0) return Form("%9.4f", x);
+    float e = TMath::Floor(TMath::Log10(TMath::Abs(x)));
+    if (TMath::Abs(e) < 4) {
+      return Form("%9.4f", x);
+    }
+    float f = x / TMath::Power(10,e);
+    return Form("%4.2f#times10^{%d}", f, int(e));
+  }
+  //__________________________________________________________________
+  void ShowFit(Double_t x1, Double_t y1, const char* title, 
+              TF1* f, Double_t dx=0, Double_t dy=0.05)
+  {
+    Double_t x = x1, y = y1;
+    TLatex* latex = new TLatex(x, y, title);
+    latex->SetTextFont(132);
+    latex->SetTextSize(0.8*dy);
+    latex->SetNDC();
+    latex->Draw();
+    x -= dx;
+    y -= dy;
+    const Double_t eqDx=0.1;
+    Double_t chi2 = f->GetChisquare();
+    Int_t    ndf  = f->GetNDF();
+    Double_t prob = f->GetProb();
+    latex->DrawLatex(x, y, "#chi^{2}/NDF");
+    latex->DrawLatex(x+eqDx, y, Form("= %7.4f/%3d=%5.2f (%7.3f%%)", 
+                                    chi2, ndf, chi2/ndf, 100*prob));
+    Int_t     n = f->GetNpar();
+    Double_t* p = f->GetParameters();
+    Double_t* e = f->GetParErrors();
+    for (int i = 0; i < n; i++) { 
+      x -= dx;
+      y -= dy;
+      latex->DrawLatex(x, y, f->GetParName(i));
+      latex->DrawLatex(x+eqDx, y, Form("= %s", PrettyFloat(p[i])));
+      latex->DrawLatex(x+2.2*eqDx, y, Form("#pm %s", PrettyFloat(e[i])));
+    }
+  }      
+  //__________________________________________________________________
   Bool_t Finish()
   {
     Info("Finish", "Will draw results");
@@ -127,20 +220,245 @@ public:
     
     if (fMult->GetEntries() <= 0) return kFALSE;
     
-    TCanvas* c = new TCanvas("c", "C");
+    TCanvas* c = new TCanvas("c", "C", 1200, 800);
+    fCleanup.Add(c);
     c->cd();
-    // c->SetLogy();
-    fMult->GetXaxis()->SetRangeUser(0,4);
-    fMult->Scale(1. / fMult->GetEntries());
-    fMult->SetStats(kFALSE);
-    fMult->SetFillColor(2);
-    fMult->SetFillStyle(3001);
-    fMult->Draw("hist e");
+    c->SetLogy();
+    c->SetTopMargin(0.05);
+    c->SetRightMargin(0.05);
+    c->SetFillColor(0);
+    c->SetBorderMode(0);
+
+    TLegend* leg = new TLegend(.1, .1, .4, .2);
+    leg->SetFillColor(0);
+    leg->SetBorderSize(1);
+
+    DrawMult(c, leg);
+
+    Double_t xmax=0, xmin=0, ymax=0;
+    FindMinMax(xmin, xmax, ymax);
+
+    FitLandau(xmin, xmax, ymax, leg);
+    DrawResponse(xmax, ymax, leg);
     
+    // TF1* f = FitMultiLandau(xmin, xmax, leg);
+    // DrawLandaus(f);
+
+    c->cd();
+    leg->Draw();
+
+    c->Modified();
+    c->Update();
+    c->cd();
+    c->SaveAs("esd_eloss.png");
+
+    fCleanup.Add(leg);
+
     return kTRUE;
     
-    Double_t mean, rms;
-    MaxInRange(fMult, 0.2, mean, rms);
+  }
+  //__________________________________________________________________
+  /** 
+   * Draw the multiplicity distribution
+   * 
+   */
+  void DrawMult(TCanvas* /* c */, TLegend* leg) 
+  {
+    fMult->GetXaxis()->SetRangeUser(0.2,20);
+    Double_t integral = fMult->Integral();
+    Info("DrawMult", "Integral in range [0.2,20]=%f (%f entries)", 
+        integral, fMult->GetEntries());
+    fMult->Scale(1. / integral);
+    Info("DrawMult", "Integral in range [0.2,20]=%f (%f entries)", 
+        fMult->Integral(), fMult->GetEntries());
+    Double_t max = 1.5 * fMult->GetMaximum();
+    fMult->GetXaxis()->SetRangeUser(0,4);
+    fMult->SetMaximum(max);
+    fMult->SetStats(kFALSE);
+    fMult->Draw("e same");
+    leg->AddEntry(fMult, "Strip signal", "lf");
+
+  }
+  //__________________________________________________________________
+  /** 
+   * Find the minimum and maximum values in range
+   * 
+   * @param xmin 
+   * @param xmax 
+   * @param ymax 
+   */
+  void FindMinMax(Double_t& xmin, Double_t& xmax, Double_t& ymax)
+  {
+    fMult->GetXaxis()->SetRangeUser(0.1,4); 
+    TSpectrum    s(10);
+    Int_t        nPeak  = s.Search(fMult);
+    fMult->GetXaxis()->SetRangeUser(0.1, 4);
+    Int_t        bmax   = fMult->GetMaximumBin();
+    xmax                = fMult->GetBinCenter(bmax);
+    fMult->GetXaxis()->SetRangeUser(0.1, xmax);
+    Double_t     bmin   = fMult->GetMinimumBin();
+    xmin                = fMult->GetBinCenter(bmin);
+    ymax                = fMult->GetBinContent(bmax);
+    Info("Finish", "Simple peak finding found x_max=%f, x_min=%f y_max=%g", 
+        xmax, xmin, ymax);
+
+    if (nPeak > 1) {
+      TPolyMarker* pm = static_cast<TPolyMarker*>(fMult->GetListOfFunctions()
+                                                 ->FindObject("TPolyMarker"));
+
+      // Peaks are ordered by size 
+      Double_t* peakX = pm->GetX();
+      Double_t* peakY = pm->GetY();
+      Double_t  xlmax = peakX[1];
+      xmax            = peakX[0];
+      ymax            = peakY[0];
+      if (xmax < xlmax) { 
+       xmax  = xlmax; 
+       xlmax = peakX[0]; 
+       ymax  = peakY[1];
+      }
+
+      fMult->GetXaxis()->SetRangeUser(xlmax, xmax);
+      bmin  = fMult->GetMinimumBin();
+      xmin  = fMult->GetBinCenter(bmin);
+      Info("Finish", "Spectrum peak finding found x_max=%f, x_min=%f y_max=%g", 
+          xmax, xmin, ymax);
+    }
+  }
+  //__________________________________________________________________
+  /** 
+   * Fit a landau and a landau+gaussian to the data
+   * 
+   * @param xmin  Minimum of peak range
+   * @param xmax  Maximum of peak range
+   * @param ymax  Y value in MIP peak
+   * @param leg   Legend
+   */
+  void FitLandau(Double_t xmin, Double_t xmax, Double_t& ymax, TLegend* leg)
+  {
+    fMult->GetXaxis()->SetRangeUser(xmin, xmax+(xmax-xmin));
+    Double_t mean = fMult->GetMean();
+    Double_t var  = fMult->GetRMS();
+    Info("Finish", "Peak range [%f,%f] mean=%f, var=%f", 
+        xmin, xmax+(xmax-xmin), mean, var);
+    Double_t lowCut = mean-var;
+    Double_t hiCut  = 4; // 2*mean;
+    Info("Finish", "Low cut set to %f", lowCut);
+    fMult->GetXaxis()->SetRangeUser(0, hiCut);
+    
+    TF1* pl = MakeLandau(lowCut, hiCut, xmax, var/2);
+    fMult->Fit(pl, "NEM", "", lowCut, hiCut);
+    ymax = pl->GetMaximum();
+    Info("Finish", "Maximum of landau is at %f (y_max=%g)", 
+        pl->GetMaximumX(), ymax);
+
+    TF1* gl  = MakeFoldLandau(lowCut, hiCut, pl, xmax, var/2);
+    gl->SetLineColor(kRed+1);
+    // gl->SetParLimits(1, xmax-var, xmax+var);
+    fMult->Fit(gl, "NEM", "", lowCut, hiCut);
+    TF1* l  = MakeLandau(lowCut, hiCut, xmax, var/2);
+    l->SetLineColor(kGreen+1);
+    l->SetParameters(gl->GetParameter(0),
+                    gl->GetParameter(1),
+                    gl->GetParameter(2));
+    TF1* g  = new TF1("g", "gaus", lowCut, hiCut);
+    g->SetParNames("A", "#mu", "#sigma");
+    g->SetLineColor(kBlue+1);
+    g->SetParameters(gl->GetParameter(3)*gl->GetParameter(0),
+                    gl->GetParameter(1),
+                    gl->GetParameter(2));
+    fMult->GetXaxis()->SetRangeUser(0,4);
+    fMult->DrawCopy("E");
+    fMult->DrawCopy("HIST SAME");
+    pl->Draw("same");
+    gl->Draw("same");
+    g->Draw("Same");
+    l->Draw("Same");
+    fCleanup.Add(pl);
+    fCleanup.Add(l);
+    fCleanup.Add(g);
+    fCleanup.Add(gl);
+
+    ShowFit(.5, .90, "Landau", pl);
+    ShowFit(.5, .65, "Landau+Gaussian", gl);
+
+    leg->AddEntry(pl, Form("Landau fit - #chi^{2}/NDF=%f", 
+                          pl->GetChisquare()/pl->GetNDF()), "l");
+    leg->AddEntry(gl, Form("Landau+Gaussian fit - #chi^{2}/NDF=%f", 
+                          gl->GetChisquare()/gl->GetNDF()), "l");
+    leg->AddEntry(l, "Landau part", "l");
+    leg->AddEntry(g, "Gaussian part", "l");
+  }
+  
+  //__________________________________________________________________
+  /** 
+   * Superimpose the response graph from the RPP 
+   * 
+   * @param ymax Y value of multiplicity spectra in the landau peak
+   * @param leg  Legend
+   */  
+  void DrawResponse(Double_t xmax, Double_t ymax, TLegend* leg) 
+  {
+    TGraph*   resp = GetResp();
+    // TGraph*   corr = GetCorr();
+    Double_t* x    = resp->GetX();
+    Double_t* y    = resp->GetY(); // [MeV cm^2/g]
+    TGraph*   gr   = new TGraph(resp->GetN());
+    gr->SetName(Form("%sCorr", resp->GetName()));
+    gr->SetTitle(resp->GetTitle());
+    gr->SetLineStyle(resp->GetLineStyle());
+    gr->SetLineColor(kMagenta+1);
+    gr->SetLineWidth(2);
+    TGraph*   gr2   = new TGraph(resp->GetN());
+    gr2->SetName(Form("%sCorr", resp->GetName()));
+    gr2->SetTitle(resp->GetTitle());
+    gr2->SetLineStyle(resp->GetLineStyle());
+    gr2->SetLineColor(kCyan+1);
+    gr2->SetLineWidth(2);
+    // const Double_t rho = 2.33;   // [g/cm^3] 
+    // const Double_t thk = 0.320;  // [cm]
+    const Double_t mip = 1.664;  // [MeV cm^2/g]
+    // const Double_t bgm = 3.4601; // beta*gamma of a MIP
+    // Double_t  xs2 = corr->Eval(bgm); // [1]
+    // Double_t  xss = 1.1;
+    Double_t  xs  = 1/mip;
+    Double_t xxmax = 0;
+    Double_t yymax = 0;
+    for (Int_t i = 0; i < gr->GetN(); i++) {
+      if (y[i] > yymax) { 
+       yymax = y[i];
+       xxmax = xs * x[i];
+      }
+      gr->SetPoint(i, x[i] * xs, y[i] * ymax);
+    }
+    Info("DrawResponse", "Maximum at x=%f (xmax=%f)", xxmax, xmax);
+    Double_t xs2 = xmax / xxmax;
+    Info("DrawResponse", "Correction factor: %f", xs2);
+    for (Int_t i = 0; i < gr->GetN(); i++) {
+      gr2->SetPoint(i, x[i] * xs * xs2, y[i] * ymax);
+    }
+    gr->Draw("C same");
+    gr2->Draw("C same");
+
+    leg->AddEntry(gr, "Response", "l");
+    leg->AddEntry(gr2, "Response", "l");
+  }
+  //__________________________________________________________________
+  /** 
+   * Fit sum of landaus to the multiplicity distribution
+   * 
+   * @param xmin Minimum of range 
+   * @param xmax Maximum of range 
+   * @param leg  Legend
+   * 
+   * @return Fitted function 
+   */    
+  TF1* FitMultiLandau(Double_t xmin, Double_t xmax, TLegend* leg)
+  {
+    fMult->GetXaxis()->SetRangeUser(xmin, xmax+(xmax-xmin));
+    Double_t mean = fMult->GetMean();
+    Double_t rms  = fMult->GetRMS();
+
     Double_t x1   = mean-rms/2; // .75; // .8;  // .65 / fCorr;
     Double_t x2   = mean+rms/2; // 1.3; // 1.7; // fCorr;
     TF1*     l1   = FitPeak(1, fMult, x1, x2);
@@ -174,9 +492,8 @@ public:
     std::cout << "Range: " << x1 << "-" << x3 << std::endl;
     
     fMult->GetXaxis()->SetRangeUser(0, 4);
-    fMult->Fit(f, "0QR", "", x1, x3);
-    fMult->Fit(f, "ME0R", "E1", x1, x3);
-    fMult->DrawClone("same hist");
+    fMult->Fit(f, "NQR", "", x1, x3);
+    fMult->Fit(f, "NMER", "E1", x1, x3);
 
     l1->SetLineColor(3);
     l1->SetLineWidth(2);
@@ -192,16 +509,26 @@ public:
     f->SetRange(0, 4);
     f->Draw("same");
 
-    TLegend* l = new TLegend(0.6, 0.6, .89, .89);
-    l->AddEntry(l1, "1 particle Landau", "l");
-    if (l2) l->AddEntry(l2, "2 particle Landau", "l");
-    l->AddEntry(f,  "1+2 particle Landau", "l");
-    l->SetFillColor(0);
-    l->Draw("same");
+    fCleanup.Add(l1);
+    if (l2) fCleanup.Add(l2);
+    fCleanup.Add(f);
 
 
-#if 0
-    c = new TCanvas("c2", "Landaus");
+    leg->AddEntry(l1, "1 particle Landau", "l");
+    if (l2) leg->AddEntry(l2, "2 particle Landau", "l");
+    leg->AddEntry(f,  "1+2 particle Landau", "l");
+
+    return f;
+  }
+  //__________________________________________________________________
+  /** 
+   * Draw landau functions in a separate canvas 
+   * 
+   * @param f Multi-landau function 
+   */  
+  void DrawLandaus(TF1* f)
+  {
+    TCanvas* c = new TCanvas("c2", "Landaus");
     // c->SetLogy();
     fMult->DrawClone("axis");
     f->Draw("same");
@@ -216,11 +543,12 @@ public:
     ll2->SetLineColor(4);
     ll2->Draw("same");
 
-    Double_t y1  = fMult->GetMinimum() * 1.1;
-    Double_t y2  = fMult->GetMaximum() * .9;
-    Double_t xc1 = p1[1]-3*p1[2];
-    Double_t xc2 = p2[1]-2*p2[2];
-    Double_t xc3 = p2[1]-2*p2[2]+diff;
+    Double_t diff = ll2->GetParameter(1)-ll1->GetParameter(1);
+    Double_t y1   = fMult->GetMinimum() * 1.1;
+    Double_t y2   = fMult->GetMaximum() * .9;
+    Double_t xc1  = p1[1]-3*p1[2];
+    Double_t xc2  = p2[1]-2*p2[2];
+    Double_t xc3  = p2[1]-2*p2[2]+diff;
     TLine* c1 = new TLine(xc1, y1, xc1, y2);
     c1->Draw("same");
     TLine* c2 = new TLine(xc2, y1, xc2, y2);
@@ -228,16 +556,194 @@ public:
     TLine* c3 = new TLine(xc3, y1, xc3, y2);
     c3->Draw("same");
 
-    l = new TLegend(0.6, 0.6, .89, .89);
+    TLegend* l = new TLegend(0.6, 0.6, .89, .89);
     l->AddEntry(ll1, "1 particle Landau", "l");
     l->AddEntry(ll2, "2 particle Landau", "l");
     l->AddEntry(f,  "1+2 particle Landau", "l");
     l->SetFillColor(0);
     l->Draw("same");
-#endif
 
+    c->Modified();
+    c->Update();
+    c->cd();
+  }
 
-    return kTRUE;
+  //__________________________________________________________________
+  /** 
+   * Get the response functin @f$ f(\Delta_p/x)@f$ from Review of
+   * Particle Physics (fig. 27.7).  It is scaled to the value at MPV.
+   *
+   * @return Graph of response 
+   */ 
+  TGraph* GetResp()
+  {
+    static TGraph*  graph = 0;
+    if (!graph) {
+      graph = new TGraph;
+      graph->SetName("si_resp");
+      graph->SetTitle("f(#Delta/x) scaled to the MPV value ");
+      graph->GetHistogram()->SetXTitle("#Delta/x (MeVcm^{2}/g)");
+      graph->GetHistogram()->SetYTitle("f(#Delta/x)");
+      graph->SetLineColor(kBlue+1);
+      graph->SetLineWidth(2);
+      graph->SetFillColor(kBlue+1);
+      graph->SetMarkerStyle(21);
+      graph->SetMarkerSize(0.6);
+#if 0
+      // Figure 27.7 or Review of Particle physics - Straggeling function in 
+      // silicon of 500 MeV Pions, normalised to unity at the most probable 
+      // value.   
+      graph->SetPoint(0,0.808094,0.00377358);
+      graph->SetPoint(1,0.860313,0.0566038);
+      graph->SetPoint(2,0.891645,0.116981);
+      graph->SetPoint(3,0.912533,0.181132);
+      graph->SetPoint(4,0.928198,0.260377);
+      graph->SetPoint(5,0.938642,0.320755);
+      graph->SetPoint(6,0.954308,0.377358);
+      graph->SetPoint(7,0.964752,0.433962);
+      graph->SetPoint(8,0.975196,0.490566);
+      graph->SetPoint(9,0.98564,0.550943);
+      graph->SetPoint(10,0.996084,0.611321);
+      graph->SetPoint(11,1.00653,0.667925);
+      graph->SetPoint(12,1.02219,0.732075);
+      graph->SetPoint(13,1.03264,0.784906);
+      graph->SetPoint(14,1.0483,0.845283);
+      graph->SetPoint(15,1.06397,0.901887);
+      graph->SetPoint(16,1.09008,0.958491);
+      graph->SetPoint(17,1.10574,0.984906);
+      graph->SetPoint(18,1.13708,1);
+      graph->SetPoint(19,1.13708,1);
+      graph->SetPoint(20,1.15796,0.988679);
+      graph->SetPoint(21,1.17363,0.966038);
+      graph->SetPoint(22,1.19974,0.916981);
+      graph->SetPoint(23,1.2154,0.89434);
+      graph->SetPoint(24,1.23629,0.837736);
+      graph->SetPoint(25,1.2624,0.784906);
+      graph->SetPoint(26,1.28329,0.724528);
+      graph->SetPoint(27,1.3094,0.664151);
+      graph->SetPoint(28,1.32507,0.611321);
+      graph->SetPoint(29,1.3564,0.550943);
+      graph->SetPoint(30,1.41384,0.445283);
+      graph->SetPoint(31,1.44517,0.392453);
+      graph->SetPoint(32,1.48695,0.335849);
+      graph->SetPoint(33,1.52872,0.286792);
+      graph->SetPoint(34,1.58094,0.237736);
+      graph->SetPoint(35,1.63838,0.196226);
+      graph->SetPoint(36,1.68016,0.169811);
+      graph->SetPoint(37,1.75326,0.135849);
+      graph->SetPoint(38,1.81593,0.113208);
+      graph->SetPoint(39,1.89426,0.0981132);
+      graph->SetPoint(40,1.96214,0.0830189);
+      graph->SetPoint(41,2.0718,0.0641509);
+      graph->SetPoint(42,2.19191,0.0490566);
+      graph->SetPoint(43,2.31723,0.0415094);
+      graph->SetPoint(44,2.453,0.0301887);
+      graph->SetPoint(45,2.53133,0.0264151);
+      graph->SetPoint(46,2.57833,0.0264151);
+#else
+      graph->SetPoint(0,0.8115124,0.009771987);
+      graph->SetPoint(1,0.9198646,0.228013);
+      graph->SetPoint(2,0.996614,0.5895765);
+      graph->SetPoint(3,1.041761,0.8241042);
+      graph->SetPoint(4,1.059819,0.8794788);
+      graph->SetPoint(5,1.077878,0.9348534);
+      graph->SetPoint(6,1.100451,0.980456);
+      graph->SetPoint(7,1.141084,0.9967427);
+      graph->SetPoint(8,1.204289,0.9153094);
+      graph->SetPoint(9,1.276524,0.742671);
+      graph->SetPoint(10,1.402935,0.465798);
+      graph->SetPoint(11,1.515801,0.3029316);
+      graph->SetPoint(12,1.73702,0.1465798);
+      graph->SetPoint(13,1.985327,0.08143322);
+      graph->SetPoint(14,2.301354,0.04234528);
+      graph->SetPoint(15,2.56772,0.02931596);
+#endif
+    }
+    return graph;
+  }
+  //__________________________________________________________________
+  /** 
+   * Get the correction to Bethe-Bloc from Review of Particle Physics
+   * (fig 27.8).
+   *
+   * @return correction graph
+   */
+  TGraph* GetCorr() 
+  {
+    static TGraph* graph = 0;
+    if (!graph) {
+      graph = new TGraph(14);
+      graph->SetName("graph");
+      graph->SetTitle("(#Delta_{p}/x)/(dE/dx)|_{mip} for 320#mu Si");
+      graph->GetHistogram()->SetXTitle("#beta#gamma = p/m");
+      graph->GetHistogram()->SetYTitle("#frac{#Delta_{p}/x)}{(dE/dx)|_{mip}}");
+      graph->SetFillColor(1);
+      graph->SetLineColor(7);
+      graph->SetMarkerStyle(21);
+      graph->SetMarkerSize(0.6);
+      graph->SetPoint(0,1.196058,0.9944915);
+      graph->SetPoint(1,1.28502,0.9411017);
+      graph->SetPoint(2,1.484334,0.8559322);
+      graph->SetPoint(3,1.984617,0.7491525);
+      graph->SetPoint(4,2.658367,0.6983051);
+      graph->SetPoint(5,3.780227,0.6779661);
+      graph->SetPoint(6,4.997358,0.6741525);
+      graph->SetPoint(7,8.611026,0.684322);
+      graph->SetPoint(8,15.28296,0.6995763);
+      graph->SetPoint(9,41.54516,0.7186441);
+      graph->SetPoint(10,98.91461,0.7288136);
+      graph->SetPoint(11,203.2734,0.7326271);
+      graph->SetPoint(12,505.6421,0.7338983);
+      graph->SetPoint(13,896.973,0.7338983);
+    }
+    return graph;
+  }
+  //__________________________________________________________________
+  /** 
+   * Make a Landau function object. 
+   * 
+   * @param min  Minimum of fit range 
+   * @param max  Maximum of fit range
+   * @param p    Peak position
+   * @param v    Variance around peak
+   * 
+   * @return Landau function object
+   */
+  TF1* MakeLandau(Double_t min, Double_t max, Double_t p, Double_t v)
+  {
+    TF1* f = new TF1("l", "landau", min, max);
+    f->SetParNames("A", "#delta", "#xi");
+    f->SetParameters(1, p, p/10);
+    if      (false) f->FixParameter(1,p);
+    else if (false) f->SetParLimits(1, p-v, p+v);
+    return f;
+  }
+  //__________________________________________________________________
+  /** 
+   * Make a Landau, folded with a gaussian, function object
+   * 
+   * @param min  Minimum of fit range 
+   * @param max  Maximum of fit range
+   * @param l    Seed Landau function object
+   * @param p    Peak position
+   * @param v    Variance around peak
+   * 
+   * @return Landau+Gaus function object
+   */
+  TF1* MakeFoldLandau(Double_t min, Double_t max, TF1* l, 
+                     Double_t p, Double_t v)
+  {
+    TF1* f = new TF1("gl", &foldLandau, min, max, 5);
+    f->SetParNames("A", "#delta", "#xi", "B", "#sigma");
+    f->SetParameters(l->GetParameter(0), 
+                    l->GetParameter(1),
+                    l->GetParameter(2),
+                    l->GetParameter(0)/1000,
+                    l->GetParameter(2));
+    f->SetParLimits(3, 1e-7, 1e-1);
+    if      (false) f->FixParameter(1,p);
+    else if (true)  f->SetParLimits(1, p-2*v, p+2*v);
+    return f;
   }
 
   ClassDef(DrawESD,0);
index eb5a7a3..d58bc85 100644 (file)
@@ -14,6 +14,7 @@
 #include <TH2D.h>
 #include <AliFMDHit.h>
 #include <AliFMDInput.h>
+#include <AliStack.h>
 #include <iostream>
 #include <TStyle.h>
 #include <TArrayF.h>
@@ -42,17 +43,21 @@ class DrawHits : public AliFMDInput
 private:
   TH2D* fElossVsPMQ; // Histogram 
   TH1D* fEloss;
+  TH1D* fBetaGamma;
   TParticlePDG* fPdg;
   const Double_t fRho;
   const Double_t fBetaGammaMip;
 public:
   //__________________________________________________________________
   DrawHits(const char* pdgName="pi+",
-          Int_t m=1000, Double_t emin=1, Double_t emax=1000, 
+          Int_t m=500, Double_t emin=1, Double_t emax=1000, 
           Int_t n=900, Double_t tmin=1e-2, Double_t tmax=1e3) 
     : AliFMDInput("galice.root"),
+      fElossVsPMQ(0),
+      fEloss(0),
+      fBetaGamma(0),
       fPdg(0), 
-      fRho(2.33),
+      fRho(2.33), // 2.33),
       fBetaGammaMip(3.4601)
   { 
     AddLoad(kKinematics);
@@ -61,8 +66,9 @@ public:
     fPdg                = pdgDB->GetParticle(pdgName);
     if (!fPdg) Warning("DrawHits", "Particle %s not found", pdgName);
 
-    TArrayF tkine(MakeLogScale(n, tmin, tmax));
-    TArrayF eloss(MakeLogScale(m, emin, emax));
+    TArrayF tkine(MakeLogScale(n,    tmin, tmax));
+    TArrayF betag(MakeLogScale(n/10, tmin, tmax));
+    TArrayF eloss(MakeLogScale(m,    emin, emax));
     TString name("elossVsPMQ");
     TString title(Form("#Delta E/#Delta x / q^{2} vs. p/m, %s", 
                       (pdgName ? pdgName : "")));
@@ -74,10 +80,16 @@ public:
     fElossVsPMQ->Sumw2();
     fEloss = new TH1D("eloss", "#Delta E/#Delta x / q^{2}", 
                      eloss.fN-1, eloss.fArray);
-    fEloss->SetFillColor(2);
+    fEloss->SetFillColor(kRed);
     fEloss->SetFillStyle(3001);
     fEloss->SetXTitle("#Delta E/#Delta x / q^{2} [MeV/cm]");
     fEloss->Sumw2();
+
+    fBetaGamma = new TH1D("betaGamma", "#beta#gamma of particles", 
+                         betag.fN-1, betag.fArray);
+    fBetaGamma->SetXTitle("#beta#gamma");
+    fBetaGamma->SetFillColor(kBlue+1);
+    fBetaGamma->SetFillStyle(3001);
   }
   //__________________________________________________________________
   Bool_t ProcessHit(AliFMDHit* hit, TParticle* p) 
@@ -93,22 +105,76 @@ public:
     }
     // if (!p->IsPrimary()) return kTRUE;
     if (hit->IsStop()) return kTRUE;
-    if (hit->Length() == 0) return kTRUE;
-
-    Float_t x = hit->P();
-    Float_t y = hit->Edep()/hit->Length();
+    if (hit->Length() == 0) { 
+      Warning("ProcessHit", "Hit in %s has 0 length", hit->GetName());
+      return kTRUE;
+    }
+    
     Float_t q = hit->Q() / 3.;
     Float_t m = hit->M();
     if (m == 0 || q == 0) return kTRUE;
 
+    TLorentzVector pp;
+    p->Momentum(pp);
+    Double_t betagamma = 0;
+    Info("ProcessHit", "%s (%s) beta=%f", p->GetPDG()->GetName(), 
+        fStack->IsPhysicalPrimary(hit->Track()) ? "primary" : "secondary", 
+        pp.Beta());
+    if (pp.Beta() <= 1 && pp.Beta() >= 0) 
+      betagamma = pp.Beta() * pp.Gamma();
+    fBetaGamma->Fill(betagamma);
+#if 0
+    if (betagamma < 10) { 
+      Info("ProcessHit", "%s (%s) beta=%f gamma=%f beta*gamma=%f", 
+          p->GetPDG()->GetName(), 
+          fStack->IsPhysicalPrimary(hit->Track()) ? "primary" : 
+          "secondary", 
+          pp.Beta(), pp.Gamma(), betagamma);
+      return kTRUE;
+    }
+#endif
+
+    Float_t x = hit->P();
+    Float_t y = hit->Edep()/hit->Length();
+
     x /= hit->M();
-    y /= q * q;
+    // y /= q * q;
     fElossVsPMQ->Fill(x, y);
     fEloss->Fill(y);
     // fNHits++;
     return kTRUE;
   }
   //__________________________________________________________________
+  void ShowFit(Double_t x1, Double_t y1, const char* title, 
+              TF1* f, Double_t dx=0, Double_t dy=0.05)
+  {
+    Double_t x = x1, y = y1;
+    TLatex* latex = new TLatex(x, y, title);
+    latex->SetTextFont(132);
+    latex->SetTextSize(0.8*dy);
+    latex->SetNDC();
+    latex->Draw();
+    x -= dx;
+    y -= dy;
+    const Double_t eqDx=0.1;
+    Double_t chi2 = f->GetChisquare();
+    Int_t    ndf  = f->GetNDF();
+    Double_t prob = f->GetProb();
+    latex->DrawLatex(x, y, "#chi^{2}/NDF");
+    latex->DrawLatex(x+eqDx, y, Form("= %7.4f/%3d=%5.2f (%3d%%)", 
+                                    chi2, ndf, chi2/ndf, int(100*prob)));
+    Int_t     n = f->GetNpar();
+    Double_t* p = f->GetParameters();
+    Double_t* e = f->GetParErrors();
+    for (int i = 0; i < n; i++) { 
+      x -= dx;
+      y -= dy;
+      latex->DrawLatex(x, y, f->GetParName(i));
+      latex->DrawLatex(x+eqDx, y, Form("= %7.4f", p[i]));
+      latex->DrawLatex(x+2*eqDx, y, Form("#pm %7.4f", e[i]));
+    }
+  }      
+  //__________________________________________________________________
   Bool_t Finish()
   {
     gStyle->SetPalette(1);
@@ -120,7 +186,9 @@ public:
     gStyle->SetCanvasBorderSize(0);
     gStyle->SetPadColor(0);
     gStyle->SetPadBorderSize(0);
-    TCanvas* c = new TCanvas("elossVsP", "Energy loss versus momentum");
+    gStyle->SetOptStat(0);
+    TCanvas* c = new TCanvas("elossVsP", "Energy loss versus momentum", 
+                            1200, 800);
     c->SetLogy();
     c->SetLogx();
 
@@ -137,6 +205,7 @@ public:
     // fElossVsPMQ->Draw("ACOLZ same");
     TLegend* l       = new TLegend(.5, .6, .89, .89);
     // l->AddEntry(fElossVsPMQ, fElossVsPMQ->GetTitle(), "pf");
+    l->SetFillColor(0);
     l->AddEntry(mate,        mate->GetTitle(),    "lf");
     l->AddEntry(bb,          bb->GetTitle(),      "l");
     l->AddEntry(nodelta,     nodelta->GetTitle(), "l");
@@ -154,36 +223,114 @@ public:
     c->Modified();
     c->Update();
     c->cd();
+    c->SaveAs("eloss_bethe.png");
 
-    c = new TCanvas("eloss", "Energy loss per unit material");
+    c = new TCanvas("cEloss", "Energy loss per unit material",
+                   1200, 800);
+    c->SetRightMargin(0.05);
+    c->SetTopMargin(0.05);
+    c->SetLeftMargin(0.05);
+    fEloss->SetStats(kFALSE);
     // c->SetLogx();
+    TF1* land     = new TF1("land", "landau");
+    land->SetParNames("A", "MPV", "width");
+    land->SetLineWidth(2);
+    land->SetLineColor(kGreen+1);
+
+    TF1* landgaus = new TF1("landgaus", "gaus(0)+landau(3)");
+    landgaus->SetParNames("A", "#mu", "#sigma", "B", "MPV", "width");
+    landgaus->SetLineWidth(2);
+    landgaus->SetLineColor(kMagenta+1);
+    TGraph*  corr  = GetCorr();
+    TGraph*  resp  = GetResp();
     if (fEloss->GetEntries() != 0) { 
       c->SetLogy();
       fEloss->Scale(1. / fEloss->GetEntries());
       fEloss->GetXaxis()->SetRangeUser(1, 10);
-      fEloss->Fit("landau", "", "", 1, 10);
-      TF1* land = fEloss->GetFunction("landau");
-      land->SetLineWidth(2);
-      Double_t max = fEloss->GetMaximum();
-      TGraph* resp  = GetResp();
-      Double_t* x   = resp->GetX();
-      Double_t* y   = resp->GetY();
-      TGraph*   g   = new TGraph(resp->GetN());
-      TGraph*   co  = GetCorr();
-      std::cout << "Correction factor: " << co->Eval(fBetaGammaMip) 
-               << std::endl;
-      Double_t  xs  = fRho; // * 1.19; // / 
-      for (Int_t i = 0; i < g->GetN(); i++) 
-       g->SetPoint(i, x[i] * xs, y[i] * max);
-      g->Draw("C same");
+      fEloss->Fit(land, "+", "", 2, 10);
+      landgaus->SetParameters(land->GetParameter(0) / 100, 
+                             land->GetParameter(1), 
+                             land->GetParameter(2), 
+                             land->GetParameter(0),
+                             land->GetParameter(1), 
+                             land->GetParameter(2));
+      fEloss->Fit(landgaus, "+", "", 1, 10);
+      fEloss->DrawCopy("HIST SAME");
     }
     
-    l = new TLegend(.6, .6, .89, .89);
+    // fEloss->DrawCopy("E SAME");
+    // land->Draw("same");
+    // landgaus->Draw("same");
+    Double_t max = fEloss->GetMaximum();
+    Double_t* x   = resp->GetX();
+    Double_t* y   = resp->GetY();
+    TGraph*   g   = new TGraph(resp->GetN());
+    g->SetName(Form("%sCorr", resp->GetName()));
+    g->SetTitle(resp->GetTitle());
+    g->SetLineStyle(resp->GetLineStyle());
+    g->SetLineColor(resp->GetLineColor());
+    g->SetLineWidth(resp->GetLineWidth());
+    Double_t  xs2 = corr->Eval(fBetaGammaMip);
+    Double_t xss   = 1.1;
+    Double_t  xs  = fRho * xss;
+    std::cout << "Correction factor: " << xs2 << std::endl;
+    for (Int_t i = 0; i < g->GetN(); i++) 
+      g->SetPoint(i, x[i] * xs, y[i] * max);
+    g->Draw("C same");
+    
+    l = new TLegend(.05, .6, .4, .95);
+    l->SetFillColor(0);
+    l->SetBorderSize(1);
     l->AddEntry(fEloss, fEloss->GetTitle(), "lf");
-    // l->AddEntry(land,   "Landau fit", "l");
-    // l->AddEntry(resp,   "f(#Delta_{p}/x) [RPP fig 27.8]", "l");
+    if (land) 
+      l->AddEntry(land,   Form("Landau fit\t- #chi^{2}/NDF=%7.5f", 
+                              land->GetChisquare()/land->GetNDF()), "l");
+    if (landgaus) 
+      l->AddEntry(landgaus,   
+                 Form("Landau+Gauss fit\t- #chi^{2}/NDF=%7.5f", 
+                      landgaus->GetChisquare()/landgaus->GetNDF()), "l");
+    l->AddEntry(resp,   Form("f(%s#Delta/x) 320#mum Si [RPP fig 27.7]",
+                            xss != 1 ? Form("%4.2f#times", xss) : ""), 
+               "l");
     l->Draw("same");
+    
+    fEloss->GetYaxis()->SetRangeUser(1e-4, 100);
+    ShowFit(0.45,.9, "Landau+Gaus", landgaus, 0, 0.02);
+    ShowFit(0.70,.9, "Landau", land, 0, 0.02);
+
+    c->Modified();
+    c->Update();
+    c->cd();
+    c->SaveAs("eloss_landau.png");
 
+
+    c = new TCanvas("cBetaGamma", "beta gamma", 1200, 800);
+    c->SetLogx();
+    fBetaGamma->Draw();
+    Int_t mipbin = fBetaGamma->FindBin(fBetaGammaMip) + 1;
+    Int_t maxbin = fBetaGamma->GetNbinsX();
+    Int_t total  = fBetaGamma->Integral();
+    Int_t over   = fBetaGamma->Integral(mipbin,maxbin);
+    TH1*  res    = (TH1*)fBetaGamma->Clone("overMip");
+    res->SetFillColor(kRed+1);
+    for (int i = 0; i < mipbin; i++) 
+      res->SetBinContent(i, 0);
+    res->Draw("same");
+    std::cout << "Percentage over MIP : " << float(over) / total << std::endl;
+
+    Double_t yy = fBetaGamma->GetBinContent(mipbin) * 1.1; 
+    a = new TArrow(fBetaGammaMip, 0, fBetaGammaMip, yy, 3, "<|");
+    a->SetAngle(30);
+    a->Draw("same");
+    t = new TLatex(fBetaGammaMip, yy, "Minimum Ionising");
+    t->SetTextSize(0.04);
+    t->SetTextAlign(21);
+    t->Draw("same");
+    c->Modified();
+    c->Update();
+    c->cd();
+    c->SaveAs("eloss_betagamma.png");
+    
     return kTRUE;
   }
 
@@ -212,9 +359,9 @@ public:
       graph->GetHistogram()->SetXTitle("#beta#gamma");
       graph->GetHistogram()->SetYTitle("#Delta E/#Delta x [MeV/cm]");
       graph->SetFillColor(0);
-      graph->SetLineColor(2);
-      graph->SetLineStyle(1);
-      graph->SetLineWidth(1);
+      graph->SetLineColor(kRed+1);
+      graph->SetLineStyle(2);
+      graph->SetLineWidth(2);
       graph->SetName("full_stop");
       graph->SetTitle("Stopping (MeVcm^{2}/g) [RPP fig 27.1]");
       graph->SetPoint(0,0.001461622,40.17542);
@@ -256,9 +403,9 @@ public:
       graph->GetHistogram()->SetYTitle("(MeVcm^{2}/g)");
       graph->GetHistogram()->SetXTitle("#beta#gamma");
       graph->SetFillColor(0);
-      graph->SetLineColor(3);
-      graph->SetLineStyle(1);
-      graph->SetLineWidth(1);
+      graph->SetLineColor(kGreen+1);
+      graph->SetLineStyle(2);
+      graph->SetLineWidth(2);
       graph->SetPoint(0,0.001461622,40.17542);
       graph->SetPoint(1,0.003775053,91.28429);
       graph->SetPoint(2,0.01178769,202.7359);
@@ -298,9 +445,9 @@ public:
       graph->GetHistogram()->SetYTitle("(MeVcm^{2}/g)");
       graph->GetHistogram()->SetXTitle("#beta#gamma");
       graph->SetFillColor(0);
-      graph->SetLineStyle(1);
-      graph->SetLineColor(4);
-      graph->SetLineWidth(1);
+      graph->SetLineColor(kBlue+1);
+      graph->SetLineWidth(2);
+      graph->SetLineStyle(2);
       graph->SetPoint(0,0.001,24.3298);
       graph->SetPoint(1,0.003117649,74.35105);
       graph->SetPoint(2,0.008675042,172.8318);
@@ -337,10 +484,10 @@ public:
                      "electronic only  [RPP fig. 27.6]");
       graph->GetHistogram()->SetYTitle("(MeVcm^{2}/g)");
       graph->GetHistogram()->SetXTitle("#mu E_{kin} (GeV)");
-      graph->SetFillColor(1);
-      graph->SetLineStyle(1);
-      graph->SetLineWidth(1);
-      graph->SetLineColor(1);
+      graph->SetFillColor(0);
+      graph->SetLineStyle(2);
+      graph->SetLineWidth(2);
+      graph->SetLineColor(kMagenta+1);
       graph->SetMarkerStyle(21);
       graph->SetMarkerSize(0.6);
       graph->SetPoint(0,0.1,1.346561);
@@ -373,11 +520,11 @@ public:
       gre->SetTitle("Energy loss 300#mu Si [GFMATE]");
       gre->GetHistogram()->SetXTitle("#beta#gamma");
       gre->GetHistogram()->SetYTitle("#Delta E/#Delta x [MeV/cm]");
-      gre->SetFillColor(6);
-      gre->SetFillStyle(3002);
-      gre->SetLineColor(1);
+      gre->SetFillColor(kGray);
+      gre->SetFillStyle(3001); // 0 /* 3002 */);
+      gre->SetLineColor(kGray+1);
       gre->SetLineStyle(1);
-      gre->SetLineWidth(1);
+      gre->SetLineWidth(2);
       gre->SetPoint(0,7.16486e-05,1218.84);
       gre->SetPoint(1,9.25378e-05,1221.38);
       gre->SetPoint(2,0.000119517,1180.12);
@@ -474,25 +621,79 @@ public:
       for (Int_t i = 0; i < gre->GetN(); i++) 
        gre->SetPointError(i, 0, 2 * 0.1 * y[i]); // ! 1 sigma
     }
-    gre->DrawClone("c3 same");
+    gre->Draw("c3 same");
     return gre;
   }
 
   /** Get the response functin @f$ f(\Delta_p/x)@f$ from Review of
-      Particle Physics (fig. 27.2).  It is scaled to the value at
+      Particle Physics (fig. 27.7).  It is scaled to the value at
       MPV. */ 
   TGraph* GetResp()
   {
     static TGraph*  graph = 0;
     if (!graph) {
-      graph = new TGraph(16);
+      graph = new TGraph;
       graph->SetName("si_resp");
-      graph->SetTitle("f(#Delta_{p}/x) scaled to the MPV value ");
-      graph->GetHistogram()->SetXTitle("#Delta_{p}/x (MeVcm^{2}/g)");
-      graph->GetHistogram()->SetYTitle("f(#Delta_{p}/x)");
-      graph->SetFillColor(1);
+      graph->SetTitle("f(#Delta/x) scaled to the MPV value ");
+      graph->GetHistogram()->SetXTitle("#Delta/x (MeVcm^{2}/g)");
+      graph->GetHistogram()->SetYTitle("f(#Delta/x)");
+      graph->SetLineColor(kBlue+1);
+      graph->SetLineWidth(2);
+      graph->SetFillColor(kBlue+1);
       graph->SetMarkerStyle(21);
       graph->SetMarkerSize(0.6);
+#if 0
+      // Figure 27.7 or Review of Particle physics - Straggeling function in 
+      // silicon of 500 MeV Pions, normalised to unity at the most probable 
+      // value.   
+      graph->SetPoint(0,0.808094,0.00377358);
+      graph->SetPoint(1,0.860313,0.0566038);
+      graph->SetPoint(2,0.891645,0.116981);
+      graph->SetPoint(3,0.912533,0.181132);
+      graph->SetPoint(4,0.928198,0.260377);
+      graph->SetPoint(5,0.938642,0.320755);
+      graph->SetPoint(6,0.954308,0.377358);
+      graph->SetPoint(7,0.964752,0.433962);
+      graph->SetPoint(8,0.975196,0.490566);
+      graph->SetPoint(9,0.98564,0.550943);
+      graph->SetPoint(10,0.996084,0.611321);
+      graph->SetPoint(11,1.00653,0.667925);
+      graph->SetPoint(12,1.02219,0.732075);
+      graph->SetPoint(13,1.03264,0.784906);
+      graph->SetPoint(14,1.0483,0.845283);
+      graph->SetPoint(15,1.06397,0.901887);
+      graph->SetPoint(16,1.09008,0.958491);
+      graph->SetPoint(17,1.10574,0.984906);
+      graph->SetPoint(18,1.13708,1);
+      graph->SetPoint(19,1.13708,1);
+      graph->SetPoint(20,1.15796,0.988679);
+      graph->SetPoint(21,1.17363,0.966038);
+      graph->SetPoint(22,1.19974,0.916981);
+      graph->SetPoint(23,1.2154,0.89434);
+      graph->SetPoint(24,1.23629,0.837736);
+      graph->SetPoint(25,1.2624,0.784906);
+      graph->SetPoint(26,1.28329,0.724528);
+      graph->SetPoint(27,1.3094,0.664151);
+      graph->SetPoint(28,1.32507,0.611321);
+      graph->SetPoint(29,1.3564,0.550943);
+      graph->SetPoint(30,1.41384,0.445283);
+      graph->SetPoint(31,1.44517,0.392453);
+      graph->SetPoint(32,1.48695,0.335849);
+      graph->SetPoint(33,1.52872,0.286792);
+      graph->SetPoint(34,1.58094,0.237736);
+      graph->SetPoint(35,1.63838,0.196226);
+      graph->SetPoint(36,1.68016,0.169811);
+      graph->SetPoint(37,1.75326,0.135849);
+      graph->SetPoint(38,1.81593,0.113208);
+      graph->SetPoint(39,1.89426,0.0981132);
+      graph->SetPoint(40,1.96214,0.0830189);
+      graph->SetPoint(41,2.0718,0.0641509);
+      graph->SetPoint(42,2.19191,0.0490566);
+      graph->SetPoint(43,2.31723,0.0415094);
+      graph->SetPoint(44,2.453,0.0301887);
+      graph->SetPoint(45,2.53133,0.0264151);
+      graph->SetPoint(46,2.57833,0.0264151);
+#else
       graph->SetPoint(0,0.8115124,0.009771987);
       graph->SetPoint(1,0.9198646,0.228013);
       graph->SetPoint(2,0.996614,0.5895765);
@@ -509,6 +710,7 @@ public:
       graph->SetPoint(13,1.985327,0.08143322);
       graph->SetPoint(14,2.301354,0.04234528);
       graph->SetPoint(15,2.56772,0.02931596);
+#endif
     }
     return graph;
   }