FairRoot/PandaRoot
Functions
radiation.cxx File Reference
#include <iostream>
#include <sstream>
#include <string>
#include <TString.h>
#include <TChain.h>
#include <TClonesArray.h>
#include <TROOT.h>
#include <TVector3.h>
#include <TH1.h>
#include <TH2.h>
#include <TApplication.h>
#include <TCanvas.h>
#include <FairRadMapPoint.h>
#include <PndLmdDim.h>

Go to the source code of this file.

Functions

double last_percent (0)
 
void DrawProgressBar (int len, double percent)
 
int main (int nargs, char **args)
 

Function Documentation

void DrawProgressBar ( int  len,
double  percent 
)

Definition at line 58 of file radiation.cxx.

References i, and last_percent().

58  {
59  if ((int)(last_percent*100) == (int)(percent*100)) return;
60  //cout << " drawing " << endl;
61  cout << "\x1B[2K"; // Erase the entire current line.
62  cout << "\x1B[0E"; // Move to the beginning of the current line.
63  string progress;
64  for (int i = 0; i < len; ++i) {
65  if (i < static_cast<int>(len * percent)) {
66  progress += "=";
67  } else {
68  progress += " ";
69  }
70  }
71  cout << "[" << progress << "] " << (static_cast<int>(100 * percent)) << "%";
72  flush(cout); // Required.
73  last_percent = percent;
74 }
double last_percent(-1.)
Int_t i
Definition: run_full.C:25
double last_percent ( )
int main ( int  nargs,
char **  args 
)

Definition at line 76 of file radiation.cxx.

References DrawProgressBar(), PndLmdDim::Get_instance(), MCFile, nEvents, PndLmdDim::Read_transformation_matrices(), startEvent, storePath, and TString.

