FMD/scripts/DrawXsection.C

   1 //____________________________________________________________________
   2 //
   3 // Script to draw a X-section, LOSS, or range made with MakeXsection
   4 //
   5 /** @ingroup xsec_script
   6     @param scale
   7     @param filename
   8     @param var
   9     @param medName
  10     @param thick
  11     @param pdgName
  12 */
  13 void
  14 DrawXsection(Bool_t scale=kFALSE,
  15              const char* filename="xsec.root",
  16              const char* var="LOSS",
  17              const char* medName="FMD_Si$",
  18              Double_t thick=.03,
  19              const char* pdgName="pi+")
  20 {
  21   TFile*   file = TFile::Open(filename, "READ");
  22   TTree*   tree = static_cast<TTree*>(file->Get(Form("%s_%s",medName,
  23                                                      pdgName)));
  24   TLeaf* tb   = tree->GetLeaf("T");
  25   TLeaf* vb   = tree->GetLeaf(var);
  26   if (!vb) {
  27     std::cerr << "Leaf " << var << " not found" << std::endl;
  28     return;
  29   }
  30   Float_t tkine, value;
  31   tb->SetAddress(&tkine);
  32   vb->SetAddress(&value);
  33   Int_t n = tree->GetEntries();
  34
  35   Float_t xscale = 1;
  36   Float_t yscale = 1;
  37   if (scale) {
  38     TDatabasePDG* pdgDb = TDatabasePDG::Instance();
  39     TParticlePDG* pdgP  = pdgDb->GetParticle(pdgName);
  40     if (!pdgP) {
  41       std::cerr << "Couldn't find particle " << pdgName << std::endl;
  42       return;
  43     }
  44     Double_t m = pdgP->Mass();
  45     Double_t q = pdgP->Charge() / 3;
  46     if (m == 0 || q == 0) {
  47       std::cerr  << "Mass is 0" << std::endl;
  48       return;
  49     }
  50     xscale = 1 / m;
  51     yscale = 1 / (q * q);
  52   }
  53
  54   TGraphErrors* graph = new TGraphErrors(n);
  55   for (Int_t i = 0; i < n; i++) {
  56     tree->GetEntry(i);
  57     Double_t x = tkine*xscale;
  58     Double_t y = value*yscale;
  59     graph->SetPoint(i, x, y);
  60     // 5 sigma
  61     graph->SetPointError(i, 0, 5 * .1 * y);
  62   }
  63   TCanvas* c = new TCanvas("c","c");
  64   c->SetLogx();
  65   c->SetLogy();
  66   graph->SetLineWidth(2);
  67   graph->SetFillStyle(3001);
  68   graph->SetFillColor(6);
  69   graph->Draw("L");
  70   graph->DrawClone("AL3");
  71   c->Modified();
  72   c->Update();
  73   c->cd();
  74   c->SaveAs("xsec.C");
  75
  76 }
  77
  78 //____________________________________________________________________
  79 //
  80 // EOF
  81 //