EICMCAnalysis/analysis_eicsmear/hera_dis_fillbins.C

0001 int
0002 hera_dis_fillbins()
0003 {
0004   gSystem->Load("libeicsmear");
0005
0006   /* Open input files. */
0007   TFile *file_mc = new TFile("/gpfs/mnt/gpfs04/sphenix/user/nfeege/data/TEST/pythia.ep.27.5x820.1Mevents.RadCor=0.Q2gt1.root");
0008
0009   /* Get trees from files. */
0010   TTree *tree = (TTree*)file_mc->Get("EICTree");
0011
0012   /* Output file. */
0013   TFile *fout = new TFile("hera_dis_histo.root", "RECREATE");
0014
0015   erhic::EventPythia *event  = NULL;
0016
0017   tree->SetBranchAddress("event", &event);
0018
0019   /* Create n-dimensional histogram to collect statistic for SIDIS events in kinematics bins:
0020    *
0021    * Event properties:
0022    * - x
0023    * - Q2
0024    */
0025   const int nbins_x  = 10;
0026   const int nbins_Q2 = 10;
0027
0028   const int hn_dis_ndim = 2;
0029   int hn_dis_nbins[] = { nbins_x,
0030                          nbins_Q2 };
0031   double hn_dis_xmin[] = {0., 0. };
0032   double hn_dis_xmax[] = {0., 0. };
0033
0034   THnSparse* hn_dis = new THnSparseF("hn_dis",
0035                                      "DIS Kinematis; x; Q2;",
0036                                      hn_dis_ndim,
0037                                      hn_dis_nbins,
0038                                      hn_dis_xmin,
0039                                      hn_dis_xmax );
0040
0041   hn_dis->GetAxis(0)->SetName("x");
0042   hn_dis->GetAxis(1)->SetName("Q2");
0043
0044   double min_x = -5;
0045   double max_x = 0;
0046   double width_x = (max_x - min_x) / nbins_x;
0047   double bins_x[nbins_x+1];
0048   for( int i =0; i <= nbins_x; i++)
0049     {
0050       cout << TMath::Power( 10, min_x + i*width_x) << endl;
0051       bins_x[i] = TMath::Power( 10, min_x + i*width_x);
0052     }
0053
0054   double min_Q2 = 0;
0055   double max_Q2 = 5;
0056   double width_Q2 = (max_Q2 - min_Q2) / nbins_Q2;
0057   double bins_Q2[nbins_Q2+1];
0058   for( int i =0; i <= nbins_Q2; i++)
0059     bins_Q2[i] = TMath::Power( 10, min_Q2 + i*width_Q2);
0060
0061   hn_dis->SetBinEdges(0,bins_x);
0062   hn_dis->SetBinEdges(1,bins_Q2);
0063
0064   /* print all bin centers */
0065   TH2F* hxQ2 = (TH2F*)hn_dis->Projection(1,0);
0066   for ( int bin_x = 1; bin_x <= hxQ2->GetNbinsX(); bin_x++ )
0067     {
0068
0069       for ( int bin_y = 1; bin_y <= hxQ2->GetNbinsY(); bin_y++ )
0070         {
0071       cout << "x = " << TMath::Power(10, 0.5 * ( ( TMath::Log10( hxQ2->GetXaxis()->GetBinLowEdge(bin_x) ) )
0072                              + ( TMath::Log10( hxQ2->GetXaxis()->GetBinLowEdge(bin_x) + hxQ2->GetXaxis()->GetBinWidth(bin_x) ) ) ) ) << endl;
0073
0074       cout << "Q2 = " << TMath::Power(10, 0.5 * ( ( TMath::Log10( hxQ2->GetYaxis()->GetBinLowEdge(bin_y) ) )
0075                              + ( TMath::Log10( hxQ2->GetYaxis()->GetBinLowEdge(bin_y) + hxQ2->GetYaxis()->GetBinWidth(bin_y) ) ) ) ) << endl;
0076     }
0077     }
0078
0079   /* Loop over all events in tree. */
0080   for ( unsigned ievent = 0; ievent < tree->GetEntries(); ievent++ )
0081     {
0082       if ( ievent%1000 == 0 )
0083         cout << "Processing event " << ievent << endl;
0084
0085       /* load event */
0086       tree->GetEntry(ievent);
0087
0088       /* Cut on scattered lepton properties */
0089       float Emin = 5;
0090       if ( event->ScatteredLepton()->GetE() < Emin )
0091     continue;
0092
0093       if ( event->ScatteredLepton()->GetTheta() >= ( 177./ 180. * TMath::Pi() ) )
0094     continue;
0095
0096       /* Cut on kinematics */
0097       float x = event->GetTrueX();
0098       float Q2 = event->GetTrueQ2();
0099       float y = event->GetTrueY();
0100
0101       if ( Q2 < 1 )
0102     continue;
0103
0104       if ( y >= 0.9 || y <= 0.01 )
0105     continue;
0106
0107
0108       double fill_hn_dis[] = {x, Q2};
0109       hn_dis->Fill( fill_hn_dis );
0110     }
0111
0112   /* Write histogram. */
0113   hn_dis->Write();
0114
0115   /* Close output file. */
0116   fout->Close();
0117
0118   return 0;
0119 }