FairRoot/PandaRoot
Functions
NeutronAnalysis_COSY_edit.C File Reference
#include "TH1D.h"
#include "TH2D.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 "TSystem.h"
#include "TStyle.h"
#include <stdio.h>
#include <iostream>

Go to the source code of this file.

Functions

int NeutronAnalysis_COSY_edit (TString Filename_ext="hypGeCOSYBeamTime2014_protons_mom2.7799999999999998GeV_1000000Evts_GausianSmearedBeam_9.root", Int_t StartEvent=0, Int_t NoOfEvents=100000)
 

Function Documentation

int NeutronAnalysis_COSY_edit ( TString  Filename_ext = "hypGeCOSYBeamTime2014_protons_mom2.7799999999999998GeV_1000000Evts_GausianSmearedBeam_9.root",
Int_t  StartEvent = 0,
Int_t  NoOfEvents = 100000 
)

Definition at line 20 of file NeutronAnalysis_COSY_edit.C.

References b, Bool_t, ctime, Double_t, PndHypGePoint::GetDetectorID(), PndHypGePoint::GetEnergyLoss(), PndHypGePoint::GetpdgCode(), PndHypGePoint::GetPx(), PndHypGePoint::GetPy(), PndHypGePoint::GetPz(), PndMCTrack::GetStartVertex(), PndHypGePoint::GetTrackID(), PndHypGePoint::GetX(), PndHypGePoint::GetY(), PndHypGePoint::GetZ(), gGeoManager, i, outfile, Pi, rtime, timer, TString, and verbose.

