FairRoot/PandaRoot
Functions
scattered_particles.cxx File Reference
#include <TString.h>
#include <TH1.h>
#include <TH2.h>
#include <sstream>
#include <iostream>
#include <vector>
#include <map>
#include <TVectorT.h>
#include <TCanvas.h>
#include <TApplication.h>
#include <TChain.h>
#include <TROOT.h>
#include <TSystem.h>
#include <TClonesArray.h>
#include <TRotation.h>
#include <TFile.h>
#include <TGaxis.h>
#include <PndMCTrack.h>
#include <PndSdsMCPoint.h>
#include <FairRunAna.h>
#include <FairRtdbRun.h>
#include <FairRuntimeDb.h>
#include <FairParRootFileIo.h>
#include <TGeoVolume.h>
#include <TGraph.h>
#include <TF1.h>
#include <algorithm>
#include <fstream>
#include <TStyle.h>
#include <TColor.h>
#include <TImage.h>
#include <TLine.h>
#include <stdlib.h>
#include <TPolyLine.h>
#include <PndLmdDim.h>

Go to the source code of this file.

Functions

double last_percent (-1.)
 
void DrawProgressBar (int len, double percent)
 
void scattered_particles ()
 
int main ()
 

Function Documentation

void DrawProgressBar ( int  len,
double  percent 
)
double last_percent ( 1.)
int main ( void  )

Definition at line 493 of file scattered_particles.cxx.

References scattered_particles().

493  {
494  //TApplication myapp("myapp", 0, 0);
496  //myapp.Run();
497  return 0;
498 }
void scattered_particles()
void scattered_particles ( )

Definition at line 52 of file scattered_particles.cxx.

References DigiFile, Double_t, DrawProgressBar(), fRun, PndLmdDim::Get_sensor_by_id(), PndMCTrack::GetMomentum(), PndMCTrack::GetPdgCode(), PndSdsMCPoint::GetPosition(), PndSdsMCPoint::GetSensorID(), PndMCTrack::GetStartVertex(), PndLmdDim::Instance(), PndMCTrack::IsGeneratorCreated(), PndLmdDim::n_planes, name, nEvents, nHits, nplanes(), out, parFile, parInput1, r, PndLmdDim::Read_transformation_matrices(), rtdb, sqrt(), startEvent, storePath, PndLmdDim::Transform_global_to_lmd_local(), TString, x, and y.

Referenced by main().

