<file_sep>{


  #include <cstdlib>
//   #include <iostream>
  #include <fstream>
  #include <stdio.h>
  #include <string.h>
  #include <vector>
  #include <math.h>
  #include <TH1F.h>

  using namespace std;

  gROOT->Reset();
  gROOT->SetStyle("Plain");

  // DEFINING FILE AND TREE VARIABLES

  const int NENERGIES= 1;
  const int nATT = 10;
  
  TFile * RunFile[NENERGIES];
  

  RunFile[0] = new TFile("data/1GeV_test_Attenuation_curves.root","READ");   //
//   RunFile[1] = new TFile("data/50GeV_full_sim_1mind.root","READ");   //
//   RunFile[2] = new TFile("data/50GeV_full_sim_2mind.root","READ");   //
//   RunFile[3] = new TFile("data/50GeV_full_sim_3mind.root","READ");   //
//   RunFile[4] = new TFile("data/50GeV_full_sim_5mind.root","READ");   //  
//   RunFile[5] = new TFile("data/50GeV_full_sim_7mind.root","READ");   //
  
  double input_energy[NENERGIES];
  input_energy[0] = 1;
//   input_energy[1] = 20;
//   input_energy[2] = 50;
//   input_energy[3] = 100;
//   input_energy[4] = 120;
//   input_energy[5] = 150;
  
  double input_att[nATT];
  input_att[0] = 0.117681;
  input_att[1] = 0.0637243;
  input_att[2] = 0.0439325;
  input_att[3] = 0.0321204;
  input_att[4] = 0.0192074;
  input_att[5] = 0.0124216;
  input_att[6] = 0.00660631;
  input_att[7] = 0.00431397;
  input_att[8] = 0.00229227;
  input_att[9] = 0.000864042;
  
//   input_att[10] = 0.001;
  
  cout << "almeno fin qui ..." << endl;
  
  TGraphErrors * grResolution = new TGraphErrors ();
  TGraphErrors * grLinearity = new TGraphErrors ();
  
  TGraphErrors * grResolutionAtt [nATT];
  TGraphErrors * grLinearityAtt [nATT];
  
  for (int iAtt = 0; iAtt < nATT; iAtt++) {
    grResolutionAtt[iAtt] = new TGraphErrors ();
    grLinearityAtt[iAtt] = new TGraphErrors ();	
    
  }
  
  cout << "almeno fin qui ..." << endl;
  
  TH1F * hDepoEn[NENERGIES];
  TH1F * hAttEn[NENERGIES][nATT];
  TH1F * hAttCalib_En[NENERGIES][nATT];
  
  TH1F * hInputAbsorption[NENERGIES];
  
  int NBINS = 1000;
  char name[50];
  
  for (int iEn = 0; iEn < NENERGIES; iEn ++) {  
    
    sprintf (name, "hDepoEn_%d", iEn);
    hDepoEn [iEn] = new TH1F (name, name, NBINS, 0.3*input_energy[iEn], input_energy[iEn]);
    
    for (int iAtt = 0; iAtt < nATT; iAtt++){
      sprintf (name, "hAttEn_%d_ATT_%d", iEn, iAtt);
//       hAttEn[iEn][iAtt] = new TH1F (name, name, NBINS, 0.3*input_energy[iEn]*input_att[iAtt], input_energy[iEn]*input_att[iAtt]);
      hAttEn[iEn][iAtt] = new TH1F (name, name, NBINS, 0, input_energy[iEn]*input_att[iAtt]*2);
      
      sprintf (name, "hAttCalib_En_%d_ATT_%d", iEn, iAtt);
      hAttCalib_En[iEn][iAtt] = new TH1F (name, name, NBINS, 0.3*input_energy[iEn], input_energy[iEn]*2);
    }
  }
  
  /// defining tree variables
  
    Float_t InitialPositionX;
    Float_t InitialPositionY;
    Float_t InitialPositionZ;
  
    Float_t Total_energy ;
    Float_t Total_ion_energy ;
    Float_t Total_ion_energy_att[nATT] ;


  ///******************************************///
  /// Run over different beam position
  ///******************************************///

  for (int iEn = 0; iEn < NENERGIES; iEn++) {


    TTree* TreeRun = (TTree*) RunFile[iEn].Get("tree");
    hInputAbsorption[iEn] = (TH1F*) RunFile[iEn].Get("hAttenuation");

    TreeRun->SetBranchAddress("InitialPositionX",&InitialPositionX);
    TreeRun->SetBranchAddress("InitialPositionY",&InitialPositionY);
    TreeRun->SetBranchAddress("InitialPositionZ",&InitialPositionZ);

    TreeRun->SetBranchAddress("Total_energy",&Total_energy);
    TreeRun->SetBranchAddress("Total_ion_energy",&Total_ion_energy);
    TreeRun->SetBranchAddress("Total_ion_energy_att",&Total_ion_energy_att);


      ///******************************************///
      ///		 Run over events	    ///
      ///******************************************///

      int NEVENTS = TreeRun->GetEntries();
      cout << "iEn = " << iEn << " :: nEvents = " << NEVENTS << endl;
      
      for (Int_t iEvt= 0; iEvt < NEVENTS; iEvt++) {//&& iEvt < 5

	TreeRun->GetEntry(iEvt);
	cout << " iEvt = " << iEvt << endl;
	hDepoEn[iEn]->Fill(Total_ion_energy);
	
	for (int iAtt = 0; iAtt < nATT; iAtt++){
	  cout << "input_att[" << iAtt << "] = " << Total_ion_energy_att[iAtt]/0.84 <<";" << endl;
	  hAttEn[iEn][iAtt]->Fill(Total_ion_energy_att[iAtt]);
	}
	
      }	   //end of events loop
      
      
/*
    TF1 * fEnergies = new TF1 ("fEnergies", " [0]*exp((x-220)/[1]) + [2]", 0, 220);
    
    fEnergies->SetParameters(0.1, 20, 0.01);
//     fEnergies->FixParameter(1,220);
    fEnergies->SetLineColor(iEn+1);
    for (int j=0; j< 10; j++) hEfficiency[iEn]->Fit("fEnergies", "QR");
    
    
    TF1 * fDecayTime = new TF1 ("fDecayTime", " [0]*exp(-x/[1])", 0, 100);
    fDecayTime->SetParameters(100, 20);
    fDecayTime->SetLineColor(iEn+1);
    for (int j = 0; j< 10; j++) hPulseShape[iEn]->Fit("fDecayTime", "QR");
  
    */
    
    hInputAbsorption[iEn]->SetLineColor(iEn+1);
    hDepoEn[iEn]->SetLineColor(iEn+1);
    for (int iAtt = 0; iAtt < nATT; iAtt++) hAttEn[iEn][iAtt]->SetLineColor(iEn+1);
    
    
//     if (iEn == 0) double initialLY = hLightOutput[0]->GetMean();
    grResolution->SetPoint(iEn, input_energy[iEn], hDepoEn[iEn]->GetRMS()/hDepoEn[iEn]->GetMean());
    grLinearity->SetPoint(iEn, input_energy[iEn], hDepoEn[iEn]->GetMean()/input_energy[iEn]);
    
    for (int iAtt = 0; iAtt < nATT; iAtt++){
      grResolutionAtt[iAtt]->SetPoint(iEn, input_energy[iEn], hAttEn[iEn][iAtt]->GetRMS()/hAttEn[iEn][iAtt]->GetMean());
//       grResolution[iAtt]->SetPointError(iEn, 0, hLightOutput[iEn]->GetMeanError());
      grLinearityAtt[iAtt]->SetPoint(iEn, input_energy[iEn], hAttEn[iEn][iAtt]->GetMean()/input_energy[iEn]);
    }
    
   }  //end of energy loop
    

  /// ********* DRAWING ***************** ///

    TCanvas * cInputAbsorption = new TCanvas ("cInputAbsorption", "cInputAbsorption", 600, 600);
    hInputAbsorption[0]->Draw();
    hInputAbsorption[0]->GetXaxis()->SetTitle("Wavelength (nm)");
    hInputAbsorption[0]->GetYaxis()->SetTitle("Absorption Length (mm)");
    for (int iEn = 1; iEn < NENERGIES; iEn++){
      hInputAbsorption[iEn]->Draw("same");
    }
    gPad->SetGrid();
    
    
    TCanvas * cDepoEn = new TCanvas ("cDepoEn", "cDepoEn", 600, 600);
    hDepoEn[0]->Draw();
    hDepoEn[0]->GetXaxis()->SetTitle("Deposited Energy (GeV)");
    hDepoEn[0]->GetYaxis()->SetTitle("Counts");
    for (int iEn = 1; iEn < NENERGIES; iEn++){
      hDepoEn[iEn]->Draw("same");
    }
    gPad->SetGrid();
    
    
    int selEn = 0;
    TCanvas * cAttEn = new TCanvas ("cAttEn", "cAttEn", 600, 600);
//     cAttEn->Divide(2,3);
    hAttEn[selEn][0]->Draw();
    hAttEn[selEn][0]->GetXaxis()->SetTitle("Measured Energy (GeV)");
    hAttEn[selEn][0]->GetYaxis()->SetTitle("Counts");
    
    for (int iAtt = 1; iAtt < nATT; iAtt++){
      hAttEn[selEn][iAtt]->Draw("same");
    }
    gPad->SetGrid();
    
    
    TCanvas * cResolution = new TCanvas ("cResolution", "cResolution", 600, 600);
    grResolution->SetMarkerStyle (21);
    grResolution->GetXaxis()->SetTitle("Energy (GeV)");
    grResolution->GetYaxis()->SetTitle("#sigma_{E}/E");
    grResolution->Draw("ALPE");
    gPad->SetGrid();

    TCanvas * cLinearity = new TCanvas ("cLinearity", "cLinearity", 600, 600);
    grLinearity->SetMarkerStyle (21);
    grLinearity->GetXaxis()->SetTitle("Energy (GeV)");
    grLinearity->GetYaxis()->SetTitle("Linearity");
    grLinearity->GetXaxis()->SetRangeUser(0, 12);
    grLinearity->GetYaxis()->SetRangeUser(0, 11);
    grLinearity->Draw("ALPE");
    gPad->SetGrid();

}


