Bug fix for reading raw data - thanks Frederic YERMIA, Diego, and Gines MARTINEZ
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 14 Sep 2007 13:06:51 +0000 (13:06 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 14 Sep 2007 13:06:51 +0000 (13:06 +0000)
13 files changed:
FMD/AliFMDDisplay.h
FMD/AliFMDPattern.h
FMD/AliFMDRawReader.cxx
FMD/AliFMDRawStream.cxx
FMD/Config.C
FMD/Simulate.C
FMD/scripts/DisplayDigits.C
FMD/scripts/DrawESD.C [new file with mode: 0644]
FMD/scripts/DrawHits.C
FMD/scripts/DrawHitsRecs.C
FMD/scripts/GetXsection.C
FMD/scripts/XSection.C [new file with mode: 0644]
FMD/scripts/checkSizes.sh [new file with mode: 0755]

index 159a04e..3c1f1eb 100644 (file)
@@ -128,6 +128,7 @@ protected:
       fHits(0),
       fCanvas(0),
       fPad(0),
+      fButtons(0),
       fSlider(0),
       fZoomMode(0),
       fX0(0),
@@ -140,7 +141,11 @@ protected:
       fYPixel(0),
       fOldXPixel(0),
       fOldYPixel(0),
-      fLineDrawn(0)
+      fLineDrawn(0), 
+      fOnlyFMD(kTRUE),
+      fSpec(0), 
+      fSpecCut(0),
+      fAux(0)
   { } 
   /** Assignment operator 
       @return Reference to this object */
index b29c5e0..bca5489 100644 (file)
@@ -79,7 +79,7 @@ public:
   private:
     /** Copy constructor 
        - Not implemented. */ 
-    AliFMDPatternDetector(const AliFMDPattern&);
+    AliFMDPatternDetector(const AliFMDPatternDetector&);
     /** Assignement operator 
        -- Not implemented */
     AliFMDPatternDetector& operator=(const AliFMDPatternDetector&);
index a34175e..dcf488e 100644 (file)
@@ -122,7 +122,11 @@ AliFMDRawReader::ReadAdcs(TClonesArray* array)
   UInt_t hwaddr = 0;
   // Data array is approx twice the size needed. 
   UShort_t data[2048];
-  while (input.ReadChannel(ddl, hwaddr, last, data)) {
+
+  Bool_t isGood = kTRUE;
+  while (isGood) {
+    isGood = input.ReadChannel(ddl, hwaddr, last, data);
+
     AliFMDDebug(5, ("Read channel 0x%x of size %d", hwaddr, last));
     UShort_t det, sec, str;
     Char_t   ring;
index 4ad4ff6..89b4b83 100644 (file)
@@ -64,7 +64,13 @@ AliFMDRawStream::ReadChannel(UInt_t& ddl, UInt_t& addr,
     if (last > 0x3FF) {
       AliFMDDebug(30, ("Last is 0x%x, so reading a new word", last));
       next   = Next();
-      if (!next) break;
+      if(!next){
+       addr = GetPrevHWAddress();
+       ddl  = GetPrevDDLNumber();
+       len  = l+1; // Need to add one - l points to last valid index
+       last = signal;
+       break;
+      }
       signal = GetSignal();
       if (GetHWAddress() != GetPrevHWAddress() && GetPrevHWAddress() >= 0) {
        AliFMDDebug(15, ("New hardware address, was 0x%x, now 0x%x", 
index 6317aeb..240d320 100644 (file)
@@ -386,7 +386,7 @@ Config()
   gMC->SetProcess("MUNU",1);
   gMC->SetProcess("CKOV",1);
   gMC->SetProcess("HADR",1);
-  gMC->SetProcess("LOSS",2);
+  gMC->SetProcess("LOSS",2); // 0:none 1,3:dray 2:nodray 4:nofluct (def:2)
   gMC->SetProcess("MULS",1);
   gMC->SetProcess("RAYL",1);
 
index 5273269..9d38d89 100644 (file)
@@ -24,7 +24,7 @@ void
 Simulate()
 {
   AliSimulation sim;
-  AliLog::SetModuleDebugLevel("FMD", 1);
+  // AliLog::SetModuleDebugLevel("FMD", 1);
   sim.SetConfigFile("$(ALICE_ROOT)/FMD/Config.C");
   // sim.SetMakeSDigits("FMD");
   sim.SetMakeDigits("FMD"); 
@@ -32,7 +32,7 @@ Simulate()
   // sim.SetMakeDigitsFromHits("FMD"); 
   TStopwatch w; 
   w.Start(); 
-  sim.Run(1);  
+  sim.Run(10);  
   w.Stop(); 
   w.Print(); 
 }
index 5b293a4..5f162b4 100644 (file)
@@ -12,6 +12,7 @@ DisplayDigits()
 {
   AliCDBManager* cdb = AliCDBManager::Instance();
   cdb->SetDefaultStorage("local://$ALICE_ROOT");
+  cdb->SetRun(0);
   gSystem->Load("libFMDutil.so");
   AliFMDDisplay* d = new AliFMDDisplay;
   d->AddLoad(AliFMDInput::kDigits);
diff --git a/FMD/scripts/DrawESD.C b/FMD/scripts/DrawESD.C
new file mode 100644 (file)
index 0000000..1e24aa0
--- /dev/null
@@ -0,0 +1,223 @@
+//____________________________________________________________________
+//
+// $Id$
+//
+// Script that contains a class to draw eloss from hits, versus ADC
+// counts from digits, using the AliFMDInputHits class in the util library. 
+//
+// It draws the energy loss versus the p/(mq^2).  It can be overlayed
+// with the Bethe-Bloc curve to show how the simulation behaves
+// relative to the expected. 
+//
+// Use the script `Compile.C' to compile this class using ACLic. 
+//
+#include <TH1D.h>
+#include <AliFMDHit.h>
+#include <AliFMDDigit.h>
+#include <AliFMDInput.h>
+#include <AliFMDUShortMap.h>
+#include <AliFMDFloatMap.h>
+#include <AliFMDRecPoint.h>
+#include <AliESDFMD.h>
+#include <AliLog.h>
+#include <iostream>
+#include <TStyle.h>
+#include <TArrayF.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <TF1.h>
+#include <TSpectrum.h>
+#include <TLegend.h>
+#include <TLine.h>
+
+/** @class DrawESD
+    @brief Draw digit ADC versus Rec point mult
+    @code 
+    Root> .L Compile.C
+    Root> Compile("DrawESD.C")
+    Root> DrawESD c
+    Root> c.Run();
+    @endcode
+    @ingroup FMD_script
+ */
+class DrawESD : public AliFMDInput
+{
+private:
+  TH1D* fMult; // Histogram 
+  const Double_t fCorr;
+public:
+  //__________________________________________________________________
+  TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
+  {
+    TArrayF bins(n+1);
+    Float_t dp   = n / TMath::Log10(max / min);
+    Float_t pmin = TMath::Log10(min);
+    bins[0]      = min;
+    for (Int_t i = 1; i < n+1; i++) {
+      Float_t p = pmin + i / dp;
+      bins[i]   = TMath::Power(10, p);
+    }
+    return bins;
+  }
+  //__________________________________________________________________
+  DrawESD(Int_t n=420, Double_t mmin=-0.5, Double_t mmax=20.5) 
+    : fCorr(1) // 0.68377 / 1.1)
+  { 
+    AddLoad(kESD);
+    fMult = new TH1D("mult", " Multiplicity (strip)", n, mmin, mmax);
+    fMult->SetXTitle("Strip Multiplicity");
+  }
+  //__________________________________________________________________
+  /** Begining of event
+      @param ev Event number
+      @return @c false on error */
+  Bool_t Begin(Int_t ev) 
+  {
+    return AliFMDInput::Begin(ev);
+  }
+  //__________________________________________________________________
+  Bool_t ProcessESD(UShort_t /* det */, 
+                   Char_t   /* ring */, 
+                   UShort_t /* sec */, 
+                   UShort_t /* strip */, 
+                   Float_t  /* eta */, 
+                   Float_t  mult)
+  {
+    // Cache the energy loss 
+    fMult->Fill(mult/fCorr);
+    return kTRUE;
+  }
+  //__________________________________________________________________
+  Bool_t Finish()
+  {
+    gStyle->SetPalette(1);
+    gStyle->SetOptTitle(0);
+    gStyle->SetOptFit(1111111);
+    gStyle->SetCanvasColor(0);
+    gStyle->SetCanvasBorderSize(0);
+    gStyle->SetPadColor(0);
+    gStyle->SetPadBorderSize(0);
+    
+    TCanvas* c = new TCanvas("c", "C");
+    c->cd();
+    c->SetLogy();
+    // fMult->GetXaxis()->SetRangeUser(0.2,4);
+    // fMult->Sumw2();
+    // fMult->Scale(1. / fMult->GetMaximum());
+    fMult->SetStats(kFALSE);
+    fMult->SetFillColor(2);
+    fMult->SetFillStyle(3001);
+    fMult->Draw();
+    Double_t x1 = .75; // .8;  // .65 / fCorr;
+    Double_t x2 = 1.3; // 1.7; // fCorr;
+    Double_t x3 = 2.5; // 2.7; // fCorr;
+    Double_t x4 = 3.7; // fCorr;
+    TF1* l1 = new TF1("landau1", "landau", x1, x2);
+    TF1* l2 = new TF1("landau2", "landau", x2, x3);
+    // TF1* l3 = new TF1("landau3", "landau", x3, x4);
+    TF1* f  = new TF1("user", "landau(0)+landau(3)", x1, x3);
+
+    fMult->SetStats(kTRUE);
+    fMult->GetXaxis()->SetRangeUser(0, 4);
+    fMult->Fit(l1, "E0L", "", x1, x2);
+    fMult->Fit(l2, "E0L", "", x2, x3);
+    // fMult->Fit(l3, "E0L", "", x3, x4);
+    f->SetParNames("A_{1}", "Mpv_{1}", "#sigma_{1}",
+                  "A_{2}", "Mpv_{2}", "#sigma_{2}",
+                  "A_{3}", "Mpv_{3}", "#sigma_{3}");
+    f->SetParameters(l1->GetParameter(0), 
+                    l1->GetParameter(1), 
+                    l1->GetParameter(2), 
+                    l2->GetParameter(0), 
+                    l2->GetParameter(1), 
+                    l2->GetParameter(2),
+                    l1->GetParameter(0), 
+                    l1->GetParameter(1) * 3, 
+                    l1->GetParameter(2) * 3);
+    fMult->Fit(f, "E0L", "", x1, x3);
+    fMult->Fit(f, "MEL", "E1", x1, x3);
+
+    TH1* h = static_cast<TH1*>(fMult->DrawClone("same hist"));
+
+    // l1->SetLineStyle(2);
+    l1->SetLineColor(3);
+    l1->SetLineWidth(2);
+    l1->SetRange(0, 4);
+    l1->Draw("same");
+    // l2->SetLineStyle(3);
+    l2->SetLineColor(4);
+    l2->SetLineWidth(2);
+    l2->SetRange(0, 4);
+    l2->Draw("same");
+    // l3->SetLineStyle(3);
+    // l3->SetLineWidth(2);
+    // l3->SetRange(0, 5);
+    // l3->Draw("same");
+    f->SetLineWidth(2);
+    f->SetRange(0, 4);
+    f->Draw("same");
+
+    TLegend* l = new TLegend(0.2, 0.6, .6, .89);
+    l->AddEntry(l1, "1 particle Landau", "l");
+    l->AddEntry(l2, "2 particle Landau", "l");
+    l->AddEntry(f,  "1+2 particle Landau", "l");
+    l->SetFillColor(0);
+    l->Draw("same");
+
+
+    c = new TCanvas("c2", "Landaus");
+    c->SetLogy();
+    fMult->DrawClone("axis");
+    f->Draw("same");
+    TF1* ll1 = new TF1("ll1", "landau", 0, 4);
+    ll1->SetParameters(f->GetParameter(0), 
+                      f->GetParameter(1), 
+                      f->GetParameter(2));
+    ll1->SetLineColor(3);
+    ll1->Draw("same");
+    TF1* ll2 = new TF1("ll2", "landau", 0, 4);
+    ll2->SetParameters(f->GetParameter(3), 
+                      f->GetParameter(4), 
+                      f->GetParameter(5));
+    ll2->SetLineColor(4);
+    ll2->Draw("same");
+
+    Double_t y1  = fMult->GetMinimum() * 1.1;
+    Double_t y2  = fMult->GetMaximum() * .9;
+    Double_t xc1 = f->GetParameter(1)-3*f->GetParameter(2);
+    Double_t xc2 = f->GetParameter(4)-2*f->GetParameter(5);
+    TLine* c1 = new TLine(xc1, y1, xc1, y2);
+    c1->Draw("same");
+    TLine* c2 = new TLine(xc2, y1, xc2, y2);
+    c2->Draw("same");
+
+    l = new TLegend(0.2, 0.6, .6, .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");
+
+#if 0
+    c = new TCanvas("s", "Spectrum");
+    TSpectrum* s = new TSpectrum(16);
+    h->GetXaxis()->SetRangeUser(0.3, 20);
+    s->Search(h, .15, "", 0.01);
+    c->Update();
+    TH1* b = s->Background(h);
+    b->SetFillColor(4);
+    b->SetFillStyle(3001);
+    b->Draw("same");
+#endif
+
+    return kTRUE;
+  }
+
+  ClassDef(DrawESD,0);
+  
+};
+
+//____________________________________________________________________
+//
+// EOF
+//
index 98e3a39..6845cf4 100644 (file)
 #include <TParticle.h>
 #include <TCanvas.h>
 #include <TGraphErrors.h>
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
+#include <TLegend.h>
+#include <TArrow.h>
+#include <TLatex.h>
+#include <TF1.h>
 
 /** @class DrawHits
     @brief Draw hit energy loss
@@ -35,7 +41,10 @@ class DrawHits : public AliFMDInput
 {
 private:
   TH2D* fElossVsPMQ; // Histogram 
-  Int_t fPdg;
+  TH1D* fEloss;
+  TParticlePDG* fPdg;
+  const Double_t fRho;
+  const Double_t fBetaGammaMip;
 public:
   //__________________________________________________________________
   TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
@@ -55,23 +64,34 @@ public:
     return bins;
   }
   //__________________________________________________________________
-  DrawHits(Int_t m=1000, Double_t emin=1, Double_t emax=1000, 
-          Int_t n=900, Double_t tmin=1e-2, Double_t tmax=1e3, 
-          Int_t pdg=211) 
-    : fPdg(pdg)
+  DrawHits(const char* pdgName="pi+",
+          Int_t m=1000, Double_t emin=1, Double_t emax=1000, 
+          Int_t n=900, Double_t tmin=1e-2, Double_t tmax=1e3) 
+    : fPdg(0), 
+      fRho(2.33),
+      fBetaGammaMip(3.4601)
   { 
     AddLoad(kKinematics);
     AddLoad(kHits);
+    TDatabasePDG* pdgDB = TDatabasePDG::Instance();
+    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));
-    TString name(Form("elossVsP%s", (fPdg == 0 ? "MQ" : "")));
-    TString title(Form("#Delta E/#Delta x vs. p%s", 
-                      (fPdg == 0 ? "/(mq^{2})" : "")));
+    TString name("elossVsPMQ");
+    TString title(Form("#Delta E/#Delta x / q^{2} vs. p/m, %s", 
+                      (pdgName ? pdgName : "")));
     fElossVsPMQ = new TH2D(name.Data(), title.Data(), 
                           tkine.fN-1, tkine.fArray, 
                           eloss.fN-1, eloss.fArray);
-    fElossVsPMQ->SetXTitle(Form("p%s", (fPdg == 0 ? "/(mq^{2})" : " [GeV]")));
-    fElossVsPMQ->SetYTitle("#Delta E/#Delta x [MeV/cm]");
+    fElossVsPMQ->SetXTitle("p/(mq^{2})=#beta#gamma/q^{2}");
+    fElossVsPMQ->SetYTitle("#Delta E/#Delta x / q^{2} [MeV/cm]");
+    fEloss = new TH1D("eloss", "#Delta E/#Delta x / q^{2}", 
+                     eloss.fN-1, eloss.fArray);
+    fEloss->SetFillColor(2);
+    fEloss->SetFillStyle(3001);
+    fEloss->SetXTitle("#Delta E/#Delta x / q^{2} [MeV/cm]");
   }
   //__________________________________________________________________
   Bool_t ProcessHit(AliFMDHit* hit, TParticle* p) 
@@ -87,235 +107,452 @@ public:
     }
     if (!p->IsPrimary()) return kTRUE;
     if (hit->IsStop()) return kTRUE;
+
     Float_t x = hit->P();
     Float_t y = hit->Edep()/hit->Length();
-    if (fPdg != 0) {
-      if (TMath::Abs(hit->Pdg()) != fPdg) return kTRUE;
-    }
-    else {
-      Float_t q = hit->Q() / 3.;
-      if (hit->M() == 0 || q == 0) return kTRUE;
-      x /= hit->M();
-      y /= q * q;
-    }
+    Float_t q = hit->Q() / 3.;
+    Float_t m = hit->M();
+    if (m == 0 || q == 0) return kTRUE;
+
+    x /= hit->M();
+    y /= q * q;
     fElossVsPMQ->Fill(x, y);
+    fEloss->Fill(y);
+
     return kTRUE;
   }
   //__________________________________________________________________
   Bool_t Finish()
   {
-    TCanvas* c = new TCanvas("c", "C");
-    c->SetLogy();
-    c->SetLogx();
     gStyle->SetPalette(1);
-    gStyle->SetOptTitle(0);
+    // gStyle->SetOptTitle(0);
+    gStyle->SetTitleBorderSize(1);
+    gStyle->SetTitleFillColor(0);
+    gStyle->SetCanvasColor(0);
     gStyle->SetCanvasColor(0);
     gStyle->SetCanvasBorderSize(0);
     gStyle->SetPadColor(0);
     gStyle->SetPadBorderSize(0);
+    TCanvas* c = new TCanvas("elossVsP", "Energy loss versus momentum");
+    c->SetLogy();
+    c->SetLogx();
+
+    TString title(fElossVsPMQ->GetTitle());
+    title.Append(Form(", %d events", fEventCount));
+    fElossVsPMQ->SetTitle(title.Data());
     fElossVsPMQ->SetStats(kFALSE);
-    fElossVsPMQ->Draw("COLZ box");
+    fElossVsPMQ->Draw("AXIS");
+    fElossVsPMQ->Draw("ACOLZ same");
+    TGraph*  mate    = FromGFMATE();
+    TGraph*  bb      = FromRPPFull();
+    TGraph*  nodelta = FromRPPNoDelta();
+    TGraph*  norad   = FromRPPNoRad();
+    TGraph*  mean    = FromRPPMean();
+    // fElossVsPMQ->Draw("ACOLZ same");
+    TLegend* l       = new TLegend(.5, .6, .89, .89);
+    // l->AddEntry(fElossVsPMQ, fElossVsPMQ->GetTitle(), "pf");
+    l->AddEntry(mate,        mate->GetTitle(),    "lf");
+    l->AddEntry(bb,          bb->GetTitle(),      "l");
+    l->AddEntry(nodelta,     nodelta->GetTitle(), "l");
+    l->AddEntry(norad,       norad->GetTitle(),   "l");
+    l->AddEntry(mean,        mean->GetTitle(),    "l");
+    l->Draw("same");
+    Double_t min = fElossVsPMQ->GetYaxis()->GetFirst();
+    TArrow* a = new TArrow(fBetaGammaMip, min, fBetaGammaMip, 35*min, 03, "<|");
+    a->SetAngle(30);
+    a->Draw("same");
+    TLatex* t = new TLatex(fBetaGammaMip, 40*min, "Minimum Ionising");
+    t->SetTextSize(0.04);
+    t->SetTextAlign(21);
+    t->Draw("same");
     c->Modified();
     c->Update();
     c->cd();
+
+    c = new TCanvas("eloss", "Energy loss per unit material");
+    // c->SetLogx();
+    c->SetLogy();
+    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");
+
+    l = new TLegend(.6, .6, .89, .89);
+    l->AddEntry(fEloss, fEloss->GetTitle(), "lf");
+    l->AddEntry(land,   "Landau fit", "l");
+    l->AddEntry(resp,   "f(#Delta_{p}/x) [RPP fig 27.8]");
+    l->Draw("same");
+
     return kTRUE;
   }
-  void SuperImposeBetheBloc() 
+
+  /** Scale a graph by density (multiply) and mass (divide). 
+      @param graph Graph to scale 
+      @param density If @c true, scale by the Si density
+      ($\rho=2.33$/cm^3$).  The y axis is assumed to have units of
+      $MeVg^{-1}cm^2$. 
+      @param mass Mass to scale with. The x axis is assumed to be the
+      kinetic energy of the particles in units of $GeV$.  */
+  void ScaleGraph(TGraph* graph, bool density=true, double mass=1) 
+  {
+      Double_t*      x   = graph->GetX();
+      Double_t*      y   = graph->GetY();
+      const Double_t rho = (density ? fRho : 1);
+      for (Int_t i = 0; i < graph->GetN(); i++) 
+       graph->SetPoint(i, x[i] / mass, y[i] * rho); 
+  }    
+  /** Draw pure Bethe-Bloc from Review of Particle Physics, fig. 27.1 
+      @return TGraph object */ 
+  TGraph* FromRPPFull() 
+  {
+    static TGraph* graph = 0;
+    if (!graph) { 
+      graph = new TGraph(20);
+      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->SetName("full_stop");
+      graph->SetTitle("Stopping (MeVcm^{2}/g) [RPP fig 27.1]");
+      graph->SetPoint(0,0.001461622,40.17542);
+      graph->SetPoint(1,0.003775053,91.28429);
+      graph->SetPoint(2,0.01178769,202.7359);
+      graph->SetPoint(3,0.01722915,212.1938);
+      graph->SetPoint(4,0.03162278,172.8318);
+      graph->SetPoint(5,0.06028646,91.28429);
+      graph->SetPoint(6,0.09506529,51.62633);
+      graph->SetPoint(7,0.433873,5.281682);
+      graph->SetPoint(8,1.255744,1.808947);
+      graph->SetPoint(9,2.393982,1.440177);
+      graph->SetPoint(10,3.499097,1.407715);
+      graph->SetPoint(11,10.92601,1.542122);
+      graph->SetPoint(12,60.28646,1.85066);
+      graph->SetPoint(13,236.3885,2.121938);
+      graph->SetPoint(14,468.0903,2.324538);
+      graph->SetPoint(15,1208.976,2.987085);
+      graph->SetPoint(16,6670.768,7.961412);
+      graph->SetPoint(17,23341.67,24.3298);
+      graph->SetPoint(18,110651.2,104.6651);
+      graph->SetPoint(19,264896.9,260.5203);
+      ScaleGraph(graph);
+    }
+    graph->Draw("C same");
+    return graph;
+  }
+
+  /** Draw pure Bethe-Bloc from Review of Particle Physics, fig. 27.1,
+      but without delta electrons  
+      @return TGraph object */ 
+  TGraph* FromRPPNoDelta() 
   {
-    // This is for pi+, made with MakeXsection.C and DrawXsection.C
-    TGraphErrors *gre = new TGraphErrors(91);
-    gre->SetName("BetheBlocPi");
-    gre->SetTitle("BetheBlocPi");
-    gre->SetFillColor(6);
-    gre->SetFillStyle(3001);
-    gre->SetLineWidth(2);
-    gre->SetPoint(0,7.16486e-05,1218.84);
-    gre->SetPointError(0,0,609.418);
-    gre->SetPoint(1,9.25378e-05,1221.38);
-    gre->SetPointError(1,0,610.689);
-    gre->SetPoint(2,0.000119517,1180.12);
-    gre->SetPointError(2,0,590.058);
-    gre->SetPoint(3,0.000154362,1100.31);
-    gre->SetPointError(3,0,550.156);
-    gre->SetPoint(4,0.000199367,996.621);
-    gre->SetPointError(4,0,498.31);
-    gre->SetPoint(5,0.000257492,886.005);
-    gre->SetPointError(5,0,443.003);
-    gre->SetPoint(6,0.000332563,780.483);
-    gre->SetPointError(6,0,390.241);
-    gre->SetPoint(7,0.000429522,684.927);
-    gre->SetPointError(7,0,342.463);
-    gre->SetPoint(8,0.000554749,599.407);
-    gre->SetPointError(8,0,299.703);
-    gre->SetPoint(9,0.000716486,522.375);
-    gre->SetPointError(9,0,261.187);
-    gre->SetPoint(10,0.000925378,452.497);
-    gre->SetPointError(10,0,226.249);
-    gre->SetPoint(11,0.00119517,389.101);
-    gre->SetPointError(11,0,194.551);
-    gre->SetPoint(12,0.00154362,331.974);
-    gre->SetPointError(12,0,165.987);
-    gre->SetPoint(13,0.00199367,280.969);
-    gre->SetPointError(13,0,140.485);
-    gre->SetPoint(14,0.00257492,235.689);
-    gre->SetPointError(14,0,117.844);
-    gre->SetPoint(15,0.00332564,196.156);
-    gre->SetPointError(15,0,98.078);
-    gre->SetPoint(16,0.00429522,162.402);
-    gre->SetPointError(16,0,81.2012);
-    gre->SetPoint(17,0.00554749,133.87);
-    gre->SetPointError(17,0,66.9351);
-    gre->SetPoint(18,0.00716486,109.959);
-    gre->SetPointError(18,0,54.9797);
-    gre->SetPoint(19,0.00925378,90.2035);
-    gre->SetPointError(19,0,45.1017);
-    gre->SetPoint(20,0.0119517,74.1317);
-    gre->SetPointError(20,0,37.0658);
-    gre->SetPoint(21,0.0154362,60.8988);
-    gre->SetPointError(21,0,30.4494);
-    gre->SetPoint(22,0.0199367,49.9915);
-    gre->SetPointError(22,0,24.9957);
-    gre->SetPoint(23,0.0257492,40.9812);
-    gre->SetPointError(23,0,20.4906);
-    gre->SetPoint(24,0.0332564,33.5739);
-    gre->SetPointError(24,0,16.7869);
-    gre->SetPoint(25,0.0429522,27.5127);
-    gre->SetPointError(25,0,13.7563);
-    gre->SetPoint(26,0.0554749,22.5744);
-    gre->SetPointError(26,0,11.2872);
-    gre->SetPoint(27,0.0716486,18.5674);
-    gre->SetPointError(27,0,9.28372);
-    gre->SetPoint(28,0.0925378,15.3292);
-    gre->SetPointError(28,0,7.66462);
-    gre->SetPoint(29,0.119517,12.7231);
-    gre->SetPointError(29,0,6.36156);
-    gre->SetPoint(30,0.154362,10.6352);
-    gre->SetPointError(30,0,5.31759);
-    gre->SetPoint(31,0.199367,8.97115);
-    gre->SetPointError(31,0,4.48558);
-    gre->SetPoint(32,0.257492,7.65358);
-    gre->SetPointError(32,0,3.82679);
-    gre->SetPoint(33,0.332564,6.61909);
-    gre->SetPointError(33,0,3.30955);
-    gre->SetPoint(34,0.429522,5.81614);
-    gre->SetPointError(34,0,2.90807);
-    gre->SetPoint(35,0.554749,5.20286);
-    gre->SetPointError(35,0,2.60143);
-    gre->SetPoint(36,0.716486,4.74533);
-    gre->SetPointError(36,0,2.37267);
-    gre->SetPoint(37,0.925378,4.40987);
-    gre->SetPointError(37,0,2.20494);
-    gre->SetPoint(38,1.19517,4.17077);
-    gre->SetPointError(38,0,2.08538);
-    gre->SetPoint(39,1.54362,4.014);
-    gre->SetPointError(39,0,2.007);
-    gre->SetPoint(40,1.99367,3.92577);
-    gre->SetPointError(40,0,1.96288);
-    gre->SetPoint(41,2.57492,3.89199);
-    gre->SetPointError(41,0,1.946);
-    gre->SetPoint(42,3.32564,3.90063);
-    gre->SetPointError(42,0,1.95032);
-    gre->SetPoint(43,4.29522,3.94146);
-    gre->SetPointError(43,0,1.97073);
-    gre->SetPoint(44,5.54749,4.00597);
-    gre->SetPointError(44,0,2.00299);
-    gre->SetPoint(45,7.16486,4.08725);
-    gre->SetPointError(45,0,2.04362);
-    gre->SetPoint(46,9.25378,4.17985);
-    gre->SetPointError(46,0,2.08993);
-    gre->SetPoint(47,11.9517,4.27962);
-    gre->SetPointError(47,0,2.13981);
-    gre->SetPoint(48,15.4362,4.38347);
-    gre->SetPointError(48,0,2.19174);
-    gre->SetPoint(49,19.9367,4.48919);
-    gre->SetPointError(49,0,2.2446);
-    gre->SetPoint(50,25.7492,4.59523);
-    gre->SetPointError(50,0,2.29762);
-    gre->SetPoint(51,33.2564,4.70052);
-    gre->SetPointError(51,0,2.35026);
-    gre->SetPoint(52,42.9522,4.80435);
-    gre->SetPointError(52,0,2.40218);
-    gre->SetPoint(53,55.4749,4.90625);
-    gre->SetPointError(53,0,2.45312);
-    gre->SetPoint(54,71.6486,5.00589);
-    gre->SetPointError(54,0,2.50295);
-    gre->SetPoint(55,92.5378,5.10279);
-    gre->SetPointError(55,0,2.55139);
-    gre->SetPoint(56,119.517,5.19654);
-    gre->SetPointError(56,0,2.59827);
-    gre->SetPoint(57,154.362,5.28758);
-    gre->SetPointError(57,0,2.64379);
-    gre->SetPoint(58,199.367,5.37581);
-    gre->SetPointError(58,0,2.6879);
-    gre->SetPoint(59,257.492,5.46109);
-    gre->SetPointError(59,0,2.73055);
-    gre->SetPoint(60,332.564,5.54335);
-    gre->SetPointError(60,0,2.77167);
-    gre->SetPoint(61,429.522,5.62248);
-    gre->SetPointError(61,0,2.81124);
-    gre->SetPoint(62,554.749,5.69843);
-    gre->SetPointError(62,0,2.84922);
-    gre->SetPoint(63,716.486,5.77122);
-    gre->SetPointError(63,0,2.88561);
-    gre->SetPoint(64,925.378,5.84093);
-    gre->SetPointError(64,0,2.92046);
-    gre->SetPoint(65,1195.17,5.9077);
-    gre->SetPointError(65,0,2.95385);
-    gre->SetPoint(66,1543.62,5.97165);
-    gre->SetPointError(66,0,2.98582);
-    gre->SetPoint(67,1993.67,6.03292);
-    gre->SetPointError(67,0,3.01646);
-    gre->SetPoint(68,2574.92,6.09171);
-    gre->SetPointError(68,0,3.04586);
-    gre->SetPoint(69,3325.64,6.14827);
-    gre->SetPointError(69,0,3.07413);
-    gre->SetPoint(70,4295.22,6.20286);
-    gre->SetPointError(70,0,3.10143);
-    gre->SetPoint(71,5547.49,6.25577);
-    gre->SetPointError(71,0,3.12788);
-    gre->SetPoint(72,7164.86,6.30725);
-    gre->SetPointError(72,0,3.15363);
-    gre->SetPoint(73,9253.78,6.35757);
-    gre->SetPointError(73,0,3.17878);
-    gre->SetPoint(74,11951.7,6.38446);
-    gre->SetPointError(74,0,3.19223);
-    gre->SetPoint(75,15436.2,6.38446);
-    gre->SetPointError(75,0,3.19223);
-    gre->SetPoint(76,19936.7,6.38446);
-    gre->SetPointError(76,0,3.19223);
-    gre->SetPoint(77,25749.2,6.38446);
-    gre->SetPointError(77,0,3.19223);
-    gre->SetPoint(78,33256.4,6.38446);
-    gre->SetPointError(78,0,3.19223);
-    gre->SetPoint(79,42952.2,6.38446);
-    gre->SetPointError(79,0,3.19223);
-    gre->SetPoint(80,55474.9,6.38446);
-    gre->SetPointError(80,0,3.19223);
-    gre->SetPoint(81,71648.6,6.38446);
-    gre->SetPointError(81,0,3.19223);
-    gre->SetPoint(82,92537.8,6.38446);
-    gre->SetPointError(82,0,3.19223);
-    gre->SetPoint(83,119517,6.38446);
-    gre->SetPointError(83,0,3.19223);
-    gre->SetPoint(84,154362,6.38446);
-    gre->SetPointError(84,0,3.19223);
-    gre->SetPoint(85,199367,6.38446);
-    gre->SetPointError(85,0,3.19223);
-    gre->SetPoint(86,257492,6.38446);
-    gre->SetPointError(86,0,3.19223);
-    gre->SetPoint(87,332563,6.38446);
-    gre->SetPointError(87,0,3.19223);
-    gre->SetPoint(88,429522,6.38446);
-    gre->SetPointError(88,0,3.19223);
-    gre->SetPoint(89,554749,6.38446);
-    gre->SetPointError(89,0,3.19223);
-    gre->SetPoint(90,716486,6.38446);
-    gre->SetPointError(90,0,3.19223);
-    gre->Draw("l same");
-    gre->DrawClone("l3 same");
-    gPad->Modified();
-    gPad->Update();
-    gPad->cd();
+    static TGraph* graph = 0;
+    if (!graph) { 
+      graph = new TGraph(20);
+      graph->SetName("stop_nodelta");
+      graph->SetTitle("Stopping w/o #delta's [RPP fig 27.1]");
+      graph->GetHistogram()->SetYTitle("(MeVcm^{2}/g)");
+      graph->GetHistogram()->SetXTitle("#beta#gamma");
+      graph->SetFillColor(0);
+      graph->SetLineColor(3);
+      graph->SetLineStyle(1);
+      graph->SetLineWidth(1);
+      graph->SetPoint(0,0.001461622,40.17542);
+      graph->SetPoint(1,0.003775053,91.28429);
+      graph->SetPoint(2,0.01178769,202.7359);
+      graph->SetPoint(3,0.01722915,212.1938);
+      graph->SetPoint(4,0.03162278,172.8318);
+      graph->SetPoint(5,0.06028646,91.28429);
+      graph->SetPoint(6,0.09506529,51.62633);
+      graph->SetPoint(7,0.433873,5.281682);
+      graph->SetPoint(8,1.255744,1.808947);
+      graph->SetPoint(9,2.304822,1.473387);
+      graph->SetPoint(10,3.921088,1.473387);
+      graph->SetPoint(11,8.064796,1.614064);
+      graph->SetPoint(12,26.15667,1.936996);
+      graph->SetPoint(13,264.8969,2.489084);
+      graph->SetPoint(14,544.8334,2.665278);
+      graph->SetPoint(15,1163.949,2.853945);
+      graph->SetPoint(16,5312.204,3.19853);
+      graph->SetPoint(17,15374.93,3.424944);
+      graph->SetPoint(18,49865.73,3.667384);
+      graph->SetPoint(19,634158.5,4.110185);
+      ScaleGraph(graph);
+    }
+    graph->Draw("C same");
+    return graph;
+  }
+
+  /** Draw pure Bethe-Bloc from Review of Particle Physics, fig. 27.1,
+      but without delta electrons  
+      @return TGraph object */ 
+  TGraph* FromRPPNoRad() 
+  {
+    static TGraph* graph = 0;
+    if (!graph) { 
+      graph = new TGraph(18);
+      graph->SetName("norad_stop");
+      graph->SetTitle("Stopping w/o radiative loss [RPP fig. 27.1]");
+      graph->GetHistogram()->SetYTitle("(MeVcm^{2}/g)");
+      graph->GetHistogram()->SetXTitle("#beta#gamma");
+      graph->SetFillColor(0);
+      graph->SetLineStyle(1);
+      graph->SetLineColor(4);
+      graph->SetLineWidth(1);
+      graph->SetPoint(0,0.001,24.3298);
+      graph->SetPoint(1,0.003117649,74.35105);
+      graph->SetPoint(2,0.008675042,172.8318);
+      graph->SetPoint(3,0.01782497,212.1938);
+      graph->SetPoint(4,0.02704573,189.3336);
+      graph->SetPoint(5,0.07481082,70.29816);
+      graph->SetPoint(6,0.3300035,8.524974);
+      graph->SetPoint(7,0.819559,2.489084);
+      graph->SetPoint(8,1.447084,1.651284);
+      graph->SetPoint(9,2.555097,1.440177);
+      graph->SetPoint(10,4.026598,1.407715);
+      graph->SetPoint(11,32.38084,1.728318);
+      graph->SetPoint(12,97.19733,1.893336);
+      graph->SetPoint(13,1732.539,2.170869);
+      graph->SetPoint(14,11098.58,2.324538);
+      graph->SetPoint(15,32075.46,2.378141);
+      graph->SetPoint(16,221655.8,2.546482);
+      graph->SetPoint(17,593830.6,2.605203);
+      ScaleGraph(graph);
+    }
+    graph->Draw("C same");
+    return graph;
+  }
+
+  /** Draw pure Bethe-Bloc from Review of Particle Physics, fig. 27.6 
+      @return TGraph object */ 
+  TGraph* FromRPPMean() 
+  {
+    static TGraph* graph = 0;
+    if (!graph) { 
+      graph = new TGraph(12);
+      graph->SetName("mean_eloss");
+      graph->SetTitle("Mean #Delta E/#Delta x - 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->SetMarkerStyle(21);
+      graph->SetMarkerSize(0.6);
+      graph->SetPoint(0,0.1,1.346561);
+      graph->SetPoint(1,0.1435819,1.230159);
+      graph->SetPoint(2,0.2061576,1.156085);
+      graph->SetPoint(3,0.3698076,1.124339);
+      graph->SetPoint(4,0.4620113,1.124339);
+      graph->SetPoint(5,0.8521452,1.145503);
+      graph->SetPoint(6,1.909707,1.177249);
+      graph->SetPoint(7,4.048096,1.198413);
+      graph->SetPoint(8,12.66832,1.219577);
+      graph->SetPoint(9,48.17031,1.230159);
+      graph->SetPoint(10,285.8863,1.230159);
+      graph->SetPoint(11,894.6674,1.230159);
+      const Double_t m   = 0.10566; // Muon 
+      ScaleGraph(graph, true, m);
+    }
+    graph->Draw("C same");
+    return graph;
+  }
+
+  /** Draw energy loss as obtained from GEANT 3.21 GFMATE. 
+      @return TGraph object */
+  TGraph* FromGFMATE() 
+  {
+    static TGraphErrors* gre = 0;
+    if (!gre) {
+      gre = new TGraphErrors(91);
+      gre->SetName("ELOSS");
+      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->SetLineStyle(1);
+      gre->SetLineWidth(1);
+      gre->SetPoint(0,7.16486e-05,1218.84);
+      gre->SetPoint(1,9.25378e-05,1221.38);
+      gre->SetPoint(2,0.000119517,1180.12);
+      gre->SetPoint(3,0.000154362,1100.31);
+      gre->SetPoint(4,0.000199367,996.621);
+      gre->SetPoint(5,0.000257492,886.005);
+      gre->SetPoint(6,0.000332563,780.483);
+      gre->SetPoint(7,0.000429522,684.927);
+      gre->SetPoint(8,0.000554749,599.407);
+      gre->SetPoint(9,0.000716486,522.375);
+      gre->SetPoint(10,0.000925378,452.497);
+      gre->SetPoint(11,0.00119517,389.101);
+      gre->SetPoint(12,0.00154362,331.974);
+      gre->SetPoint(13,0.00199367,280.969);
+      gre->SetPoint(14,0.00257492,235.689);
+      gre->SetPoint(15,0.00332564,196.156);
+      gre->SetPoint(16,0.00429522,162.402);
+      gre->SetPoint(17,0.00554749,133.87);
+      gre->SetPoint(18,0.00716486,109.959);
+      gre->SetPoint(19,0.00925378,90.2035);
+      gre->SetPoint(20,0.0119517,74.1317);
+      gre->SetPoint(21,0.0154362,60.8988);
+      gre->SetPoint(22,0.0199367,49.9915);
+      gre->SetPoint(23,0.0257492,40.9812);
+      gre->SetPoint(24,0.0332564,33.5739);
+      gre->SetPoint(25,0.0429522,27.5127);
+      gre->SetPoint(26,0.0554749,22.5744);
+      gre->SetPoint(27,0.0716486,18.5674);
+      gre->SetPoint(28,0.0925378,15.3292);
+      gre->SetPoint(29,0.119517,12.7231);
+      gre->SetPoint(30,0.154362,10.6352);
+      gre->SetPoint(31,0.199367,8.97115);
+      gre->SetPoint(32,0.257492,7.65358);
+      gre->SetPoint(33,0.332564,6.61909);
+      gre->SetPoint(34,0.429522,5.79837);
+      gre->SetPoint(35,0.554749,5.148);
+      gre->SetPoint(36,0.716486,4.65024);
+      gre->SetPoint(37,0.925378,4.27671);
+      gre->SetPoint(38,1.19517,3.99831);
+      gre->SetPoint(39,1.54362,3.79877);
+      gre->SetPoint(40,1.99367,3.6629);
+      gre->SetPoint(41,2.57492,3.57594);
+      gre->SetPoint(42,3.32564,3.52565);
+      gre->SetPoint(43,4.29522,3.50206);
+      gre->SetPoint(44,5.54749,3.49715);
+      gre->SetPoint(45,7.16486,3.50467);
+      gre->SetPoint(46,9.25378,3.51988);
+      gre->SetPoint(47,11.9517,3.53932);
+      gre->SetPoint(48,15.4362,3.56054);
+      gre->SetPoint(49,19.9367,3.58189);
+      gre->SetPoint(50,25.7492,3.60231);
+      gre->SetPoint(51,33.2564,3.62113);
+      gre->SetPoint(52,42.9522,3.638);
+      gre->SetPoint(53,55.4749,3.65275);
+      gre->SetPoint(54,71.6486,3.66537);
+      gre->SetPoint(55,92.5378,3.67586);
+      gre->SetPoint(56,119.517,3.68433);
+      gre->SetPoint(57,154.362,3.69105);
+      gre->SetPoint(58,199.367,3.6962);
+      gre->SetPoint(59,257.492,3.69997);
+      gre->SetPoint(60,332.564,3.70257);
+      gre->SetPoint(61,429.522,3.70421);
+      gre->SetPoint(62,554.749,3.70511);
+      gre->SetPoint(63,716.486,3.7055);
+      gre->SetPoint(64,925.378,3.70559);
+      gre->SetPoint(65,1195.17,3.70558);
+      gre->SetPoint(66,1543.62,3.70557);
+      gre->SetPoint(67,1993.67,3.70555);
+      gre->SetPoint(68,2574.92,3.70553);
+      gre->SetPoint(69,3325.64,3.70552);
+      gre->SetPoint(70,4295.22,3.7055);
+      gre->SetPoint(71,5547.49,3.70548);
+      gre->SetPoint(72,7164.86,3.70547);
+      gre->SetPoint(73,9253.78,3.70545);
+      gre->SetPoint(74,11951.7,3.70544);
+      gre->SetPoint(75,15436.2,3.70544);
+      gre->SetPoint(76,19936.7,3.70544);
+      gre->SetPoint(77,25749.2,3.70544);
+      gre->SetPoint(78,33256.4,3.70544);
+      gre->SetPoint(79,42952.2,3.70544);
+      gre->SetPoint(80,55474.9,3.70544);
+      gre->SetPoint(81,71648.6,3.70544);
+      gre->SetPoint(82,92537.8,3.70544);
+      gre->SetPoint(83,119517,3.70544);
+      gre->SetPoint(84,154362,3.70544);
+      gre->SetPoint(85,199367,3.70544);
+      gre->SetPoint(86,257492,3.70544);
+      gre->SetPoint(87,332563,3.70544);
+      gre->SetPoint(88,429522,3.70544);
+      gre->SetPoint(89,554749,3.70544);
+      gre->SetPoint(90,716486,3.70544);
+      // Double_t* x = gre->GetX();
+      Double_t* y = gre->GetY();
+      for (Int_t i = 0; i < gre->GetN(); i++) 
+       gre->SetPointError(i, 0, 2 * 0.1 * y[i]); // ! 1 sigma
+    }
+    gre->DrawClone("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
+      MPV. */ 
+  TGraph* GetResp()
+  {
+    static TGraph*  graph = 0;
+    if (!graph) {
+      graph = new TGraph(16);
+      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->SetMarkerStyle(21);
+      graph->SetMarkerSize(0.6);
+      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);
+    }
+    return graph;
+  }
+
+  /** Get the correction to Bethe-Bloc from Review of Particle Physics
+      (fig 27.8). 
+  */
+  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->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;
   }
   
   ClassDef(DrawHits,0);
index 46b3bc1..5cbd4e1 100644 (file)
@@ -27,6 +27,7 @@
 #include <TTree.h>
 #include <AliStack.h>
 #include <AliLog.h>
+#include <TF1.h>
 
 //____________________________________________________________________
 /** @class DrawHitsRecs
@@ -48,6 +49,8 @@ private:
   TH1D* fDiffE;       // Histogram 
   TH2D* fHitsVsRecM;  // Histogram 
   TH2D* fDiffM;       // Histogram 
+  TH1*  fHitEloss;
+  TH1*  fRecEloss;
   AliFMDEdepMap  fMap;
   AliFMDFloatMap fEta;
   AliFMDFloatMap fPhi;
@@ -117,7 +120,16 @@ public:
     fDiffM->SetXTitle("M_{sim} - M_{rec}");
     fDiffM->SetYTitle("|#eta|");
     // fDiffM->SetYTitle("Detector");
+
+    fHitEloss = new TH1D("hitEloss", "#frac{#Delta E_{sim}}{#Delta x} (MeV/cm)", 
+                        100, 0, 10);
+    fHitEloss->SetFillColor(2);
+    fHitEloss->SetFillStyle(3001);
     
+    fRecEloss = new TH1D("recEloss", "#frac{#Delta E_{rec}}{#Delta x} (MeV/cm)", 
+                        100, 0, 10);
+    fRecEloss->SetFillColor(4);
+    fRecEloss->SetFillStyle(3001);
   }
   //__________________________________________________________________
   /** Begining of event
@@ -147,6 +159,7 @@ public:
       if (!kine->IsPrimary()) return kTRUE;
     }
     
+    if (hit->Edep()/hit->Length() > 0.1) fHitEloss->Fill(hit->Edep() / hit->Length());
     fMap(det, rng, sec, str).fEdep += hit->Edep();
     fMap(det, rng, sec, str).fN++;
     return kTRUE;
@@ -192,6 +205,7 @@ public:
     }
     if (nhit > 0) fHitsVsRecM->Fill(nhit, single->Particles());
     fDiffM->Fill(nhit - single->Particles(), TMath::Abs(single->Eta()));
+    if (single->Edep()/.03 > 0.1) fRecEloss->Fill(single->Edep() / 0.0300);
     return kTRUE;
   }
   //__________________________________________________________________
@@ -232,6 +246,20 @@ public:
     c->SetLogz();
     fDiffM->Draw("colz");
 
+    c = new TCanvas("c6", "Hit Eloss, Reco Eloss");
+    fRecEloss->Scale(1./fRecEloss->GetMaximum());
+    fRecEloss->Draw();
+    fRecEloss->Fit("landau", "", "SAME", 2, 4);
+    TF1* recResp = new TF1(*fRecEloss->GetFunction("landau"));
+    fHitEloss->Scale(1./fHitEloss->GetMaximum());
+    fHitEloss->Draw("same");
+    fHitEloss->Fit("landau", "", "SAME", 2, 10);
+    TF1* hitResp = new TF1(*fHitEloss->GetFunction("landau"));
+    std::cout << "Hit MPV,width: " << hitResp->GetParameter(1) << ","
+             << hitResp->GetParameter(2) << "\n"
+             << "Rec MPV,width: " << recResp->GetParameter(1) << ","
+             << recResp->GetParameter(2) << std::endl;
+    c->SetLogy();
 
     return kTRUE;
   }
index 3298339..9b5a635 100644 (file)
 //
 // Note, that VMC _must_ be the TGeant3TGeo VMC. 
 //
+#ifdef __CINT__
+XSection()
+{
+  gROOT->ProcessLine(".x Compile.C(\"$(ALICE_ROOT)/FMD/scripts/XSection.C\"");
+  gAlice->InitMC("$(ALICE_ROOT)/FMD/Config.C");
+  TFile* file = TFile::Open("xsec.root", "RECREATE");
+  GetXsection("FMD_Si$", "pi+");
+  file->Close();
+}
+#else  
 #include <TArrayF.h>
 #include <TTree.h>
 #include <TMath.h>
 #include <TGeoManager.h>
 #include <TGeoMedium.h>
 #include <TGeant3.h>
+#include <TGraph.h>
+#include <TAxis.h>
+#include <cmath>
+#include <string>
+#include <vector>
+#include <sstream>
+#include <iomanip>
 /** @defgroup xsec_script X-section script 
     @ingroup FMD_script 
 */
  */
 struct Mech 
 {
-  char*    name;
-  char*    title;
-  char*    unit;
-  TArrayF  values;
-  int      status;
+  char* fName;
+  char* fTitle;
+  char* fUnit;
 };
 
+Mech fgMechs[] = 
+  {{ "HADF","total hadronic x-section according to FLUKA",      "cm^{1}"},
+   { "INEF","hadronic inelastic x-section according to FLUKA",  "cm^{1}"},
+   { "ELAF","hadronic elastic x-section according to FLUKA",    "cm^{1}"},
+   { "HADG","total hadronic x-section according to GHEISHA",    "cm^{1}"},
+   { "INEG","hadronic inelastic x-section according to GHEISHA","cm^{1}"},
+   { "ELAG","hadronic elastic x-section according to GHEISHA",  "cm^{1}"},
+   { "FISG","nuclear fission x-section according to GHEISHA",   "cm^{1}"},
+   { "CAPG","neutron capture x-section according to GHEISHA",   "cm^{1}"},
+   { "LOSS","stopping power",                                   "cm^{1}"},
+   { "PHOT","photoelectric x-section",                          "MeV/cm"},
+   { "ANNI","positron annihilation x-section",                  "cm^{1}"},
+   { "COMP","Compton effect x-section",                         "cm^{1}"},
+   { "BREM","bremsstrahlung x-section",                         "cm^{1}"},
+   { "PAIR","photon and muon direct- pair x-section",           "cm^{1}"},
+   { "DRAY","delta-rays x-section",                             "cm^{1}"},
+   { "PFIS","photo-fission x-section",                          "cm^{1}"},
+   { "RAYL","Rayleigh scattering x-section",                    "cm^{1}"},
+   { "MUNU","muon-nuclear interaction x-section",               "cm^{1}"},
+   { "RANG","range",                                            "cm"},
+   { "STEP","maximum step",                                     "cm"},
+   { 0, 0, 0,}};
+
+//____________________________________________________________________
+/** @ingroup xsec_script
+ */
+struct MechValue 
+{
+  MechValue(const Mech& m, size_t n) 
+    : fMech(m), fValues(n), fStatus(0)
+  {}
+  bool Get(TGeant3* mc, std::vector<float>& t, 
+          std::vector<float>& c, int medNo, int pidNo) 
+  {
+    TGraph g(t.size());
+    g.SetName(fMech.fName);
+    g.SetTitle(fMech.fTitle);
+    g.GetXaxis()->SetTitle("p [GeV]");
+    g.GetYaxis()->SetTitle(fMech.fUnit);
+    int ixst;
+    mc->Gftmat(medNo, pidNo, fMech.fName, t.size(), 
+              &(t[0]), &(fValues[0]), &(c[0]), ixst);
+    fStatus = ixst;
+    if (!fStatus) return false;
+    for (size_t i = 0; i < t.size(); i++) g.SetPoint(i, t[i], fValues[i]);
+    g.Write();
+    return true;
+  }
+  const Mech&        fMech;
+  std::vector<float> fValues;
+  bool               fStatus;
+};
+
+//____________________________________________________________________
+/** @ingroup xsec_script
+ */
+struct XSection 
+{
+  XSection(size_t n=91, float emin=1e-5, float emax=1e4) 
+    : fTKine(n), 
+      fCuts(5)
+  {
+    float dp   = 1. / log10(emax/emin);
+    float pmin = log10(emin);
+    fTKine[0]  = emin;
+    for (size_t i = 1; i < fTKine.size(); i++) {
+      float el  = pmin + i * dp;
+      fTKine[i] = pow(10, el);
+    }
+    for (float_array::iterator i = fCuts.begin(); i != fCuts.end(); ++i) 
+      *i = 1e-4;
+    Mech* mech = &(fgMechs[0]);
+    size_t i = 0;
+    while (mech->fName) {
+      fMechs.push_back(new MechValue(*mech, n));
+      mech++;
+    }
+  }
+  void Run(const std::string& medName, const std::string& pdgName) 
+  {
+    TGeant3* mc = static_cast<TGeant3*>(gMC);
+    if (!mc) {
+      std::cerr << "Couldn't get VMC" << std::endl;
+      return;
+    }
+    TGeoMedium* medium = gGeoManager->GetMedium(medName.c_str());
+    if (!medium) {
+      std::cerr << "Couldn't find medium " << medName << std::endl;
+      return;
+    }
+    int medNo = medium->GetMaterial()->GetUniqueID();
+    TDatabasePDG* pdgDb = TDatabasePDG::Instance();
+    TParticlePDG* pdgP  = pdgDb->GetParticle(pdgName.c_str());
+    if (!pdgP) {
+      std::cerr << "Couldn't find particle " << pdgName << std::endl;
+      return;
+    }
+    int pdgNo = pdgP->PdgCode();
+    int pidNo = mc->IdFromPDG(pdgNo);
+
+    std::stringstream vars;
+    vars << "T/F";
+
+    size_t nOk   = 0;
+    // Loop over defined mechanisms 
+    for (mech_array::iterator i = fMechs.begin(); i != fMechs.end(); ++i) {
+      if (!(*i)->Get(mc, fTKine, fCuts, medNo, pidNo))continue;
+      vars << ":" << (*i)->fMech.fName;
+      nOk ++;
+    }
+    
+    std::stringstream tName;
+    tName << medName << "_" << pdgName;
+    TTree* tree = new TTree(tName.str().c_str(), tName.str().c_str());
+
+    float_array cache(nOk+1);
+    tree->Branch("xsec", &(cache[0]), vars.str().c_str());
+    for (size_t i = 0; i < fTKine.size(); i++) {
+      cache[0] = fTKine[i];
+      int k = 0;
+      for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
+       if (!(*j)->fStatus) continue;
+       cache[++k] = (*j)->fValues[i];
+      }
+      tree->Fill();
+    }
+    tree->Write();
+  }
+protected: 
+  typedef std::vector<float> float_array;
+  typedef std::vector<MechValue*> mech_array;
+  float_array  fTKine;
+  float_array  fCuts;
+  mech_array   fMechs;
+};
 
+#if 0
 //____________________________________________________________________
 /** @ingroup xsec_script
     @param medName 
@@ -148,6 +299,8 @@ GetXsection(const char* medName, const char* pdgName,
   }
   tree->Write();
 }
+#endif
+#endif
 //____________________________________________________________________
 //
 // EOF
diff --git a/FMD/scripts/XSection.C b/FMD/scripts/XSection.C
new file mode 100644 (file)
index 0000000..90e5095
--- /dev/null
@@ -0,0 +1,280 @@
+//____________________________________________________________________
+//
+// $Id$
+//
+// Script to get the various cross sections, energy loss, and ranges
+// of a particular particle type in a particular medium. 
+//
+// This script should be compiled to speed it up.  
+// 
+// It creates a tree on the current output file, with the relevant
+// information. 
+//
+// Note, that VMC _must_ be the TGeant3TGeo VMC. 
+//
+#include <TArrayF.h>
+#include <TTree.h>
+#include <TMath.h>
+#include <iostream>
+#include <TVirtualMC.h>
+#include <TDatabasePDG.h>
+#include <TString.h>
+#include <TGeoManager.h>
+#include <TGeoMedium.h>
+#include <TGeant3.h>
+#include <TGraph.h>
+#include <TGraphErrors.h>
+#include <TAxis.h>
+#include <TFile.h>
+#include <TH2.h>
+#include <TLegend.h>
+#include <AliRun.h>
+#include <TCanvas.h>
+#include <cmath>
+#include <string>
+#include <vector>
+#include <sstream>
+#include <iomanip>
+/** @defgroup xsec_script X-section script 
+    @ingroup FMD_script 
+*/
+//____________________________________________________________________
+/** @ingroup xsec_script
+ */
+struct Mech 
+{
+  char* fName;
+  char* fTitle;
+  char* fUnit;
+  int   fColor;
+  int   fStyle;
+};
+
+Mech fgMechs[] = 
+  {{ "HADF","FLUKA total hadronic x-section",        "cm^{-1}",1,1},
+   { "INEF","FLUKA hadronic inelastic x-section",    "cm^{-1}",2,1},
+   { "ELAF","FLUKA hadronic elastic x-section",      "cm^{-1}",3,1},
+   { "HADG","GHEISHA total hadronic x-section",      "cm^{-1}",4,1},
+   { "INEG","GHEISHA hadronic inelastic x-section",  "cm^{-1}",6,1},
+   { "ELAG","GHEISHA hadronic elastic x-section",    "cm^{-1}",7,1},
+   { "FISG","GHEISHA nuclear fission x-section",     "cm^{-1}",1,2},
+   { "CAPG","GHEISHA neutron capture x-section",     "cm^{-1}",2,2},
+   { "LOSS","stopping power",                        "MeV/cm", 3,2},
+   { "PHOT","photoelectric x-section",               "cm^{-1}",4,2},
+   { "ANNI","positron annihilation x-section",       "cm^{-1}",6,2},
+   { "COMP","Compton effect x-section",              "cm^{-1}",7,2},
+   { "BREM","bremsstrahlung x-section",              "cm^{-1}",1,3},
+   { "PAIR","photon and muon direct- pair x-section","cm^{-1}",2,3},
+   { "DRAY","delta-rays x-section",                  "cm^{-1}",3,3},
+   { "PFIS","photo-fission x-section",               "cm^{-1}",4,3},
+   { "RAYL","Rayleigh scattering x-section",         "cm^{-1}",6,3},
+   { "MUNU","muon-nuclear interaction x-section",    "cm^{-1}",7,3},
+   { "RANG","range",                                 "cm",     1,4},
+   { "STEP","maximum step",                          "cm",     2,4},
+   { 0, 0, 0, 0, 0}};
+
+//____________________________________________________________________
+/** @ingroup xsec_script
+ */
+struct MechValue 
+{
+  MechValue(const Mech& m, size_t n) 
+    : fMech(m), fValues(n), fStatus(0), fELoss(0)
+  {
+    fGraph.SetName(fMech.fName);
+    fGraph.SetTitle(fMech.fTitle);
+    fGraph.GetXaxis()->SetTitle("#beta#gamma");
+    fGraph.GetYaxis()->SetTitle(fMech.fUnit);
+    fGraph.SetLineColor(fMech.fColor);
+    fGraph.SetLineStyle(fMech.fStyle);
+    fGraph.SetLineWidth(2);
+    std::string name(fMech.fName);
+    if (name != "LOSS") return;
+    
+    fELoss = new TGraphErrors(n);
+    fELoss->SetName("ELOSS");
+    fELoss->SetTitle(fMech.fTitle);
+    fELoss->GetXaxis()->SetTitle("#beta#gamma");
+    fELoss->GetYaxis()->SetTitle(fMech.fUnit);
+    fELoss->SetLineColor(fMech.fColor);
+    fELoss->SetLineStyle(fMech.fStyle);
+    fELoss->SetLineWidth(2);
+    fELoss->SetFillColor(fMech.fColor+1);
+    fELoss->SetFillStyle(3001);
+  }
+  bool Get(TGeant3* mc, std::vector<float>& t, 
+          std::vector<float>& c, int medNo, int pidNo, 
+          float mass) 
+  {
+    int ixst;
+    mc->Gftmat(medNo, pidNo, fMech.fName, t.size(), 
+              &(t[0]), &(fValues[0]), &(c[0]), ixst);
+    fStatus = ixst;
+    if (!fStatus) return false;
+    fGraph.Set(t.size());
+    for (size_t i = 0; i < t.size(); i++) {
+      fGraph.SetPoint(i, t[i]/mass, fValues[i]);
+      if (!fELoss) continue;
+      fELoss->SetPoint(i, t[i]/mass, fValues[i]);
+      // ~ 5 sigma
+      fELoss->SetPointError(i, 0, .5 * fValues[i]);
+    }
+    fGraph.Write();
+    if (fELoss) fELoss->Write();
+    return true;
+  }
+  TGraph& Draw()
+  {
+    if (fELoss) fELoss->DrawClone("4 same");
+    fGraph.DrawClone("l same");
+    return fGraph;
+  }
+  const Mech&        fMech;
+  std::vector<float> fValues;
+  bool               fStatus;
+  TGraph             fGraph;
+  TGraphErrors*      fELoss;
+};
+
+//____________________________________________________________________
+/** @ingroup xsec_script
+ */
+struct XSections 
+{
+  XSections(size_t n=91, float emin=1e-5, float emax=1e4) 
+    : fTKine(n), 
+      fCuts(5)
+  {
+    float dp   = 1. / log10(emax/emin);
+    float pmin = log10(emin);
+    fTKine[0]  = emin;
+    for (size_t i = 1; i < fTKine.size(); i++) {
+      float el  = pmin + i * dp;
+      fTKine[i] = pow(10, el);
+    }
+    for (float_array::iterator i = fCuts.begin(); i != fCuts.end(); ++i) 
+      *i = 1e-4;
+    Mech* mech = &(fgMechs[0]);
+    while (mech->fName) {
+      fMechs.push_back(new MechValue(*mech, n));
+      mech++;
+    }
+  }
+  void Run(const std::string& medName, const std::string& pdgName) 
+  {
+    TGeant3* mc = static_cast<TGeant3*>(gMC);
+    if (!mc) {
+      std::cerr << "Couldn't get VMC" << std::endl;
+      return;
+    }
+    TGeoMedium* medium = gGeoManager->GetMedium(medName.c_str());
+    if (!medium) {
+      std::cerr << "Couldn't find medium " << medName << std::endl;
+      return;
+    }
+    int medNo = medium->GetMaterial()->GetUniqueID();
+    TDatabasePDG* pdgDb = TDatabasePDG::Instance();
+    fPDG                = pdgDb->GetParticle(pdgName.c_str());
+    if (!fPDG) {
+      std::cerr << "Couldn't find particle " << pdgName << std::endl;
+      return;
+    }
+    int pdgNo = fPDG->PdgCode();
+    int pidNo = mc->IdFromPDG(pdgNo);
+
+    std::stringstream vars;
+    vars << "betagamma/F";
+
+    size_t nOk   = 0;
+    // Loop over defined mechanisms 
+    for (mech_array::iterator i = fMechs.begin(); i != fMechs.end(); ++i) {
+      if (!(*i)->Get(mc, fTKine, fCuts, medNo, pidNo, fPDG->Mass()))continue;
+      vars << ":" << (*i)->fMech.fName;
+      nOk ++;
+    }
+    
+    std::stringstream tName;
+    tName << medName << "_" << pdgName;
+    TTree* tree = new TTree(tName.str().c_str(), tName.str().c_str());
+
+    float_array cache(nOk+1);
+    tree->Branch("xsec", &(cache[0]), vars.str().c_str());
+    for (size_t i = 0; i < fTKine.size(); i++) {
+      cache[0] = fTKine[i] / fPDG->Mass();
+      int k = 0;
+      for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
+       if (!(*j)->fStatus) continue;
+       cache[++k] = (*j)->fValues[i];
+      }
+      tree->Fill();
+    }
+    tree->Write();
+  }
+  void Draw() 
+  {
+    float min = 100000;
+    float max = 0;
+    std::vector<TGraph*> gs;
+    float_array bg(fTKine.size());
+    for (size_t i = 0; i < fTKine.size(); i++) bg[i] = fTKine[i]/fPDG->Mass();
+    for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
+      if (!(*j)->fStatus) continue;
+      for (size_t i = 0; i < fTKine.size(); i++) {
+       if ((*j)->fValues[i] == 0) continue;
+       min = std::min(min, (*j)->fValues[i]);
+       max = std::max(max, (*j)->fValues[i]);
+      }
+    }
+
+    TCanvas* c = new TCanvas("c", "C", 700, 700);
+    c->SetFillColor(0);
+    c->SetLogy();
+    c->SetLogx();
+    c->SetGridy();
+    c->SetGridx();
+    float_array y(101);
+    float ymin = log10(min);
+    float dy   = (log10(max)+.5 - log10(min)) / y.size();
+    for (size_t i = 1; i < y.size(); i++)  y[i] = pow(10, ymin + i * dy);
+    TH2* f = new TH2F("x", "X-sec",bg.size()-1,&(bg[0]),y.size()-1,&(y[0]));
+    f->SetXTitle("#beta#gamma");
+    f->SetDirectory(0);
+    f->SetStats(kFALSE);
+    f->Draw();
+    TLegend* l = new TLegend(0.45, 0.125, 0.90, 0.45);
+    l->SetFillColor(0);
+    // l->SetFillStyle(0);
+    for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
+      if (!(*j)->fStatus) continue;
+      TGraph& g = (*j)->Draw();
+      l->AddEntry(&g, g.GetTitle(), "l");
+    }
+    l->Draw("same");
+  }
+protected: 
+  typedef std::vector<float> float_array;
+  typedef std::vector<MechValue*> mech_array;
+  float_array   fTKine;
+  float_array   fCuts;
+  mech_array    fMechs;
+  TParticlePDG* fPDG;
+};
+
+bool init = false;
+
+void XSection(const char* pdgName="pi-")
+{
+  if (!init) {
+    gAlice->InitMC("$(ALICE_ROOT)/FMD/Config.C");
+    init = true;
+  }
+  TFile* file = TFile::Open("xsec.root", "RECREATE");
+  XSections xs;
+  xs.Run("FMD_Si$", pdgName);
+  xs.Draw();
+  file->Close();
+}
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/FMD/scripts/checkSizes.sh b/FMD/scripts/checkSizes.sh
new file mode 100755 (executable)
index 0000000..a04a16f
--- /dev/null
@@ -0,0 +1,96 @@
+#!/bin/bash
+
+extra=" DPMJET         \
+        TDPMjet                \
+        EPEMGEN                \
+        TEPEMGEN       \
+        HBTP           \
+        THbtp          \
+        HERWIG         \
+        THerwig                \
+        HIJING         \
+        THijing                \
+        ISAJET         \
+        TIsajet                \
+        LHAPDF         \
+        MEVSIM         \
+        TMEVSIM                \
+        MICROCERN      \
+        PDF            \
+        PYTHIA6                \
+        TPHIC"
+base="  ALIFAST                \
+       ALIROOT         \
+       ANALYSIS        \
+       CONTAINERS      \
+        CRT            \
+       DISPLAY         \
+       EMCAL           \
+       EVE             \
+        EVGEN          \
+       FASTSIM         \
+       FLOW            \
+       FMD             \
+        HBTAN          \
+       HLT             \
+       ITS             \
+       JETAN           \
+        LHC            \
+       MONITOR         \
+       MUON            \
+       PHOS            \
+        PMD            \
+       PWG0            \
+       PWG2            \
+       PWG3            \
+        RALICE         \
+       RAW             \
+       RICH            \
+       SHUTTLE         \
+        START          \
+       STEER           \
+       STRUCT          \
+       TOF             \
+        TPC            \
+       TRD             \
+       VZERO           \
+       ZDC"
+
+cat <<EOF > exclude
+*/tgt_*/*
+*/html/*
+.#*
+*/CVS*
+*~ 
+*.root
+*.so
+*.o
+EOF
+
+get_size()
+{
+    s=`du -X exclude -kc $1 | tail -n 1 | awk 'BEGIN {FS=" "}{print $1}'`
+    printf "\t%-30s\t%10d kB\n" $1 $s 
+    total=`echo ${total} + ${s} | bc`
+}
+    
+echo "Extras:"
+total=0
+for e in $extra ; do 
+    get_size $e
+done 
+for i in `seq 1 56` ; do echo -n "-" ; done
+mb=`echo $total / 1024 | bc` 
+printf "\n\t%-30s\t%10d kB = %10d MB\n" "Total" $total $mb
+
+echo "Base:"
+total=0
+for b in $base ; do 
+    get_size $b 
+done 
+for i in `seq 1 56` ; do echo -n "-" ; done
+mb=`echo $total / 1024 | bc` 
+printf "\n\t%-30s\t%10d kB = %10d MB\n" "Total" $total $mb
+
+        
+rm -f exclude 
\ No newline at end of file