Skip to content
Snippets Groups Projects
makeplot_input.C 3.53 KiB
Newer Older
  • Learn to ignore specific revisions
  • Jihee Kim's avatar
    Jihee Kim committed
    ////////////////////////////////////////
    // Read simulation ROOT input file
    // Write ASCII file
    ////////////////////////////////////////
    int makeplot_input(const char* input_fname = "sim_output/sim_emcal_electrons_input.root", const char *fout = "results/sim_1000k_info_1GeV.txt")
    {
      // Setting figures
      gROOT->SetStyle("Plain");
      gStyle->SetLineWidth(3);
      gStyle->SetOptStat("nem");
      gStyle->SetPadTickX(1);
      gStyle->SetPadTickY(1);
      gStyle->SetPadGridX(1);
      gStyle->SetPadGridY(1);
      gStyle->SetPadLeftMargin(0.14);
      gStyle->SetPadRightMargin(0.14);
    
      // Output ASCII file
      FILE *fp = fopen(fout,"a");
      if(!fp)
            fprintf(stderr,"error: can't open %s\n", fout);
    
      // Input ROOT file
      TFile *f = new TFile(input_fname,"READ");
      TTree *t = (TTree *)f->Get("events");
    
      // Set Branch status and addressed
      t->SetMakeClass(1);
      t->SetBranchStatus("*", 0);
    
      Int_t CrystalEcalHits_;
      t->SetBranchStatus("CrystalEcalHits", 1);
      t->SetBranchAddress("CrystalEcalHits", &CrystalEcalHits_);
    
      Int_t mcparticles_;
      t->SetBranchStatus("mcparticles", 1);
      t->SetBranchAddress("mcparticles", &mcparticles_);
    
      const Int_t kMaxCrystalEcalHits = 100000;
      Double_t truth_deposit[kMaxCrystalEcalHits];
      Double_t energyDeposit[kMaxCrystalEcalHits];
      t->SetBranchStatus("CrystalEcalHits.truth.deposit",1);
      t->SetBranchStatus("CrystalEcalHits.energyDeposit",1);
      t->SetBranchAddress("CrystalEcalHits.truth.deposit",truth_deposit);
      t->SetBranchAddress("CrystalEcalHits.energyDeposit",energyDeposit);
    
      const Int_t kMaxmcparticles = 100000;
      Double_t px[kMaxmcparticles];
      Double_t py[kMaxmcparticles];
      Double_t pz[kMaxmcparticles];
      Double_t mass[kMaxmcparticles];
      t->SetBranchStatus("mcparticles.psx",1);
      t->SetBranchStatus("mcparticles.psy",1);
      t->SetBranchStatus("mcparticles.psz",1);
      t->SetBranchStatus("mcparticles.mass",1);
      t->SetBranchAddress("mcparticles.psx",px);
      t->SetBranchAddress("mcparticles.psy",py);
      t->SetBranchAddress("mcparticles.psz",pz);
      t->SetBranchAddress("mcparticles.mass",mass);
    
      // Total number of entries
      Int_t nentries = t->GetEntries();
    
      // Variables are used in calculation
      Double_t total_thr_e;       // Thrown energy [GeV]
      Double_t truth_sim_e;       // truth simulated deposit energy [GeV]
      Double_t total_truth_sim_e; // Add up truth simulated deposit energy [GeV]
      Double_t sim_e;             // simulated deposit energy [GeV]
      Double_t total_sim_e;       // Add up truth simulated deposit energy [GeV]
      Double_t momentum2;         // momentum squared
      Double_t mass2;             // mass squared
    
      // Loop over event by event
      for (int ievent = 0; ievent < nentries; ievent++)
      {
    	t->GetEntry(ievent);
    
            Int_t nCrystalEcalHits = CrystalEcalHits_;
            Int_t nmcparticle      = 2;
    
            total_thr_e       = 0.0;
            total_truth_sim_e = 0.0;
            total_sim_e       = 0.0;
    
            momentum2 = px[nmcparticle]*px[nmcparticle]+py[nmcparticle]*py[nmcparticle]+pz[nmcparticle]*pz[nmcparticle];
            mass2     = mass[nmcparticle]*mass[nmcparticle];
    
            total_thr_e = sqrt(momentum2 + mass2)/1.e+3;
    
            // Loop over cluster by cluster
            for (int isimhit=0; isimhit < nCrystalEcalHits; isimhit++)
            {
                    truth_sim_e        = truth_deposit[isimhit];
                    sim_e              = energyDeposit[isimhit];
                    total_truth_sim_e += truth_sim_e;
                    total_sim_e       += energyDeposit[isimhit]/1.e+3;
            }
            fprintf(fp,"%d\t%2.4f\t%2.4f\t%2.4f\n",ievent,total_thr_e,total_truth_sim_e,total_sim_e);
      }
      fclose(fp);
      return 0;
    }