FairRoot/PandaRoot
Functions
rad_dose_studies.C File Reference
#include <iostream>
#include <TCanvas.h>
#include <TApplication.h>
#include <TROOT.h>
#include <TChain.h>
#include <TStyle.h>
#include <TGaxis.h>
#include <TColor.h>
#include <TClonesArray.h>
#include <TH1.h>
#include <TH2.h>
#include <PndLmdDim.h>
#include <PndMCTrack.h>
#include <PndSdsMCPoint.h>

Go to the source code of this file.

Functions

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

Function Documentation

void DrawProgressBar ( int  len,
double  percent 
)

Definition at line 203 of file rad_dose_studies.C.

References i, and last_percent().

Referenced by rad_dose_studies().

203  {
204  if ((int) (last_percent * 100) == (int) (percent * 100))
205  return;
206  //cout << " drawing " << endl;
207  cout << "\x1B[2K"; // Erase the entire current line.
208  cout << "\x1B[0E"; // Move to the beginning of the current line.
209  string progress;
210  for (int i = 0; i < len; ++i) {
211  if (i < static_cast<int> (len * percent)) {
212  progress += "=";
213  } else {
214  progress += " ";
215  }
216  }
217  cout << "[" << progress << "] " << (static_cast<int> (100 * percent))
218  << "%";
219  flush( cout); // Required.
220  last_percent = percent;
221 }
i
Int_t i
Definition: run_full.C:25
last_percent
double last_percent(-1.)
double last_percent ( -  1.)

Referenced by DrawProgressBar().

int main ( void  )

Definition at line 32 of file rad_dose_studies.C.

References rad_dose_studies(), and Root_Appearance().

32  {
33  Root_Appearance();
34  TApplication myapp("myapp", 0, 0);
35  rad_dose_studies();
36  myapp.Run();
37  return 0;
38 }
rad_dose_studies
int rad_dose_studies()
Definition: rad_dose_studies.C:40
Root_Appearance
void Root_Appearance()
Definition: rad_dose_studies.C:171
int rad_dose_studies ( )

Definition at line 40 of file rad_dose_studies.C.

References DrawProgressBar(), PndLmdDim::Get_histogram_Plane(), PndLmdDim::Get_instance(), PndLmdDim::Get_sensor_by_id(), PndSdsMCPoint::GetPosition(), PndSdsMCPoint::GetSensorID(), nEvents, nHits, PndLmdDim::Transform_global_to_lmd_local(), and tsim.

Referenced by main().

