EICMCAnalysis/analysis_eicsmear/eic_sphenix_fill_xq2.C

0001 /* STL includes */
0002 #include <iostream>
0003
0004 /* ROOT includes */
0005 #include <TROOT.h>
0006 //#include <gSystem.h>
0007 #include <TFile.h>
0008 #include <TKey.h>
0009 #include <TTree.h>
0010 #include <TH2F.h>
0011 #include <THnSparse.h>
0012
0013 /* eicsmear includes */
0014 #include <eicsmear/smear/EventSmear.h>
0015 #include <eicsmear/erhic/EventPythia.h>
0016 #include <eicsmear/erhic/ParticleMC.h>
0017
0018 using namespace std;
0019
0020 int
0021 eic_sphenix_fill_xq2( TString filename_output,
0022               TString filename_mc,
0023               TString filename_mc_smeared,
0024               bool debug = false )
0025 {
0026   /* Uncomment this line when running without compilation */
0027   //  gSystem->Load("libeicsmear");
0028
0029   /* Selection cuts for scattered electron (after smearing) */
0030   /* Minimum momentum:  this makes sure the momentum was smeared, i.e. is within tracking acceptance */
0031   double electron_pmin = 0.01; // GeV
0032
0033   /* Minimumenergy: Acceptance cut for good electron/pion separation */
0034   double electron_emin = 3.0; // GeV
0035
0036   /* Open input files. */
0037   TFile *file_mc = new TFile(filename_mc, "OPEN");
0038   TFile *file_mc_smeared = new TFile(filename_mc_smeared, "OPEN");
0039
0040   /* Get trees from files. */
0041   TTree *tree = (TTree*)file_mc->Get("EICTree");
0042   TTree *tree_smeared = (TTree*)file_mc_smeared->Get("Smeared");
0043
0044
0045   /* Get event generator parameters (cross section, number of events, ...) from file. */
0046   /* Input file may include multiple cross sections because it was merged from multiple files, so loop over all crossSection
0047      objects and calculate average; same for sum of nEvents. */
0048   double xsec = 0;
0049   unsigned nxsec = 0;
0050   int nevent = 0;
0051   TIter next(file_mc->GetListOfKeys());
0052   TKey *key;
0053   while ((key = (TKey*)next()))
0054     {
0055       if ( TString(key->GetName()) == "crossSection" )
0056     {
0057       xsec += (TString(((TObjString*)key->ReadObj())->GetName())).Atof();
0058       nxsec++;
0059       continue;
0060     }
0061
0062       if ( TString(key->GetName()) == "nEvents" )
0063     {
0064       nevent += (TString(((TObjString*)key->ReadObj())->GetName())).Atof();
0065       continue;
0066     }
0067     }
0068
0069   /* Get average cross section per file */
0070   if ( nxsec > 0 )
0071     xsec /= nxsec;
0072
0073   TObjString* gen_crossSection = (TObjString*)file_mc->Get("crossSection");
0074   TObjString* gen_nEvents = (TObjString*)file_mc->Get("nEvents");
0075
0076   /* Update strings with total values */
0077   gen_crossSection->SetString( TString::Format("%f", xsec) );
0078   gen_nEvents->SetString( TString::Format("%d", nevent) );
0079
0080   /* Convert cross section from microbarn (10^-6) to femtobarn (10^-15) */
0081   float xsec_fb = xsec * 1e9;
0082
0083   /* Calculate integrated luminosity (inverse femtobarn) */
0084   float lumi_ifb = (float)nevent / (float)xsec_fb;
0085
0086   /* Add Luminosity strings */
0087   TObjString* gen_lumi_ifb = new TObjString();
0088   gen_lumi_ifb->SetString( TString::Format("%f", lumi_ifb) );
0089
0090   /* Print parameters */
0091   cout << "crossSection: " << gen_crossSection->GetName() << endl;
0092   cout << "nEvents: " << gen_nEvents->GetName() << endl;
0093   cout << "luminosity: " << gen_lumi_ifb->GetName() << endl;
0094
0095   /* Output file. */
0096   TFile *file_out = new TFile(filename_output, "RECREATE");
0097
0098   /* Add friend to match branches in trees. */
0099   tree->AddFriend(tree_smeared);
0100
0101   erhic::EventPythia *event  = NULL;
0102   Smear::Event       *eventS = NULL;
0103
0104   tree->SetBranchAddress("event", &event);
0105   tree->SetBranchAddress("eventS", &eventS);
0106
0107   /* Create histogram for eta of particles to check acceptance of detectors. */
0108   TH1F* h_eta = new TH1F("h_eta", ";#eta;dN/d#eta", 100, -5, 5);
0109   TH1F* h_eta_accept = (TH1F*)h_eta->Clone("h_eta_accept");
0110   h_eta_accept->SetLineColor(kGreen);
0111
0112   /* Create n-dimensional histogram to collect statistic for DIS and SIDIS events in kinematics bins:
0113    *
0114    * Event properties:
0115    * - x
0116    * - Q2
0117    *
0118    * Particle properties:
0119    * - z
0120    * - pT
0121    */
0122
0123   /* Define bin ranges for DIS and SIDIS histograms */
0124
0125   /* x */
0126   const int nbins_x  = 25; // 5 per decade
0127
0128   double min_x = -5;
0129   double max_x = 0;
0130   double width_x = (max_x - min_x) / nbins_x;
0131   double bins_x[nbins_x+1];
0132   for( int i =0; i <= nbins_x; i++)
0133     {
0134       //cout << TMath::Power( 10, min_x + i*width_x) << endl;
0135       bins_x[i] = TMath::Power( 10, min_x + i*width_x);
0136     }
0137
0138   /* Q2 */
0139   const int nbins_Q2 = 20; // 4 per decade
0140
0141   double min_Q2 = -1;
0142   double max_Q2 = 4;
0143   double width_Q2 = (max_Q2 - min_Q2) / nbins_Q2;
0144   double bins_Q2[nbins_Q2+1];
0145   for( int i =0; i <= nbins_Q2; i++)
0146     bins_Q2[i] = TMath::Power( 10, min_Q2 + i*width_Q2);
0147
0148   /* z */
0149   const int nbins_z  = 20;
0150
0151   double min_z = 0.0;
0152   double max_z = 1.0;
0153   double width_z = (max_z - min_z) / nbins_z;
0154   double bins_z[nbins_z+1];
0155   for( int i =0; i <= nbins_z; i++)
0156     {
0157       //cout << min_z + i*width_z << endl;
0158       bins_z[i] = min_z + i*width_z;
0159     }
0160
0161   /* pT */
0162   const int nbins_pT1 = 5;
0163   const int nbins_pT2 = 4;
0164   const int nbins_pT = nbins_pT1 + nbins_pT2;
0165
0166   double min_pT = 0.0;
0167   double width_pT1 = 0.2;
0168   double width_pT2 = 1.0;
0169
0170   double bins_pT[nbins_pT+1];
0171   for( int i =0; i <= nbins_pT1; i++)
0172     {
0173       //cout << min_pT + i*width_pT1 << endl;
0174       bins_pT[i] = min_pT + i*width_pT1;
0175     }
0176   for( int i =1; i <= nbins_pT2; i++)
0177     {
0178       //cout << min_pT + nbins_pT1*width_pT1 + i*width_pT2<< endl;
0179       bins_pT[nbins_pT1 + i] = min_pT + nbins_pT1*width_pT1 + i*width_pT2;
0180     }
0181
0182
0183   /* Create DIS histogram- one entry per event */
0184   const int hn_dis_ndim = 2;
0185   int hn_dis_nbins[] = { nbins_x,
0186                          nbins_Q2 };
0187   double hn_dis_xmin[] = {0., 0. };
0188   double hn_dis_xmax[] = {0., 0. };
0189
0190   THnSparse* hn_dis = new THnSparseF("hn_dis_electron",
0191                                      "DIS Kinematis Per Event (Electron); x; Q2;",
0192                                      hn_dis_ndim,
0193                                      hn_dis_nbins,
0194                                      hn_dis_xmin,
0195                                      hn_dis_xmax );
0196
0197   hn_dis->GetAxis(0)->SetName("x");
0198   hn_dis->GetAxis(1)->SetName("Q2");
0199
0200   hn_dis->SetBinEdges(0,bins_x);
0201   hn_dis->SetBinEdges(1,bins_Q2);
0202
0203   /* clone histogram for ACCEPTED events */
0204   THnSparse* hn_dis_accept = (THnSparse*)hn_dis->Clone("hn_dis_electron_accept");
0205   hn_dis_accept->SetTitle("DIS Kinematis Per Event (Accepted)");
0206
0207   /* Create SIDIS histogram- one entry per particle */
0208   const int hn_sidis_ndim = 4;
0209   int hn_sidis_nbins[] = { nbins_x,
0210                            nbins_Q2,
0211                            nbins_z,
0212                nbins_pT };
0213   double hn_sidis_xmin[] = {0., 0., 0., 0. };
0214   double hn_sidis_xmax[] = {0., 0., 0., 0. };
0215
0216   THnSparse* hn_sidis_pion_plus = new THnSparseF("hn_sidis_pion_plus",
0217                                                  "SIDIS Kinematis Per Particle (Postitively Charged Pion);x;Q2;z;pT;",
0218                                                  hn_sidis_ndim,
0219                                                  hn_sidis_nbins,
0220                                                  hn_sidis_xmin,
0221                                                  hn_sidis_xmax );
0222
0223   hn_sidis_pion_plus->GetAxis(0)->SetName("x");
0224   hn_sidis_pion_plus->GetAxis(1)->SetName("Q2");
0225   hn_sidis_pion_plus->GetAxis(2)->SetName("z");
0226   hn_sidis_pion_plus->GetAxis(3)->SetName("pT");
0227
0228   /* Set logarithmic bin ranges. */
0229   hn_sidis_pion_plus->SetBinEdges(0,bins_x);
0230   hn_sidis_pion_plus->SetBinEdges(1,bins_Q2);
0231   hn_sidis_pion_plus->SetBinEdges(2,bins_z);
0232   hn_sidis_pion_plus->SetBinEdges(3,bins_pT);
0233
0234   /* clone SIDIS histogram for other particle identities */
0235   THnSparse* hn_sidis_pion_minus = (THnSparseF*)hn_sidis_pion_plus->Clone("hn_sidis_pion_minus");
0236   hn_sidis_pion_minus->SetTitle("SIDIS Kinematis Per Particle (Negatively Charged Pion)");
0237
0238   THnSparse* hn_sidis_kaon_plus  = (THnSparseF*)hn_sidis_pion_plus->Clone("hn_sidis_kaon_plus");
0239   hn_sidis_kaon_plus->SetTitle("SIDIS Kinematis Per Particle (Positively Charged Kaon)");
0240
0241   THnSparse* hn_sidis_kaon_minus = (THnSparseF*)hn_sidis_pion_plus->Clone("hn_sidis_kaon_minus");
0242   hn_sidis_kaon_minus->SetTitle("SIDIS Kinematis Per Particle (Negatively Charged Kaon)");
0243
0244   /* clone histogram for ACCEPTED events */
0245   THnSparse* hn_sidis_pion_plus_accept  = (THnSparse*)hn_sidis_pion_plus->Clone("hn_sidis_pion_plus_accept");
0246   hn_sidis_pion_plus_accept->SetTitle("SIDIS Kinematis Per Particle (Positively Charged Pion, Accepted)");
0247
0248   THnSparse* hn_sidis_pion_minus_accept = (THnSparseF*)hn_sidis_pion_plus->Clone("hn_sidis_pion_minus_accept");
0249   hn_sidis_pion_minus_accept->SetTitle("SIDIS Kinematis Per Particle (Negatively Charged Pion, Accepted)");
0250
0251   THnSparse* hn_sidis_kaon_plus_accept  = (THnSparseF*)hn_sidis_pion_plus->Clone("hn_sidis_kaon_plus_accept");
0252   hn_sidis_kaon_plus_accept->SetTitle("SIDIS Kinematis Per Particle (Positively Charged Kaon, Accepted)");
0253
0254   THnSparse* hn_sidis_kaon_minus_accept = (THnSparseF*)hn_sidis_pion_plus->Clone("hn_sidis_kaon_minus_accept");
0255   hn_sidis_kaon_minus_accept->SetTitle("SIDIS Kinematis Per Particle (Negatively Charged Kaon, Accepted)");
0256
0257   /* print all bin centers for x bins */
0258   TH2F* hxQ2 = (TH2F*)hn_dis->Projection(1,0);
0259
0260   /* Two alternative to determine bin center- which one is the CORRECT on a log scale axis?? */
0261   //for ( int bin_x = 1; bin_x <= hxQ2->GetNbinsX(); bin_x++ )
0262   //  {
0263   //    /* Option 1 matches eRHIC high energy document Fig 11 */
0264   //    fprintf( stdout, "xC = %.2e\n",
0265   //           hxQ2->GetXaxis()->GetBinCenter(bin_x) );
0266   //  }
0267
0268   /* Option 2 matches HERA table */
0269   if ( debug )
0270     {
0271       for ( int bin_x = 1; bin_x <= hxQ2->GetNbinsX(); bin_x++ )
0272     {
0273       fprintf( stdout, "x = %.2e\n",
0274            TMath::Power(10, 0.5 * ( ( TMath::Log10( hxQ2->GetXaxis()->GetBinLowEdge(bin_x) ) )
0275                         + ( TMath::Log10( hxQ2->GetXaxis()->GetBinLowEdge(bin_x) + hxQ2->GetXaxis()->GetBinWidth(bin_x) ) ) ) ) );
0276     }
0277     }
0278
0279   /* print bin centers for Q2 bins */
0280   if ( debug )
0281     {
0282       for ( int bin_Q2 = 1; bin_Q2 <= hxQ2->GetNbinsY(); bin_Q2++ )
0283     {
0284
0285       fprintf( stdout, "Q2 = %.2e\n",
0286            TMath::Power(10, 0.5 * ( ( TMath::Log10( hxQ2->GetYaxis()->GetBinLowEdge(bin_Q2) ) )
0287                         + ( TMath::Log10( hxQ2->GetYaxis()->GetBinLowEdge(bin_Q2) + hxQ2->GetYaxis()->GetBinWidth(bin_Q2) ) ) ) ) );
0288     }
0289     }
0290
0291   /* Loop over all events in tree. */
0292   unsigned max_event = tree->GetEntries();
0293
0294   if ( debug )
0295     max_event = 10000;
0296
0297   for ( unsigned ievent = 0; ievent < max_event; ievent++ )
0298     {
0299       if ( ievent%1000 == 0 )
0300         cout << "Processing event " << ievent << endl;
0301
0302       /* load event */
0303       tree->GetEntry(ievent);
0304
0305       /* Cut on EVENT kinematics */
0306       float y = event->GetTrueY();
0307       if ( y > 0.95 || y < 0.01 )
0308         continue;
0309
0310       float x = event->GetTrueX();
0311       float Q2 = event->GetTrueQ2();
0312
0313       double fill_hn_dis[] = {x, Q2};
0314       hn_dis->Fill( fill_hn_dis );
0315
0316       /* Scattered lepton found within acceptance of traacking system (= valid smeared momentum P)? */
0317       bool electron_detected = false;
0318       if ( eventS->ScatteredLepton() && eventS->ScatteredLepton()->GetP() > electron_pmin && eventS->ScatteredLepton()->GetE() > electron_emin )
0319     electron_detected = true;
0320
0321       /* Continue if electron not detected */
0322       if ( !electron_detected )
0323     continue;
0324
0325       /* execute the part below only if scattered electron within acceptance */
0326       hn_dis_accept->Fill( fill_hn_dis );
0327
0328
0329       /* For SIDIS: Loop over all final state particles in this event */
0330       unsigned ntracks = event->GetNTracks();
0331
0332       for ( unsigned itrack = 0; itrack < ntracks; itrack++ )
0333         {
0334           erhic::ParticleMC  * iparticle = event->GetTrack( itrack );
0335           Smear::ParticleMCS * ismeared  = eventS->GetTrack( itrack );
0336
0337           /* Check status */
0338           if ( iparticle->GetStatus() != 1 )
0339             continue;
0340
0341           /* Get z of particle */
0342           float z = iparticle->GetZ();
0343
0344           /* skip low z paricles- soft physics and beam remnants */
0345           if ( z <= 0.15 )
0346             continue;
0347
0348       /* skip particles outside +/- 4 pseudorapidity */
0349       if ( iparticle->GetEta() < -4 || iparticle->GetEta() > 4 )
0350         continue;
0351
0352       /* Get pT of particle w.r.t. exchange boson of interaction */
0353       float pT = iparticle->GetPtVsGamma();
0354
0355           /* Prepare array to fill histogram */
0356           double fill_hn_sidis[] = {x, Q2, z, pT};
0357
0358           /* get TRUE pid */
0359           int pid = iparticle->Id().Code();
0360           //TParticlePDG *pid_pdg =  iparticle->Id().Info();
0361
0362           /* Use true PID to choose which histogram to fill */
0363           /* Pi+ */
0364           if ( pid == 211 )
0365         {
0366           hn_sidis_pion_plus->Fill( fill_hn_sidis );
0367
0368           if ( ismeared->GetP() > 0.1 && ismeared->Id().Code() == 211 )
0369         hn_sidis_pion_plus_accept->Fill( fill_hn_sidis );
0370         }
0371
0372           /* Pi - */
0373           else if ( pid == -211 )
0374         {
0375           hn_sidis_pion_minus->Fill( fill_hn_sidis );
0376
0377           if ( ismeared->GetP() > 0.1 && ismeared->Id().Code() == -211 )
0378         hn_sidis_pion_minus_accept->Fill( fill_hn_sidis );
0379         }
0380
0381           /* K+ */
0382           else if ( pid == 321 )
0383         {
0384           hn_sidis_kaon_plus->Fill( fill_hn_sidis );
0385
0386           if ( ismeared->GetP() > 0.1 && ismeared->Id().Code() == 321 )
0387         hn_sidis_kaon_plus_accept->Fill( fill_hn_sidis );
0388         }
0389
0390           /* K- */
0391           else if ( pid == -321 )
0392         {
0393           hn_sidis_kaon_minus->Fill( fill_hn_sidis );
0394
0395           if ( ismeared->GetP() > 0.1 && ismeared->Id().Code() == -321 )
0396         hn_sidis_kaon_minus_accept->Fill( fill_hn_sidis );
0397         }
0398     } // end loop over particles
0399
0400     } // end loop over events
0401
0402   /* Write histogram. */
0403   hn_dis->Write();
0404
0405   hn_dis_accept->Write();
0406
0407   hn_sidis_pion_plus->Write();
0408   hn_sidis_pion_minus->Write();
0409   hn_sidis_kaon_plus->Write();
0410   hn_sidis_kaon_minus->Write();
0411
0412   hn_sidis_pion_plus_accept->Write();
0413   hn_sidis_pion_minus_accept->Write();
0414   hn_sidis_kaon_plus_accept->Write();
0415   hn_sidis_kaon_minus_accept->Write();
0416
0417   h_eta->Write();
0418   h_eta_accept->Write();
0419
0420   gen_crossSection->Write("crossSection");
0421   gen_nEvents->Write("nEvents");
0422   //gen_nTrials->Write("nTrials");
0423   gen_lumi_ifb->Write("luminosity");
0424
0425   /* Close output file. */
0426   file_out->Close();
0427
0428   return 0;
0429 }
0430
0431
0432 /* MAIN function */
0433 int main( int argc , char* argv[] )
0434 {
0435   if ( argc < 4 || argc > 5 )
0436     {
0437       cout << "Usage: " << argv[0] << " <filename_output> <filename_EICTree> <filename_EICTree_smeared> <optional: 'debug' for debug mode" << endl;
0438       return 1;
0439     }
0440
0441   cout << "Running eic_sphenix_fill_xq2 with: \n" << endl;
0442   cout << " - Output file:            " << argv[1] << endl;
0443   cout << " - EICTree input file:     " << argv[2] << endl;
0444
0445   if ( argc == 4 )
0446     {
0447       cout << " - EICTree (smeared) file: " << argv[3] << endl;
0448       eic_sphenix_fill_xq2( argv[1], argv[2], argv[3] );
0449     }
0450   else if ( argc == 5 )
0451     {
0452       cout << " - EICTree (smeared) file: " << argv[3] << endl;
0453       cout << " ==== DEBUG MODE ==== " << endl;
0454       eic_sphenix_fill_xq2( argv[1], argv[2], argv[3], true );
0455     }
0456
0457   return 0;
0458 }