10 #include<TApplication.h>
14 #include<TClonesArray.h>
20 #include<FairRunAna.h>
22 #include<FairRtdbRun.h>
23 #include<FairRuntimeDb.h>
24 #include<FairParRootFileIo.h>
25 #include<TGeoVolume.h>
43 #include <TPolyLine.h>
53 cout <<
" analysis tool vers. 1.0 " << endl;
55 gROOT->Macro(
"$VMCWORKDIR/gconfig/rootlogon.C");
56 gSystem->Load(
"libSds");
57 gSystem->Load(
"libSdsReco");
58 gSystem->Load(
"libLmd");
59 gSystem->Load(
"libLmdReco");
60 gSystem->Load(
"libLmdTrk");
64 gROOT->SetStyle(
"Plain");
65 const Int_t NRGBs = 5;
66 const Int_t NCont = 255;
67 Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
68 Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
69 Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
70 Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
71 TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
72 gStyle->SetNumberContours(NCont);
73 gStyle->SetTitleFont(10*13+2,
"xyz");
74 gStyle->SetTitleSize(0.06,
"xyz");
75 gStyle->SetTitleOffset(0.8,
"y");
76 gStyle->SetTitleOffset(1.3,
"z");
77 gStyle->SetLabelFont(10*13+2,
"xyz");
78 gStyle->SetLabelSize(0.06,
"xyz");
79 gStyle->SetLabelOffset(0.009,
"xyz");
80 gStyle->SetPadBottomMargin(0.16);
81 gStyle->SetPadTopMargin(0.16);
82 gStyle->SetPadLeftMargin(0.10);
83 gStyle->SetPadRightMargin(0.10);
84 gStyle->SetOptTitle(1);
85 gStyle->SetOptStat(1);
87 gStyle->SetFrameFillColor(0);
88 gStyle->SetFrameFillStyle(0);
89 TGaxis::SetMaxDigits(3);
97 TString simMC = storePath +
"/Lumi_MC";
100 TChain tMC(
"pndsim");
106 TChain tdigiHits(
"pndsim");
107 tdigiHits.Add(DigiFile);
110 TString fakefile = storePath +
"/fake";
115 TString out = storePath +
"/analysis_results";
118 TFile *output_histfile =
new TFile(out,
"RECREATE");
127 FairRunAna *
fRun =
new FairRunAna();
128 fRun->SetInputFile(simMC);
129 fRun->AddFriend(DigiFile);
133 fRun->SetOutputFile(fakefile);
137 FairRuntimeDb*
rtdb = fRun->GetRuntimeDb();
138 FairParRootFileIo*
parInput1 =
new FairParRootFileIo(kTRUE);
139 parInput1->open(parFile.Data());
140 rtdb->setFirstInput(parInput1);
143 TClonesArray* true_tracks =
new TClonesArray(
"PndMCTrack");
144 tMC.SetBranchAddress(
"MCTrack", &true_tracks);
146 TClonesArray* true_points =
new TClonesArray(
"PndSdsMCPoint");
147 tMC.SetBranchAddress(
"LMDPoint", &true_points);
158 TH2D* hist_gen_events =
new TH2D(
"hist_gen_events",
"generated events", 100, 0., 10, 100, -3.1416, 3.1416);
159 TH2D* hist_acc_events =
new TH2D(
"hist_acc_events",
"accepted events", 100, 0., 10, 100, -3.1416, 3.1416);
161 TH2D* hist_sec_hit_theta[
nplanes];
164 TH2D* hist_sec_hit_edep[
nplanes];
165 TH2D* hist_prim_hit_edep[
nplanes];
167 TH2D* hist_angle_vs_origin[
nplanes];
168 for (
unsigned int iplane = 0; iplane <
nplanes; iplane++){
171 name <<
"hist_all_hits_plane_"<< iplane;
172 title <<
"Hit distribution of all hits at plane " << iplane;
173 hist_all_hits[iplane] =
new TH2D(name.str().c_str(), title.str().c_str(), 100, -10, 10, 100, -10, 10);
176 name <<
"hist_sec_hits_plane_"<< iplane;
177 title <<
"Hit distribution of secondary hits at plane " << iplane;
178 hist_sec_hits[iplane] =
new TH2D(name.str().c_str(), title.str().c_str(), 100, -10, 10, 100, -10, 10);
181 name <<
"hist_prim_hits_plane_"<< iplane;
182 title <<
"Hit distribution of primary hits at plane " << iplane;
183 hist_prim_hits[iplane] =
new TH2D(name.str().c_str(), title.str().c_str(), 100, -10, 10, 100, -10, 10);
186 name <<
"hist_sec_hit_theta_plane_"<< iplane;
187 title <<
"track angle distribution of secondary hits at plane " << iplane;
188 hist_sec_hit_theta[iplane] =
new TH2D(name.str().c_str(), title.str().c_str(), 100, 0, 1, 50, 0, 10);
191 name <<
"hist_sec_hit_origin_plane_"<< iplane;
192 title <<
"track origin z distribution of secondary hits at plane " << iplane;
193 hist_sec_hit_z[iplane] =
new TH2D(name.str().c_str(), title.str().c_str(), 50, -50, 50, 50, 0, 10);
197 name <<
"hist_sec_hit_edep_plane_"<< iplane;
198 title <<
"secondary edep at plane " << iplane;
199 hist_sec_hit_edep[iplane] =
new TH2D(name.str().c_str(), title.str().c_str(), 100, 0, 30, 50, 0, 10);
203 name <<
"hist_prim_hit_edep_plane_"<< iplane;
204 title <<
"primary edep at plane " << iplane;
205 hist_prim_hit_edep[iplane] =
new TH2D(name.str().c_str(), title.str().c_str(), 100, 0, 30, 50, 0, 10);
209 name <<
"hist_angle_vs_origin_plane_"<< iplane;
210 title <<
"hit angle vs track origin z " << iplane;
211 hist_angle_vs_origin[iplane] =
new TH2D(name.str().c_str(), title.str().c_str(), 200, -20, 50, 100, 0, 1);
215 int nEvents = tMC.GetEntries();
216 cout <<
" reading " << nEvents <<
" Events " << endl;
218 for (Int_t j = 0; j <
nEvents ; j++) {
221 tdigiHits.GetEntry(j);
222 const int nParticles = true_tracks->GetEntriesFast();
227 for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
233 primary_track = mctrk;
238 cout <<
"found " << npbar <<
" anti-protons" << endl;
240 int nHits = true_points->GetEntriesFast();
242 for (Int_t iHit = 0; iHit <
nHits; iHit++) {
244 if (mcpoint->GetTrackID() < 0)
continue;
246 mcpoint->GetTrackID());
249 TVector3 _mctrack(mcpoint->GetPx(), mcpoint->GetPy(),
259 int ihalf, iplane, imodule, iside, idie, isensor;
260 lmddim->
Get_sensor_by_id(sensID, ihalf, iplane, imodule, iside, idie, isensor);
263 double x = mcpoint_lmd.X();
264 double y = mcpoint_lmd.Y();
265 double r =
sqrt(x*x+y*y);
266 double Eloss = mcpoint->GetEnergyLoss()*1e6;
268 double startpos = posMC_lmd.Z();
269 double startr =
sqrt(posMC_lmd.X()*posMC_lmd.X()+posMC_lmd.Y()*posMC_lmd.Y());
270 double trackangle = mctrack_lmd.Theta();
273 if (startr < 3.6) hist_all_hits[iplane]->Fill(x, y);
274 if (startr < 3.6) hist_angle_vs_origin[iplane]->Fill(startpos, trackangle);
277 hist_prim_hit_edep[iplane]->Fill(Eloss,r);
279 if (startr < 3.6) hist_prim_hits[iplane]->Fill(x, y);
284 hist_sec_hit_edep[iplane]->Fill(Eloss,r);
286 if (startr < 3.6) hist_sec_hits[iplane]->Fill(x, y);
287 if (startr < 3.6) hist_sec_hit_theta[iplane]->Fill(trackangle, r);
288 if (startr < 3.6) hist_sec_hit_z[iplane]->Fill(startpos, r);
294 cout <<
"found " << nHits <<
" hit(s) with " << nSdsHits
295 <<
" sds hit(s)" << endl;
298 output_histfile->cd();
300 TCanvas canvas(
"canvas",
"Results", 800, 800);
301 canvas.SetMargin(0.15, 0.15, 0.15, 0.15);
303 hist_gen_events->Draw(
"COLZ");
305 hist_gen_events->GetXaxis()->SetTitle(
"#theta / mrad");
306 hist_gen_events->GetYaxis()->SetTitle(
"#phi / phi");
307 canvas.Print(
"Secondary_results.ps(");
308 hist_acc_events->Divide(hist_gen_events);
309 hist_acc_events->Draw(
"COLZ");
311 hist_acc_events->GetXaxis()->SetTitle(
"#theta / mrad");
312 hist_acc_events->GetYaxis()->SetTitle(
"#phi / rad");
313 canvas.Print(
"Secondary_results.ps(");
315 for (
unsigned int iplane = 0; iplane <
nplanes; iplane++){
318 hist_all_hits[iplane]->Draw(
"COLZ");
320 hist_all_hits[iplane]->GetXaxis()->SetTitle(
"x[cm]");
321 hist_all_hits[iplane]->GetYaxis()->SetTitle(
"y[cm]");
322 canvas.Print(
"Secondary_results.ps(");
323 hist_all_hits[iplane]->Write();
325 for (
unsigned int iplane = 0; iplane <
nplanes; iplane++){
328 hist_sec_hits[iplane]->Draw(
"COLZ");
330 hist_sec_hits[iplane]->GetXaxis()->SetTitle(
"x[cm]");
331 hist_sec_hits[iplane]->GetYaxis()->SetTitle(
"y[cm]");
332 hist_sec_hits[iplane]->Write();
333 canvas.Print(
"Secondary_results.ps(");
335 for (
unsigned int iplane = 0; iplane <
nplanes; iplane++){
338 hist_prim_hits[iplane]->Draw(
"COLZ");
340 hist_prim_hits[iplane]->GetXaxis()->SetTitle(
"x[cm]");
341 hist_prim_hits[iplane]->GetYaxis()->SetTitle(
"y[cm]");
342 hist_prim_hits[iplane]->Write();
343 canvas.Print(
"Secondary_results.ps(");
345 for (
unsigned int iplane = 0; iplane <
nplanes; iplane++){
348 hist_sec_hit_theta[iplane]->Draw(
"COLZ");
350 hist_sec_hit_theta[iplane]->GetXaxis()->SetTitle(
"#theta[rad]");
351 hist_sec_hit_theta[iplane]->GetYaxis()->SetTitle(
"r[cm]");
352 hist_sec_hit_theta[iplane]->Write();
353 canvas.Print(
"Secondary_results.ps(");
355 for (
unsigned int iplane = 0; iplane <
nplanes; iplane++){
358 hist_sec_hit_z[iplane]->Draw(
"COLZ");
360 hist_sec_hit_z[iplane]->GetXaxis()->SetTitle(
"z[cm]");
361 hist_sec_hit_z[iplane]->GetYaxis()->SetTitle(
"r[cm]");
362 hist_sec_hit_z[iplane]->Write();
363 canvas.Print(
"Secondary_results.ps(");
365 for (
unsigned int iplane = 0; iplane <
nplanes; iplane++){
368 hist_sec_hit_edep[iplane]->Draw(
"COLZ");
370 hist_sec_hit_edep[iplane]->GetXaxis()->SetTitle(
"Eloss[keV]");
371 hist_sec_hit_edep[iplane]->GetYaxis()->SetTitle(
"r[cm]");
372 hist_sec_hit_edep[iplane]->Write();
373 canvas.Print(
"Secondary_results.ps(");
375 for (
unsigned int iplane = 0; iplane <
nplanes; iplane++){
378 hist_prim_hit_edep[iplane]->Draw(
"COLZ");
380 hist_prim_hit_edep[iplane]->GetXaxis()->SetTitle(
"Eloss[keV]");
381 hist_prim_hit_edep[iplane]->GetYaxis()->SetTitle(
"r[cm]");
382 hist_prim_hit_edep[iplane]->Write();
383 canvas.Print(
"Secondary_results.ps(");
385 for (
unsigned int iplane = 0; iplane <
nplanes; iplane++){
388 stringstream name_prim;
389 name_prim <<
"hist_prim_hit_edep_projX_plane_" << iplane;
390 TH1D* hist_prim = hist_prim_hit_edep[iplane]->ProjectionX(name_prim.str().c_str());
391 stringstream name_sec;
392 name_sec <<
"hist_sec_hit_edep_projX_plane_" << iplane;
393 TH1D* hist_sec = hist_sec_hit_edep[iplane]->ProjectionX(name_sec.str().c_str());
394 hist_prim->DrawNormalized();
395 hist_sec->DrawNormalized(
"same");
399 canvas.Print(
"Secondary_results.ps(");
401 for (
unsigned int iplane = 0; iplane <
nplanes; iplane++){
404 hist_angle_vs_origin[iplane]->Draw(
"COLZ");
406 hist_angle_vs_origin[iplane]->GetXaxis()->SetTitle(
"z[cm]");
407 hist_angle_vs_origin[iplane]->GetYaxis()->SetTitle(
"#theta[rad]");
408 hist_angle_vs_origin[iplane]->Write();
409 canvas.Print(
"Secondary_results.ps(");
412 for (
unsigned int iplane = 0; iplane <
nplanes; iplane++){
417 name <<
"hist_sec_hist_norm_prim_plane_" << iplane;
418 title <<
"Secondary hit distribution normalized to primaries at plane " << iplane;
419 TH2D* _hist =
new TH2D(*hist_sec_hits[iplane]);
420 _hist->SetNameTitle(name.str().c_str(), title.str().c_str());
421 _hist->Divide(hist_prim_hits[iplane]);
423 _hist->SetMaximum(.2);
426 canvas.Print(
"Secondary_results.ps(");
428 for (
unsigned int iplane = 0; iplane <
nplanes; iplane++){
433 name <<
"hist_sec_hist_norm_all_plane_" << iplane;
434 title <<
"Secondary hit distribution normalized to all hits at plane " << iplane;
435 TH2D* _hist =
new TH2D(*hist_sec_hits[iplane]);
436 _hist->SetNameTitle(name.str().c_str(), title.str().c_str());
437 _hist->Divide(hist_all_hits[iplane]);
439 _hist->SetMaximum(.2);
442 canvas.Print(
"Secondary_results.ps(");
445 for (
unsigned int iplane = 1; iplane <
nplanes; iplane++){
450 name <<
"hist_sec_hist_norm_plane0_plane_" << iplane;
451 title <<
"Secondary hit distribution normalized to first plane at plane " << iplane;
452 TH2D* _hist =
new TH2D(*hist_sec_hits[iplane]);
453 _hist->SetNameTitle(name.str().c_str(), title.str().c_str());
454 _hist->Divide(hist_sec_hits[0]);
459 canvas.Print(
"Secondary_results.ps(");
462 canvas.Print(
"Secondary_results.ps)");
464 cout <<
" converting ps to pdf " << endl;
465 system(
"ps2pdf Secondary_results.ps");
466 cout <<
" done " << endl;
470 output_histfile->Close();
474 if ((
int) (
last_percent * 100) == (
int) (percent * 100))
480 for (
int i = 0;
i < len; ++
i) {
481 if (
i < static_cast<int> (len * percent)) {
487 cout <<
"[" << progress <<
"] " << (
static_cast<int> (100 * percent))
int main(int argc, char **argv)
void scattered_particles()
friend F32vec4 sqrt(const F32vec4 &a)
static PndLmdDim * Instance()
TVector3 Transform_global_to_lmd_local(const TVector3 &point, bool isvector=false, bool aligned=true)
TVector3 GetMomentum() const
void DrawProgressBar(int len, double percent)
Int_t GetSensorID() const
Bool_t IsGeneratorCreated(void) const
void Get_sensor_by_id(const int sensor_id, int &ihalf, int &iplane, int &imodule, int &iside, int &idie, int &isensor)
FairParRootFileIo * parInput1
void Read_transformation_matrices(string filename="", bool aligned=true, int version_number=geometry_version)
TVector3 GetPosition() const
TVector3 GetStartVertex() const