FairRoot/PandaRoot
Functions
AllNeutronAnalysis_job_save.C File Reference

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 1 of file AllNeutronAnalysis_job_save.C.

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

2 {
3  bool verbose = false;
4  // ----- Load libraries ------------------------------------------------
5  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
6  gSystem->Load("libHypGe");
7  // ----- Timer --------------------------------------------------------
8  TStopwatch timer;
9  timer.Start();
10  // ------------------------------------------------------------------------
11 
12 
13  //Files
14  //TString Path = getenv("SIMDATADIR");
15  TString Filename = "$SIMDATADIR/Neutron/"+Filename_ext;
16  //TString Filename = "$SIMDATADIR/TripleBall40Offset10_urqmd_pbarC_1_5000Evts.root";
17  TString outfile= "$SIMDATADIR/Neutron/Ana/Ana" + Filename_ext;
18  TFile* InputFile = new TFile(Filename);
19  TFile* OutputFile = new TFile(outfile,"RECREATE");
20 
21  //getting simulation branches from input file
22  TTree *b=(TTree *) InputFile->Get("pndsim") ;
23  TClonesArray* hit_bar=new TClonesArray("PndHypGePoint");
24  b->SetBranchAddress("HypGePoint",&hit_bar);//Branch names
25  TClonesArray* mc_bar=new TClonesArray("PndMCTrack");
26  b->SetBranchAddress("MCTrack",&mc_bar);//Branch names
27 
28  //declaring histograms
29  //histogram with polar angle of detector interactions of primary neutrons
30  TH1D* hNHits = new TH1D("hNHits","Polar angle of primary neutrons interacting with the crystals",180, 90,180);
31  hNHits->SetXTitle("#Theta [°]");
32  hNHits->SetYTitle("Counts");
33  //histogram with the polar momentum of primary neutrons interacting with the crystals
34  TH1D* hNHits_p = new TH1D("hNHits_p","Direction of polar momentum of primary neutrons interacting with the crystals",180, 90,180);
35  hNHits_p->SetXTitle("#Theta [°]");
36  hNHits_p->SetYTitle("Counts");
37  //histogram with the polar momentum of primary neutrons in the solid angle of the crystal (not necessarily interacting with the crystal)
38  TH1D* hNAll = new TH1D("hNAll","All Neutrons",180, 90,180);
39  hNAll->SetXTitle("#Theta [°]");
40  hNAll->SetYTitle("Counts");
41  //histogram with E_kin of primary neutrons
42  TH1D* hNMom = new TH1D("hNMom","E_{kin} of neutrons",1000, 0,0.5);
43  hNMom->SetXTitle("E_{kin} of neutrons [GeV]");
44  hNMom->SetYTitle("Counts");
45  //histogram to see which detector is hit
46  TH1D* hCrystalHit = new TH1D("hCrystalHit","Hits per Crystal",1700,1,1700);
47  hCrystalHit->SetXTitle("Crystal number");
48  hCrystalHit->SetYTitle("Counts");
49 
50 
51  Int_t *ActualTrackID ;
52  Int_t nEvents = b->GetEntriesFast();
53  cout<< "Number of Simulated Events: "<<nEvents<<endl;
54 
55  for (Int_t k=0; k<nEvents; k++)
56  {
57  b->GetEntry(k);
58  if (!((k*100)% nEvents))
59  {
60  cout << "Event number :\t" << k << endl;
61  }
62  //if(verbose) cout<<"Event No "<<j<<endl;
63  for (Int_t i=0; i<hit_bar->GetEntriesFast(); i++)
64  {
65  //cout <<"Hit "<< i << endl;
66  if (i == 0) // because of trackID
67  {
68  //cout << hit_bar->GetEntriesFast()<<endl;
69  PndHypGePoint *hitgam=(PndHypGePoint*)hit_bar->At(i);
70  ActualTrackID=-1;
71  PndMCTrack *mcgam = (PndMCTrack*)mc_bar->At(hitgam->GetTrackID());
72  if (hitgam->GetpdgCode()==2112 && mcgam->GetMotherID()==-1)
73  {
74  ActualTrackID= hitgam->GetTrackID();
75  TVector3 NHit;
76  NHit.SetX(hitgam->GetX());
77  NHit.SetY(hitgam->GetY());
78  NHit.SetZ(hitgam->GetZ()+55);
79  hNHits->Fill(180/TMath::Pi()*NHit.Theta()); //fill histogram with polar angle of detector interactions of primary neutrons
80 
81  TVector3 NHit_p;
82  NHit_p.SetX(hitgam->GetPx());
83  NHit_p.SetY(hitgam->GetPy());
84  NHit_p.SetZ(hitgam->GetPz());
85  hNHits_p->Fill(180/TMath::Pi()*NHit_p.Theta()); //fill histogram with the polar momentum of primary neutrons interacting with the crystals
86  hCrystalHit->Fill(hitgam->GetDetectorID()); //fill histogram to see which detector is hit
87  //hNHits->Fill(180/TMath::Pi()*TMath::ACos(hitgam->GetZ()/TMath::Sqrt(hitgam->GetX()*hitgam->GetX()+hitgam->GetY()*hitgam->GetY()+hitgam->GetZ()* hitgam->GetZ()))); //Fill spectrum
88  }
89  if (mcgam->GetPdgCode()==2112 && mcgam->GetMotherID()==-1)
90  {
91  TVector3 NAllMom = mcgam->GetMomentum();
92  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)
93  hNMom->Fill(1./2.*NAllMom.Mag2()/0.939); //fill histogram with E_kin of primary neutrons
94  //cout << 1./2.*NAllMom.Mag2()/0.939<<endl;
95  }
96  }
97  else
98  {
99  PndHypGePoint *hitgam=(PndHypGePoint*)hit_bar->At(i);
100  PndMCTrack *mcgam = (PndMCTrack*)mc_bar->At(hitgam->GetTrackID());
101  //rande= gRandom->Gaus(hitgam->GetEnergyLoss(),0.000003);
102  if (hitgam->GetpdgCode()==2112 && mcgam->GetMotherID()==-1)
103  {
104  if (ActualTrackID != hitgam->GetTrackID())
105  {
106  TVector3 NHit;
107  NHit.SetX(hitgam->GetX());
108  NHit.SetY(hitgam->GetY());
109  NHit.SetZ(hitgam->GetZ()+55);
110  hNHits->Fill(180/TMath::Pi()*NHit.Theta());
111  TVector3 NHit_p;
112  NHit_p.SetX(hitgam->GetPx());
113  NHit_p.SetY(hitgam->GetPy());
114  NHit_p.SetZ(hitgam->GetPz());
115  hNHits_p->Fill(180/TMath::Pi()*NHit_p.Theta());
116  ActualTrackID = hitgam->GetTrackID();
117  hCrystalHit->Fill(hitgam->GetDetectorID());
118  }
119  }
120  if (mcgam->GetPdgCode()==2112 && mcgam->GetMotherID()==-1)
121  {
122  TVector3 NAllMom = mcgam->GetMomentum();
123  hNAll->Fill(180/TMath::Pi()*NAllMom.Theta());
124  hNMom->Fill(1./2.*NAllMom.Mag2()/0.939);
125  //cout << 1./2.*NAllMom.Mag2()/0.939<<endl;
126  }
127 
128  }
129  }
130 
131 
132  }// end for j (events)
133 
134 
135 
136  hNHits->Draw();
137 
138  //Analysis of spectrum
139  hNHits->Write();
140  hNHits_p->Write();
141  hNAll->Write();
142  hNMom->Write();
143  hCrystalHit->Write();
144 
145  OutputFile->Close();
146 
147 
148 
149 
150  // ----- Finish -------------------------------------------------------
151  timer.Stop();
152  Double_t rtime = timer.RealTime();
153  Double_t ctime = timer.CpuTime();
154  cout << endl << endl;
155  cout << "Macro finished succesfully." << endl;
156 
157  cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
158  cout << endl;
159  // ------------------------------------------------------------------------
160 
161  return 0;
162 }
i
Int_t i
Definition: run_full.C:25
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::GetPz
Double_t GetPz() const
Definition: PndHypGePoint.h:59
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
PndHypGePoint::GetPy
Double_t GetPy() const
Definition: PndHypGePoint.h:58
nEvents
Int_t nEvents
Definition: hit_dirc.C:11
timer
TStopwatch timer
Definition: hit_dirc.C:51
PndHypGePoint::GetPx
Double_t GetPx() const
Definition: PndHypGePoint.h:57
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::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