52  {
53  cout << " analysis tool vers. 1.0 " << endl;
54 
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");
61 
62  TPad foo; // never remove this line :-)))
63  if(1){
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);
86  gROOT->ForceStyle();
87  gStyle->SetFrameFillColor(0);
88  gStyle->SetFrameFillStyle(0);
89  TGaxis::SetMaxDigits(3);
90  }
91 
92  // ------------------------------------------------------------------------
93 
94  TString storePath = ".";
95  TString startEvent = "";
96  // ---- Input files --------------------------------------------------------
97  TString simMC = storePath + "/Lumi_MC";
98  simMC += startEvent;
99  simMC += ".root";
100  TChain tMC("pndsim");
101  tMC.Add(simMC);
102 
103  TString DigiFile = storePath + "/Lumi_digi";
104  DigiFile += startEvent;
105  DigiFile += ".root";
106  TChain tdigiHits("pndsim");
107  tdigiHits.Add(DigiFile);
108 
109  // ---- Output file ----------------------------------------------------------------
110  TString fakefile = storePath + "/fake";
111  fakefile += startEvent;
112  fakefile += ".root";
113 
114  // ---- Output file ----------------------------------------------------------------
115  TString out = storePath + "/analysis_results";
116  out += startEvent;
117  out += ".root";
118  TFile *output_histfile = new TFile(out, "RECREATE");
119  // ---------------------------------------------------------------------------------
120 
121  // Parameter file << needed for geane back tracking
122  TString parFile = storePath+"/Lumi_Params";
123  parFile += startEvent;
124  parFile += ".root";
125 
126  // ----- Reconstruction run -------------------------------------------
127  FairRunAna *fRun = new FairRunAna();
128  fRun->SetInputFile(simMC);
129  fRun->AddFriend(DigiFile);
130  //fRun->AddFriend(RecoFile);
131  //fRun->AddFriend(CandFile);
132  //fRun->AddFriend(TrkFile);
133  fRun->SetOutputFile(fakefile);
134  // ------------------------------------------------------------------------
135 
136  // ----- Parameter database --------------------------------------------
137  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
138  FairParRootFileIo* parInput1 = new FairParRootFileIo(kTRUE);
139  parInput1->open(parFile.Data());
140  rtdb->setFirstInput(parInput1);
141 
142  //--- MC info -----------------------------------------------------------------
143  TClonesArray* true_tracks = new TClonesArray("PndMCTrack");
144  tMC.SetBranchAddress("MCTrack", &true_tracks); //True Track to compare
145 
146  TClonesArray* true_points = new TClonesArray("PndSdsMCPoint");
147  tMC.SetBranchAddress("LMDPoint", &true_points); //True Points to compare
148 
149  //----------------------------------------------------------------------------------
150 
151 
152  PndLmdDim* lmddim = PndLmdDim::Instance();
154  const unsigned int nplanes = lmddim->n_planes;
155  TH2D* hist_all_hits[nplanes];
156  TH2D* hist_sec_hits[nplanes];
157  TH2D* hist_prim_hits[nplanes];
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);
160 
161  TH2D* hist_sec_hit_theta[nplanes];
162  TH2D* hist_sec_hit_z[nplanes];
163 
164  TH2D* hist_sec_hit_edep[nplanes];
165  TH2D* hist_prim_hit_edep[nplanes];
166 
167  TH2D* hist_angle_vs_origin[nplanes];
168  for (unsigned int iplane = 0; iplane < nplanes; iplane++){
169  stringstream name;
170  stringstream title;
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);
174  name.str("");
175  title.str("");
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);
179  name.str("");
180  title.str("");
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);
184  name.str("");
185  title.str("");
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);
189  name.str("");
190  title.str("");
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);
194 
195  name.str("");
196  title.str("");
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);
200 
201  name.str("");
202  title.str("");
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);
206 
207  name.str("");
208  title.str("");
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);
212  }
213 
214 
215  int nEvents = tMC.GetEntries();
216  cout << " reading " << nEvents << " Events " << endl;
217 
218  for (Int_t j = 0; j < nEvents ; j++) {
219  DrawProgressBar(50, (j + 1) / ((double) nEvents));
220  tMC.GetEntry(j);
221  tdigiHits.GetEntry(j);
222  const int nParticles = true_tracks->GetEntriesFast();
223 
224  PndMCTrack* primary_track = NULL;
225 
226  int npbar = 0;
227  for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
228  // if(nParticles!=4) continue;
229  PndMCTrack *mctrk = (PndMCTrack*) true_tracks->At(iParticle);
230  Int_t mcID = mctrk->GetPdgCode();
231  if (mctrk->IsGeneratorCreated()) {//mcID == -2212) {
232  npbar++;
233  primary_track = mctrk;
234  hist_gen_events->Fill(mctrk->GetMomentum().Theta()*1e3, mctrk->GetMomentum().Phi());
235  }
236  }
237  if (0 && npbar != 1)
238  cout << "found " << npbar << " anti-protons" << endl;
239 
240  int nHits = true_points->GetEntriesFast();
241  int nSdsHits = 0;
242  for (Int_t iHit = 0; iHit < nHits; iHit++) {
243  PndSdsMCPoint* mcpoint = (PndSdsMCPoint*) true_points->At(iHit);
244  if (mcpoint->GetTrackID() < 0) continue;
245  PndMCTrack *mctrk = (PndMCTrack*) true_tracks->At(
246  mcpoint->GetTrackID());
247  if (mcpoint) {
248  TVector3 _mcpoint((mcpoint->GetPosition())); // position at the sensor entrance
249  TVector3 _mctrack(mcpoint->GetPx(), mcpoint->GetPy(),
250  mcpoint->GetPz()); // momentum of the track at the entrance
251 
252  TVector3 momMC = mctrk->GetMomentum(); // momentum in the primary vertex
253  TVector3 posMC = mctrk->GetStartVertex(); // position of the primary vertex
254 
255  const TVector3& mcpoint_lmd = lmddim->Transform_global_to_lmd_local(_mcpoint);
256  const TVector3& mctrack_lmd = lmddim->Transform_global_to_lmd_local(_mctrack, true);
257 
258  int sensID = mcpoint->GetSensorID();
259  int ihalf, iplane, imodule, iside, idie, isensor;
260  lmddim->Get_sensor_by_id(sensID, ihalf, iplane, imodule, iside, idie, isensor);
261  nSdsHits++;
262 
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;
267  const TVector3& posMC_lmd = lmddim->Transform_global_to_lmd_local(posMC);
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();
271 
272  if (Eloss > 5.5){
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);
275  }
276  if (mctrk->IsGeneratorCreated()) { // take only anti-protons from the generator
277  hist_prim_hit_edep[iplane]->Fill(Eloss,r);
278  if (Eloss > 5.5){ // minimal E deposit for protons to be detected in HV-MAPS (by T.Weber)
279  if (startr < 3.6) hist_prim_hits[iplane]->Fill(x, y);
280  if (primary_track)
281  hist_acc_events->Fill(primary_track->GetMomentum().Theta()*1e3, primary_track->GetMomentum().Phi());
282  }
283  } else { // it is a secondary
284  hist_sec_hit_edep[iplane]->Fill(Eloss,r);
285  if (Eloss > 5.5){
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);
289  }
290  }
291  }
292  }
293  if (0 && nHits > 0)
294  cout << "found " << nHits << " hit(s) with " << nSdsHits
295  << " sds hit(s)" << endl;
296  }
297 
298  output_histfile->cd();
299 
300  TCanvas canvas("canvas", "Results", 800, 800);
301  canvas.SetMargin(0.15, 0.15, 0.15, 0.15);
302  canvas.cd();
303  hist_gen_events->Draw("COLZ");
304  gPad->Update();
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");
310  gPad->Update();
311  hist_acc_events->GetXaxis()->SetTitle("#theta / mrad");
312  hist_acc_events->GetYaxis()->SetTitle("#phi / rad");
313  canvas.Print("Secondary_results.ps(");
314 
315  for (unsigned int iplane = 0; iplane < nplanes; iplane++){
316  canvas.cd();
317  gPad->Clear();
318  hist_all_hits[iplane]->Draw("COLZ");
319  gPad->Update();
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();
324  }
325  for (unsigned int iplane = 0; iplane < nplanes; iplane++){
326  canvas.cd();
327  gPad->Clear();
328  hist_sec_hits[iplane]->Draw("COLZ");
329  gPad->Update();
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(");
334  }
335  for (unsigned int iplane = 0; iplane < nplanes; iplane++){
336  canvas.cd();
337  gPad->Clear();
338  hist_prim_hits[iplane]->Draw("COLZ");
339  gPad->Update();
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(");
344  }
345  for (unsigned int iplane = 0; iplane < nplanes; iplane++){
346  canvas.cd();
347  gPad->Clear();
348  hist_sec_hit_theta[iplane]->Draw("COLZ");
349  gPad->Update();
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(");
354  }
355  for (unsigned int iplane = 0; iplane < nplanes; iplane++){
356  canvas.cd();
357  gPad->Clear();
358  hist_sec_hit_z[iplane]->Draw("COLZ");
359  gPad->Update();
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(");
364  }
365  for (unsigned int iplane = 0; iplane < nplanes; iplane++){
366  canvas.cd();
367  gPad->Clear();
368  hist_sec_hit_edep[iplane]->Draw("COLZ");
369  gPad->Update();
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(");
374  }
375  for (unsigned int iplane = 0; iplane < nplanes; iplane++){
376  canvas.cd();
377  gPad->Clear();
378  hist_prim_hit_edep[iplane]->Draw("COLZ");
379  gPad->Update();
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(");
384  }
385  for (unsigned int iplane = 0; iplane < nplanes; iplane++){
386  canvas.cd();
387  gPad->Clear();
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");
396  gPad->Update();
397  hist_prim->Write();
398  hist_sec->Write();
399  canvas.Print("Secondary_results.ps(");
400  }
401  for (unsigned int iplane = 0; iplane < nplanes; iplane++){
402  canvas.cd();
403  gPad->Clear();
404  hist_angle_vs_origin[iplane]->Draw("COLZ");
405  gPad->Update();
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(");
410  }
411  // normalize
412  for (unsigned int iplane = 0; iplane < nplanes; iplane++){
413  canvas.cd();
414  gPad->Clear();
415  stringstream name;
416  stringstream title;
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]);
422  _hist->Draw("COLZ");
423  _hist->SetMaximum(.2);
424  gPad->Update();
425  _hist->Write();
426  canvas.Print("Secondary_results.ps(");
427  }
428  for (unsigned int iplane = 0; iplane < nplanes; iplane++){
429  canvas.cd();
430  gPad->Clear();
431  stringstream name;
432  stringstream title;
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]);
438  _hist->Draw("COLZ");
439  _hist->SetMaximum(.2);
440  gPad->Update();
441  _hist->Write();
442  canvas.Print("Secondary_results.ps(");
443  }
444  // normalize
445  for (unsigned int iplane = 1; iplane < nplanes; iplane++){
446  canvas.cd();
447  gPad->Clear();
448  stringstream name;
449  stringstream title;
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]);
455  _hist->Draw("COLZ");
456  //_hist->SetMaximum(.2);
457  gPad->Update();
458  _hist->Write();
459  canvas.Print("Secondary_results.ps(");
460  }
461  gPad->Clear();
462  canvas.Print("Secondary_results.ps)");
463 
464  cout << " converting ps to pdf " << endl;
465  system("ps2pdf Secondary_results.ps");
466  cout << " done " << endl;
467 
468  //output_histfile->Write();
469 
470  output_histfile->Close();
471 }
const int nplanes(120)
double r
Definition: RiemannTest.C:14
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
static PndLmdDim * Instance()
Definition: PndLmdDim.cxx:249
TVector3 Transform_global_to_lmd_local(const TVector3 &point, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3113
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
void DrawProgressBar(int len, double percent)
Int_t startEvent
TString storePath
Int_t GetSensorID() const
Definition: PndSdsMCPoint.h:89
Bool_t IsGeneratorCreated(void) const
Definition: PndMCTrack.h:84
FairRunAna * fRun
Definition: hit_dirc.C:58
TString DigiFile
int nHits
Definition: RiemannTest.C:16
Double_t
TString parFile
Definition: hit_dirc.C:14
void Get_sensor_by_id(const int sensor_id, int &ihalf, int &iplane, int &imodule, int &iside, int &idie, int &isensor)
Definition: PndLmdDim.h:526
Int_t nEvents
Definition: hit_dirc.C:11
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
TFile * out
Definition: reco_muo.C:20
TString name
FairParRootFileIo * parInput1
Definition: hit_dirc.C:67
void Read_transformation_matrices(string filename="", bool aligned=true, int version_number=geometry_version)
Definition: PndLmdDim.cxx:1515
TVector3 GetPosition() const
Definition: PndSdsMCPoint.h:90
Double_t x
Double_t y
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
unsigned int n_planes
Definition: PndLmdDim.h:262