<file_sep>#include "CreateTree.hh"



CreateTree* CreateTree::fInstance = NULL;



CreateTree::CreateTree(TString name, bool energy_fiber, bool init_data, bool pos_fiber, bool opPhotons, bool timing, bool double_ro)
{
  if( fInstance )
  {
    return;
  }
  
  this -> ENERGY_FIBER 	= energy_fiber;
  this -> INIT_DATA    	= init_data;
  this -> POS_FIBER    	= pos_fiber;
  this -> OPPHOTONS     = opPhotons;
  this -> TIMING        = timing;
  this -> DOUBLE        = double_ro;
  
  this -> fInstance = this;
  this -> fname     = name;
  this -> ftree     = new TTree(name,name);
  
  this->GetTree()->Branch("Event",&this->Event,"Event/I");
  
//   this->GetTree()->Branch("Tot_phot_cer",		   &this->Tot_phot_cer,			 "Tot_phot_cer/F");
//   this->GetTree()->Branch("Tot_phot_scint",		   &this->Tot_phot_scint,			 "Tot_phot_scint/F");
  
//   this->GetTree()->Branch("Total_delta_world",        &this->Total_delta_world,        "Total_delta_world/F");
//   this->GetTree()->Branch("Total_energy_world",       &this->Total_energy_world,       "Total_energy_world/F");
//   this->GetTree()->Branch("Total_ion_energy_world",   &this->Total_ion_energy_world,   "Total_ion_energy_world/F");
//   this->GetTree()->Branch("Total_nonion_energy_world",&this->Total_nonion_energy_world,"Total_nonion_energy_world/F");
  
//   this->GetTree()->Branch("Total_em_energy",&this->Total_em_energy,"Total_em_energy/F");
  
  if( this -> INIT_DATA )
  {
 //   this->GetTree()->Branch("InitialEnergy",&this->InitialEnergy,"InitialEnergy/F");
    
    this->GetTree()->Branch("InitialPositionX",&this->InitialPositionX,"InitialPositionX/F");
    this->GetTree()->Branch("InitialPositionY",&this->InitialPositionY,"InitialPositionY/F");
    this->GetTree()->Branch("InitialPositionZ",&this->InitialPositionZ,"InitialPositionZ/F");
    
//    this->GetTree()->Branch("InitalMomentumDirectionX",&this->InitalMomentumDirectionX,"InitalMomentumDirectionX/F");
//    this->GetTree()->Branch("InitalMomentumDirectionY",&this->InitalMomentumDirectionY,"InitalMomentumDirectionY/F");
//    this->GetTree()->Branch("InitalMomentumDirectionZ",&this->InitalMomentumDirectionZ,"InitalMomentumDirectionZ/F");
  }
    
  if( this -> ENERGY_FIBER)
  {
    this->GetTree()->Branch("Total_energy_dead_material",&this->Total_energy_dead_material,"Total_energy_dead_material/F");
    this->GetTree()->Branch("Plastic_tag_ID",&this->Plastic_tag_ID,"Plastic_tag_ID/I");
    this->GetTree()->Branch("Total_energy",&this->Total_energy,"Total_energy/F");
    this->GetTree()->Branch("Total_ion_energy",&this->Total_ion_energy,"Total_ion_energy/F");	  
    this->GetTree()->Branch("Total_ion_energy_att",&this->Total_ion_energy_att,"Total_ion_energy_att[13]/F");
    this->GetTree()->Branch("Total_ion_energy_att_rear",&this->Total_ion_energy_att_rear,"Total_ion_energy_att_rear[13]/F");
    this->GetTree()->Branch("Total_ion_energy_att_front",&this->Total_ion_energy_att_front,"Total_ion_energy_att_front[13]/F");
  }
  
  if( this -> POS_FIBER)
  { 
    this->GetTree()->Branch("depositionX",&depositionX);
    this->GetTree()->Branch("depositionY",&depositionY);
    this->GetTree()->Branch("depositionZ",&depositionZ);  
    this->GetTree()->Branch("Energy_deposited",&Energy_deposited);
    this->GetTree()->Branch("Time_deposit",&Time_deposit);
  }
  
  if( this -> OPPHOTONS)
  {
    this->GetTree()->Branch("opPhoton_n",&this->opPhoton_n,"opPhoton_n/I");
//    this->GetTree()->Branch("opPhoton_n_type",&this->opPhoton_n_type);
    
//    this->GetTree()->Branch("opPhoton_n_ext",&this->opPhoton_n_ext,"opPhoton_n_ext/I");
//    this->GetTree()->Branch("opPhoton_n_ext_type",&this->opPhoton_n_ext_type);
//    this->GetTree()->Branch("opPhoton_n_ext_end",&this->opPhoton_n_ext_end);
    
    this->GetTree()->Branch("opPhoton_n_det",&this->opPhoton_n_det,"opPhoton_n_det/I");
//    this->GetTree()->Branch("opPhoton_n_det_type",&this->opPhoton_n_det_type);
    this->GetTree()->Branch("opPhoton_n_det_end",&this->opPhoton_n_det_end);
       
//    this->GetTree()->Branch("opPhoton_waveLength_det",&opPhoton_waveLength_det);
//    this->GetTree()->Branch("opPhoton_time_det",&opPhoton_time_det);
//    this->GetTree()->Branch("opPhoton_trackLength_det",&opPhoton_trackLength_det);
  }
  
  this->Clear();
}



