FairRoot/PandaRoot
AnalyseThetaRadiusCorrelation.C
Go to the documentation of this file.
1 // root macro to analyze the simulation output
2 //void convertMCPoints()
3 
4 #include "TROOT.h"
5 #include "TSystem.h"
6 #include "TFile.h"
7 #include "TStopwatch.h"
8 #include "TString.h"
9 #include "TH1D.h"
10 #include "TH2D.h"
11 #include "TF1.h"
12 #include "TClonesArray.h"
13 #include "TMath.h"
14 #include "TTree.h"
15 #include <fstream>
16 #include <string>
17 #include "TVector3.h"
18 
19 int AnalyseThetaRadiusCorrelation()
20 {
21 
22  // ----- Load libraries ------------------------------------------------
23 //``gSystem->Load("fstream.h");
24  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
25 
26  gSystem->Load("libHyp");
27  gSystem->Load("libHypGe");
28 
29  // ----- Timer --------------------------------------------------------
30  TStopwatch timer;
31  timer.Start();
32  // ------------------------------------------------------------------------
33  TString CompleteFilename = "/data/work/kpha1/steinen/CrystalsOnly/TripleBall40Offset20STTCrystalsOnly_0.01MeV_10000000Evts.root";
34  TFile* g = new TFile(CompleteFilename);
35 
36 //Output Files
37  TString Path = getenv("SIMDATADIR");
38  TString outfile= Path+"/CrystalsOnly/Ana/AnaTripleBall40Offset20STTCrystalsOnly_0.01MeV_10000000Evts";
39 
40  TString txtfileName = outfile;
41  outfile +=".root";
42  txtfileName += ".txt";
43  TFile* fi = new TFile(outfile,"RECREATE");
44 
45  ofstream txtfile;
46  txtfile.open(txtfileName);
47  txtfile << "File read:" << CompleteFilename << endl;
48  //photons from hyp electromag. decay
49  TTree *b=(TTree *) g->Get("pndsim") ;
50  TClonesArray* hit_bar=new TClonesArray("PndHypGePoint");
51  b->SetBranchAddress("HypGePoint",&hit_bar);//Branch names
52  TClonesArray* mc_bar=new TClonesArray("PndMCTrack");
53  b->SetBranchAddress("MCTrack",&mc_bar);//Branch names
54 
55 
56  TString Name = "#Theta detector distance correlation;#Theta [#circ];Distance [cm]";
57  TH2D* hThetaR = new TH2D("hThetaR",Name.Data(),180,90,180,250,15,40);
58 
59  TH2D* hRing1 = new TH2D("hRing1","#Theta detector ring 1 distance correlation;#Theta [#circ];Distance [cm]",180,90,180,250,15,40);
60  TH2D* hRing2 = new TH2D("hRing2","#Theta detector ring 2 distance correlation;#Theta [#circ];Distance [cm]",180,90,180,250,15,40);
61  TH2D* hRing3 = new TH2D("hRing3","#Theta detector ring 3 distance correlation;#Theta [#circ];Distance [cm]",180,90,180,250,15,40);
62  TH2D* hRing4 = new TH2D("hRing4","#Theta detector ring 4 distance correlation;#Theta [#circ];Distance [cm]",180,90,180,250,15,40);
63  TH2D* hRing5 = new TH2D("hRing5","#Theta detector ring 5 distance correlation;#Theta [#circ];Distance [cm]",180,90,180,250,15,40);
64  TH2D* hRing6 = new TH2D("hRing6","#Theta detector ring 6 distance correlation;#Theta [#circ];Distance [cm]",180,90,180,250,15,40);
65  TH2D* hRing7 = new TH2D("hRing7","#Theta detector ring 7 distance correlation;#Theta [#circ];Distance [cm]",180,90,180,250,15,40);
66  TH2D* hRing8 = new TH2D("hRing8","#Theta detector ring 8 distance correlation;#Theta [#circ];Distance [cm]",180,90,180,250,15,40);
67 
68  TH1D* hRing1Theta = new TH1D("hRing1Theta","#Theta distribution of ring 1 ;#Theta [#circ];Counts",180,90,180);
69  TH1D* hRing2Theta = new TH1D("hRing2Theta","#Theta distribution of ring 2 ;#Theta [#circ];Counts",180,90,180);
70  TH1D* hRing3Theta = new TH1D("hRing3Theta","#Theta distribution of ring 3 ;#Theta [#circ];Counts",180,90,180);
71  TH1D* hRing4Theta = new TH1D("hRing4Theta","#Theta distribution of ring 4 ;#Theta [#circ];Counts",180,90,180);
72  TH1D* hRing5Theta = new TH1D("hRing5Theta","#Theta distribution of ring 5 ;#Theta [#circ];Counts",180,90,180);
73  TH1D* hRing6Theta = new TH1D("hRing6Theta","#Theta distribution of ring 6 ;#Theta [#circ];Counts",180,90,180);
74  TH1D* hRing7Theta = new TH1D("hRing7Theta","#Theta distribution of ring 7 ;#Theta [#circ];Counts",180,90,180);
75  TH1D* hRing8Theta = new TH1D("hRing8Theta","#Theta distribution of ring 8 ;#Theta [#circ];Counts",180,90,180);
76 
77  TH1D* hRing1Distance = new TH1D("hRing1Distance","Crystal distance of ring 1 ;#Theta [#circ];Counts",250,15,40);
78  TH1D* hRing2Distance = new TH1D("hRing2Distance","Crystal distance of ring 2 ;#Theta [#circ];Counts",250,15,40);
79  TH1D* hRing3Distance = new TH1D("hRing3Distance","Crystal distance of ring 3 ;#Theta [#circ];Counts",250,15,40);
80  TH1D* hRing4Distance = new TH1D("hRing4Distance","Crystal distance of ring 4 ;#Theta [#circ];Counts",250,15,40);
81  TH1D* hRing5Distance = new TH1D("hRing5Distance","Crystal distance of ring 5 ;#Theta [#circ];Counts",250,15,40);
82  TH1D* hRing6Distance = new TH1D("hRing6Distance","Crystal distance of ring 6 ;#Theta [#circ];Counts",250,15,40);
83  TH1D* hRing7Distance = new TH1D("hRing7Distance","Crystal distance of ring 7 ;#Theta [#circ];Counts",250,15,40);
84  TH1D* hRing8Distance = new TH1D("hRing8Distance","Crystal distance of ring 8 ;#Theta [#circ];Counts",250,15,40);
85 
86  bool verbose = false;
87  Int_t MotherId,Motherpdg;
88 
89 
90  TVector3 vecs , pos;
91  int mcpdg = -1,ev;
92  Double_t mult,En,Eng,Enth;
93 
94  //vector<int> event;
95  int count;
96  set<int> SetOfCrystalHit;
97  set<int>::iterator it;
98 
99  Int_t nEvents = b->GetEntriesFast()/100;
100  cout<< "Number of Simulated Events: "<<nEvents<<endl;
101  txtfile<< "Number of Simulated Events: "<<nEvents<<endl;
102  TVector3 Hit;
103 
104  for (Int_t k=0; k<nEvents; k++)
105  {
106  b->GetEntry(k);
107  if (!((k*1000)% nEvents))
108  {
109  cout << k << endl;
110  }
111 
112  for (Int_t i=0; i<hit_bar->GetEntriesFast(); i++)
113  {
114  PndHypGePoint *hitgam=(PndHypGePoint*)hit_bar->At(i);
115  PndMCTrack *mcgam = (PndMCTrack*)mc_bar->At(hitgam->GetTrackID());
116 
117  Hit.SetX(hitgam->GetX());
118  Hit.SetY(hitgam->GetY());
119  Hit.SetZ(hitgam->GetZ()+55);
120  Double_t theta = 180/TMath::Pi()*Hit.Theta();
121  //cout << 180/TMath::Pi()*Hit.Theta()<< ", " << Hit.Mag() << endl;
122  hThetaR->Fill(theta, Hit.Mag());
123 
124  switch (hitgam->GetDetectorID())
125  {
126  case 101:
127  //iRing6+=Content;
128  hRing6->Fill(theta, Hit.Mag()); hRing6Theta->Fill(theta); hRing6Distance->Fill(Hit.Mag());
129  break;
130  case 102:
131  //iRing7+=Content;
132  hRing7->Fill(theta, Hit.Mag()); hRing7Theta->Fill(theta); hRing7Distance->Fill(Hit.Mag());
133  break;
134  case 103:
135  //iRing6+=Content;
136  hRing6->Fill(theta, Hit.Mag()); hRing6Theta->Fill(theta); hRing6Distance->Fill(Hit.Mag());
137  break;
138  case 204:
139  //iRing7+=Content;
140  hRing7->Fill(theta, Hit.Mag()); hRing7Theta->Fill(theta); hRing7Distance->Fill(Hit.Mag());
141  break;
142  case 205:
143  //iRing6+=Content;
144  hRing6->Fill(theta, Hit.Mag()); hRing6Theta->Fill(theta); hRing6Distance->Fill(Hit.Mag());
145  break;
146  case 206:
147  //iRing6+=Content;
148  hRing6->Fill(theta, Hit.Mag()); hRing6Theta->Fill(theta); hRing6Distance->Fill(Hit.Mag());
149  break;
150  case 307:
151  //iRing6+=Content;
152  hRing6->Fill(theta, Hit.Mag()); hRing6Theta->Fill(theta); hRing6Distance->Fill(Hit.Mag());
153  break;
154  case 308:
155  //iRing6+=Content;
156  hRing6->Fill(theta, Hit.Mag()); hRing6Theta->Fill(theta); hRing6Distance->Fill(Hit.Mag());
157  break;
158  case 309:
159  //iRing7+=Content;
160  hRing7->Fill(theta, Hit.Mag()); hRing7Theta->Fill(theta); hRing7Distance->Fill(Hit.Mag());
161  break;
162  case 410:
163  //iRing4+=Content;
164  hRing4->Fill(theta, Hit.Mag()); hRing4Theta->Fill(theta); hRing4Distance->Fill(Hit.Mag());
165  break;
166  case 411:
167  //iRing1+=Content;
168  hRing1->Fill(theta, Hit.Mag()); hRing1Theta->Fill(theta); hRing1Distance->Fill(Hit.Mag());
169  break;
170  case 412:
171  //iRing5+=Content;
172  hRing5->Fill(theta, Hit.Mag()); hRing5Theta->Fill(theta); hRing5Distance->Fill(Hit.Mag());
173  break;
174  case 513:
175  //iRing1+=Content;
176  hRing1->Fill(theta, Hit.Mag()); hRing1Theta->Fill(theta); hRing1Distance->Fill(Hit.Mag());
177  break;
178  case 514:
179  //iRing5+=Content;
180  hRing5->Fill(theta, Hit.Mag()); hRing5Theta->Fill(theta); hRing5Distance->Fill(Hit.Mag());
181  break;
182  case 515:
183  //iRing3+=Content;
184  hRing3->Fill(theta, Hit.Mag()); hRing3Theta->Fill(theta); hRing3Distance->Fill(Hit.Mag());
185  break;
186  case 616:
187  //iRing5+=Content;
188  hRing5->Fill(theta, Hit.Mag()); hRing5Theta->Fill(theta); hRing5Distance->Fill(Hit.Mag());
189  break;
190  case 617:
191  //iRing2+=Content;
192  hRing2->Fill(theta, Hit.Mag()); hRing2Theta->Fill(theta); hRing2Distance->Fill(Hit.Mag());
193  break;
194  case 618:
195  //iRing2+=Content;
196  hRing2->Fill(theta, Hit.Mag()); hRing2Theta->Fill(theta); hRing2Distance->Fill(Hit.Mag());
197  break;
198  case 719:
199  //iRing1+=Content;
200  hRing1->Fill(theta, Hit.Mag()); hRing1Theta->Fill(theta); hRing1Distance->Fill(Hit.Mag());
201  break;
202  case 720:
203  //iRing3+=Content;
204  hRing3->Fill(theta, Hit.Mag()); hRing3Theta->Fill(theta); hRing3Distance->Fill(Hit.Mag());
205  break;
206  case 721:
207  //iRing5+=Content;
208  hRing5->Fill(theta, Hit.Mag()); hRing5Theta->Fill(theta); hRing5Distance->Fill(Hit.Mag());
209  break;
210  case 822:
211  //iRing4+=Content;
212  hRing4->Fill(theta, Hit.Mag()); hRing4Theta->Fill(theta); hRing4Distance->Fill(Hit.Mag());
213  break;
214  case 823:
215  //iRing5+=Content;
216  hRing5->Fill(theta, Hit.Mag()); hRing5Theta->Fill(theta); hRing5Distance->Fill(Hit.Mag());
217  break;
218  case 824:
219  //iRing1+=Content;
220  hRing1->Fill(theta, Hit.Mag()); hRing1Theta->Fill(theta); hRing1Distance->Fill(Hit.Mag());
221  break;
222  case 925:
223  //iRing6+=Content;
224  hRing6->Fill(theta, Hit.Mag()); hRing6Theta->Fill(theta); hRing6Distance->Fill(Hit.Mag());
225  break;
226  case 926:
227  //iRing6+=Content;
228  hRing6->Fill(theta, Hit.Mag()); hRing6Theta->Fill(theta); hRing6Distance->Fill(Hit.Mag());
229  break;
230  case 927:
231  //iRing7+=Content;
232  hRing7->Fill(theta, Hit.Mag()); hRing7Theta->Fill(theta); hRing7Distance->Fill(Hit.Mag());
233  break;
234  case 1028:
235  //iRing7+=Content;
236  hRing7->Fill(theta, Hit.Mag()); hRing7Theta->Fill(theta); hRing7Distance->Fill(Hit.Mag());
237  break;
238  case 1029:
239  //iRing6+=Content;
240  hRing6->Fill(theta, Hit.Mag()); hRing6Theta->Fill(theta); hRing6Distance->Fill(Hit.Mag());
241  break;
242  case 1030:
243  //iRing6+=Content;
244  hRing6->Fill(theta, Hit.Mag()); hRing6Theta->Fill(theta); hRing6Distance->Fill(Hit.Mag());
245  break;
246  case 1131:
247  //iRing6+=Content;
248  hRing6->Fill(theta, Hit.Mag()); hRing6Theta->Fill(theta); hRing6Distance->Fill(Hit.Mag());
249  break;
250  case 1132:
251  //iRing7+=Content;
252  hRing7->Fill(theta, Hit.Mag()); hRing7Theta->Fill(theta); hRing7Distance->Fill(Hit.Mag());
253  break;
254  case 1133:
255  //iRing6+=Content;
256  hRing6->Fill(theta, Hit.Mag()); hRing6Theta->Fill(theta); hRing6Distance->Fill(Hit.Mag());
257  break;
258  case 1234:
259  //iRing4+=Content;
260  hRing4->Fill(theta, Hit.Mag()); hRing4Theta->Fill(theta); hRing4Distance->Fill(Hit.Mag());
261  break;
262  case 1235:
263  //iRing5+=Content;
264  hRing5->Fill(theta, Hit.Mag()); hRing5Theta->Fill(theta); hRing5Distance->Fill(Hit.Mag());
265  break;
266  case 1236:
267  //iRing1+=Content;
268  hRing1->Fill(theta, Hit.Mag()); hRing1Theta->Fill(theta); hRing1Distance->Fill(Hit.Mag());
269  break;
270  case 1337:
271  //iRing1+=Content;
272  hRing1->Fill(theta, Hit.Mag()); hRing1Theta->Fill(theta); hRing1Distance->Fill(Hit.Mag());
273  break;
274  case 1338:
275  //iRing3+=Content;
276  hRing3->Fill(theta, Hit.Mag()); hRing3Theta->Fill(theta); hRing3Distance->Fill(Hit.Mag());
277  break;
278  case 1339:
279  //iRing5+=Content;
280  hRing5->Fill(theta, Hit.Mag()); hRing5Theta->Fill(theta); hRing5Distance->Fill(Hit.Mag());
281  break;
282  case 1440:
283  //iRing5+=Content;
284  hRing5->Fill(theta, Hit.Mag()); hRing5Theta->Fill(theta); hRing5Distance->Fill(Hit.Mag());
285  break;
286  case 1441:
287  //iRing2+=Content;
288  hRing2->Fill(theta, Hit.Mag()); hRing2Theta->Fill(theta); hRing2Distance->Fill(Hit.Mag());
289  break;
290  case 1442:
291  //iRing2+=Content;
292  hRing2->Fill(theta, Hit.Mag()); hRing2Theta->Fill(theta); hRing2Distance->Fill(Hit.Mag());
293  break;
294  case 1543:
295  //iRing1+=Content;
296  hRing1->Fill(theta, Hit.Mag()); hRing1Theta->Fill(theta); hRing1Distance->Fill(Hit.Mag());
297  break;
298  case 1544:
299  //iRing5+=Content;
300  hRing5->Fill(theta, Hit.Mag()); hRing5Theta->Fill(theta); hRing5Distance->Fill(Hit.Mag());
301  break;
302  case 1545:
303  //iRing3+=Content;
304  hRing3->Fill(theta, Hit.Mag()); hRing3Theta->Fill(theta); hRing3Distance->Fill(Hit.Mag());
305  break;
306  case 1646:
307  //iRing4+=Content;
308  hRing4->Fill(theta, Hit.Mag()); hRing4Theta->Fill(theta); hRing4Distance->Fill(Hit.Mag());
309  break;
310  case 1647:
311  //iRing1+=Content;
312  hRing1->Fill(theta, Hit.Mag()); hRing1Theta->Fill(theta); hRing1Distance->Fill(Hit.Mag());
313  break;
314  case 1648:
315  //iRing5+=Content;
316  hRing5->Fill(theta, Hit.Mag()); hRing5Theta->Fill(theta); hRing5Distance->Fill(Hit.Mag());
317  break;
318  default:
319  break;
320  }
321 
322  }//end for i (points in event)
323 
324  }// end for j (events)
325  cout << "event loop finished"<< endl;
326 
327 
328  //TCanvas* can3 = new TCanvas("can3","germanium detector",0,0,1000,1000);
329  //can3->cd();
330 
331  hThetaR->Draw("");
332 
333  hThetaR->Write();
334 
335  hRing1->Write();
336  hRing2->Write();
337  hRing3->Write();
338  hRing4->Write();
339  hRing5->Write();
340  hRing6->Write();
341  hRing7->Write();
342  hRing8->Write();
343 
344  hRing1Theta->Write();
345  hRing2Theta->Write();
346  hRing3Theta->Write();
347  hRing4Theta->Write();
348  hRing5Theta->Write();
349  hRing6Theta->Write();
350  hRing7Theta->Write();
351  hRing8Theta->Write();
352 
353  hRing1Distance->Write();
354  hRing2Distance->Write();
355  hRing3Distance->Write();
356  hRing4Distance->Write();
357  hRing5Distance->Write();
358  hRing6Distance->Write();
359  hRing7Distance->Write();
360  hRing8Distance->Write();
361 
362  fi->Close();
363  txtfile.close();
364 
365 
366 
367  // ----- Finish -------------------------------------------------------
368  timer.Stop();
369  Double_t rtime = timer.RealTime();
370  Double_t ctime = timer.CpuTime();
371  cout << endl << endl;
372  cout << "Macro finished succesfully." << endl;
373  cout << "Output file is " << outfile << endl;
374  cout << "Parameter file is " << txtfileName << endl;
375  cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
376  cout << endl;
377  // ------------------------------------------------------------------------
378 
379  return 0;
380 }
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
theta
double theta
Definition: macro/detectors/emc/dedicated/reco_analys.C:24
Motherpdg
Int_t Motherpdg
Definition: GammaSpectraAnalysis.C:39
i
Int_t i
Definition: run_full.C:25
b
TTree * b
Definition: GammaSpectraAnalysis.C:24
verbose
#define verbose
Definition: PndRadMapBoxMesh.cxx:15
ev
int ev
Definition: GammaSpectraAnalysis.C:44
g
TFile * g
Definition: GammaSpectraAnalysis.C:22
PndHypGePoint::GetTrackID
Int_t GetTrackID() const
Definition: PndHypGePoint.h:51
Enth
Double_t Enth
Definition: GammaSpectraAnalysis.C:45
hit_bar
TClonesArray * hit_bar
Definition: GammaSpectraAnalysis.C:25
vecs
TVector3 vecs
Definition: GammaSpectraAnalysis.C:43
En
Double_t En
Definition: GammaSpectraAnalysis.C:45
fi
TFile * fi
Definition: createRootGeoFile.C:238
CompleteFilename
TString CompleteFilename(TString mu, TString Q, TString Psf)
Definition: OmegaQuadrupoleFourEnergiesCombination.C:1
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
nEvents
Int_t nEvents
Definition: hit_dirc.C:11
timer
TStopwatch timer
Definition: hit_dirc.C:51
Eng
Double_t Eng
Definition: GammaSpectraAnalysis.C:45
mcpdg
int mcpdg
Definition: GammaSpectraAnalysis.C:44
ctime
Double_t ctime
Definition: hit_dirc.C:114
PndHypGePoint::GetZ
Double_t GetZ() const
Definition: PndHypGePoint.h:56
PndHypGePoint::GetX
Double_t GetX() const
Definition: PndHypGePoint.h:54
PndHypGePoint::GetY
Double_t GetY() const
Definition: PndHypGePoint.h:55
mc_bar
TClonesArray * mc_bar
Definition: GammaSpectraAnalysis.C:27
count
int count
Definition: GammaSpectraAnalysis.C:48
rtime
Double_t rtime
Definition: hit_dirc.C:113
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
PndHypGePoint::GetDetectorID
Int_t GetDetectorID() const
Definition: PndHypGePoint.h:53
AnalyseThetaRadiusCorrelation
int AnalyseThetaRadiusCorrelation()
Definition: AnalyseThetaRadiusCorrelation.C:19
MotherId
Int_t MotherId
Definition: GammaSpectraAnalysis.C:39
mult
Double_t mult
Definition: GammaSpectraAnalysis.C:45
outfile
TString outfile
Definition: createRootFscGeometryFile.C:57
Generated on Sun Aug 18 2019 03:01:06 for FairRoot/PandaRoot by   doxygen 1.8.8