40  {
41  cout << " lmd radiation dose studies based on energy deposit in sensors " << endl;
42  TCanvas* canvas = new TCanvas("canvas", "canvas", 600, 600);
43  canvas->Divide(2,2);
44 
45  // load the MC File
46  TChain tMC("pndsim");
47  tMC.Add("./Lumi_MC_0.root");
48 
49  //--- assign MC info -----------------------------------------------------
50  TClonesArray* true_tracks = new TClonesArray("PndMCTrack");
51  tMC.SetBranchAddress("MCTrack", &true_tracks); //True Track to compare
52 
53  TClonesArray* true_points = new TClonesArray("PndSdsMCPoint");
54  tMC.SetBranchAddress("LMDPoint", &true_points); //True Points to compare
55 
56  // get the lmd dim help class
57  PndLmdDim& lmddim = PndLmdDim::Get_instance();
58 
59  // setup histograms
60  TH1F* hist_eloss_total = new TH1F("hist_eloss_total", "energy deposit in 50#mum x 2cm x2cm sensors", 200, 0, 60);
61  TH1F* hist_eloss = new TH1F("hist_eloss", "energy deposit scaled to 1#mum", 200, 0, 1200);
62  hist_eloss_total->SetXTitle("E_{dep} [keV]");
63  hist_eloss->SetXTitle("E_{dep} [eV]");
64 
65  TH2* plane3_hits = new TH2F("plane3_hits", "hits on plane 3", 200, -10, 10, 200, -10, 10); //lmddim.Get_histogram_Plane(0, 0);
66  plane3_hits->SetXTitle("X [cm]");
67  plane3_hits->SetYTitle("Y [cm]");
68  plane3_hits->SetZTitle("# hits");
69  TH2Poly* plane3_Edep = lmddim.Get_histogram_Plane(3, 0);
70  plane3_Edep->SetNameTitle("plane3_Edep", "E dep on plane 3");
71  plane3_Edep->SetZTitle("E_{dep} [J]");
72  TH2Poly* plane3_rate = lmddim.Get_histogram_Plane(3, 0);
73  plane3_rate->SetNameTitle("plane3_rate", "rate on plane 3");
74  plane3_rate->SetZTitle("rate [kHz]");
75  TH2Poly* plane3_rad = lmddim.Get_histogram_Plane(3, 0);
76  plane3_rad->SetNameTitle("plane3_rad", "IEL dose on plane 3");
77  plane3_rad->SetZTitle("IEL dose [Gy/a]");
78 
79  // read events and analyze
80  int nEvents = tMC.GetEntries();
81  int nEvents_counted(0);
82  int nhitstotal(0);
83 
84  cout << " reading " << nEvents << " Events " << endl;
85  for (Int_t j = 0; j < nEvents ; j++) {
86  nEvents_counted++;
87  DrawProgressBar(50, (j + 1) / ((double) nEvents));
88  tMC.GetEntry(j);
89 
90  int nHits = true_points->GetEntriesFast();
91  int nSdsHits = 0;
92  for (Int_t iHit = 0; iHit < nHits; iHit++) {
93  PndSdsMCPoint* mcpoint = (PndSdsMCPoint*) true_points->At(iHit);
94  //if (mcpoint->GetTrackID() < 0) continue;
95  // cout << mcpoint->GetEnergyLoss() << endl;
96  hist_eloss_total->Fill(mcpoint->GetEnergyLoss()*1.e6); // in keV
97  hist_eloss->Fill(mcpoint->GetEnergyLoss()*1.e9/50.); // in eV
98 
99  int sensID = mcpoint->GetSensorID();
100  int ihalf, iplane, imodule, iside, idie, isensor;
101  lmddim.Get_sensor_by_id(sensID, ihalf, iplane, imodule, iside, idie, isensor);
102  TVector3 hit_pos = mcpoint->GetPosition();
103  TVector3 hit_pos_lmd = lmddim.Transform_global_to_lmd_local(hit_pos, false);
104 
105  nhitstotal++;
106 
107  if (iplane == 3 && iside == 0){
108  plane3_hits->Fill(hit_pos_lmd.X(), hit_pos_lmd.Y());
109  plane3_rate->Fill(hit_pos_lmd.X(), hit_pos_lmd.Y());
110  plane3_Edep->Fill(hit_pos_lmd.X(), hit_pos_lmd.Y(), mcpoint->GetEnergyLoss()*1.e9*1.602e-19); // in J
111  plane3_rad->Fill(hit_pos_lmd.X(), hit_pos_lmd.Y(), mcpoint->GetEnergyLoss()*1.e9*1.602e-19); // in J
112  }
113 /*
114  if (mcpoint->GetTrackID() < 0) continue;
115  PndMCTrack *mctrk = (PndMCTrack*) true_tracks->At(
116  mcpoint->GetTrackID());
117  if (mctrk->IsGeneratorCreated()) { // take only anti-protons from the generator
118  if (mcpoint) {
119  TVector3 _mcpoint((mcpoint->GetPosition())); // position at the sensor entrance
120  TVector3 _mctrack(mcpoint->GetPx(), mcpoint->GetPy(),
121  mcpoint->GetPz()); // momentum of the track at the entrance
122 
123  TVector3 momMC = mctrk->GetMomentum(); // momentum in the primary vertex
124  TVector3 posMC = mctrk->GetStartVertex(); // position of the primary vertex
125  }
126  } */
127  }
128  }
129  // cross section for N generated events at 1.5 GeV/c momentum was
130  // sigm_coul = 36.683 mbarn
131  // sigm_had_el = 25.3021 mbarn
132  // sigm_inel = 62.71 mbarn according to Anastasias note
133  // -> sig_total = 124.63 mbarn or 124.63 x 10-27 cm^-2
134  double sig_total = 124.63e-27;//62e-27;//124.63e-27; // cm^2
135  // the mean instantanious luminosity at high luminosity mode is 1.6 x 10^32
136  // but we use here for reference 1 x 10^32 which can be scaled easily to any other number
137  double lumi_ref = 1e32; // cm^-2 s^-1
138  // therefore we simulated an interaction rate of
139  double rate_total = sig_total * lumi_ref; // 1/s or Hz
140  cout << "\n The total rate at 1.5 GeV/c and " << lumi_ref << " mean luminosity to normalize to is " << rate_total << " Hz" << endl;
141  // The simulated time scale is then
142  double tsim = nEvents_counted / rate_total; // s
143  cout << nEvents_counted << " events correspond to a real time of " << tsim << " s " << endl;
144  plane3_rate->Scale(1./tsim*1e-3); // in kHz
145  //plane3_rate->SetMaximum(150);
146  cout << " total rate of " << nhitstotal << " hits is " << nhitstotal/tsim*1e-6 << "MHz" << endl;
147  plane3_rad->Scale(1./tsim*60.*60.*24.*365.*0.5/(2.33e-3*(2*2*50.e-4))); // 50% duty cycle over a year and normalized to the mass in kg of one sensor
148  //canvas->cd(1);
149  plane3_Edep->Draw("colz");
150  //hist_eloss_total->Draw();
151  //canvas->cd(2);
152  plane3_hits->Draw("colz");
153  canvas->Print("MC_hit_distr_1_5_plane_3_side_0.pdf");
154  canvas->Print("MC_hit_distr_1_5_plane_3_side_0.root");
155  hist_eloss->Draw();
156  canvas->Print("MC_IEL_distrib_1_5.pdf");
157  canvas->Print("MC_IEL_distrib_1_5.root");
158  //canvas->cd(3);
159  plane3_rate->Draw("colz");
160  cout << " max rate is " << plane3_rate->GetBinContent(plane3_rate->GetMaximumBin()) << " kHz" << endl;
161  canvas->Print("MC_hit_rate_1_5_plane_3_side_0.pdf");
162  canvas->Print("MC_hit_rate_1_5_plane_3_side_0.root");
163  //canvas->cd(4);
164  plane3_rad->Draw("colz");
165  cout << " max dose is " << plane3_rad->GetBinContent(plane3_rad->GetMaximumBin()) << " Gy" << endl;
166  canvas->Print("MC_IEL_dose_1_5_plane_3_side_0.pdf");
167  canvas->Print("MC_IEL_does_1_5_plane_3_side_0.root");
168  return 0;
169 }
tsim
TTree * tsim
Definition: anal_hit_digi_cluster_fwendcap.C:10
PndLmdDim::Get_histogram_Plane
TH2Poly * Get_histogram_Plane(int iplane, int iside, bool aligned=true, bool lmd_frame=true, bool pixel_subdivision=false)
Definition: PndLmdDim.cxx:3658
PndLmdDim::Transform_global_to_lmd_local
TVector3 Transform_global_to_lmd_local(const TVector3 &point, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3113
PndSdsMCPoint::GetSensorID
Int_t GetSensorID() const
Definition: PndSdsMCPoint.h:89
nHits
int nHits
Definition: RiemannTest.C:16
PndLmdDim::Get_sensor_by_id
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
nEvents
Int_t nEvents
Definition: hit_dirc.C:11
PndSdsMCPoint::GetPosition
TVector3 GetPosition() const
Definition: PndSdsMCPoint.h:90
DrawProgressBar
void DrawProgressBar(int len, double percent)
Definition: rad_dose_studies.C:203
PndLmdDim::Get_instance
static PndLmdDim & Get_instance()
Definition: PndLmdDim.cxx:242
void Root_Appearance ( )

Definition at line 171 of file rad_dose_studies.C.

References Double_t.

Referenced by main().

171  {
172  TPad foo; // never remove this line :-)))
173  if(1){
174  gROOT->SetStyle("Plain");
175  const Int_t NRGBs = 5;
176  const Int_t NCont = 255;
177  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
178  Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
179  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
180  Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
181  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
182  gStyle->SetNumberContours(NCont);
183  gStyle->SetTitleFont(10*13+2,"xyz");
184  gStyle->SetTitleSize(0.06, "xyz");
185  gStyle->SetTitleOffset(1,"xy");
186  gStyle->SetTitleOffset(1.3,"z");
187  gStyle->SetLabelFont(10*13+2,"xyz");
188  gStyle->SetLabelSize(0.06,"xyz");
189  gStyle->SetLabelOffset(0.009,"xyz");
190  gStyle->SetPadBottomMargin(0.2);
191  gStyle->SetPadTopMargin(0.15);
192  gStyle->SetPadLeftMargin(0.15);
193  gStyle->SetPadRightMargin(0.20);
194  gStyle->SetOptTitle(0);
195  gStyle->SetOptStat(0);
196  gROOT->ForceStyle();
197  gStyle->SetFrameFillColor(0);
198  gStyle->SetFrameFillStyle(0);
199  TGaxis::SetMaxDigits(3);
200  }
201 }
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
Generated on Sun Aug 18 2019 03:02:47 for FairRoot/PandaRoot by   doxygen 1.8.8