int AllNeutronAnalysis_job_edit(TString Filename_ext) { bool verbose = false; // ----- Load libraries ------------------------------------------------ gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C"); gSystem->Load("libHypGe"); // ----- Timer -------------------------------------------------------- TStopwatch timer; timer.Start(); // ------------------------------------------------------------------------ //Files //TString Path = getenv("SIMDATADIR"); TString Filename = "$SIMDATADIR/Neutron/"+Filename_ext; //TString Filename = "$SIMDATADIR/TripleBall40Offset10_urqmd_pbarC_1_5000Evts.root"; TString outfile= "$SIMDATADIR/Neutron/Ana/Ana" + Filename_ext; TFile* InputFile = new TFile(Filename); TFile* OutputFile = new TFile(outfile,"RECREATE"); //getting simulation branches from input file TTree *b=(TTree *) InputFile->Get("cbmsim") ; TClonesArray* hit_bar=new TClonesArray("PndHypGePoint"); b->SetBranchAddress("HypGePoint",&hit_bar);//Branch names TClonesArray* mc_bar=new TClonesArray("PndMCTrack"); b->SetBranchAddress("MCTrack",&mc_bar);//Branch names //declaring histograms //histogram with polar angle of detector interactions of primary neutrons TH1D* hNHits = new TH1D("hNHits","Polar angle of primary neutrons interacting with the crystals",180, 90,180); hNHits->SetXTitle("#Theta [#circ]"); hNHits->SetYTitle("Counts / 0.5 #circ"); //histogram with the polar momentum of primary neutrons interacting with the crystals TH1D* hNHits_p = new TH1D("hNHits_p","Direction of polar momentum of primary neutrons interacting with the crystals",180, 90,180); hNHits_p->SetXTitle("#Theta [#circ]"); hNHits_p->SetYTitle("Counts / 0.5 #circ"); //histogram with the polar momentum of primary neutrons in the solid angle of the crystal (not necessarily interacting with the crystal) TH1D* hNAll = new TH1D("hNAll","All Neutrons",180, 90,180); hNAll->SetXTitle("#Theta [#circ]"); hNAll->SetYTitle("Counts / 0.5 #circ"); //histogram with E_kin of primary neutrons TH1D* hNMom = new TH1D("hNMom","E_{kin} of neutrons",1000, 0,0.5); hNMom->SetXTitle("E_{kin} of neutrons [GeV]"); hNMom->SetYTitle("Counts / 0.5 MeV"); //histogram to see which detector is hit TH1D* hCrystalHit = new TH1D("hCrystalHit","Hits per Crystal",1700,1,1700); hCrystalHit->SetXTitle("Crystal number"); hCrystalHit->SetYTitle("Counts per crystal"); //histogram to see the energy deposited by the neutron TH1D* hNeutronEnergyDeposit = new TH1D("hNeutronEnergyDeposit","Energy deposited by Neutrons per Hit",500,0,0.005); hNeutronEnergyDeposit->SetXTitle("Energy [GeV]"); hNeutronEnergyDeposit->SetYTitle("Counts / 10 keV"); Int_t *ActualTrackID ; Int_t nEvents = b->GetEntriesFast(); cout<< "Number of Simulated Events: "<GetEntry(k); if (!((k*100)% nEvents)) // mark every % of proceeded events { cout << "Event number :\t" << k << endl; } //if(verbose) cout<<"Event No "<GetEntriesFast(); i++) { //if (i==0) //Hits++; //cout <<"Hit "<< i << endl; PndHypGePoint *hitgam=(PndHypGePoint*)hit_bar->At(i); PndMCTrack *mcgam = (PndMCTrack*)mc_bar->At(hitgam->GetTrackID()); TVector3 NMom = mcgam->GetMomentum(); if ( hitgam->GetpdgCode()==2112 && mcgam->GetMotherID() == -1 && !(DetID == hitgam->GetDetectorID() && TrackID==hitgam->GetTrackID())) // primary neutron, last part: not(same paricle and crystal) -> against multi counting of a single neutron { DetID = hitgam->GetDetectorID(); TrackID = hitgam->GetTrackID(); cout << "Event :\t" << k <<"\tHit : \t" << i << "\tDetector :\t" << DetID << "\tTrackID :\t"<< TrackID<<"\tMother :\t"<< mcgam->GetMotherID()<Fill(hitgam->GetDetectorID()); //fill histogram to see which detector is hit TVector3 NHit; // fill vector with hit coordinates to get theta of the hit NHit.SetX(hitgam->GetX()); NHit.SetY(hitgam->GetY()); NHit.SetZ(hitgam->GetZ()+55); hNHits->Fill(180/TMath::Pi()*NHit.Theta()); //fill histogram with polar angle of detector interactions of primary neutrons TVector3 NHit_p; // fill vector with hit momentum to get theta of the hit momentum NHit_p.SetX(hitgam->GetPx()); NHit_p.SetY(hitgam->GetPy()); NHit_p.SetZ(hitgam->GetPz()); hNHits_p->Fill(180/TMath::Pi()*NHit_p.Theta()); //fill histogram with the polar momentum of primary neutrons interacting with the crystals hNMom->Fill(1./2.* NMom.Mag2()/0.939); //fill histogram with E_kin of primary neutrons } if (hitgam->GetpdgCode()==2112) // neutron { hNeutronEnergyDeposit->Fill(hitgam->GetEnergyLoss()); } } // loop over all entries in mctrack cout << "mcbar length :\t" <GetEntriesFast() << endl; for (Int_t m = 0; mGetEntriesFast();m++) { //cout << "mcbar length :\t" <GetEntriesFast() << endl; PndMCTrack *mcgam = (PndMCTrack*)mc_bar->At(m); TVector3 NAllMom = mcgam->GetMomentum(); if (mcgam->GetPdgCode() == 2112 && mcgam->GetMotherID() == -1) { cout << "event :\t"<GetPdgCode() << "\tmcmotherid :\t"<< mcgam->GetMotherID() <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) } } }// end for k (events) //hNAll->DrawCopy(); //hNHits->SetLineColor(kRed); //hNHits->DrawCopy("same"); //Analysis of spectrum hNHits->Write(); hNHits_p->Write(); hNAll->Write(); hNMom->Write(); hCrystalHit->Write(); hNeutronEnergyDeposit->Write(); //hNeutronEnergyDeposit->DrawCopy(); OutputFile->Close(); //hNeutronEnergyDeposit->Draw(); //hNHits->Draw(); // ----- Finish ------------------------------------------------------- timer.Stop(); Double_t rtime = timer.RealTime(); Double_t ctime = timer.CpuTime(); cout << endl << endl; cout << "Macro finished succesfully." << endl; cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl; cout << endl; cout << "Hits\t" << Hits << endl; // ------------------------------------------------------------------------ return iNeutrons; }