CreateTree::~CreateTree()
{}



bool CreateTree::Write()
{
  TString filename = this->GetName();
  filename+=".root";
  TFile* file = new TFile(filename,"RECREATE");
  this->GetTree()->Write();
  file->Write();
  file->Close();
  return true;
}



void CreateTree::Clear()
{
  Event	= 0;
  /*
  Tot_phot_cer = 0;
  Tot_phot_scint = 0;*/
  /*
  Total_delta_world         = 0;
  Total_energy_world        = 0;
  Total_ion_energy_world    = 0;
  Total_nonion_energy_world = 0;*/
  
  if( this->INIT_DATA )
  {
    InitialEnergy = 0;
    
    InitialPositionX = 0;
    InitialPositionY = 0;
    InitialPositionZ = 0;
    
    InitalMomentumDirectionX = 0;
    InitalMomentumDirectionY = 0;
    InitalMomentumDirectionZ = 0;
  }
  
  if( this->ENERGY_FIBER )
  {

	Total_energy_dead_material = 0;
	Plastic_tag_ID = -1;
        Total_energy = 0;
        Total_ion_energy = 0;	  		    	  
	for (int j = 0; j< 13; j++)
	{
	  Total_ion_energy_att      [j] = 0;
	  Total_ion_energy_att_rear [j] = 0;
	  Total_ion_energy_att_front[j] = 0;
	}
  }
  
  if( this->POS_FIBER )
  {
    depositionX.clear();		
    depositionY.clear();		
    depositionZ.clear();
    Energy_deposited.clear();
    Time_deposit.clear();
  }
  
  if( this->OPPHOTONS )
  {

    opPhoton_n = 0;
    opPhoton_n_type.clear();
      
    opPhoton_n_ext = 0;
    opPhoton_n_ext_type.clear();
    opPhoton_n_ext_end.clear();
      
    opPhoton_n_det = 0;
    opPhoton_n_det_type.clear();
    opPhoton_n_det_end.clear();

    opPhoton_waveLength_det.clear();
    opPhoton_time_det.clear();
    opPhoton_trackLength_det.clear();
  }
}
<file_sep>ECAL Endcap tapered PbWO4 crystal
==========

Single PbWO4 crystal of CMS EE tapered geometry, read-out using a PMT/VPT like photodetector on the rear face.

Possibility to add a front photodetector to study double-side read-out configuration.

Software allows to study the light collection efficiency inside the crystal as well as the energy resolution and linearity of response to high energy particles by factorizing the processes of:
1. shower development + energy deposition
2. convolution with look-up table of light collection efficiency curves for different levels of radiation damage
 
Effects of containment and beam spot size can also be investigated.
