FairRoot/PandaRoot
Functions
tools/RadMapTool/plot_radmap.C File Reference
#include "TROOT.h"
#include "TString.h"
#include "TFile.h"
#include "TChain.h"
#include "TClonesArray.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TH3D.h"
#include "THnSparse.h"
#include "TVector3.h"
#include "FairRadMapPoint.h"
#include <fstream>
#include <string>
#include <cmath>
#include <G4PhysicalConstants.hh>
#include <G4SystemOfUnits.hh>
#include "PndRadMapBoxMesh.h"

Go to the source code of this file.

Functions

void plot_radmap (Long64_t nevreq=1000, TString inputfile="RadMap_Out_Sim.root", TString outputfile="RadMap_Out.root")
 
int main ()
 

Function Documentation

int main ( void  )

Definition at line 32 of file tools/RadMapTool/plot_radmap.C.

References plot_radmap().

32  {
33  plot_radmap();
34 }
void plot_radmap(Long64_t nevreq=1000, TString inputfile="RadMap_Out_Sim.root", TString outputfile="RadMap_Out.root")
void plot_radmap ( Long64_t  nevreq = 1000,
TString  inputfile = "RadMap_Out_Sim.root",
TString  outputfile = "RadMap_Out.root" 
)

Definition at line 37 of file tools/RadMapTool/plot_radmap.C.

References Edep, Fluence, i, npoints, p, XY, ZY, and Zz.

Referenced by main().