76  {
77  // load first all dictionaries for ROOT in order
78  // to read the files correctly
79  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
80  TApplication myapp("myapp",0,0);
81  TString storePath = ".";
82  // set the storepath to somewhere else, if provided
83  if (nargs == 2) {
84  storePath = args[1];
85  } else {
86  if (nargs == 1)
87  cout << " Hint: you may change the standard path from "
88  << storePath.Data()
89  << " to any other folder by providing the path as an argument! "
90  << endl;
91  if (nargs > 2)
92  cout << " Error: too many arguments! " << endl;
93  }
94  TString startEvent = "";
95  // Input file (MC events)
96  TString MCFile = storePath + "/Lumi_MC";
97  MCFile += startEvent;
98  MCFile += ".root";
99  TChain tMC("pndsim");
100  tMC.Add(MCFile);
101 
102  std::cout << "MCFile : " << MCFile.Data() << std::endl;
103 
104  //--- MC info -----------------------------------------------------------------
105  TClonesArray* true_tracks = new TClonesArray("PndMCTrack");
106  tMC.SetBranchAddress("MCTrack", &true_tracks); //True Tracks to compare
107 
108  TClonesArray* true_points = new TClonesArray("PndSdsMCPoint");
109  tMC.SetBranchAddress("LMDPoint", &true_points); //True Points to compare
110 
111  TClonesArray* radmap_points = new TClonesArray("FairRadMapPoint");
112  tMC.SetBranchAddress("RadMap", & radmap_points); // Radiation map points
113 
114  TClonesArray* radlen_points = new TClonesArray("FairRadLenPoint");
115  tMC.SetBranchAddress("RadLen", & radlen_points); // Radiation length points
116  //----------------------------------------------------------------------------------
117 
118  TCanvas canvas_rad_map_x_y("canvas_rad_map_x_y", "radiation map x y", 600, 600);
119  TH2D radiation_map_x_y("radiation_map_x_y", "radiation map x y", 500, 15, 45, 500, -15, 15);
120  radiation_map_x_y.Draw("COLZ");
121  radiation_map_x_y.SetXTitle("X [cm]");
122  radiation_map_x_y.SetYTitle("Y [cm]");
123  radiation_map_x_y.SetZTitle("dose");
124  canvas_rad_map_x_y.SetLogz();
125 
126  TCanvas canvas_rad_map_x_z("canvas_rad_map_x_z", "radiation map x z", 1000, 600);
127  TH2D radiation_map_x_z("radiation_map_x_z", "radiation map x z", 1000, 1100, 1200, 400, 10, 50);
128  radiation_map_x_z.Draw("COLZ");
129  radiation_map_x_z.SetXTitle("Z [cm]");
130  radiation_map_x_z.SetYTitle("X [cm]");
131  radiation_map_x_z.SetZTitle("dose");
132  canvas_rad_map_x_z.SetLogz();
133 
134  TCanvas canvas_rad_map_x_y_plane_0("canvas_rad_map_x_y_plane_0", "radiation map x y of plane 0", 600, 600);
135  TH2D radiation_map_x_y_plane_0("radiation_map_x_y_plane_0", "radiation map x y of plane 0", 150, 10, 40, 150, -15, 15);
136  radiation_map_x_y_plane_0.Draw("COLZ");
137  radiation_map_x_y_plane_0.SetXTitle("X [cm]");
138  radiation_map_x_y_plane_0.SetYTitle("Y [cm]");
139  radiation_map_x_y_plane_0.SetZTitle("dose");
140  canvas_rad_map_x_y_plane_0.SetLogz();
141 
142  // initialize the lmd dimension class
144  lmddim.Read_transformation_matrices("matrices_perfect.txt", false);
145  lmddim.Read_transformation_matrices("matrices.txt", true);
146 
147  int nEvents = tMC.GetEntries();
148  cout << " reading " << nEvents << " Events " << endl;
149  for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
150  DrawProgressBar(50, (iEvent+1)/((double)nEvents));
151  tMC.GetEntry(iEvent);
152 
153  //PndSdsMCPoint* mcpoint = (PndSdsMCPoint*) true_points->At(iHit);
154  int nradmappoints = radmap_points->GetEntriesFast();
155  for (int iradmappoint = 0; iradmappoint < nradmappoints; iradmappoint++){
156  FairRadMapPoint* radmap_point = (FairRadMapPoint*) radmap_points->At(iradmappoint);
157  const TVector3& position = radmap_point->GetPosition();
158  double dose = radmap_point->GetDose();
159 
160  int detID = radmap_point->GetDetectorID();
161 
162  //int ihalf, iplane, imodule, iside, idie, isensor;
163  //lmddim.Get_sensor_by_id(detID, ihalf, iplane, imodule, iside, idie, isensor);
164 
165  //if (iplane == 0) does not work
166  if (dose > 0){
167  if (position.Z() > 1105 && position.Z() < 1130)
168  radiation_map_x_y_plane_0.Fill(position.X(), position.Y(), dose);
169  radiation_map_x_y.Fill(position.X(), position.Y(), dose);
170  radiation_map_x_z.Fill(position.Z(), position.X(), dose);
171  }
172  }
173 
174  if ((iEvent%1000)==0){
175  canvas_rad_map_x_y.cd();
176  radiation_map_x_y.Draw("COLZ");
177  canvas_rad_map_x_y.Update();
178  canvas_rad_map_x_y_plane_0.cd();
179  radiation_map_x_y_plane_0.Draw("COLZ");
180  canvas_rad_map_x_y_plane_0.Update();
181  canvas_rad_map_x_z.cd();
182  radiation_map_x_z.Draw("COLZ");
183  canvas_rad_map_x_z.Update();
184  }
185  }
186 
187  canvas_rad_map_x_y.Print("dose_x_y.pdf");
188  canvas_rad_map_x_y_plane_0.Print("dose_x_y_plane_0.pdf");
189  canvas_rad_map_x_z.Print("dose_x_z.pdf");
190 
191 
192 
193 /*
194  TString DigiFile = storePath + "/Lumi_digi";
195  DigiFile += startEvent;
196  DigiFile += ".root";
197  // Digi file
198  TString RecoFile = storePath + "/Lumi_reco";
199  RecoFile += startEvent;
200  RecoFile += ".root";
201  // TCand file
202  TString CandFile = storePath + "/Lumi_TCand";
203  CandFile += startEvent;
204  CandFile += ".root";
205  // Parameter file << this one
206  TString parFile = storePath + "/Lumi_Params";
207  parFile += startEvent;
208  parFile += ".root";
209  // Track file
210  TString TrkFile = storePath + "/Lumi_Track";
211  TrkFile += startEvent;
212  TrkFile += ".root";
213 
214  // ---- Load libraries -------------------------------------------------
215  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
216  gSystem->Load("libSds");
217  gSystem->Load("libSdsReco");
218  gSystem->Load("libLmd");
219  gSystem->Load("libLmdReco");
220  gSystem->Load("libLmdTrk");
221  // ------------------------------------------------------------------------
222  // Output file
223  TString outFile = storePath + "/Lumi_Testing_";
224  outFile += startEvent;
225  outFile += ".root";
226 
227  // Output histogram file
228  TString outFilehist = storePath + "/Lumi_performance_histograms_";
229  outFilehist += startEvent;
230  outFilehist += ".root";
231 
232  std::cout << "DigiFile: " << DigiFile.Data() << std::endl;
233  std::cout << "RecoFile: " << RecoFile.Data() << std::endl;
234  std::cout << "TCandFile: " << CandFile.Data() << std::endl;
235  std::cout << "TrackFile: " << TrkFile.Data() << std::endl;
236  std::cout << "TestingFile: " << outFile.Data() << std::endl;
237  // --- Now choose concrete engines for the different tasks -------------
238  // ------------------------------------------------------------------------
239 
240 
241  // In general, the following parts need not be touched
242  // ========================================================================
243 
244 
245  // ----- Timer --------------------------------------------------------
246  TStopwatch timer;
247  timer.Start();
248  // ------------------------------------------------------------------------
249 
250 
251  // ----- Reconstruction run -------------------------------------------
252  FairRunAna *fRun = new FairRunAna();
253  fRun->SetInputFile(MCFile);
254  //fRun->AddFriend(DigiFile);
255  //fRun->AddFriend(RecoFile);
256  //fRun->AddFriend(CandFile);
257  //fRun->AddFriend(TrkFile);
258  fRun->SetOutputFile(outFile);
259  // ------------------------------------------------------------------------
260 
261  FairGeane *Geane = new FairGeane();
262  fRun->AddTask(Geane);
263 
264  // ----- Parameter database --------------------------------------------
265  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
266  FairParRootFileIo* parInput1 = new FairParRootFileIo(kTRUE);
267  parInput1->open(parFile.Data());
268  rtdb->setFirstInput(parInput1);
269 
270  // =========================================================================
271  // ====== Geane Back-Propagating ======
272  // =========================================================================
273 
274  // ----- SDS hit producer --------------------------------------------
275 
276  // PndEmcMapper *emcMap = PndEmcMapper::Instance(6);
277  Double_t fpBeam = 1.5;
278  TVector3 IP(0, 0, 0);
279  PndLmdPerformanceTask* lmdperformance = new PndLmdPerformanceTask();
280  lmdperformance->SetVerbose(0);
281  fRun->AddTask(lmdperformance);
282  rtdb->setOutput(parInput1);
283  //lmdperformance->SetHistFilename(outFilehist);
284  rtdb->print();
285  // ===== End of Geane =====
286  // =========================================================================
287 
288  // ----- Intialise and run --------------------------------------------
289  fRun->Init();
290  // // PndEmcMapper *emcMap = PndEmcMapper::Instance(6);
291  // PndEmcMapper *emcMap = PndEmcMapper::Instance();
292  // //Geane->SetField(fRun->GetField());
293  FairRootManager* ioman = FairRootManager::Instance();
294  if (!ioman) {
295  std::cout << "Error: "
296  << "RootManager not instantiated!" << std::endl;
297  return 1;
298  }
299  fRun->Run(0, ioman->CheckMaxEventNo()); // hmm does not work yet, seems to return 0 but then is runs over all what was intended.
300  // ------------------------------------------------------------------------
301  rtdb->saveOutput();
302  rtdb->print();
303  // ----- Finish -------------------------------------------------------
304  timer.Stop();
305  Double_t rtime = timer.RealTime();
306  Double_t ctime = timer.CpuTime();
307  cout << endl << endl;
308  cout << "Macro finished succesfully." << endl;
309  cout << "Output file is " << outFile << endl;
310  cout << "Parameter file is " << parFile << endl;
311  cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
312  cout << endl;
313  // ------------------------------------------------------------------------
314 
315  //myapp.Run();
316 */
317  return 0;
318 }
void DrawProgressBar(int len, double percent)
Int_t startEvent
TString storePath
Int_t nEvents
Definition: hit_dirc.C:11
void Read_transformation_matrices(string filename="", bool aligned=true, int version_number=geometry_version)
Definition: PndLmdDim.cxx:1515
static PndLmdDim & Get_instance()
Definition: PndLmdDim.cxx:242
TString MCFile