FairRoot/PandaRoot
PndHypGeCOSYBackgroundAna.cxx
Go to the documentation of this file.
1 /******************************************************
2 
3 Analysis Task created by M.Steinen steinen@kph.uni-mainz.de
4 Analysis of Gamma Simulation with hypGe detectors
5 *******************************************************/
6 
7 
8 #include "PndHypGeCOSYBackgroundAna.h"
9 #include <iostream>
10 #include <TMath.h>
11 #include "PndHypGePoint.h"
12 #include "PndMCTrack.h"
13 #include <time.h>
14 #include <TCanvas.h>
15 #include <TSpectrum.h>
16 #include <TStyle.h>
17 
18 //#include "RhoBase/TFactory.h"
19 
20 
21 using namespace std;
22 
23 
24 PndHypGeCOSYBackgroundAna::PndHypGeCOSYBackgroundAna()
25  : FairTask("Panda HypGeGammaAna Task")
26 {
27 
28 }
29 PndHypGeCOSYBackgroundAna::PndHypGeCOSYBackgroundAna(Int_t nEvents =0)
30  : FairTask("Panda HypGeGammaAna Task")
31 {
32  NumberOfEvents= nEvents;
33 }
34 
35 PndHypGeCOSYBackgroundAna::~PndHypGeCOSYBackgroundAna()
36 {
37  delete hNHits;
38  delete hCrystalHit;
39  delete hNeutronOrigin;
40  delete hCrystalOrigin;
41 }
42 
43 
44 InitStatus PndHypGeCOSYBackgroundAna::Init()
45 {
46  cout << "HypGe COSYBackgroundSim Ana:\t" << "Starting Init of HypGe Gamma Ana" << endl;
47 
48  //----------get FairRootManager-------------------------------------
49  FairRootManager* ioman = FairRootManager::Instance();
50  if ( ! ioman ) {
51  cout << "-E- PndEmcHitProducer::Init: "
52  << "RootManager not instantiated!" << endl;
53  return kFATAL;
54  }
55 
56  if (fVerbose == 3)
57  cout << "HypGe COSYBackgroundSim Ana:\t" << "Getting FairRootManager finished" << endl;
58 
59  //----------get input arrays---------------------------------------------
60  fMcTr = (TClonesArray*) ioman->GetObject("MCTrack");
61  //fHyp = (TClonesArray*) ioman->GetObject("HypPoint");
62  fHypGe = (TClonesArray*) ioman->GetObject("HypGePoint");
63 
64  if (fVerbose ==3)
65  {
66  cout << "HypGe COSYBackgroundSim Ana:\t" << "Getting input arrays finished" << endl;
67  }
68  //-----------extracting file name--------------------------------------
69  fName = ioman->GetInFile()->GetName();
70  if (fVerbose == 3)
71  {
72  cout << "HypGe COSYBackgroundSim Ana:\t" << "Full input file name with path: " << fName << endl;
73  cout << "HypGe COSYBackgroundSim Ana:\t" << "Chopping first " << fName.Last('/')+1 << " chars"<< endl;
74  }
75  fName.Remove(0,fName.Last('/')+1);
76  if (fVerbose == 3)
77  cout << "HypGe COSYBackgroundSim Ana:\t" << "File name: " <<fName<< endl;
78 
79  if (fVerbose == 3)
80  cout << "HypGe COSYBackgroundSim Ana:\t" << "File name extraction finished" << endl;
81 
82  //----------Get GeoManager---------------------------------------------
83 
84  fgeom = (TGeoManager*)gROOT->FindObject("FAIRGeom");
85  if (fVerbose == 3)
86  cout << "HypGe COSYBackgroundSim Ana:\t" << "Found GeoManager" << endl;
87  //----------create histograms------------------------------------------
88 
89  //histogram with polar angle of detector interactions of primary neutrons
90  hNHits = new TH1D("hNHits","Polar angle of primary neutrons interacting with the crystals",180, 90,180);
91  hNHits->SetXTitle("#Theta [#circ]");
92  hNHits->SetYTitle("Counts / 0.5 #circ");
93 
94  //histogram to see which detector is hit
95  hCrystalHit = new TH1D("hCrystalHit","Hits per Crystal",2,0,2);
96  hCrystalHit->SetXTitle("Crystal number");
97  hCrystalHit->SetYTitle("Counts per crystal");
98  hCrystalHit->GetXaxis()->SetBinLabel(1,"Crystal 30 cm"); //Crystal_101
99  hCrystalHit->GetXaxis()->SetBinLabel(2,"Crystal 90 cm"); //Crystal_202
100 
101  //histogram to see where neutron is coming from
102  hNeutronOrigin = new TH1D("hNeutronOrigin","Origin of Neutrons",5,0,5);
103  hNeutronOrigin->SetXTitle("Neutron Origin");
104  hNeutronOrigin->SetYTitle("Counts");
105  hNeutronOrigin->GetXaxis()->SetBinLabel(1,"Target");
106  hNeutronOrigin->GetXaxis()->SetBinLabel(2,"Beam dump");
107  hNeutronOrigin->GetXaxis()->SetBinLabel(3,"Cave");
108  hNeutronOrigin->GetXaxis()->SetBinLabel(4,"Capsule");
109  hNeutronOrigin->GetXaxis()->SetBinLabel(5,"Crystal");
110 
111  //histogram for correlation of neutron origin and crystal
112  hCrystalOrigin = new TH2D("hCrystalOrigin","Crystal -Neutron Origin - Correlation",2,0,2,5,0,5);
113  hCrystalOrigin->SetXTitle("Crystal number");
114  hCrystalOrigin->SetYTitle("Neutron Origin");
115  hCrystalOrigin->GetXaxis()->SetBinLabel(1,"Crystal 30 cm"); //Crystal_101
116  hCrystalOrigin->GetXaxis()->SetBinLabel(2,"Crystal 90 cm"); //Crystal_202
117  hCrystalOrigin->GetYaxis()->SetBinLabel(1,"Target");
118  hCrystalOrigin->GetYaxis()->SetBinLabel(2,"Beam dump");
119  hCrystalOrigin->GetYaxis()->SetBinLabel(3,"Cave");
120  hCrystalOrigin->GetYaxis()->SetBinLabel(4,"Capsule");
121  hCrystalOrigin->GetYaxis()->SetBinLabel(5,"Crystal");
122 
123  //histogram for kin. energy of neutrons
124  hNeutronEkin = new TH1D("hNeutronEkin","E_{kin} of neutrons;E_{kin} of neutrons [MeV]; Counts", 200,0,20);
125 
126  //histogram for correlation of E_kin and neutron origin
127  hNeutronEkinOrigin = new TH2D("hNeutronEkinOrigin","E_{kin} -Neutron Origin - Correlation;E_{kin} of neutrons [MeV]; Neutron Origin",200,0,20,5,0,5);
128 
129  hNeutronEkinOrigin->GetYaxis()->SetBinLabel(1,"Target");
130  hNeutronEkinOrigin->GetYaxis()->SetBinLabel(2,"Beam dump");
131  hNeutronEkinOrigin->GetYaxis()->SetBinLabel(3,"Cave");
132  hNeutronEkinOrigin->GetYaxis()->SetBinLabel(4,"Capsule");
133  hNeutronEkinOrigin->GetYaxis()->SetBinLabel(5,"Crystal");
134 
135  //histogram with energy deposited by neutrons in detector
136  hNeutronEnergyLoss = new TH1D("hNeutronEnergyLoss","Energy loss of neutrons inside a crystal;Energy loss of neutrons [MeV]; Counts", 2000,0,2);
137 
138 
139 
140  //histogram with all particles in both crystals created outside of the detector
141  hAllParticlesGe = new TH1D("hAllParticlesGe","Particles interaction in the crystals;PDG Code of particle; Counts", 2e4,-1e4,1e4);
142  //histogram with all particles in first crystal created outside of the detector
143  hAllParticlesCrystal1 = new TH1D("hAllParticlesCrystal1","Particles interaction in the first crystal (30 cm);PDG Code of particle; Counts", 5000,-2500,2500);
144  //histogram with all particles in second crystal created outside of the detector
145  hAllParticlesCrystal2 = new TH1D("hAllParticlesCrystal2","Particles interaction in the second crystal (90 cm);PDG Code of particle; Counts", 5000,-2500,2500);
146  //histogram with all particles in piezo created outside of the crystal
147  hAllParticlesPiezo = new TH1D("hAllParticlesPiezo","Particles interaction in the piezo;PDG Code of particle; Counts", 2e4,-1e4,1e4);
148  //histogram with all particles in first Sipm created outside of the detector
149  hAllParticlesSiPm1 = new TH1D("hAllParticlesSiPm1","Particles interaction in the first SiPm (60 cm);PDG Code of particle; Counts", 5000,-2500,2500);
150  //histogram with all particles in second Sipm created outside of the detector
151  hAllParticlesSiPm2 = new TH1D("hAllParticlesSiPm2","Particles interaction in the second SiPm (30 cm);PDG Code of particle; Counts", 5000,-2500,2500);
152 
153  // histogram for correlation of E_kin and particle PDG Code for all detectors
154  hEkinAllParticles = new TH2D("hEkinAllParticles","E_{kin} - Particle - Correlation;E_{kin} of particles [MeV]; particle",200,0,200,9,1,10);
155  hEkinAllParticles->GetYaxis()->SetBinLabel(1,"#pi^{-}");
156  hEkinAllParticles->GetYaxis()->SetBinLabel(2,"#mu^{+}");
157  hEkinAllParticles->GetYaxis()->SetBinLabel(3,"e^{+}");
158  hEkinAllParticles->GetYaxis()->SetBinLabel(4,"e^{-}");
159  hEkinAllParticles->GetYaxis()->SetBinLabel(5,"#mu^{-}");
160  hEkinAllParticles->GetYaxis()->SetBinLabel(6,"#gamma");
161  hEkinAllParticles->GetYaxis()->SetBinLabel(7,"#pi^{+}");
162  hEkinAllParticles->GetYaxis()->SetBinLabel(8,"neutron");
163  hEkinAllParticles->GetYaxis()->SetBinLabel(9,"proton");
164 
165  // histogram for correlation of E_kin and particle PDG Code for crystal1 (30 cm)
166  hEkinAllParticlesCrystal1 = new TH2D("hEkinAllParticlesCrystal1","E_{kin} - Particle - Correlation (crystal1);E_{kin} of particles [MeV]; particle",200,0,200,9,1,10);
167  hEkinAllParticlesCrystal1->GetYaxis()->SetBinLabel(1,"#pi^{-}");
168  hEkinAllParticlesCrystal1->GetYaxis()->SetBinLabel(2,"#mu^{+}");
169  hEkinAllParticlesCrystal1->GetYaxis()->SetBinLabel(3,"e^{+}");
170  hEkinAllParticlesCrystal1->GetYaxis()->SetBinLabel(4,"e^{-}");
171  hEkinAllParticlesCrystal1->GetYaxis()->SetBinLabel(5,"#mu^{-}");
172  hEkinAllParticlesCrystal1->GetYaxis()->SetBinLabel(6,"#gamma");
173  hEkinAllParticlesCrystal1->GetYaxis()->SetBinLabel(7,"#pi^{+}");
174  hEkinAllParticlesCrystal1->GetYaxis()->SetBinLabel(8,"neutron");
175  hEkinAllParticlesCrystal1->GetYaxis()->SetBinLabel(9,"proton");
176 
177  // histogram for correlation of E_kin and particle PDG Code for crystal2 (90 cm)
178  hEkinAllParticlesCrystal2 = new TH2D("hEkinAllParticlesCrystal2","E_{kin} - Particle - Correlation (Crystal2);E_{kin} of particles [MeV]; particle",200,0,200,9,1,10);
179  hEkinAllParticlesCrystal2->GetYaxis()->SetBinLabel(1,"#pi^{-}");
180  hEkinAllParticlesCrystal2->GetYaxis()->SetBinLabel(2,"#mu^{+}");
181  hEkinAllParticlesCrystal2->GetYaxis()->SetBinLabel(3,"e^{+}");
182  hEkinAllParticlesCrystal2->GetYaxis()->SetBinLabel(4,"e^{-}");
183  hEkinAllParticlesCrystal2->GetYaxis()->SetBinLabel(5,"#mu^{-}");
184  hEkinAllParticlesCrystal2->GetYaxis()->SetBinLabel(6,"#gamma");
185  hEkinAllParticlesCrystal2->GetYaxis()->SetBinLabel(7,"#pi^{+}");
186  hEkinAllParticlesCrystal2->GetYaxis()->SetBinLabel(8,"neutron");
187  hEkinAllParticlesCrystal2->GetYaxis()->SetBinLabel(9,"proton");
188 
189  // histogram for correlation of E_kin and particle PDG Code for piezo
190  hEkinAllParticlesPiezo = new TH2D("hEkinAllParticlesPiezo","E_{kin} - Particle - Correlation (Piezo);E_{kin} of particles [MeV]; particle",200,0,200,9,1,10);
191  hEkinAllParticlesPiezo->GetYaxis()->SetBinLabel(1,"#pi^{-}");
192  hEkinAllParticlesPiezo->GetYaxis()->SetBinLabel(2,"#mu^{+}");
193  hEkinAllParticlesPiezo->GetYaxis()->SetBinLabel(3,"e^{+}");
194  hEkinAllParticlesPiezo->GetYaxis()->SetBinLabel(4,"e^{-}");
195  hEkinAllParticlesPiezo->GetYaxis()->SetBinLabel(5,"#mu^{-}");
196  hEkinAllParticlesPiezo->GetYaxis()->SetBinLabel(6,"#gamma");
197  hEkinAllParticlesPiezo->GetYaxis()->SetBinLabel(7,"#pi^{+}");
198  hEkinAllParticlesPiezo->GetYaxis()->SetBinLabel(8,"neutron");
199  hEkinAllParticlesPiezo->GetYaxis()->SetBinLabel(9,"proton");
200 
201  // histogram for correlation of E_kin and particle PDG Code for SiPm1 (60 cm)
202  hEkinAllParticlesSiPm1 = new TH2D("hEkinAllParticlesSiPm1","E_{kin} - Particle - Correlation (SiPm1);E_{kin} of particles [MeV]; particle",200,0,200,9,1,10);
203  hEkinAllParticlesSiPm1->GetYaxis()->SetBinLabel(1,"#pi^{-}");
204  hEkinAllParticlesSiPm1->GetYaxis()->SetBinLabel(2,"#mu^{+}");
205  hEkinAllParticlesSiPm1->GetYaxis()->SetBinLabel(3,"e^{+}");
206  hEkinAllParticlesSiPm1->GetYaxis()->SetBinLabel(4,"e^{-}");
207  hEkinAllParticlesSiPm1->GetYaxis()->SetBinLabel(5,"#mu^{-}");
208  hEkinAllParticlesSiPm1->GetYaxis()->SetBinLabel(6,"#gamma");
209  hEkinAllParticlesSiPm1->GetYaxis()->SetBinLabel(7,"#pi^{+}");
210  hEkinAllParticlesSiPm1->GetYaxis()->SetBinLabel(8,"neutron");
211  hEkinAllParticlesSiPm1->GetYaxis()->SetBinLabel(9,"proton");
212 
213  // histogram for correlation of E_kin and particle PDG Code for SiPm2 (30 cm)
214  hEkinAllParticlesSiPm2 = new TH2D("hEkinAllParticlesSiPm2","E_{kin} - Particle - Correlation (SiPm2);E_{kin} of particles [MeV]; particle",200,0,200,9,1,10);
215  hEkinAllParticlesSiPm2->GetYaxis()->SetBinLabel(1,"#pi^{-}");
216  hEkinAllParticlesSiPm2->GetYaxis()->SetBinLabel(2,"#mu^{+}");
217  hEkinAllParticlesSiPm2->GetYaxis()->SetBinLabel(3,"e^{+}");
218  hEkinAllParticlesSiPm2->GetYaxis()->SetBinLabel(4,"e^{-}");
219  hEkinAllParticlesSiPm2->GetYaxis()->SetBinLabel(5,"#mu^{-}");
220  hEkinAllParticlesSiPm2->GetYaxis()->SetBinLabel(6,"#gamma");
221  hEkinAllParticlesSiPm2->GetYaxis()->SetBinLabel(7,"#pi^{+}");
222  hEkinAllParticlesSiPm2->GetYaxis()->SetBinLabel(8,"neutron");
223  hEkinAllParticlesSiPm2->GetYaxis()->SetBinLabel(9,"proton");
224 
225  //histogram for correlation of E_kin and gamma origin
226  hGammaEkinOrigin = new TH2D("hGammaEkinOrigin","E_{kin} -#gamma Origin - Correlation;E_{kin} of neutrons [MeV]; Gamma Origin",200,0,20,5,0,5);
227  hGammaEkinOrigin->GetYaxis()->SetBinLabel(1,"Target");
228  hGammaEkinOrigin->GetYaxis()->SetBinLabel(2,"Beam dump");
229  hGammaEkinOrigin->GetYaxis()->SetBinLabel(3,"Cave");
230  hGammaEkinOrigin->GetYaxis()->SetBinLabel(4,"Capsule");
231  hGammaEkinOrigin->GetYaxis()->SetBinLabel(5,"Crystal");
232 
234 
235  if (fVerbose == 3)
236  cout << "HypGe COSYBackgroundSim Ana:\t" << "Histogram creation finished" << endl;
237 
238  //----------initialize EvtCount----------------------------------------
239  EvtCount = 0;
240 
241  //----------Get number of events from FairRootManager-----------------------
242 
243  if (NumberOfEvents== 0)
244  NumberOfEvents= ioman->GetInChain()->GetEntries();
245  cout << "HypGe COSYBackgroundSim Ana:\tNumber of Events to analyze: "<< NumberOfEvents << endl;
246 
247 
248  cout << "HypGe COSYBackgroundSim Ana:\t" << "Init of HypGe COSYBackgroundSim Ana finished succesfully" << endl;
249 }
250 
251 void PndHypGeCOSYBackgroundAna::Exec(Option_t*)
252 {
253  Double_t NeutronEnergyLossArray[10];
254  for(Int_t i = 0; i < 10; i++)
255  NeutronEnergyLossArray[i]=0;
256  Int_t nNeutrons = -1;
257  Int_t NeutronTrackID =-1;
258  for (Int_t i=0; i < fHypGe->GetEntriesFast(); i++)
259  {
260  PndHypGePoint *hitgam=(PndHypGePoint*)fHypGe->At(i);
261  PndMCTrack *mcgam = (PndMCTrack*)fMcTr->At(hitgam->GetTrackID());
262  //only Neutrons
263  StartVertex = mcgam->GetStartVertex();
264  VertexVolumeName =fgeom->FindNode(StartVertex.X(),StartVertex.Y(),StartVertex.Z())->GetName();
265  if (hitgam->GetpdgCode() == 2112)
266  {
267 
268  //cout <<"Event " <<EvtCount<< " \t\tStartVertex Radius: " <<StartVertex.Mag()<<endl;
269 
270 
271 
272 
273  //Fill hNeutronOrigin
274 
275  //cout << VertexVolumeName << endl;
276 
277  if(hitgam->GetDetectorID()==101)
278  {
279  hCrystalHit->Fill("Crystal 30 cm",1);
280  if( VertexVolumeName.Contains("Target"))
281  {
282  hNeutronOrigin->Fill("Target",1);
283  hCrystalOrigin->Fill("Crystal 30 cm","Target",1);
284  }
285  else
286  {
287  if (VertexVolumeName.Contains("BeamDump"))
288  {
289  hNeutronOrigin->Fill("Beam dump",1);
290  hCrystalOrigin->Fill("Crystal 30 cm","Beam dump",1);
291  }
292  else
293  {
294  if (VertexVolumeName.Contains("cave"))
295  {
296  hNeutronOrigin->Fill("Cave",1);
297  hCrystalOrigin->Fill("Crystal 30 cm", "Cave",1);
298  }
299  else
300  {
301  if (VertexVolumeName.Contains ("Capsule"))
302  {
303  hNeutronOrigin->Fill("Capsule",1);
304  hCrystalOrigin->Fill("Crystal 30 cm","Capsule",1);
305  }
306  else
307  {
308  if (VertexVolumeName.Contains("Crystal"))
309  {
310  hNeutronOrigin->Fill("Crystal",1);
311  hCrystalOrigin->Fill("Crystal 30 cm","Crystal",1);
312 
313  }
314  }
315  }
316  }
317  }
318  }
319  else
320  if(hitgam->GetDetectorID()==202)
321  {
322  hCrystalHit->Fill("Crystal 90 cm",1);
323  if( VertexVolumeName.Contains("Target"))
324  {
325  hNeutronOrigin->Fill("Target",1);
326  hCrystalOrigin->Fill("Crystal 90 cm","Target",1);
327  }
328  else
329  {
330  if (VertexVolumeName.Contains("BeamDump"))
331  {
332  hNeutronOrigin->Fill("Beam dump",1);
333  hCrystalOrigin->Fill("Crystal 90 cm","Beam dump",1);
334  }
335  else
336  {
337  if (VertexVolumeName.Contains("cave"))
338  {
339  hNeutronOrigin->Fill("Cave",1);
340  hCrystalOrigin->Fill("Crystal 90 cm", "Cave",1);
341  }
342  else
343  {
344  if (VertexVolumeName.Contains ("Capsule"))
345  {
346  hNeutronOrigin->Fill("Capsule",1);
347  hCrystalOrigin->Fill("Crystal 90 cm","Capsule",1);
348  }
349  else
350  {
351  if (VertexVolumeName.Contains("Crystal"))
352  {
353  hNeutronOrigin->Fill("Crystal",1);
354  hCrystalOrigin->Fill("Crystal 90 cm","Crystal",1);
355 
356  }
357  }
358  }
359  }
360  }
361  }
362 
363  //E_kin stuff
364  NeutronMomentum.SetX(hitgam->GetPx());
365  NeutronMomentum.SetY(hitgam->GetPy());
366  NeutronMomentum.SetZ(hitgam->GetPz());
367  NeutronEkin = 1./2.*NeutronMomentum.Mag2()/0.939565*1000;
368  //cout << NeutronEkin << "\t\t\t" << NeutronMomentum.Mag2()<<endl;
369  if (!VertexVolumeName.Contains("Crystal") && ! VertexVolumeName.Contains ("Capsule")) //created outside of the detector
370  {
371  hNeutronEkin->Fill(NeutronEkin);
372  if (NeutronTrackID !=hitgam->GetTrackID()) // new neutron (different TrackID)
373  {
374  NeutronTrackID = hitgam->GetTrackID();
375  nNeutrons++;
376  }
377  NeutronEnergyLossArray[nNeutrons] += hitgam->GetEnergyLoss()*1000;
378  //cout << hitgam->GetEnergyLoss()*1000<< endl;
379  }
380  //fill E_kin - Neutron origin correlation
381  if( VertexVolumeName.Contains("Target"))
382  {
383  hNeutronEkinOrigin->Fill(NeutronEkin,"Target",1);
384  }
385  else
386  {
387  if (VertexVolumeName.Contains("BeamDump"))
388  {
389  hNeutronEkinOrigin->Fill(NeutronEkin,"Beam dump",1);
390  }
391  else
392  {
393  if (VertexVolumeName.Contains("cave"))
394  {
395  hNeutronEkinOrigin->Fill(NeutronEkin, "Cave",1);
396  }
397  else
398  {
399  if (VertexVolumeName.Contains ("Capsule"))
400  {
401  hNeutronEkinOrigin->Fill(NeutronEkin,"Capsule",1);
402  }
403  else
404  {
405  if (VertexVolumeName.Contains("Crystal"))
406  {
407  hNeutronEkinOrigin->Fill(NeutronEkin,"Crystal",1);
408 
409  }
410  }
411  }
412  }
413  }
414  } // end of only neutrons
415 
416  //all particles in Crystal
417  if (!VertexVolumeName.Contains("Crystal") && ! VertexVolumeName.Contains ("Capsule")) //created outside of the detector
418  {
419  if(hitgam->GetDetectorID()==101 || hitgam->GetDetectorID()==102)
420  hAllParticlesGe->Fill(hitgam->GetpdgCode());
421  if(hitgam->GetDetectorID()==101)
422  hAllParticlesCrystal1->Fill(hitgam->GetpdgCode());
423  if(hitgam->GetDetectorID()==202)
424  hAllParticlesCrystal2->Fill(hitgam->GetpdgCode());
425  if(hitgam->GetDetectorID()==1000)
426  hAllParticlesPiezo ->Fill(hitgam->GetpdgCode());
427  if(hitgam->GetDetectorID()==1001)
428  hAllParticlesSiPm1->Fill(hitgam->GetpdgCode());
429  if(hitgam->GetDetectorID()==1002)
430  hAllParticlesSiPm2 ->Fill(hitgam->GetpdgCode());
431  ParticleMomentum.SetX(hitgam->GetPx());
432  ParticleMomentum.SetY(hitgam->GetPy());
433  ParticleMomentum.SetZ(hitgam->GetPz());
434  Double_t ParticleMass = 0;
435  TString ParticleName = "";
436  switch ((int) hitgam->GetpdgCode())
437  {
438  case -211: ParticleMass = 0.13957; ParticleName ="#pi^{-}"; break; //pi-
439  case -13: ParticleMass = 0.10561; ParticleName ="#mu^{+}"; break; //mu+
440  case -11: ParticleMass = 0.000511; ParticleName ="e^{+}"; break; //e+
441  case 11: ParticleMass = 0.000511; ParticleName ="e^{-}"; break; //e-
442  case 13: ParticleMass = 0.10561; ParticleName ="#mu^{-}"; break; //mu-
443  case 22: ParticleMass = 0; ParticleName ="#gamma"; break; //gamma
444  case 211: ParticleMass = 0.13957; ParticleName ="#pi^{+}"; break; //pi+
445  case 2112: ParticleMass = 0.939565; ParticleName ="neutron"; break; //n
446  case 2212: ParticleMass = 0.93827; ParticleName ="proton"; break; //p
447  default: ParticleMass = -1; std::cout.precision(10); cout << "Found new particle:\t" << hitgam->GetpdgCode() <<"\t" << VertexVolumeName << endl; break;
448  }
449  if(ParticleMass != -1)
450  {
451  ParticleEkin = (sqrt(ParticleMass*ParticleMass+ParticleMomentum.Mag2())-ParticleMass)*1000;//Ekin in MeV
452  //all particles (but gammas) in Crystal - E_{kin} - correlation
453  if (!VertexVolumeName.Contains("Crystal") && ! VertexVolumeName.Contains ("Capsule")) //created outside of the detector
454  {
455  hEkinAllParticles->Fill( ParticleEkin , ParticleName.Data(),1);
456  if(hitgam->GetDetectorID()==101)
457  hEkinAllParticlesCrystal1->Fill(ParticleEkin , ParticleName.Data(),1);
458  if(hitgam->GetDetectorID()==202)
459  hEkinAllParticlesCrystal2->Fill(ParticleEkin , ParticleName.Data(),1);
460  if(hitgam->GetDetectorID()==1000)
461  hEkinAllParticlesPiezo ->Fill(ParticleEkin , ParticleName.Data(),1);
462  if(hitgam->GetDetectorID()==1001)
463  hEkinAllParticlesSiPm1->Fill(ParticleEkin , ParticleName.Data(),1);
464  if(hitgam->GetDetectorID()==1002)
465  hEkinAllParticlesSiPm2 ->Fill(ParticleEkin , ParticleName.Data(),1);
466  }
467  }
468 
469  //gammas only
470  if (hitgam->GetpdgCode() == 22)
471  {
472  if( VertexVolumeName.Contains("Target"))
473  {
474  hGammaEkinOrigin->Fill(ParticleMomentum.Mag(),"Target",1);
475  }
476  else
477  {
478  if (VertexVolumeName.Contains("BeamDump"))
479  {
480  hGammaEkinOrigin->Fill(ParticleMomentum.Mag(),"Beam dump",1);
481  }
482  else
483  {
484  if (VertexVolumeName.Contains("cave"))
485  {
486  hGammaEkinOrigin->Fill(ParticleMomentum.Mag(), "Cave",1);
487  }
488  else
489  {
490  if (VertexVolumeName.Contains ("Capsule"))
491  {
492  hGammaEkinOrigin->Fill(ParticleMomentum.Mag(),"Capsule",1);
493  }
494  else
495  {
496  if (VertexVolumeName.Contains("Crystal"))
497  {
498  hGammaEkinOrigin->Fill(ParticleMomentum.Mag(),"Crystal",1);
499 
500  }
501  }
502  }
503  }
504  }
505  }
506  }
507 
508  } // end of for loop
509 
510  //fill Neutron Engergy loss histogram
511  if (nNeutrons>-1)
512  cout << nNeutrons+1 << endl;
513  for (Int_t i = 0; i <= nNeutrons; i++)
514  {
515  cout << i << "\t\t" << NeutronEnergyLossArray[i] << endl;
516  hNeutronEnergyLoss->Fill(NeutronEnergyLossArray[i]);
517  }
518 
519  EvtCount++;
520  if (!((EvtCount*100)% NumberOfEvents)) // mark every % of proceeded events
521  {
522  cout << "HypGe COSYBackgroundSim Ana:\t" << EvtCount*100/NumberOfEvents << " % done"<< endl;
523  }
524  //cout << "HypGe COSYBackgroundSim Ana:\t" << "Event Number"<< EvtCount << endl;
525 }
526 
527 void PndHypGeCOSYBackgroundAna::Finish()
528 {
529  cNHits = new TCanvas("cNHits","Neutron Hits",800,600);
530  cCrystalHit = new TCanvas("cCrystalHit","Crystal hit", 800,600);
531  cNeutronOrigin = new TCanvas("cNeutronOrigin","Origin of Neutrons", 800,600);
532  cCrystalOrigin = new TCanvas("cCrystalOrigin","Correlation of Crystal hit and neutron origin", 800, 600);
533  cNeutronEkin = new TCanvas("cNeutronEkin","E_{kin} of neutrons", 800, 600);
534  cNeutronEkinOrigin = new TCanvas("cNeutronEkinOrigin","E_{kin} - Neutron Origin - Correlation", 800, 600);
535  cNeutronEnergyLoss = new TCanvas("cNeutronEnergyLossn","Energy loss of neutrons inside a crystal", 800, 600);
536 
537  cAllParticlesGe = new TCanvas("cAllParticles","All particles interacting in the crystal", 800, 600);
538  cAllParticlesCrystal1 = new TCanvas("cAllParticlesCrystal1","All particles interacting in the first crystal", 800, 600);
539  cAllParticlesCrystal2 = new TCanvas("cAllParticlesCrystal2","All particles interacting in the second crystal", 800, 600);
540  cEkinAllParticles = new TCanvas("cEkinAllParticles","Correlation of E_{kin} of Particles and PDG code", 800, 600);
541  cGammaEkinOrigin= new TCanvas("cGammaEkinOrigin","E_{kin} - #gamma Origin - Correlation", 800, 600);
542 
543  cAllParticlesPiezo = new TCanvas("cAllParticlesPiazo","All particles interacting in the piezo", 800, 600);
544  cAllParticlesSiPm1 = new TCanvas("cAllParticlesSiPm1","All particles interacting in the SiPm1", 800, 600);
545  cAllParticlesSiPm2 = new TCanvas("cAllParticlesSiPm2","All particles interacting in the SiPm2", 800, 600);
546  cEkinAllParticlesSplit= new TCanvas("cEkinAllParticlesSplit","Correlation of E_{kin} of Particles and PDG code", 800, 600);
547  cEkinAllParticlesSplit->Divide(3,2);
548 
549  hCrystalHit->Write();
550  hNeutronOrigin->Write();
551  hCrystalOrigin->Write();
552  hNeutronEkin->Write();
553  hNeutronEkinOrigin->Write();
554  hAllParticlesGe->Write();
555  hAllParticlesCrystal1->Write();
556  hAllParticlesCrystal2->Write();
557  hEkinAllParticles->Write();
558  hGammaEkinOrigin->Write();
559  hNeutronEnergyLoss->Write();
560 
561  hAllParticlesPiezo->Write();
562  hAllParticlesSiPm1->Write();
563  hAllParticlesSiPm2->Write();
564 
565  hEkinAllParticlesCrystal1->Write();
566  hEkinAllParticlesCrystal2->Write();
567  hEkinAllParticlesPiezo ->Write();
568  hEkinAllParticlesSiPm1->Write();
569  hEkinAllParticlesSiPm2 ->Write();
570 
571 
572 
573  cCrystalHit->cd();
574  hCrystalHit->Draw();
575 
576  cNeutronOrigin->cd();
577  hNeutronOrigin->Draw();
578 
579  cCrystalOrigin->cd();
580  gStyle->SetPalette(1);
581  hCrystalOrigin->Draw("colz");
582 
583  cNeutronEkin->cd();
584  hNeutronEkin->Draw();
585 
586  cNeutronEkinOrigin->cd();
587  hNeutronEkinOrigin->Draw("LEGO2");
588 
589  cNeutronEnergyLoss->cd();
590  hNeutronEnergyLoss->Draw();
591 
592  cAllParticlesGe->cd();
593  hAllParticlesGe->Draw();
594 
595  cAllParticlesCrystal1->cd();
596  hAllParticlesCrystal1->Draw();
597 
598  cAllParticlesCrystal2->cd();
599  hAllParticlesCrystal2->Draw();
600 
601  cEkinAllParticles->cd();
602  hEkinAllParticles->Draw("colz");
603 
604  cGammaEkinOrigin->cd();
605  hGammaEkinOrigin->Draw("colz");
606 
607  cAllParticlesPiezo ->cd();
608  hAllParticlesPiezo->Draw();
609  cAllParticlesSiPm1->cd();
610  hAllParticlesSiPm1->Draw();
611  cAllParticlesSiPm2 ->cd();
612  hAllParticlesSiPm2 ->Draw();
613 
614  cEkinAllParticlesSplit->cd(1);
615  hEkinAllParticlesCrystal1->Draw("colz");
616  cEkinAllParticlesSplit->cd(2);
617  hEkinAllParticlesCrystal2->Draw("colz");
618  cEkinAllParticlesSplit->cd(3);
619  hEkinAllParticlesPiezo ->Draw("colz");
620  cEkinAllParticlesSplit->cd(4);
621  hEkinAllParticlesSiPm1->Draw("colz");
622  cEkinAllParticlesSplit->cd(5);
623  hEkinAllParticlesSiPm2 ->Draw("colz");
624 
625 
626  for(int i = 0;i<hAllParticlesGe->GetNbinsX(); i++)
627  {
628  if (hAllParticlesGe->GetBinContent(i) != 0)
629  {
630  cout << "PDG code:\t" << hAllParticlesGe->GetBinCenter(i)-0.5 << "\tCount:\t"<<hAllParticlesGe->GetBinContent(i) << endl;
631  }
632  }
633 
634  cout << "HypGe COSYBackgroundSim Ana:\tAnalysis finished succesfully" << endl;
635 }
636 
637 ClassImp(PndHypGeCOSYBackgroundAna)
PndHypGeCOSYBackgroundAna::Exec
virtual void Exec(Option_t *opt)
Definition: PndHypGeCOSYBackgroundAna.cxx:251
fVerbose
int fVerbose
Definition: poormantracks.C:24
i
Int_t i
Definition: run_full.C:25
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
PndHypGePoint::GetPz
Double_t GetPz() const
Definition: PndHypGePoint.h:59
PndHypGePoint::GetTrackID
Int_t GetTrackID() const
Definition: PndHypGePoint.h:51
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
PndHypGePoint::GetPx
Double_t GetPx() const
Definition: PndHypGePoint.h:57
Fill
HISThit_ene Fill(sum_hit_ene)
ClassImp
ClassImp(PndAnaContFact)
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
PndHypGePoint::GetDetectorID
Int_t GetDetectorID() const
Definition: PndHypGePoint.h:53
PndHypGePoint::GetEnergyLoss
Double_t GetEnergyLoss() const
Definition: PndHypGePoint.h:62
PndHypGePoint::GetpdgCode
Double_t GetpdgCode() const
Definition: PndHypGePoint.h:63
Generated on Sun Aug 18 2019 03:01:32 for FairRoot/PandaRoot by   doxygen 1.8.8