40 {
41  // open an output file
42  TFile *fout = new TFile(outputfile,"RECREATE");
43 
44  // chain for the pndsim tree in the simulation files
45  TChain mych("pndsim");
46 
47  // ifstream ifile("filecollectionn.dat");
48  // string line;
49 
50  // mych.Add("blablabla2.root");
51  mych.Add(inputfile.Data());
52 
53  Long64_t nentries = (Long64_t)mych.GetEntries();
54  cout << "we have " << nentries << " events and will process "<< ((nevreq > nentries) ? nentries : nevreq) <<" of them." << endl;
55  TClonesArray *fRadMapPoint = new TClonesArray("FairRadMapPoint");
56  mych.SetBranchAddress("RadMap", &fRadMapPoint);
57 
58  std::vector<PndRadMapBoxMesh*> BM;
59 
60  BM.push_back(new PndRadMapBoxMesh("EmcEdepXY",
61  250,0.,250.,
62  250,0.,250.,
63  1, -40, 110));
64  BM.at(BM.size()-1)->SetQuantity(Edep);
65  BM.at(BM.size()-1)->SetOrientation(XY);
66 
67  BM.push_back(new PndRadMapBoxMesh("EmcAllFluenceZX",
68  1, 50, 50,
69  250,0.,250.,
70  150, -40, 110));
71  BM.at(BM.size()-1)->SetQuantity(Fluence);
72  BM.at(BM.size()-1)->SetOrientation(ZY, 14.3, Zz);
73  BM.at(BM.size()-1)->SetVerbosityLevel(0);
74 
75  BM.push_back(new PndRadMapBoxMesh("EmcCFluenceZX",
76  1, 50, 50,
77  250,0.,250.,
78  150, -40, 110));
79  BM.at(BM.size()-1)->SetQuantity(Fluence);
80  BM.at(BM.size()-1)->SetOrientation(ZY, 14.3, Zz);
81  BM.at(BM.size()-1)->SetFilter("(Ch!=0)");//charged
82 
83  BM.push_back(new PndRadMapBoxMesh("EmcCnepFluenceZX",
84  1, 50, 50,
85  250,0.,250.,
86  150, -40, 110));
87  BM.at(BM.size()-1)->SetQuantity(Fluence);
88  BM.at(BM.size()-1)->SetOrientation(ZY, 14.3, Zz);
89  BM.at(BM.size()-1)->SetFilter("(Ch!=0) && (Pid!=11) && (Pid!=-11)");//charged but e+-
90 
91  BM.push_back(new PndRadMapBoxMesh("EmcPFluenceZX",
92  1, 50, 50,
93  250,0.,250.,
94  150, -40, 110));
95  BM.at(BM.size()-1)->SetQuantity(Fluence);
96  BM.at(BM.size()-1)->SetOrientation(ZY, 14.3, Zz);
97  BM.at(BM.size()-1)->SetFilter("(Pid==2212)");//proton
98 
99  BM.push_back(new PndRadMapBoxMesh("EmcNFluenceZX",
100  1, 50, 50,
101  250,0.,250.,
102  150, -40, 110));
103  BM.at(BM.size()-1)->SetQuantity(Fluence);
104  BM.at(BM.size()-1)->SetOrientation(ZY, 14.3, Zz);
105  BM.at(BM.size()-1)->SetFilter("(Pid==2112)");//neutron
106 
107 
108  BM.push_back(new PndRadMapBoxMesh("EmcAllFluencePawelXY",
109  210, -105. , 105.,
110  210, -105. , 105,
111  1 , 207., 207.));
112  BM.at(BM.size()-1)->SetQuantity(Fluence);
113  BM.at(BM.size()-1)->SetOrientation(XY);
114 
115  BM.push_back(new PndRadMapBoxMesh("EmcCFluencePawelXY",
116  210, -105. , 105.,
117  210, -105. , 105,
118  1 , 207., 207.));
119  BM.at(BM.size()-1)->SetQuantity(Fluence);
120  BM.at(BM.size()-1)->SetOrientation(XY);
121  BM.at(BM.size()-1)->SetFilter("(Ch!=0)");//charged (including e+-)
122 
123  BM.push_back(new PndRadMapBoxMesh("EmcCnepFluencePawelXY",
124  210, -105. , 105.,
125  210, -105. , 105,
126  1 , 207., 207.));
127  BM.at(BM.size()-1)->SetQuantity(Fluence);
128  BM.at(BM.size()-1)->SetOrientation(XY);
129  BM.at(BM.size()-1)->SetFilter("(Ch!=0) && (Pid!=11) && (Pid!=-11) ");//charged particles but not e+-
130 
131  BM.push_back(new PndRadMapBoxMesh("EmcPFluencePawelXY",
132  210, -105. , 105.,
133  210, -105. , 105,
134  1 , 207., 207.));
135  BM.at(BM.size()-1)->SetQuantity(Fluence);
136  BM.at(BM.size()-1)->SetOrientation(XY);
137  BM.at(BM.size()-1)->SetFilter("(Pid==2212)");//proton
138 
139  BM.push_back(new PndRadMapBoxMesh("EmcNFluencePawelXY",
140  210, -105. , 105.,
141  210, -105. , 105,
142  1 , 207., 207.));
143  BM.at(BM.size()-1)->SetQuantity(Fluence);
144  BM.at(BM.size()-1)->SetOrientation(XY);
145  BM.at(BM.size()-1)->SetFilter("(Pid==2112)");//proton
146 
147 
149 
150  // BM.push_back(new PndRadMapBoxMesh("EmcEdepXY",
151  // 210, -105. , 105.,
152  // 210, -105. , 105,
153  // 1 , 207., 223.));
154  // BM.at(BM.size()-1)->SetQuantity(Edep);
155  // BM.at(BM.size()-1)->SetOrientation(XY);
156 
157  // BM.push_back(new PndRadMapBoxMesh("EmcTwosXY",
158  // 210, -105. , 105.,
159  // 210, -105. , 105,
160  // 1 , 222., 223.));
161  // BM.at(BM.size()-1)->SetQuantity(Twos);
162  // BM.at(BM.size()-1)->SetOrientation(XY);
163 
164  // BM.push_back(new PndRadMapBoxMesh("EmcDoseXY",
165  // 210, -105. , 105.,
166  // 210, -105. , 105,
167  // 1 , 222., 223.));
168  // BM.at(BM.size()-1)->SetQuantity(Dose);
169  // BM.at(BM.size()-1)->SetOrientation(XY);
170 
171  // BM.push_back(new PndRadMapBoxMesh("EmcTwoXY",
172  // 1000, 2 , 3,
173  // 1000, 23.5, 24.5,
174  // 1 , 222. , 223.));
175  // BM.at(BM.size()-1)->SetQuantity(Twos);
176  // BM.at(BM.size()-1)->SetOrientation(XY);
177 
178  // BM.push_back(new PndRadMapBoxMesh("EmcTwoXZ",
179  // 210, -105. , 105.,
180  // 210, -105. , 105,
181  // 1 , 222., 223.));
182  // BM.at(BM.size()-1)->SetQuantity(Twos);
183  // BM.at(BM.size()-1)->SetOrientation(XY);
184 
185  // BM.push_back(new PndRadMapBoxMesh("EmcTwoYZ",
186  // 1 , 2 , 3,
187  // 1000, 23.5, 24.5,
188  // 1000, 222. , 223.));
189  // BM.at(BM.size()-1)->SetQuantity(Twos);
190  // BM.at(BM.size()-1)->SetOrientation(YZ);
191  // BM.push_back(new PndRadMapBoxMesh("EmcTwoXY",
192  // 210, -105. , 105.,
193  // 210, -105. , 105,
194  // 1 , 222., 223.));
195  // BM.at(BM.size()-1)->SetQuantity(Twos);
196  // BM.at(BM.size()-1)->SetOrientation(XY);
197 
198  // BM.push_back(new PndRadMapBoxMesh("EmcMassXY",
199  // 210, -105. , 105.,
200  // 210, -105. , 105,
201  // 1 , -40., 110.));
202  // BM.at(BM.size()-1)->SetQuantity(Mass);
203  // BM.at(BM.size()-1)->SetOrientation(XY);
204  // BM.at(BM.size()-1)->SetVerbosityLevel(0);
205 
206  // BM.push_back(new PndRadMapBoxMesh("EmcDensityXY",
207  // 210, -105. , 105.,
208  // 210, -105. , 105,
209  // 1 , 222., 223.));
210  // BM.at(BM.size()-1)->SetQuantity(Density);
211  // BM.at(BM.size()-1)->SetOrientation(XY);
212 
213 
214  // BM.push_back(new PndRadMapBoxMesh("EdepZX_1",
215  // 250,0,250,
216  // 1,-0.5,0.5,
217  // 1200,-200,1000));
218  // BM.at(BM.size()-1)->SetQuantity(Edep);
219  // BM.at(BM.size()-1)->SetOrientation(ZX);
220 
221  // BM.push_back(new PndRadMapBoxMesh("MassZX_1",
222  // 250,0,250,
223  // 1,-0.5,0.5,
224  // 1200,-200,1000));
225  // BM.at(BM.size()-1)->SetQuantity(Mass);
226  // BM.at(BM.size()-1)->SetOrientation(ZX);
227 
228  // BM.push_back(new PndRadMapBoxMesh("DoseZX_1",
229  // 250,0,250,
230  // 1,-0.5,0.5,
231  // 1200,-200,1000));
232  // BM.at(BM.size()-1)->SetQuantity(Dose);
233  // BM.at(BM.size()-1)->SetOrientation(ZX);
234 
235  // BM.push_back(new PndRadMapBoxMesh("EdepZY_1",
236  // 1,-0.5,0.5,
237  // 250,0,250,
238  // 1200,-200,1000));
239  // BM.at(BM.size()-1)->SetQuantity(Edep);
240  // BM.at(BM.size()-1)->SetOrientation(ZY);
241 
242  // BM.push_back(new PndRadMapBoxMesh("DoseZY_1",
243  // 1,-0.5,0.5,
244  // 250,0,250,
245  // 1200,-200,1000));
246  // BM.at(BM.size()-1)->SetQuantity(Dose);
247  // BM.at(BM.size()-1)->SetOrientation(ZY);
248 
249  // BM.push_back(new PndRadMapBoxMesh("MassZY_1",
250  // 1,-0.5,0.5,
251  // 250,0,250,
252  // 1200,-200,1000));
253  // BM.at(BM.size()-1)->SetQuantity(Mass);
254  // BM.at(BM.size()-1)->SetOrientation(ZY);
255 
256 
257 
258 
259 
260  for (Long64_t i=0; i<nentries && i<nevreq; i++) {
261  cout << "event " << i << endl;
262  mych.GetEntry(i);
263  Int_t npoints = fRadMapPoint->GetEntries();
264 
265  cout << npoints << endl;
266  // loop over points
267  for (Int_t ii=0; ii < npoints; ii++) {
268 
269  //if (i==99)cout<<ii<<" "<<npoints<<endl;
270 
271  FairRadMapPoint *p= (FairRadMapPoint *) fRadMapPoint->At(ii);
272  for(unsigned int iii = 0; iii < BM.size(); iii++){
273  BM.at(iii)->Fill(p);
274  }
275  }
276  }
277  cout << "Save " << BM.size() << endl;
278  for(unsigned int i = 0; i < BM.size(); i++){
279  BM.at(i)->Scale(1./nentries);//normalizing it to 1 event!
280  BM.at(i)->Save(fout);
281  }
282 }
Int_t i
Definition: run_full.C:25
Double_t p
Definition: anasim.C:58