Use VMC id's rather than TGeo id's

[u/mrichter/AliRoot.git] / FMD / scripts / DrawXsection.C

diff --git a/FMD/scripts/DrawXsection.C b/FMD/scripts/DrawXsection.C
index fb53a3aac86a16d4d7f2b0dc10014b68a8535305..ca45398f97484195356febe06d3d03095e37533c 100644 (file)
--- a/FMD/scripts/DrawXsection.C
+++ b/FMD/scripts/DrawXsection.C
@@ -1,11 +1,24 @@
+//____________________________________________________________________
+//
+// Script to draw a X-section, LOSS, or range made with MakeXsection
+//
+/** @ingroup xsec_script
+    @param scale 
+    @param filename 
+    @param var 
+    @param medName 
+    @param thick 
+    @param pdgName 
+*/
 void
-DrawXsection(const char* file="xsec.root", 
+DrawXsection(Bool_t scale=kFALSE, 
+            const char* filename="xsec.root", 
             const char* var="LOSS", 
             const char* medName="FMD_Si$", 
             Double_t thick=.03,
             const char* pdgName="pi+")
 {
-  TFile*   file = TFile::Open("xsec.root", "READ");
+  TFile*   file = TFile::Open(filename, "READ");
   TTree*   tree = static_cast<TTree*>(file->Get(Form("%s_%s",medName,
                                                     pdgName)));
   TLeaf* tb   = tree->GetLeaf("T");
@@ -19,25 +32,47 @@ DrawXsection(const char* file="xsec.root",
   vb->SetAddress(&value);
   Int_t n = tree->GetEntries();
 
-  TDatabasePDG* pdgDb = TDatabasePDG::Instance();
-  TParticlePDG* pdgP  = pdgDb->GetParticle(pdgName);
-  if (!pdgP) {
-    std::cerr << "Couldn't find particle " << pdgName << std::endl;
-    return;
-  }
-  Double_t m = pdgP->Mass();
-  Double_t q = pdgP->Charge() / 3;
-  if (m == 0) {
-    std::cerr  << "Mass is 0" << std::endl;
-    return;
+  Float_t xscale = 1;
+  Float_t yscale = 1;
+  if (scale) {
+    TDatabasePDG* pdgDb = TDatabasePDG::Instance();
+    TParticlePDG* pdgP  = pdgDb->GetParticle(pdgName);
+    if (!pdgP) {
+      std::cerr << "Couldn't find particle " << pdgName << std::endl;
+      return;
+    }
+    Double_t m = pdgP->Mass();
+    Double_t q = pdgP->Charge() / 3;
+    if (m == 0 || q == 0) {
+      std::cerr  << "Mass is 0" << std::endl;
+      return;
+    }
+    xscale = 1 / m;
+    yscale = 1 / (q * q);
   }
   
-  TGraph* graph = new TGraph(n);
+  TGraphErrors* graph = new TGraphErrors(n);
   for (Int_t i = 0; i < n; i++) {
     tree->GetEntry(i);
-    graph->SetPoint(i, tkine/m/q/q, value*thick);
+    Double_t x = tkine*xscale;
+    Double_t y = value*yscale;
+    graph->SetPoint(i, x, y); 
+    // 5 sigma
+    graph->SetPointError(i, 0, 5 * .1 * y);
   }
-  graph->Draw("ALP");
+  TCanvas* c = new TCanvas("c","c");
+  c->SetLogx();
+  c->SetLogy();
+  graph->SetLineWidth(2);
+  graph->SetFillStyle(3001);
+  graph->SetFillColor(6);
+  graph->Draw("L");
+  graph->DrawClone("AL3");
+  c->Modified();
+  c->Update();
+  c->cd();
+  c->SaveAs("xsec.C");
+  
 }
 
 //____________________________________________________________________