AllNeutronAnalysis_job_edit.C

Go to the documentation of this file.

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("pndsim") ;
  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: "<<nEvents<<endl;

        Int_t Hits=0;
        Int_t DetID=-10;
        Int_t TrackID=-10;
        int iNeutrons = 0;
        //event loop
        for (Int_t k=0; k<nEvents; k++)
        {
                b->GetEntry(k);
                if (!((k*100)% nEvents))        // mark every % of proceeded events
                {
                        cout << "Event number :\t" << k << endl;
                }
            //if(verbose) cout<<"Event No "<<j<<endl;

          DetID=-10;
          TrackID=-10;
                // loop over the hits of an event
                for (Int_t i=0; i<hit_bar->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()<<endl;
                                iNeutrons++;
                                hCrystalHit->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" <<mc_bar->GetEntriesFast() << endl;
                for (Int_t m = 0; m<mc_bar->GetEntriesFast();m++)
                {
                        //cout << "mcbar length :\t" <<mc_bar->GetEntriesFast() << endl;
                        PndMCTrack *mcgam = (PndMCTrack*)mc_bar->At(m);
                        TVector3 NAllMom = mcgam->GetMomentum();
                        if (mcgam->GetPdgCode() == 2112 && mcgam->GetMotherID() == -1)
                        {
                                cout << "event :\t"<<k<<"\tmcpdg :\t" << mcgam->GetPdgCode() << "\tmcmotherid :\t"<< mcgam->GetMotherID() <<endl;

                                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)
                        }
                }
        }// 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;
}