21 {
22  bool verbose = false;
23  // ----- Load libraries ------------------------------------------------
24  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
25  gSystem->Load("libHypGe");
26  // ----- Timer --------------------------------------------------------
27  TStopwatch timer;
28  timer.Start();
29  // ------------------------------------------------------------------------
30 
31  TGeoManager *fgeom;
32  //Files
33  //TString Path = getenv("SIMDATADIR");
34  TString Filename = "$SIMDATADIR/COSY/"+Filename_ext;
35  TString ParfileName = "$SIMDATADIR/COSY/"+Filename_ext;
36  ParfileName.ReplaceAll(".root","__Simparams.root");
37  //TString Filename = "$SIMDATADIR/TripleBall40Offset10_urqmd_pbarC_1_5000Evts.root";
38  TString outfile= "$SIMDATADIR/COSY/Ana/ToCombine/" + Filename_ext;
39  outfile += "/Ana";
40  outfile += Filename_ext;
41  outfile.ReplaceAll(".root","");
42  outfile += "_Start_";
43  outfile += StartEvent;
44  outfile += "_nEvents_";
45  outfile += NoOfEvents;
46  outfile += ".root";
47  cout << outfile << endl;
48  TFile* InputFile = new TFile(Filename);
49  TFile* Parfile = new TFile(ParfileName);
50  Parfile->Get("FairBaseParSet");
51  TFile* OutputFile = new TFile(outfile,"RECREATE");
52 
53  //getting simulation branches from input file
54  TTree *b=(TTree *) InputFile->Get("pndsim") ;
55  TClonesArray* fHypGe=new TClonesArray("PndHypGePoint");
56  b->SetBranchAddress("HypGePoint",&fHypGe);//Branch names
57  TClonesArray* fMcTr=new TClonesArray("PndMCTrack");
58  b->SetBranchAddress("MCTrack",&fMcTr);//Branch names
59 
60 
61  //declaring variables from PndHypGeCOSYBackgroundAna.h (macro updated for 2014 beam time changes, file not yet, 20.08.14 steinen)
62  TGeoManager *fgeom;
63  TH1D* hNHitsAngleGe; // angular distribution of neutrons in the germanium
64  TH1D* hNeutronOriginGe; // neutrons and their origin in the germanium
65  TH1D* hNHitsGeKermitBall; // neutrons in ge, active and passive neutron detector
66 
67  TH1D* hAllParticlesGe; // all particles hitting the germanium
68  TH2D* hAllParticlesEkinGe; // 2D all E_kin distribution of particles hitting the germanium
69  TH1D* hAllParticlesPiezo; // all particles hitting the first pieo
70  TH1D* hAllParticlesPiezo2; // all particles hitting the germanium
71  TH1D* hNEkinGe; // E_kin of neutrons entering the germanium
72  TH1D* hNEkinKermit; // E_kin of neutrons entering the active neutron detector
73  TH1D* hNEkinBall; // E_kin of neutrons entering the passive neutron detector
74  TH1D* hNEnergyLossGe; // Energy loss of neutrons in the germanium
75  TH1D* hNEnergyLossKermit; // Energy loss of neutrons in the active neutron detector
76  TH1D* hNEnergyLossBall; // Energy loss of neutrons in the passive neutron detector
77  TH1D* hPEkinGe; // E_kin of protons entering the germanium
78  TH1D* hPEnergyLossGe; // Energy loss of protons in the germanium
79  TH1D* hGammaSpecGe; // spectrum of gammas in the germanium
80 
81 
82  TCanvas* cNHitsGe; //canvas for hNHitsAngleGe and hNeutronOriginGe
83  TCanvas* cNHitsGeKermitBall; //canvas for hNHitsGeKermitBall
84  TCanvas* cAllParticlesGe; //canvas for hAllParticlesGe and hAllParticlesEkinGe
85  TCanvas* cAllParticlesPiezo; //canvas for hAllParticlesPiezo and hAllParticlesPiezo2
86  TCanvas* cNeutronEkinGe; //canvas for hNEkinGe and hNEnergyLossGe
87  TCanvas* cNeutronEkinKermit; //canvas for hNEkinKermit and hNEnergyLossKermit
88  TCanvas* cNeutronEkinBall; //canvas for hNEkinBall and hNEnergyLossBall
89  TCanvas* cProtonEkinGe; //canvas for hPEkinGe and hNEnergyLossGe
90  TCanvas* cGammaEnergyGe; //canvas for hGammaSpecGe
91 
92 
93  TString VertexVolumeName;
94  TVector3 NeutronMomentum;
95  Double_t NeutronEkin;
96  TVector3 ProtonMomentum;
97  Double_t ProtonEkin;
98  TVector3 ParticleMomentum;
99  Double_t ParticleEkin;
100 
101  //----------Get GeoManager---------------------------------------------
102  fgeom = (TGeoManager*)gROOT->FindObject("FAIRGeom");
103  //----------create histograms------------------------------------------
104 
105  //histogram with polar angle of detector interactions of primary neutrons
106  hNHitsAngleGe = new TH1D("hNHitsAngleGe","Polar angle of primary neutrons interacting with the crystals",180, 90,180);
107  hNHitsAngleGe->SetXTitle("#Theta [#circ]");
108  hNHitsAngleGe->SetYTitle("Counts / 0.5 #circ");
109 
110  //histogram to see where neutron is coming from
111  hNeutronOriginGe = new TH1D("hNeutronOriginGe","Origin of Neutrons",13,0,13);
112  hNeutronOriginGe->SetXTitle("Neutron Origin");
113  hNeutronOriginGe->SetYTitle("Counts");
114  hNeutronOriginGe->GetXaxis()->SetBinLabel(1,"Target"); // neutrons created in the carbon target
115  hNeutronOriginGe->GetXaxis()->SetBinLabel(2,"Target Support"); // neutrons created in the goniometer, the brackets
116  hNeutronOriginGe->GetXaxis()->SetBinLabel(3,"Al base plate"); // neutrons created in the al support plate
117  hNeutronOriginGe->GetXaxis()->SetBinLabel(4,"Plastic blocks"); // neutrons created in the plastic blocks under the al plate
118  hNeutronOriginGe->GetXaxis()->SetBinLabel(5,"Maytec frame"); // neutrons created in the maytec frame supporting the act. neutron det.
119  hNeutronOriginGe->GetXaxis()->SetBinLabel(6,"Active neutron detector"); // neutrons created in the active neutron detector
120  hNeutronOriginGe->GetXaxis()->SetBinLabel(7,"Passive neutron detector"); // neutrons created in the passive neutron detector
121  hNeutronOriginGe->GetXaxis()->SetBinLabel(8,"Crystal"); // neutrons created in the germanium crystal
122  hNeutronOriginGe->GetXaxis()->SetBinLabel(9,"Cryostat"); // neutrons created in the germanium detector cryostat
123  hNeutronOriginGe->GetXaxis()->SetBinLabel(10,"Lumi glue test"); // neutrons created in the lumi glue test stuff
124  hNeutronOriginGe->GetXaxis()->SetBinLabel(11,"Lumi electronics test"); // neutrons created in the lumi electronics test stuff
125  hNeutronOriginGe->GetXaxis()->SetBinLabel(12,"Beam dump"); // neutrons created in the beam dump
126  hNeutronOriginGe->GetXaxis()->SetBinLabel(13,"Cave"); // neutrons created in the experimental area cave
127 
128  //histogram to see which detector is hit
129  hNHitsGeKermitBall = new TH1D("hNHitsGeKermitBall","Total number of neutron hits on the 3 detectors",3,0,1);
130  hNHitsGeKermitBall->SetXTitle("Crystal number");
131  hNHitsGeKermitBall->SetYTitle("Counts per crystal");
132  hNHitsGeKermitBall->GetXaxis()->SetBinLabel(1,"Germanium crystal"); //Crystal_101
133  hNHitsGeKermitBall->GetXaxis()->SetBinLabel(2,"Active neutron detector"); //Kermit
134  hNHitsGeKermitBall->GetXaxis()->SetBinLabel(3,"Passive neutron detector"); //Ball
135 
136  //histogram with all particles in the germanium
137  hAllParticlesGe = new TH1D("hAllParticlesGe","Particles interaction in the crystals;PDG Code of particle; Counts", 2e4,-1e4,1e4);
138 
139  // histogram for correlation of E_kin and particle PDG Code for the germanium
140  hAllParticlesEkinGe = new TH2D("hAllParticlesEkinGe","E_{kin} - Particle - Correlation;E_{kin} of particles [MeV]; particle",3000,0,3000,9,1,10);
141  hAllParticlesEkinGe->GetYaxis()->SetBinLabel(1,"#pi^{-}");
142  hAllParticlesEkinGe->GetYaxis()->SetBinLabel(2,"#mu^{+}");
143  hAllParticlesEkinGe->GetYaxis()->SetBinLabel(3,"e^{+}");
144  hAllParticlesEkinGe->GetYaxis()->SetBinLabel(4,"e^{-}");
145  hAllParticlesEkinGe->GetYaxis()->SetBinLabel(5,"#mu^{-}");
146  hAllParticlesEkinGe->GetYaxis()->SetBinLabel(6,"#gamma");
147  hAllParticlesEkinGe->GetYaxis()->SetBinLabel(7,"#pi^{+}");
148  hAllParticlesEkinGe->GetYaxis()->SetBinLabel(8,"neutron");
149  hAllParticlesEkinGe->GetYaxis()->SetBinLabel(9,"proton");
150 
151  //histogram with all particles in piezo created outside of the crystal
152  hAllParticlesPiezo = new TH1D("hAllParticlesPiezo","Particles interaction in the first piezo;PDG Code of particle; Counts", 2e4,-1e4,1e4);
153  //histogram with all particles in piezo created outside of the crystal
154  hAllParticlesPiezo2 = new TH1D("hAllParticlesPiezo2","Particles interaction in the second piezo;PDG Code of particle; Counts", 2e4,-1e4,1e4);
155 
156 
157  //histogram for kin. energy of neutrons in the germanium
158  hNEkinGe = new TH1D("hNEkinGe","E_{kin} of neutrons in the germanium;E_{kin} of neutrons [MeV]; Counts", 2000,0,200);
159  //histogram for kin. energy of neutrons in the active neutron detector
160  hNEkinKermit = new TH1D("hNEkinKermit","E_{kin} of neutrons in the active neutron detector;E_{kin} of neutrons [MeV]; Counts", 2000,0,200);
161  //histogram for kin. energy of neutrons in the passive neutron detector
162  hNEkinBall = new TH1D("hNEkinBall","E_{kin} of neutrons in the passive neutron detector;E_{kin} of neutrons [MeV]; Counts", 2000,0,200);
163 
164  //histogram with energy deposited by neutrons in the germanium
165  hNEnergyLossGe = new TH1D("hNEnergyLossGe","Energy loss of neutrons inside the germanium;Energy loss of neutrons [MeV]; Counts", 2000,0,2);
166  //histogram with energy deposited by neutrons in the active neutron detector
167  hNEnergyLossKermit = new TH1D("hNEnergyLossKermit","Energy loss of neutrons inside the active neutron detector;Energy loss of neutrons [MeV]; Counts", 2000,0,2);
168  //histogram with energy deposited by neutrons in the passive neutron detector
169  hNEnergyLossBall = new TH1D("hNEnergyLossBall","Energy loss of neutrons inside the passive neutron detecor;Energy loss of neutrons [MeV]; Counts", 2000,0,2);
170 
171  //histogram for kin. energy of protons in the germanium
172  hPEkinGe = new TH1D("hPEkinGe","E_{kin} of protons in the germanium;E_{kin} of protons [MeV]; Counts", 2000,0,200);
173  //histogram with energy deposited by protons in the germanium
174  hPEnergyLossGe = new TH1D("hPEnergyLossGe","Energy loss of protons inside the germanium;Energy loss of protons [MeV]; Counts", 20000,0,200);
175 
176  //histogram of gamma spectrum in the germanium
177  hGammaSpecGe = new TH1D("hGammaSpecGe","Energy loss of gammas inside the germanium;Energy loss of gammas [MeV]; Counts/1 keV", 10000,0,10);
178 
179 
180 
181 
182  Int_t *ActualTrackID ;
183  //Int_t nEvents = b->GetEntriesFast();
184  //cout<< "Number of Simulated Events: "<<nEvents<<endl;
185 
186  Int_t Hits=0;
187  Int_t DetID=-10;
188  Int_t TrackID=-10;
189  Int_t iNeutrons = 0;
190 
191 
192 
193  TVector3 NHit;
194  TVector3 NMom;
195  TVector3 NAllMom;
196  TVector3 StartVertex;
197  TVector3 PHit;
198  TVector3 PMom;
199 
200  Int_t EndEvent;
201  if (NoOfEvents ==0)
202  EndEvent = b->GetEntriesFast();
203  else
204  EndEvent = StartEvent + NoOfEvents;
205 
206  Double_t NeutronEnergyLossArrayGe[50]; // array to store the energy loss of all neutrons in the germanium during one event
207  Double_t NeutronEnergyLossArrayKermit[50]; // array to store the energy loss of all neutrons in the active neutron detector during one event
208  Double_t NeutronEnergyLossArrayBall[50]; // array to store the energy loss of all neutrons in the passive neutron detector during one event
209  Double_t ProtonEnergyLossArrayGe[50]; // array to store the energy loss of all protons in the germanium during one event
210  Double_t GammaEnergyLossArrayGe[50]; // array to store the energy loss of all gammas in the germanium during one event
211 
212  //event loop
213  for (Int_t k=StartEvent; k< EndEvent; k++)
214  //for (Int_t k=569494; k<569498; k++)
215  {
216 
217  for(Int_t i = 0; i < 50; i++) //init of array
218  {
219  NeutronEnergyLossArrayGe[i]=0;
220  NeutronEnergyLossArrayKermit[i]=0;
221  NeutronEnergyLossArrayBall[i]=0;
222  ProtonEnergyLossArrayGe[i] =0;
223  GammaEnergyLossArrayGe[i] = 0;
224  }
225  Int_t nNeutronsGe = -1;
226  Int_t NeutronTrackIDGe =-1;
227  Bool_t newNeutronGe = 1;
228  Int_t nNeutronsKermit = -1;
229  Int_t NeutronTrackIDKermit =-1;
230  Bool_t newNeutronKermit = 1;
231  Int_t nNeutronsBall = -1;
232  Int_t NeutronTrackIDBall =-1;
233  Bool_t newNeutronBall = 1;
234  Int_t nProtonsGe = -1;
235  Int_t ProtonTrackIDGe =-1;
236  Bool_t newProton = 1;
237 
238  Int_t nGammasGe = -1;
239  Int_t GammaTrackIDGe =-1;
240 
241  b->GetEntry(k);
242  if (!((k*100)% NoOfEvents)) // mark every % of proceeded events
243  {
244  cout << "Event number :\t" << k << endl;
245  }
246  //if(verbose) cout<<"Event No "<<j<<endl;
247 
248  DetID=-10;
249  TrackID=-10;
250 
251  // loop over the hits of the event
252  for (Int_t i=0; i < fHypGe->GetEntriesFast(); i++)
253  {
254  PndHypGePoint *hitgam=(PndHypGePoint*)fHypGe->At(i);
255  PndMCTrack *mcgam = (PndMCTrack*)fMcTr->At(hitgam->GetTrackID());
256 
257  StartVertex = mcgam->GetStartVertex();
258  VertexVolumeName =gGeoManager->FindNode(StartVertex.X(),StartVertex.Y(),StartVertex.Z())->GetName();
259  //only Neutrons
260  if (hitgam->GetpdgCode() == 2112)
261  {
262  //cout <<"Event " <<EvtCount<< " \t\tStartVertex Radius: " <<StartVertex.Mag()<<endl;
263  //Fill hNeutronOriginGe
264 
265  //cout << VertexVolumeName << endl;
266 
267  //calculate E_kin of neutron
268  NeutronMomentum.SetX(hitgam->GetPx());
269  NeutronMomentum.SetY(hitgam->GetPy());
270  NeutronMomentum.SetZ(hitgam->GetPz());
271  NeutronEkin = 1./2.*NeutronMomentum.Mag2()/0.939565*1000; // in MeV
272  //cout << NeutronEkin << "\t\t\t" << NeutronMomentum.Mag2()<<endl;
273 
274 
275 
276  if(hitgam->GetDetectorID()==101)
277  {
278 
279  //cut on neutrons created outside of the germanium detecto
280  if (!VertexVolumeName.Contains("Crystal_101") && ! VertexVolumeName.Contains ("Capsule") && !VertexVolumeName.Contains ("ColdFinger"))
281  {
282  if (newNeutronGe)
283  {
284  // fill Ekin histogram
285  hNHitsGeKermitBall->Fill("Germanium crystal",1);
286  hNEkinGe->Fill(NeutronEkin);
287  newNeutronGe = 0;
288  }
289  // calculate energy loss of neutron in germanium
290  if (NeutronTrackIDGe !=hitgam->GetTrackID()) // new neutron (different TrackID)
291  {
292  NeutronTrackIDGe = hitgam->GetTrackID();
293  nNeutronsGe++;
294  newNeutronGe = 1;
295  //cout << nNeutronsGe << endl;
296  }
297  NeutronEnergyLossArrayGe[nNeutronsGe] += hitgam->GetEnergyLoss()*1000;
298  //cout << hitgam->GetEnergyLoss()*1000<< endl;
299 
300  // calculate theta of neutron
301  NHit.SetX(hitgam->GetX());
302  NHit.SetY(hitgam->GetY());
303  NHit.SetZ(hitgam->GetZ()); // no shift like in PANDA!!!!
304  hNHitsAngleGe->Fill(180/TMath::Pi()*NHit.Theta()); //fill histogram with polar angle of detector interactions of primary neutrons
305  // fill Ekin histogram
306  }
307 
308  if( VertexVolumeName.Contains("Target") || VertexVolumeName.Contains("Piezo"))
309  hNeutronOriginGe->Fill("Target",1);
310  if( VertexVolumeName.Contains("ScrewHead") || VertexVolumeName.Contains("Goniometer")|| VertexVolumeName.Contains("Bracket")|| VertexVolumeName.Contains("AlPlate"))
311  hNeutronOriginGe->Fill("Target Support",1);
312  if( VertexVolumeName.Contains("AlBasePlate"))
313  hNeutronOriginGe->Fill("Al base plate",1);
314  if( VertexVolumeName.Contains("PlasticBlock"))
315  hNeutronOriginGe->Fill("Plastic blocks",1);
316  if( VertexVolumeName.Contains("Beam") || VertexVolumeName.Contains("Pillar"))
317  hNeutronOriginGe->Fill("Maytec frame",1);
318  if( VertexVolumeName.Contains("KermitCrystal"))
319  hNeutronOriginGe->Fill("Active neutron detector",1);
320  if( VertexVolumeName.Contains("KermitCrystal"))
321  hNeutronOriginGe->Fill("Active neutron detector",1);
322  if( VertexVolumeName.Contains("NeutronBallCrystal"))
323  hNeutronOriginGe->Fill("Passive neutron detector",1);
324  if (VertexVolumeName.Contains("Crystal_101"))
325  hNeutronOriginGe->Fill("Crystal",1);
326  if (VertexVolumeName.Contains ("Capsule") || VertexVolumeName.Contains ("ColdFinger") || VertexVolumeName.Contains ("Cryostat"))
327  hNeutronOriginGe->Fill("Cryostat",1);
328  if (VertexVolumeName.Contains("Glue"))
329  hNeutronOriginGe->Fill("Lumi glue test",1);
330  if (VertexVolumeName.Contains("LumiElectronics"))
331  hNeutronOriginGe->Fill("Lumi electronics test",1);
332  if (VertexVolumeName.Contains("BeamDump"))
333  hNeutronOriginGe->Fill("Beam dump",1);
334  if (VertexVolumeName.Contains("cave"))
335  hNeutronOriginGe->Fill("Cave",1);
336 
337  }
338  if(hitgam->GetDetectorID()==2000)
339  {
340 
341  //cut on neutrons created outside of the active neutron detector
342  if (!VertexVolumeName.Contains("Kermit") )
343  {
344  if (newNeutronKermit)
345  {
346  // fill Ekin histogram
347  hNHitsGeKermitBall->Fill("Active neutron detector",1);
348  hNEkinKermit->Fill(NeutronEkin);
349  newNeutronKermit = 0;
350  }
351  // calculate energy loss of neutron in active neutron detektor
352  if (NeutronTrackIDKermit !=hitgam->GetTrackID()) // new neutron (different TrackID)
353  {
354  NeutronTrackIDKermit = hitgam->GetTrackID();
355  nNeutronsKermit++;
356  newNeutronKermit = 1;
357  //cout << nNeutronsKermit << endl;
358  }
359  NeutronEnergyLossArrayKermit[nNeutronsKermit] += hitgam->GetEnergyLoss()*1000;
360  //cout << hitgam->GetEnergyLoss()*1000<< endl;
361 
362  hNHitsGeKermitBall->Fill("Active neutron detector",1);
363  // fill histogram with Ekin of neutrons in the active neutron detector
364  hNEkinKermit->Fill(NeutronEkin);
365  }
366 
367 
368  }
369  if(hitgam->GetDetectorID()==2001)
370  {
371 
372  // cut on neutrons created outside of the passive neutron detector
373  if (!VertexVolumeName.Contains("NeutronBall") )
374  {
375  if (newNeutronBall)
376  {
377  // fill Ekin histogram
378  hNHitsGeKermitBall->Fill("Passive neutron detector",1);
379  hNEkinBall->Fill(NeutronEkin);
380  newNeutronBall = 0;
381  }
382  // calculate energy loss of a neutron in passive neutron detector
383  if (NeutronTrackIDBall !=hitgam->GetTrackID()) // new neutron (different TrackID)
384  {
385  NeutronTrackIDBall = hitgam->GetTrackID();
386  nNeutronsBall++;
387  newNeutronBall = 1;
388  //cout << nNeutronsBall << endl;
389  }
390  NeutronEnergyLossArrayBall[nNeutronsBall] += hitgam->GetEnergyLoss()*1000;
391  //cout << hitgam->GetEnergyLoss()*1000<< endl;
392 
393 
394  }
395 
396 
397  }
398  } // end of only neutrons
399  // only protons
400  /*if (hitgam->GetpdgCode() == 2212)
401  {
402  if(hitgam->GetDetectorID()==101)
403  {
404  //calculate E_kin of proton
405  ProtonMomentum.SetX(hitgam->GetPx());
406  ProtonMomentum.SetY(hitgam->GetPy());
407  ProtonMomentum.SetZ(hitgam->GetPz());
408  ProtonEkin = 1./2.*ProtonMomentum.Mag2()/0.938272*1000; // in MeV
409  //cout << ProtonEkin << "\t\t\t" << ProtonMomentum.Mag2()<<endl;
410 
411 
412  //cut on protons created outside of the germanium detector
413  if (!VertexVolumeName.Contains("Crystal_101") && ! VertexVolumeName.Contains ("Capsule") && ! VertexVolumeName.Contains ("ColdFinger"))
414  {
415  // if first occurence of the proton in the detecor fill the E_kin
416  if (newProton)
417  {
418  // fill Ekin histogram
419  hPEkinGe->Fill(ProtonEkin);
420  newProton = 0;
421  }
422  // calculate energy loss of proton in germanium
423  if (ProtonTrackIDGe !=hitgam->GetTrackID()) // new proton (different TrackID)
424  {
425  ProtonTrackIDGe = hitgam->GetTrackID();
426  nProtonsGe++;
427  newProton = 1;
428  //cout << nProtonsGe << endl;
429  }
430  ProtonEnergyLossArrayGe[nProtonsGe] += hitgam->GetEnergyLoss()*1000;
431  cout << hitgam->GetEnergyLoss()*1000<< endl;
432 
433  }
434  }
435  } // end of only protons
436 
437  // only gammas
438  if (hitgam->GetpdgCode() == 22 )
439  {
440  if(hitgam->GetDetectorID()==101)
441  {
442  // calculate energy loss of gamma in the germanium
443  if (GammaTrackIDGe !=hitgam->GetTrackID()) // new gamma (different TrackID)
444  {
445  GammaTrackIDGe = hitgam->GetTrackID();
446  nGammasGe++;
447  }
448  GammaEnergyLossArrayGe[nGammasGe] += hitgam->GetEnergyLoss()*1000;
449  cout << nGammasGe << "\t" << GammaEnergyLossArrayGe[nGammasGe]<< "\t" << hitgam->GetEnergyLoss()*1000 << endl;
450  }
451  } // end of only gammas
452 
453  //all particles in Crystal
454  if (!VertexVolumeName.Contains("Crystal_101") && ! VertexVolumeName.Contains ("Capsule") && !VertexVolumeName.Contains ("ColdFinger")) //created outside of the detector
455  {
456  if(hitgam->GetDetectorID()==101 || hitgam->GetDetectorID()==102)
457  hAllParticlesGe->Fill(hitgam->GetpdgCode());
458 
459  ParticleMomentum.SetX(hitgam->GetPx());
460  ParticleMomentum.SetY(hitgam->GetPy());
461  ParticleMomentum.SetZ(hitgam->GetPz());
462  Double_t ParticleMass = 0;
463  TString ParticleName = "";
464  switch ((int) hitgam->GetpdgCode())
465  {
466  case -211: ParticleMass = 0.13957; ParticleName ="#pi^{-}"; break; //pi-
467  case -13: ParticleMass = 0.10561; ParticleName ="#mu^{+}"; break; //mu+
468  case -11: ParticleMass = 0.000511; ParticleName ="e^{+}"; break; //e+
469  case 11: ParticleMass = 0.000511; ParticleName ="e^{-}"; break; //e-
470  case 13: ParticleMass = 0.10561; ParticleName ="#mu^{-}"; break; //mu-
471  case 22: ParticleMass = 0; ParticleName ="#gamma"; break; //gamma
472  case 211: ParticleMass = 0.13957; ParticleName ="#pi^{+}"; break; //pi+
473  case 2112: ParticleMass = 0.939565; ParticleName ="neutron"; break; //n
474  case 2212: ParticleMass = 0.93827; ParticleName ="proton"; break; //p
475  default: ParticleMass = -1; std::cout.precision(10); cout << "Found new particle:\t" << hitgam->GetpdgCode() <<"\t" << VertexVolumeName << endl; break;
476  }
477  if(ParticleMass != -1)
478  {
479  ParticleEkin = (sqrt(ParticleMass*ParticleMass+ParticleMomentum.Mag2())-ParticleMass)*1000;//Ekin in MeV
480  //all particles (but gammas) in Crystal - E_{kin} - correlation
481  if (!VertexVolumeName.Contains("Crystal_101") && ! VertexVolumeName.Contains ("Capsule")) //created outside of the detector
482  {
483  if(hitgam->GetDetectorID()==101)
484  hAllParticlesEkinGe->Fill( ParticleEkin , ParticleName.Data(),1);
485  }
486  }
487 
488  } // end of all particles part
489  //piezos
490  if(hitgam->GetDetectorID()==1000)
491  hAllParticlesPiezo ->Fill(hitgam->GetpdgCode());
492  if(hitgam->GetDetectorID()==1001)
493  hAllParticlesPiezo2 ->Fill(hitgam->GetpdgCode());*/
494  } // end of for loop over the hits
495 
496  //fill Neutron Engergy loss histograms
497  //germanium
498  if (nNeutronsGe>-1)
499  {
500  // cout << nNeutronsGe+1 << endl;
501  for (Int_t i = 0; i <= nNeutronsGe; i++)
502  {
503  // cout << i << "\t\t" << NeutronEnergyLossArrayGe[i] << endl;
504  hNEnergyLossGe->Fill(NeutronEnergyLossArrayGe[i]);
505  }
506  }
507  //active neutron detector
508  if (nNeutronsKermit>-1)
509  {
510  // cout << nNeutronsKermit+1 << endl;
511  for (Int_t i = 0; i <= nNeutronsKermit; i++)
512  {
513  // cout << i << "\t\t" << NeutronEnergyLossArrayKermit[i] << endl;
514  hNEnergyLossKermit->Fill(NeutronEnergyLossArrayKermit[i]);
515  }
516  }
517  //passive neutron detector
518  if (nNeutronsBall>-1)
519  {
520  // cout << nNeutronsBall+1 << endl;
521  for (Int_t i = 0; i <= nNeutronsBall; i++)
522  {
523  // cout << i << "\t\t" << NeutronEnergyLossArrayBall[i] << endl;
524  hNEnergyLossBall->Fill(NeutronEnergyLossArrayBall[i]);
525  }
526  }
527  if (nProtonsGe>-1)
528  {
529  // cout << nProtonsGe+1 << endl;
530  for (Int_t i = 0; i <= nProtonsGe; i++)
531  {
532  // cout << i << "\t\t" << ProtonEnergyLossArrayGe[i] << endl;
533  hPEnergyLossGe->Fill(ProtonEnergyLossArrayGe[i]);
534  }
535  }
536  if (nGammasGe>-1)
537  {
538  // cout << nGammasGe+1 << endl;
539  for (Int_t i = 0; i <= nGammasGe; i++)
540  {
541  // cout << i << "\t\t" << GammaEnergyLossArrayGe[i] << endl;
542  hGammaSpecGe->Fill(GammaEnergyLossArrayGe[i]);
543  }
544  }
545  }// end for k (events)
546 
547 
548 
549 
550  cNHitsGe = new TCanvas("cNHitsGe","Neutron Hits",1600,600);
551  cNHitsGe->Divide(2,1);
552 
553  cNHitsGeKermitBall = new TCanvas("cNHitsGeKermitBall","Neutron Hits on all 3 detectors",800,600);
554 
555  cAllParticlesGe = new TCanvas("cAllParticles","All particles interacting in the germanium crystal", 1600, 600);
556  cAllParticlesGe->Divide(2,1);
557 
558  cAllParticlesPiezo = new TCanvas("cAllParticlesPiazo","All particles interacting in the piezo", 1600, 600);
559  cAllParticlesPiezo->Divide(2,1);
560 
561  cNeutronEkinGe = new TCanvas("cNeutronEkinGe","E_{kin} of neutrons", 1600, 600);
562  cNeutronEkinGe->Divide(2,1);
563 
564  cNeutronEkinKermit = new TCanvas("cNeutronEkinKermit","E_{kin} of neutrons", 1600, 600);
565  cNeutronEkinKermit->Divide(2,1);
566 
567  cNeutronEkinBall = new TCanvas("cNeutronEkinBall","E_{kin} of neutrons", 1600, 600);
568  cNeutronEkinBall->Divide(2,1);
569 
570  cProtonEkinGe = new TCanvas("cProtonEkinGe","E_{kin} of protons", 1600, 600);
571  cProtonEkinGe->Divide(2,1);
572 
573  cGammaEnergyGe= new TCanvas("cGammaEnergyGe","#gamma spectrum", 800, 600);
574 
575  cNHitsGe->cd(1);
576  hNHitsAngleGe->DrawCopy(); // angular distribution of neutrons in the germanium
577  cNHitsGe->cd(2);
578  hNeutronOriginGe->DrawCopy();// neutrons and their origin in the germanium
579 
580  cNHitsGeKermitBall->cd();
581  hNHitsGeKermitBall->DrawCopy();// neutrons in ge, active and passive neutron detector
582 
583  cAllParticlesGe->cd(1);
584  hAllParticlesGe->DrawCopy(); // all particles hitting the germanium
585  cAllParticlesGe->cd(2);
586  hAllParticlesEkinGe->DrawCopy("colz");// 2D all E_kin distribution of particles hitting the germanium
587 
588  cAllParticlesPiezo->cd(1);
589  hAllParticlesPiezo->DrawCopy(); // all particles hitting the first pieo
590  cAllParticlesPiezo->cd(2);
591  hAllParticlesPiezo2->DrawCopy();// all particles hitting the germanium
592 
593  cNeutronEkinGe->cd(1);
594  hNEkinGe->DrawCopy(); // E_kin of neutrons entering the germanium
595  cNeutronEkinKermit->cd(1);
596  hNEkinKermit->DrawCopy(); // E_kin of neutrons entering the active neutron detector
597  cNeutronEkinBall->cd(1);
598  hNEkinBall->DrawCopy(); // E_kin of neutrons entering the passive neutron detector
599 
600  cNeutronEkinGe->cd(2);
601  hNEnergyLossGe->DrawCopy(); // Energy loss of neutrons in the germanium
602  cNeutronEkinKermit->cd(2);
603  hNEnergyLossKermit->DrawCopy(); // Energy loss of neutrons in the active neutron detector
604  cNeutronEkinBall->cd(2);
605  hNEnergyLossBall->DrawCopy(); // Energy loss of neutrons in the passive neutron detector
606 
607  cProtonEkinGe->cd(1);
608  hPEkinGe->DrawCopy(); // E_kin of protons entering the germanium
609  cProtonEkinGe->cd(2);
610  hPEnergyLossGe->DrawCopy(); // Energy loss of protons in the germanium
611 
612  cGammaEnergyGe->cd();
613  hGammaSpecGe->DrawCopy(); // spectrum of gammas in the germanium
614 
615 
616 
617  hNHitsAngleGe->Write(); // angular distribution of neutrons in the germanium
618  hNeutronOriginGe->Write();// neutrons and their origin in the germanium
619  hNHitsGeKermitBall->Write();// neutrons in ge, active and passive neutron detector
620 
621 // hAllParticlesGe->Write(); // all particles hitting the germanium
622 // hAllParticlesEkinGe->Write();// 2D all E_kin distribution of particles hitting the germanium
623 //
624 // hAllParticlesPiezo->Write(); // all particles hitting the first pieo
625 // hAllParticlesPiezo2->Write();// all particles hitting the germanium
626 //
627 //
628  hNEkinGe->Write(); // E_kin of neutrons entering the germanium
629  hNEkinKermit->Write(); // E_kin of neutrons entering the active neutron detector
630  hNEkinBall->Write(); // E_kin of neutrons entering the passive neutron detector
631 //
632 // hNEnergyLossGe->Write(); // Energy loss of neutrons in the germanium
633 // hNEnergyLossKermit->Write(); // Energy loss of neutrons in the active neutron detector
634 // hNEnergyLossBall->Write(); // Energy loss of neutrons in the passive neutron detector
635 //
636 // hPEkinGe->Write(); // E_kin of protons entering the germanium
637 // hPEnergyLossGe->Write(); // Energy loss of protons in the germanium
638 //
639 // hGammaSpecGe->Write(); // spectrum of gammas in the germanium
640 
641 
642  for(int i = 0;i<hAllParticlesGe->GetNbinsX(); i++)
643  {
644  if (hAllParticlesGe->GetBinContent(i) != 0)
645  {
646  cout << "PDG code:\t" << hAllParticlesGe->GetBinCenter(i)-0.5 << "\tCount:\t"<<hAllParticlesGe->GetBinContent(i) << endl;
647  }
648  }
649 
650  cout << "HypGe COSYBackgroundSim Ana:\tAnalysis finished succesfully" << endl;
651 
652 
653  OutputFile->Close();
654 
655 
656 
657  // ----- Finish -------------------------------------------------------
658  timer.Stop();
659  Double_t rtime = timer.RealTime();
660  Double_t ctime = timer.CpuTime();
661  cout << endl << endl;
662  cout << "Macro finished succesfully." << endl;
663 
664  cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
665  cout << endl;
666 
667  cout << "Hits\t" << Hits << endl;
668  // ------------------------------------------------------------------------
669  return iNeutrons;
670 }
i
Int_t i
Definition: run_full.C:25
b
TTree * b
Definition: GammaSpectraAnalysis.C:24
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
gGeoManager
TGeoManager * gGeoManager
Definition: createRootGeoFile.C:97
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PndHypGePoint::GetPy
Double_t GetPy() const
Definition: PndHypGePoint.h:58
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
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
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
PndHypGePoint::GetEnergyLoss
Double_t GetEnergyLoss() const
Definition: PndHypGePoint.h:62
outfile
TString outfile
Definition: createRootFscGeometryFile.C:57
PndHypGePoint::GetpdgCode
Double_t GetpdgCode() const
Definition: PndHypGePoint.h:63
Generated on Wed Apr 3 2019 10:03:14 for FairRoot/PandaRoot by   doxygen 1.8.5