Calorimetry/test/calculateECALresponceScales.C

0001 #include<iomanip> // provides setprecision
0002 #include<iostream>
0003 #include<sstream>
0004 #include<string>
0005
0006 // ROOT
0007 #include<TChain.h>
0008 #include<TFile.h>
0009 #include<TH3F.h>
0010 #include<TMath.h>
0011 #include<TProfile.h>
0012 #include<TROOT.h>
0013
0014 using namespace std;
0015
0016 std::map<int,std::vector<float>> getVectorsFromTree( TChain& tree, const TH3F& h, float minR ) {
0017   /* map.first:  global bin in TH3F h, belonging to e,eta,-1
0018    * map.second: vector of r in this bin (with this e,eta)
0019    */
0020
0021   std::map<int,std::vector<float>> vectors;
0022
0023   float r,e,eta;
0024   tree.SetBranchAddress("r",&r);
0025   tree.SetBranchAddress("e",&e);
0026   tree.SetBranchAddress("eta",&eta);
0027   for( int i=0; i<tree.GetEntries(); i++ ) {
0028     tree.GetEntry(i);
0029     if( r < minR ) continue;
0030     int bin = h.FindFixBin( e, eta, -1 );
0031
0032     // check if this bin already exists in map
0033     if( !vectors.count( bin ) ) {
0034       vectors[bin] = std::vector<float>();
0035     }
0036
0037     vectors.at(bin).push_back( r );
0038   }
0039
0040   for( auto& m : vectors ) {
0041     std::sort( m.second.begin(), m.second.end() );
0042   }
0043
0044   return vectors;
0045 }
0046
0047 void fillEmptyBinsWithUnity( TH3F& h3 ) {
0048   for( int x=0; x<h3.GetNbinsX()+2; x++ ) {
0049     for( int y=0; y<h3.GetNbinsY()+2; y++ ) {
0050       for( int z=0; z<h3.GetNbinsZ()+2; z++ ) {
0051         if( h3.GetBinContent(x,y,z) < 1e-5 ) {
0052           h3.SetBinContent(x,y,z,1.);
0053         }
0054       }
0055     }
0056   }
0057 }
0058
0059 class KKFactorsFactory {
0060   public:
0061
0062   KKFactorsFactory( const TH3F& h3_, const string& fileNameFast, const string& fileNameFull, const string& treeName ):
0063       h3(h3_),
0064       chFast(treeName.c_str()),
0065       chFull(treeName.c_str())
0066     {
0067       chFast.AddFile( fileNameFast.c_str() );
0068       chFull.AddFile( fileNameFull.c_str() );
0069     }
0070
0071   void calculate(){
0072     float minR = h3.GetZaxis()->GetXmin();
0073
0074     auto fastvectors = getVectorsFromTree( chFast, h3, minR );
0075     auto fullvectors = getVectorsFromTree( chFull, h3, minR );
0076
0077     for( auto m : fastvectors ) {
0078
0079       if( ! fullvectors.count( m.first ) ) {
0080         std::cout << "no compatible vector found for fullsim" << std::endl;
0081         continue;
0082       }
0083
0084       auto fastvector = m.second;
0085       auto fullvector = fullvectors.at( m.first );
0086
0087       TProfile profile("profile", "title", h3.GetZaxis()->GetNbins(), h3.GetZaxis()->GetXmin(), h3.GetZaxis()->GetXmax(), "s" );
0088       for( unsigned i=0; i<fastvector.size(); i++ ) {
0089         auto fa = fastvector[i];
0090         auto jRel = 1.*i/fastvector.size()*fullvector.size();
0091         int j = (int) jRel;
0092         auto fu = fullvector[j]*(jRel-j)+fullvector[j+1]*(j-jRel+1);
0093         if( fa ) profile.Fill( fa, fu/fa );
0094       }
0095
0096       int xbin, ybin, zbin;
0097       h3.GetBinXYZ( m.first, xbin, ybin, zbin );
0098
0099       for( int i=0; i<profile.GetNbinsX()+2; i++ ) {
0100         h3.SetBinContent( xbin, ybin, i, profile.GetBinContent(i) );
0101       }
0102
0103     }
0104
0105     fillEmptyBinsWithUnity( h3 );
0106
0107   }
0108
0109   TH3F GetH3() const { return h3; }
0110
0111   private:
0112   TH3F h3;
0113   TChain chFast;
0114   TChain chFull;
0115 };
0116
0117
0118