FairRoot/PandaRoot
Functions
AllNeutronAnalysis_job.C File Reference
#include "TH1D.h"
#include "TString.h"
#include "TFile.h"
#include "TTree.h"
#include "TStopwatch.h"
#include "TVector3.h"
#include "TROOT.h"
#include "TMath.h"
#include "TClonesArray.h"
#include <stdio.h>
#include <iostream>

Go to the source code of this file.

Functions

int AllNeutronAnalysis_job (TString Filename_ext)
 

Function Documentation

int AllNeutronAnalysis_job ( TString  Filename_ext)

Definition at line 15 of file AllNeutronAnalysis_job.C.

References b, ctime, Double_t, PndHypGePoint::GetDetectorID(), PndMCTrack::GetMomentum(), PndMCTrack::GetMotherID(), PndHypGePoint::GetpdgCode(), PndMCTrack::GetPdgCode(), PndMCTrack::GetStartVertex(), PndHypGePoint::GetTrackID(), PndHypGePoint::GetX(), PndHypGePoint::GetY(), PndHypGePoint::GetZ(), hit_bar, i, m, mc_bar, nEvents, outfile, Pi, rtime, theta, timer, TString, and verbose.

16 {
17  bool verbose = false;
18  // ----- Load libraries ------------------------------------------------
19  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
20  gSystem->Load("libHypGe");
21  // ----- Timer --------------------------------------------------------
22  TStopwatch timer;
23  timer.Start();
24  // ------------------------------------------------------------------------
25 
26 
27  //Files
28  //TString Path = getenv("SIMDATADIR");
29  TString Filename = "$SIMDATADIR/Neutron/"+Filename_ext;
30  //TString Filename = "$SIMDATADIR/TripleBall40Offset10_urqmd_pbarC_1_5000Evts.root";
31  TString outfile= "$SIMDATADIR/Neutron/Ana/Ana" + Filename_ext;
32  TFile* InputFile = new TFile(Filename);
33  TFile* OutputFile = new TFile(outfile,"RECREATE");
34 
35  //getting simulation branches from input file
36  TTree *b=(TTree *) InputFile->Get("pndsim") ;
37  TClonesArray* hit_bar=new TClonesArray("PndHypGePoint");
38  b->SetBranchAddress("HypGePoint",&hit_bar);//Branch names
39  TClonesArray* mc_bar=new TClonesArray("PndMCTrack");
40  b->SetBranchAddress("MCTrack",&mc_bar);//Branch names
41 
42  //declaring histograms
43  //histogram with polar angle of detector interactions of primary neutrons
44  TH1D* hNHits = new TH1D("hNHits","Polar angle of primary neutrons interacting with the crystals",180, 90,180);
45  hNHits->SetXTitle("#Theta [#circ]");
46  hNHits->SetYTitle("Counts / 0.5 #circ");
47 
48  //histogram with polar angle of crystal ring 1 interactions of primary neutrons
49  TH1D* hRing1 = new TH1D("hRing1","Polar angle of primary neutrons interacting with the crystals of ring 1",180, 90,180);
50  hRing1->SetXTitle("#Theta [#circ]");
51  hRing1->SetYTitle("Counts / 0.5 #circ");
52 
53  //histogram with polar angle of crystal ring 2 interactions of primary neutrons
54  TH1D* hRing2 = new TH1D("hRing2","Polar angle of primary neutrons interacting with the crystals of ring 2",180, 90,180);
55  hRing2->SetXTitle("#Theta [#circ]");
56  hRing2->SetYTitle("Counts / 0.5 #circ");
57 
58  //histogram with polar angle of crystal ring 3 interactions of primary neutrons
59  TH1D* hRing3 = new TH1D("hRing3","Polar angle of primary neutrons interacting with the crystals of ring 3",180, 90,180);
60  hRing3->SetXTitle("#Theta [#circ]");
61  hRing3->SetYTitle("Counts / 0.5 #circ");
62 
63  //histogram with polar angle of crystal ring 4 interactions of primary neutrons
64  TH1D* hRing4 = new TH1D("hRing4","Polar angle of primary neutrons interacting with the crystals of ring 4",180, 90,180);
65  hRing4->SetXTitle("#Theta [#circ]");
66  hRing4->SetYTitle("Counts / 0.5 #circ");
67 
68  //histogram with polar angle of crystal ring 5 interactions of primary neutrons
69  TH1D* hRing5 = new TH1D("hRing5","Polar angle of primary neutrons interacting with the crystals of ring 5",180, 90,180);
70  hRing5->SetXTitle("#Theta [#circ]");
71  hRing5->SetYTitle("Counts / 0.5 #circ");
72 
73  //histogram with polar angle of crystal ring 6 interactions of primary neutrons
74  TH1D* hRing6 = new TH1D("hRing6","Polar angle of primary neutrons interacting with the crystals of ring 6",180, 90,180);
75  hRing6->SetXTitle("#Theta [#circ]");
76  hRing6->SetYTitle("Counts / 0.5 #circ");
77 
78  //histogram with polar angle of crystal ring 7 interactions of primary neutrons
79  TH1D* hRing7 = new TH1D("hRing7","Polar angle of primary neutrons interacting with the crystals of ring 7",180, 90,180);
80  hRing7->SetXTitle("#Theta [#circ]");
81  hRing7->SetYTitle("Counts / 0.5 #circ");
82 
83  //histogram with polar angle of crystal ring 8 interactions of primary neutrons
84  TH1D* hRing8 = new TH1D("hRing8","Polar angle of primary neutrons interacting with the crystals of ring 8",180, 90,180);
85  hRing8->SetXTitle("#Theta [#circ]");
86  hRing8->SetYTitle("Counts / 0.5 #circ");
87 
88  //histogram with the polar momentum of primary neutrons interacting with the crystals
89  TH1D* hNHits_p = new TH1D("hNHits_p","Direction of polar momentum of primary neutrons interacting with the crystals",180, 90,180);
90  hNHits_p->SetXTitle("#Theta [#circ]");
91  hNHits_p->SetYTitle("Counts / 0.5 #circ");
92  //histogram with the polar momentum of primary neutrons in the solid angle of the crystal (not necessarily interacting with the crystal)
93  TH1D* hNAll = new TH1D("hNAll","All Neutrons",360, 0,180);
94  hNAll->SetXTitle("#Theta [#circ]");
95  hNAll->SetYTitle("Counts / 0.5 #circ");
96 
97  //histogram with the polar momentum of primary neutrons in the solid angle of the crystal (not necessarily interacting with the crystal)
98  TH1D* hPAll = new TH1D("hPAll","All Protons",360, 0,180);
99  hPAll->SetXTitle("#Theta [#circ]");
100  hPAll->SetYTitle("Counts / 0.5 #circ");
101 
102  //histogram with E_kin of primary neutrons
103  TH1D* hNMom = new TH1D("hNMom","E_{kin} of neutrons",1000, 0,0.5);
104  hNMom->SetXTitle("E_{kin} of neutrons [GeV]");
105  hNMom->SetYTitle("Counts / 0.5 MeV");
106  //histogram to see which detector is hit
107  TH1D* hCrystalHit = new TH1D("hCrystalHit","Hits per Crystal",2100,1,2100);
108  hCrystalHit->SetXTitle("Crystal number");
109  hCrystalHit->SetYTitle("Counts per crystal");
110  //histogram to see the energy deposited by the neutron
111  TH1D* hNeutronEnergyDeposit = new TH1D("hNeutronEnergyDeposit","Energy deposited by Neutrons per Hit",500,0,0.005);
112  hNeutronEnergyDeposit->SetXTitle("Energy [GeV]");
113  hNeutronEnergyDeposit->SetYTitle("Counts / 10 keV");
114 
115  //histogram to visualize the x-y-distribution of neutrons
116  TH2D* hNeutronXYDistribution = new TH2D("hNeutronXYDistribution","Distribution of Neutrons",40, -40,40,40,-40,40);
117  hNeutronXYDistribution->SetXTitle("X-postion [cm]");
118  hNeutronXYDistribution->SetYTitle("Y-postion [cm]");
119  gStyle->SetPalette(1);
120 
121  Int_t *ActualTrackID ;
122  Int_t nEvents = b->GetEntriesFast();
123  cout<< "Number of Simulated Events: "<<nEvents<<endl;
124 
125  Int_t Hits=0;
126  Int_t DetID=-10;
127  Int_t TrackID=-10;
128  int iNeutrons = 0;
129  //event loop
130 
131  TVector3 NHit;
132  TVector3 NHit_p;
133  TVector3 NMom;
134  TVector3 NAllMom;
135  TVector3 StartVertex;
136  for (Int_t k=0; k<nEvents; k++)
137  {
138  b->GetEntry(k);
139  if (!((k*100)% nEvents)) // mark every % of proceeded events
140  {
141  cout << "Event number :\t" << k << endl;
142  }
143  //if(verbose) cout<<"Event No "<<j<<endl;
144 
145  DetID=-10;
146  TrackID=-10;
147  // loop over the hits of an event
148  for (Int_t i=0; i<hit_bar->GetEntriesFast(); i++)
149  {
150  //if (i==0)
151  //Hits++;
152  //cout <<"Hit "<< i << endl;
153 
154  PndHypGePoint *hitgam=(PndHypGePoint*)hit_bar->At(i);
155  if (!hitgam)
156  {
157  continue;
158  }
159  PndMCTrack *mcgam = (PndMCTrack*)mc_bar->At(hitgam->GetTrackID());
160  if (!mcgam)
161  {
162  continue;
163  }
164  if(hitgam->GetpdgCode()==2112)
165  {
166  //NMom = mcgam->GetMomentum();
167  StartVertex = mcgam->GetStartVertex();
168  cout << "Sv" << StartVertex.Mag() << endl;
169  if (StartVertex.Mag() <60)
170  {
171  if ( hitgam->GetpdgCode()==2112 && mcgam->GetMotherID() == -1 && !(DetID == hitgam->GetDetectorID() && TrackID==hitgam->GetTrackID())) // primary neutron, last part: not(same particle and crystal) -> against multi counting of a single neutron
172  {
173  DetID = hitgam->GetDetectorID();
174  TrackID = hitgam->GetTrackID();
175  //cout << "Event :\t" << k <<"\tHit : \t" << i << "\tDetector :\t" << DetID << "\tTrackID :\t"<< TrackID<<"\tMother :\t"<< mcgam->GetMotherID()<<endl;
176  iNeutrons++;
177  hCrystalHit->Fill(hitgam->GetDetectorID()); //fill histogram to see which detector is hit
178 
179 
180  // fill vector with hit coordinates to get theta of the hit
181  NHit.SetX(hitgam->GetX());
182  NHit.SetY(hitgam->GetY());
183  NHit.SetZ(hitgam->GetZ()+55);
184  Double_t theta = 180/TMath::Pi()*NHit.Theta();
185  hNHits->Fill(theta); //fill histogram with polar angle of detector interactions of primary neutrons
186  cout << DetID << endl;
187  switch (DetID)
188  {
189  case 101:
190  //iRing3+=Content;
191  hRing3->Fill(theta);
192  break;
193  case 102:
194  //iRing4+=Content;
195  hRing4->Fill(theta);
196  break;
197  case 103:
198  //iRing3+=Content;
199  hRing3->Fill(theta);
200  break;
201  case 204:
202  //iRing4+=Content;
203  hRing4->Fill(theta);
204  break;
205  case 205:
206  //iRing3+=Content;
207  hRing3->Fill(theta);
208  break;
209  case 206:
210  //iRing3+=Content;
211  hRing3->Fill(theta);
212  break;
213  case 307:
214  //iRing3+=Content;
215  hRing3->Fill(theta);
216  break;
217  case 308:
218  //iRing3+=Content;
219  hRing3->Fill(theta);
220  break;
221  case 309:
222  //iRing4+=Content;
223  hRing4->Fill(theta);
224  break;
225  case 410:
226  //iRing2+=Content;
227  hRing2->Fill(theta);
228  break;
229  case 411:
230  //iRing1+=Content;
231  hRing1->Fill(theta);
232  break;
233  case 412:
234  //iRing2+=Content;
235  hRing2->Fill(theta);
236  break;
237  case 513:
238  //iRing1+=Content;
239  hRing1->Fill(theta);
240  break;
241  case 514:
242  //iRing2+=Content;
243  hRing2->Fill(theta);
244  break;
245  case 515:
246  //iRing2+=Content;
247  hRing2->Fill(theta);
248  break;
249  case 616:
250  //iRing2+=Content;
251  hRing2->Fill(theta);
252  break;
253  case 617:
254  //iRing1+=Content;
255  hRing1->Fill(theta);
256  break;
257  case 618:
258  //iRing1+=Content;
259  hRing1->Fill(theta);
260  break;
261  case 719:
262  //iRing1+=Content;
263  hRing1->Fill(theta);
264  break;
265  case 720:
266  //iRing2+=Content;
267  hRing2->Fill(theta);
268  break;
269  case 721:
270  //iRing2+=Content;
271  hRing2->Fill(theta);
272  break;
273  case 822:
274  //iRing2+=Content;
275  hRing2->Fill(theta);
276  break;
277  case 823:
278  //iRing2+=Content;
279  hRing2->Fill(theta);
280  break;
281  case 824:
282  //iRing1+=Content;
283  hRing1->Fill(theta);
284  break;
285  case 925:
286  //iRing3+=Content;
287  hRing3->Fill(theta);
288  break;
289  case 926:
290  //iRing3+=Content;
291  hRing3->Fill(theta);
292  break;
293  case 927:
294  //iRing4+=Content;
295  hRing4->Fill(theta);
296  break;
297  case 1028:
298  //iRing4+=Content;
299  hRing4->Fill(theta);
300  break;
301  case 1029:
302  //iRing3+=Content;
303  hRing3->Fill(theta);
304  break;
305  case 1030:
306  //iRing3+=Content;
307  hRing3->Fill(theta);
308  break;
309  case 1131:
310  //iRing3+=Content;
311  hRing3->Fill(theta);
312  break;
313  case 1132:
314  //iRing4+=Content;
315  hRing4->Fill(theta);
316  break;
317  case 1133:
318  //iRing3+=Content;
319  hRing3->Fill(theta);
320  break;
321  case 1234:
322  //iRing2+=Content;
323  hRing2->Fill(theta);
324  break;
325  case 1235:
326  //iRing2+=Content;
327  hRing2->Fill(theta);
328  break;
329  case 1236:
330  //iRing1+=Content;
331  hRing1->Fill(theta);
332  break;
333  case 1337:
334  //iRing1+=Content;
335  hRing1->Fill(theta);
336  break;
337  case 1338:
338  //iRing2+=Content;
339  hRing2->Fill(theta);
340  break;
341  case 1339:
342  //iRing2+=Content;
343  hRing2->Fill(theta);
344  break;
345  case 1440:
346  //iRing2+=Content;
347  hRing2->Fill(theta);
348  break;
349  case 1441:
350  //iRing1+=Content;
351  hRing1->Fill(theta);
352  break;
353  case 1442:
354  //iRing1+=Content;
355  hRing1->Fill(theta);
356  break;
357  case 1543:
358  //iRing1+=Content;
359  hRing1->Fill(theta);
360  break;
361  case 1544:
362  //iRing2+=Content;
363  hRing2->Fill(theta);
364  break;
365  case 1545:
366  //iRing2+=Content;
367  hRing2->Fill(theta);
368  break;
369  case 1646:
370  //iRing2+=Content;
371  hRing2->Fill(theta);
372  break;
373  case 1647:
374  //iRing1+=Content;
375  hRing1->Fill(theta);
376  break;
377  case 1648:
378  //iRing2+=Content;
379  hRing2->Fill(theta);
380  break;
381  default:
382  break;
383  }
384 
385 
386  // fill vector with hit momentum to get theta of the hit momentum
387  //NHit_p.SetX(hitgam->GetPx());
388  //NHit_p.SetY(hitgam->GetPy());
389  //NHit_p.SetZ(hitgam->GetPz());
390  //hNHits_p->Fill(180/TMath::Pi()*NHit_p.Theta()); //fill histogram with the polar momentum of primary neutrons interacting with the crystals
391 
392  hNeutronXYDistribution->Fill(hitgam->GetX(),hitgam->GetY());
393  //hNMom->Fill(1./2.* NMom.Mag2()/0.939); //fill histogram with E_kin of primary neutrons
394 
395  }
396 
397  //if (hitgam->GetpdgCode()==2112) // neutron
398  //{
399  //hNeutronEnergyDeposit->Fill(hitgam->GetEnergyLoss());
400  //}
401  }
402  }
403  }
404 
405  // loop over all entries in mctrack
406  //cout << "mcbar length :\t" <<mc_bar->GetEntriesFast() << endl;
407  for (Int_t m = 0; m<mc_bar->GetEntriesFast();m++)
408  {
409  //cout << "mcbar length :\t" <<mc_bar->GetEntriesFast() << endl;
410  PndMCTrack *mcgam = (PndMCTrack*)mc_bar->At(m);
411 
412  NAllMom = mcgam->GetMomentum();
413  if (mcgam->GetPdgCode() == 2112 && mcgam->GetMotherID() == -1)
414  {
415  //cout << "event :\t"<<k<<"\tmcpdg :\t" << mcgam->GetPdgCode() << "\tmcmotherid :\t"<< mcgam->GetMotherID() <<endl;
416 
417  hNAll->Fill(180/TMath::Pi()*NAllMom.Theta()); //fill histogram with the polar momentum of primary neutrons in the solid angle of the crystal (not necessarily interacting with the crystal)
418  }
419  if (mcgam->GetPdgCode() == 2212 && mcgam->GetMotherID() == -1)
420  {
421  //cout << "event :\t"<<k<<"\tmcpdg :\t" << mcgam->GetPdgCode() << "\tmcmotherid :\t"<< mcgam->GetMotherID() <<endl;
422 
423  hPAll->Fill(180/TMath::Pi()*NAllMom.Theta()); //fill histogram with the polar momentum of primary neutrons in the solid angle of the crystal (not necessarily interacting with the crystal)
424  }
425  }
426  }// end for k (events)
427 
428  TString PathFirstPart = getenv("SIMDATADIR"); //always needed for i/ofstream files, TFile can direcly use $SIMDATADIR
429  TString TxTOutFilename =PathFirstPart+"/Neutron/Ana/ana" + Filename_ext + ".txt";
430  ofstream TxTOutfile(TxTOutFilename.Data());
431  Int_t CrystalNumber = 1;
432  cout << "Bin\tCrystal\tCluster\tNeutron hits"<<endl;
433  TxTOutfile << "Bin\tCrystal\tCluster\tNeutron hits"<<endl;
434  for (Int_t iBin=0;iBin<2100;iBin++)
435  {
436 
437  if(hCrystalHit->GetBinContent(iBin))
438  {
439  cout << iBin << "\t"<< CrystalNumber <<"\t"<< iBin/100 <<"\t"<< hCrystalHit->GetBinContent(iBin)<< endl;
440 
441  TxTOutfile << iBin << "\t"<< CrystalNumber <<"\t"<< iBin/100 <<"\t"<< hCrystalHit->GetBinContent(iBin)<< endl;
442  CrystalNumber++;
443 
444  }
445  }
446  TxTOutfile.close();
447  hNeutronXYDistribution->DrawCopy("colz");
448  //hNAll->DrawCopy();
449  //hNHits->SetLineColor(kRed);
450  //hNHits->DrawCopy("same");
451 
452  //building the average over the crystals
453  hRing1->Scale(12);
454  hRing2->Scale(18);
455  hRing3->Scale(12);
456  hRing4->Scale(6);
457 
458 
459  hNHits->Write();
460  hRing1->Write();
461  hRing2->Write();
462  hRing3->Write();
463  hRing4->Write();
464  hNHits_p->Write();
465  hNAll->Write();
466  hPAll->Write();
467  hNMom->Write();
468  hCrystalHit->Write();
469  hNeutronEnergyDeposit->Write();
470  hNeutronXYDistribution->Write();
471 
472 
473 
474  //hNeutronEnergyDeposit->DrawCopy();
475 
476  OutputFile->Close();
477 
478  //hNeutronEnergyDeposit->Draw();
479  //hNHits->Draw();
480 
481 
482  // ----- Finish -------------------------------------------------------
483  timer.Stop();
484  Double_t rtime = timer.RealTime();
485  Double_t ctime = timer.CpuTime();
486  cout << endl << endl;
487  cout << "Macro finished succesfully." << endl;
488 
489  cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
490  cout << endl;
491 
492  cout << "Hits\t" << Hits << endl;
493  // ------------------------------------------------------------------------
494  return iNeutrons;
495 }
theta
double theta
Definition: macro/detectors/emc/dedicated/reco_analys.C:24
i
Int_t i
Definition: run_full.C:25
m
__m128 m
Definition: P4_F32vec4.h:28
b
TTree * b
Definition: GammaSpectraAnalysis.C:24
PndMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
verbose
#define verbose
Definition: PndRadMapBoxMesh.cxx:15
PndHypGePoint::GetTrackID
Int_t GetTrackID() const
Definition: PndHypGePoint.h:51
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
hit_bar
TClonesArray * hit_bar
Definition: GammaSpectraAnalysis.C:25
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
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
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
PndMCTrack::GetMotherID
Int_t GetMotherID() const
Definition: PndMCTrack.h:74
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
outfile
TString outfile
Definition: createRootFscGeometryFile.C:57
PndHypGePoint::GetpdgCode
Double_t GetpdgCode() const
Definition: PndHypGePoint.h:63
Generated on Sun Aug 18 2019 03:02:20 for FairRoot/PandaRoot by   doxygen 1.8.8