From 2b893216a60932c634827f02a8ba890f816fdc38 Mon Sep 17 00:00:00 2001 From: cholm Date: Fri, 14 Sep 2007 13:06:51 +0000 Subject: [PATCH] Bug fix for reading raw data - thanks Frederic YERMIA, Diego, and Gines MARTINEZ --- FMD/AliFMDDisplay.h | 7 +- FMD/AliFMDPattern.h | 2 +- FMD/AliFMDRawReader.cxx | 6 +- FMD/AliFMDRawStream.cxx | 8 +- FMD/Config.C | 2 +- FMD/Simulate.C | 4 +- FMD/scripts/DisplayDigits.C | 1 + FMD/scripts/DrawESD.C | 223 ++++++++++++ FMD/scripts/DrawHits.C | 675 ++++++++++++++++++++++++------------ FMD/scripts/DrawHitsRecs.C | 28 ++ FMD/scripts/GetXsection.C | 163 ++++++++- FMD/scripts/XSection.C | 280 +++++++++++++++ FMD/scripts/checkSizes.sh | 96 +++++ 13 files changed, 1264 insertions(+), 231 deletions(-) create mode 100644 FMD/scripts/DrawESD.C create mode 100644 FMD/scripts/XSection.C create mode 100755 FMD/scripts/checkSizes.sh diff --git a/FMD/AliFMDDisplay.h b/FMD/AliFMDDisplay.h index 159a04e3a0a..3c1f1eb6ad3 100644 --- a/FMD/AliFMDDisplay.h +++ b/FMD/AliFMDDisplay.h @@ -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 */ diff --git a/FMD/AliFMDPattern.h b/FMD/AliFMDPattern.h index b29c5e037aa..bca5489b5cf 100644 --- a/FMD/AliFMDPattern.h +++ b/FMD/AliFMDPattern.h @@ -79,7 +79,7 @@ public: private: /** Copy constructor - Not implemented. */ - AliFMDPatternDetector(const AliFMDPattern&); + AliFMDPatternDetector(const AliFMDPatternDetector&); /** Assignement operator -- Not implemented */ AliFMDPatternDetector& operator=(const AliFMDPatternDetector&); diff --git a/FMD/AliFMDRawReader.cxx b/FMD/AliFMDRawReader.cxx index a34175ed048..dcf488e0b87 100644 --- a/FMD/AliFMDRawReader.cxx +++ b/FMD/AliFMDRawReader.cxx @@ -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; diff --git a/FMD/AliFMDRawStream.cxx b/FMD/AliFMDRawStream.cxx index 4ad4ff6dbf6..89b4b83c1ad 100644 --- a/FMD/AliFMDRawStream.cxx +++ b/FMD/AliFMDRawStream.cxx @@ -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", diff --git a/FMD/Config.C b/FMD/Config.C index 6317aeb6aaa..240d3209fcf 100644 --- a/FMD/Config.C +++ b/FMD/Config.C @@ -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); diff --git a/FMD/Simulate.C b/FMD/Simulate.C index 5273269c72d..9d38d8928f5 100644 --- a/FMD/Simulate.C +++ b/FMD/Simulate.C @@ -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(); } diff --git a/FMD/scripts/DisplayDigits.C b/FMD/scripts/DisplayDigits.C index 5b293a4a368..5f162b4e988 100644 --- a/FMD/scripts/DisplayDigits.C +++ b/FMD/scripts/DisplayDigits.C @@ -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 index 00000000000..1e24aa03684 --- /dev/null +++ b/FMD/scripts/DrawESD.C @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/** @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(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 +// diff --git a/FMD/scripts/DrawHits.C b/FMD/scripts/DrawHits.C index 98e3a39afb1..6845cf4474f 100644 --- a/FMD/scripts/DrawHits.C +++ b/FMD/scripts/DrawHits.C @@ -20,6 +20,12 @@ #include #include #include +#include +#include +#include +#include +#include +#include /** @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); diff --git a/FMD/scripts/DrawHitsRecs.C b/FMD/scripts/DrawHitsRecs.C index 46b3bc13b45..5cbd4e16f29 100644 --- a/FMD/scripts/DrawHitsRecs.C +++ b/FMD/scripts/DrawHitsRecs.C @@ -27,6 +27,7 @@ #include #include #include +#include //____________________________________________________________________ /** @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; } diff --git a/FMD/scripts/GetXsection.C b/FMD/scripts/GetXsection.C index 32983398f6b..9b5a6350d5c 100644 --- a/FMD/scripts/GetXsection.C +++ b/FMD/scripts/GetXsection.C @@ -12,6 +12,16 @@ // // 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 #include #include @@ -22,6 +32,13 @@ #include #include #include +#include +#include +#include +#include +#include +#include +#include /** @defgroup xsec_script X-section script @ingroup FMD_script */ @@ -30,14 +47,148 @@ */ 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& t, + std::vector& 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 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(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_array; + typedef std::vector 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 index 00000000000..90e5095d5d3 --- /dev/null +++ b/FMD/scripts/XSection.C @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +/** @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& t, + std::vector& 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 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(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 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_array; + typedef std::vector 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 index 00000000000..a04a16ff261 --- /dev/null +++ b/FMD/scripts/checkSizes.sh @@ -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 < 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 -- 2.39.3