source/DEMPgen/ThetaCheck.C

0001 // Stephen JD Kay, University of York, 05/03/25
0002 // A short script to read in DEMPgen output and find the scattered electron theta ranges that (roughly) correspond to fixed Q2 ranges.
0003 // Run this on the .root output file DEMPgen produces for a file with a fully open e' theta range (i.e. the full range in generated Q2 covered).
0004 // It will identify the points in theta which correspond to Q2 ranges of 3-10, 10-20 and 20-35. Customise as needed.
0005 #include <string>
0006
0007 void ThetaCheck(TString infile=""){
0008
0009   // If no input file provide as argument, promot for one
0010   if(infile == ""){
0011     cout << "Enter a filename to analyse: ";
0012     cin >> infile;
0013   }
0014
0015   // Check input file exists, exit if not
0016   if(gSystem->AccessPathName(infile) == kTRUE){
0017     cerr << "!!!!! ERROR !!!!!" << endl << infile << " not found" << endl << "!!!!! ERROR !!!!!" << endl;
0018     exit(0);
0019   }
0020
0021   TString ofile_name, ofile_tmp, BeamE;
0022   TObjArray *tmp_Name_arr;
0023
0024   // Check the input file is a .root file as we would expect
0025   if(infile.Contains(".root") == false){
0026     cerr << "!!!!!!!!!!! - ERROR - !!!!!!!!!!!!" << endl;
0027     cerr << "Input files should be a root file!" << endl;
0028     cerr << "!!!!!!!!!!! - ERROR - !!!!!!!!!!!!" << endl;
0029     exit(1);
0030   }
0031   else{
0032     cout << "Opening and finding Q2 boundaries " << infile << endl;
0033   }
0034
0035   // If the input file name is a file path containing /, extract only the actual file name for further use. Assign the temporary name of the output to be the input, minus its .root extension
0036   if(infile.Contains("/")){
0037     tmp_Name_arr = infile.Tokenize("/");
0038     ofile_tmp = (((TObjString *)(tmp_Name_arr->At(tmp_Name_arr->GetLast())))->String()).ReplaceAll(".root","");
0039   }
0040   else{
0041     ofile_tmp = infile;
0042     ofile_tmp.ReplaceAll(".root","");
0043   }
0044
0045   TChain *mychain = new TChain("Events");
0046   mychain->Add(infile);
0047
0048   // Initialise reader
0049   TTreeReader tree_reader(mychain);
0050
0051   // Get particle info
0052   TTreeReaderArray<double> ScatElec_pX(tree_reader, "scat_e_px");
0053   TTreeReaderArray<double> ScatElec_pY(tree_reader, "scat_e_py");
0054   TTreeReaderArray<double> ScatElec_pZ(tree_reader, "scat_e_pz");
0055   //TTreeReaderArray<double> ScatElec_E(tree_reader, "scat_e_E");
0056   TTreeReaderArray<double> KinQ2(tree_reader, "Q2");
0057
0058   double Theta, Q2;
0059   TVector3 ScatElec;
0060
0061   // Set output file name
0062   ofile_name = Form("%s_ThetaQ2.root", ofile_tmp.Data());
0063   TH1D* H1_Theta = new TH1D("H1_Theta", "#theta_{e'}; #theta_{e'} (Deg)", 240, 120, 180);
0064   TH2D* H2_Q2_Theta = new TH2D("H2_Q2_Theta", "#theta_{e'}(Q^{2});Q^{2} (GeV^{2}); #theta_{e'} (Deg)", 400, 0, 40, 240, 120, 180);
0065   TH2D* H2_Q2_Theta_R1 = new TH2D("H2_Q2_Theta_R1", "#theta_{e'}(Q^{2}), 3 < Q^{2} < 10;Q^{2} (GeV^{2}); #theta_{e'} (Deg)", 400, 0, 40, 240, 120, 180);
0066   TH2D* H2_Q2_Theta_R2 = new TH2D("H2_Q2_Theta_R2", "#theta_{e'}(Q^{2}) 10 < Q^{2} < 20;Q^{2} (GeV^{2}); #theta_{e'} (Deg)", 400, 0, 40, 240, 120, 180);
0067   TH2D* H2_Q2_Theta_R3 = new TH2D("H2_Q2_Theta_R3", "#theta_{e'}(Q^{2}) 20 < Q^{2} < 35;Q^{2} (GeV^{2}); #theta_{e'} (Deg)", 400, 0, 40, 240, 120, 180);
0068   // Set and open output file for the histograms
0069   TFile *ofile = TFile::Open(ofile_name,"RECREATE");
0070
0071   while (tree_reader.Next()){
0072     ScatElec.SetXYZ(ScatElec_pX[0], ScatElec_pY[0], ScatElec_pZ[0]);
0073     Theta=ScatElec.Theta()*TMath::RadToDeg();
0074     Q2 = KinQ2[0];
0075     H1_Theta->Fill(Theta);
0076     H2_Q2_Theta->Fill(Q2, Theta);
0077     if (Q2 > 3 && Q2 < 10){
0078       H2_Q2_Theta_R1->Fill(Q2, Theta);
0079     }
0080     else if (Q2 > 10 && Q2 < 20){
0081       H2_Q2_Theta_R2->Fill(Q2, Theta);
0082     }
0083     if (Q2 > 20 && Q2 < 35){
0084       H2_Q2_Theta_R3->Fill(Q2, Theta);
0085     }
0086
0087   }
0088
0089   // Write histograms
0090   H1_Theta->Write();
0091   H2_Q2_Theta->Write();
0092   H2_Q2_Theta_R1->Write();
0093   H2_Q2_Theta_R2->Write();
0094   H2_Q2_Theta_R3->Write();
0095
0096   // Print to screen the angular ranges corresponding to the theta ranges in use
0097   cout << "Min elec theta = " << H2_Q2_Theta_R1->GetYaxis()->GetBinLowEdge(H2_Q2_Theta_R1->FindFirstBinAbove(0., 2)) << " deg, max elec theta = " << H2_Q2_Theta_R1->GetYaxis()->GetBinLowEdge(H2_Q2_Theta_R1->FindLastBinAbove(0., 2)) << " deg for 3 < Q2 < 10" << endl;
0098   cout << "Min elec theta = " << H2_Q2_Theta_R2->GetYaxis()->GetBinLowEdge(H2_Q2_Theta_R2->FindFirstBinAbove(0., 2)) << " deg, max elec theta = " << H2_Q2_Theta_R2->GetYaxis()->GetBinLowEdge(H2_Q2_Theta_R2->FindLastBinAbove(0., 2)) << " deg for 10 < Q2 < 20" << endl;
0099   cout << "Min elec theta = " << H2_Q2_Theta_R3->GetYaxis()->GetBinLowEdge(H2_Q2_Theta_R3->FindFirstBinAbove(0., 2)) << " deg, max elec theta = " << H2_Q2_Theta_R3->GetYaxis()->GetBinLowEdge(H2_Q2_Theta_R3->FindLastBinAbove(0., 2)) << " deg for 20 < Q2 < 35" << endl;
0100
0101   cout << H2_Q2_Theta_R1->GetEntries()/H2_Q2_Theta_R2->GetEntries() << " times more events in range 3 < Q2 < 10 than 10 < Q2 < 20" << endl;
0102   cout << H2_Q2_Theta_R1->GetEntries()/H2_Q2_Theta_R3->GetEntries() << " times more events in range 3 < Q2 < 10 than 20 < Q2 < 35" << endl;
0103
0104   ofile->Write();
0105   ofile->Close();
0106   delete ofile;
0107
0108 }