File indexing completed on 2023-03-17 11:14:51
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include "TROOT.h"
0015 #include "TStyle.h"
0016
0017 #include <string>
0018 #include <vector>
0019 #include <map>
0020 #include <sstream>
0021 #include "TChain.h"
0022 #include "TFile.h"
0023 #include "TH1.h"
0024 #include "TH2F.h"
0025 #include "TProfile.h"
0026 #include "TTree.h"
0027 #include "TKey.h"
0028 #include "Riostream.h"
0029 #include "THStack.h"
0030 #include "TCanvas.h"
0031 #include "TF1.h"
0032
0033
0034
0035
0036
0037 TString mainNamePt("hResolPtGenVSMu");
0038 TString mainNameCotgTheta("hResolCotgThetaGenVSMu");
0039 TString mainNamePhi("hResolPhiGenVSMu");
0040
0041 void draw( TDirectory *target, TList *sourcelist, const std::vector<TString> & vecNames,
0042 const bool doHalfEta, const int minEntries = 100, const int rebinX = 0, const int rebinY = 0 );
0043
0044 void ResolDraw(const TString numString = "0", const bool doHalfEta = false,
0045 const int minEntries = 100, const int rebinX = 0, const int rebinY = 0)
0046 {
0047
0048
0049
0050
0051
0052 TList *FileList;
0053
0054 std::vector<TString> vecNames;
0055
0056 TFile *Target;
0057
0058 vecNames.push_back("DeltaMassOverGenMassVsPt");
0059 vecNames.push_back("DeltaMassOverGenMassVsEta");
0060
0061 vecNames.push_back(mainNamePt + "_ResoVSPt");
0062 vecNames.push_back(mainNamePt + "_ResoVSPt_Bar");
0063 vecNames.push_back(mainNamePt + "_ResoVSPt_Endc_1.7");
0064 vecNames.push_back(mainNamePt + "_ResoVSPt_Endc_2.0");
0065 vecNames.push_back(mainNamePt + "_ResoVSPt_Endc_2.4");
0066 vecNames.push_back(mainNamePt + "_ResoVSPt_Ovlap");
0067 vecNames.push_back(mainNamePt + "_ResoVSEta");
0068 vecNames.push_back(mainNamePt + "_ResoVSPhiMinus");
0069 vecNames.push_back(mainNamePt + "_ResoVSPhiPlus");
0070
0071 vecNames.push_back(mainNameCotgTheta + "_ResoVSPt");
0072 vecNames.push_back(mainNameCotgTheta + "_ResoVSEta");
0073 vecNames.push_back(mainNameCotgTheta + "_ResoVSPhiMinus");
0074 vecNames.push_back(mainNameCotgTheta + "_ResoVSPhiPlus");
0075
0076 vecNames.push_back(mainNamePhi + "_ResoVSPt");
0077 vecNames.push_back(mainNamePhi + "_ResoVSEta");
0078 vecNames.push_back(mainNamePhi + "_ResoVSPhiMinus");
0079 vecNames.push_back(mainNamePhi + "_ResoVSPhiPlus");
0080
0081 gROOT->SetBatch(true);
0082
0083
0084
0085
0086
0087
0088
0089
0090 Target = TFile::Open( "redrawed_"+numString+".root", "RECREATE" );
0091
0092 FileList = new TList();
0093
0094
0095
0096 FileList->Add( TFile::Open(numString+"_MuScleFit.root") );
0097
0098 draw( Target, FileList, vecNames, doHalfEta, minEntries, rebinX, rebinY );
0099
0100 Target->Close();
0101 }
0102
0103 void draw( TDirectory *target, TList *sourcelist, const std::vector<TString> & vecNames,
0104 const bool doHalfEta, const int minEntries, const int rebinX, const int rebinY )
0105 {
0106
0107 TString path( (char*)strstr( target->GetPath(), ":" ) );
0108 path.Remove( 0, 2 );
0109
0110 TFile *first_source = (TFile*)sourcelist->First();
0111
0112
0113 std::map<TString, std::vector<TH1D*> > fitHistograms;
0114
0115 first_source->cd( path );
0116 TDirectory *current_sourcedir = gDirectory;
0117
0118 Bool_t status = TH1::AddDirectoryStatus();
0119 TH1::AddDirectory(kFALSE);
0120
0121
0122 TIter nextkey( current_sourcedir->GetListOfKeys() );
0123 TKey *key, *oldkey=0;
0124 while ( (key = (TKey*)nextkey())) {
0125
0126
0127 if (oldkey && !strcmp(oldkey->GetName(),key->GetName())) continue;
0128
0129
0130 first_source->cd( path );
0131 TObject *obj = key->ReadObj();
0132
0133 if ( obj->IsA()->InheritsFrom( "TH1" ) ) {
0134
0135
0136 TString objName = obj->GetName();
0137 std::vector<TString>::const_iterator namesIt = vecNames.begin();
0138 for( ; namesIt != vecNames.end(); ++namesIt ) {
0139 if( *namesIt == objName ) {
0140 std::cout << "found histogram: " << *namesIt << std::endl;
0141
0142 TDirectory * fits = (TDirectory*) target->Get("fits");
0143 if( fits == 0 ) fits = target->mkdir("fits");
0144
0145
0146 TH2F *h2 = (TH2F*)obj;
0147
0148
0149 int xBins = h2->GetNbinsX();
0150 if( rebinX != 0 ) {
0151 h2->RebinX(rebinX);
0152 xBins /= rebinX;
0153 }
0154 else if( xBins > 500 ) {
0155 h2->RebinX(2);
0156 xBins /= 2;
0157 }
0158
0159 if( rebinY != 0 ) {
0160 h2->RebinY(rebinY);
0161 }
0162
0163 TH1D * h1 = h2->ProjectionX();
0164 TH1D * h1RMS = h2->ProjectionX();
0165
0166 h1->Reset();
0167 h1->SetName(*namesIt+"_resol");
0168 h1RMS->Reset();
0169 h1RMS->SetName(*namesIt+"_resolRMS");
0170 TProfile * profile = h2->ProfileX();
0171
0172
0173
0174
0175 if( *namesIt == mainNamePt+"_ResoVSEta" && doHalfEta ) {
0176 std::cout << mainNamePt+"_ResoVSEta" << std::endl;
0177 std::cout << "bins%2 = " << xBins%2 << std::endl;
0178 for( int x=1; x<=xBins/2; ++x ) {
0179 std::stringstream fitNum;
0180 fitNum << x;
0181 TString fitName(*namesIt);
0182 fitName += "_fit_";
0183 TH1D * temp = h2->ProjectionY(fitName+fitNum.str(),x,x);
0184 TH1D * temp2 = h2->ProjectionY(fitName+fitNum.str(),xBins+1-x,xBins+1-x);
0185 temp->Add(temp2);
0186 if( temp->GetEntries() > minEntries ) {
0187 fitHistograms[*namesIt].push_back(temp);
0188 temp->Fit("gaus");
0189 double sigma = temp->GetFunction("gaus")->GetParameter(2);
0190 double sigmaError = temp->GetFunction("gaus")->GetParError(2);
0191 double rms = temp->GetRMS();
0192 double rmsError = temp->GetRMSError();
0193
0194
0195
0196
0197
0198 int xToFill = x;
0199
0200 if( *namesIt == mainNamePt+"_ResoVSEta" ) {
0201 std::cout << mainNamePt+"_ResoVSEta" << std::endl;
0202 if( x<xBins/2+1 ) xToFill = xBins+1 - x;
0203 }
0204
0205
0206 h1->SetBinContent(x, sigma);
0207 h1->SetBinError(x, sigmaError);
0208 h1->SetBinContent(xBins+1-x, 0);
0209 h1->SetBinError(xBins+1-x, 0);
0210
0211 h1RMS->SetBinContent(x, rms);
0212 h1RMS->SetBinError(x, rmsError);
0213 h1RMS->SetBinContent(xBins+1-x, 0);
0214 h1RMS->SetBinError(xBins+1-x, 0);
0215
0216 fits->cd();
0217 TCanvas * canvas = new TCanvas(temp->GetName()+TString("_canvas"), temp->GetTitle(), 1000, 800);
0218 temp->Draw();
0219 TF1 * gaussian = temp->GetFunction("gaus");
0220 gaussian->SetLineColor(kRed);
0221 gaussian->Draw("same");
0222
0223 temp->Write();
0224 }
0225 }
0226 target->cd();
0227 profile->Write();
0228
0229
0230
0231 h1->Write();
0232 h1RMS->Write();
0233 }
0234 else {
0235 for( int x=1; x<=xBins; ++x ) {
0236 std::stringstream fitNum;
0237 fitNum << x;
0238 TString fitName(*namesIt);
0239 fitName += "_fit_";
0240 TH1D * temp = h2->ProjectionY(fitName+fitNum.str(),x,x);
0241 if( temp->GetEntries() > minEntries ) {
0242 fitHistograms[*namesIt].push_back(temp);
0243 temp->Fit("gaus");
0244
0245
0246
0247
0248 double sigma = temp->GetFunction("gaus")->GetParameter(2);
0249 double sigmaError = temp->GetFunction("gaus")->GetParError(2);
0250 double rms = temp->GetRMS();
0251 double rmsError = temp->GetRMSError();
0252 if( sigma != sigma ) std::cout << "value is NaN: rms = " << rms << std::endl;
0253 if( sigma == sigma ) {
0254
0255 std::cout << "rms = " << rms << std::endl;
0256 std::cout << "rms error = " << rmsError << std::endl;
0257
0258
0259
0260 h1->SetBinContent(x, sigma);
0261 h1->SetBinError(x, sigmaError);
0262 h1RMS->SetBinContent(x, rms);
0263 h1RMS->SetBinError(x, rmsError);
0264 }
0265
0266 fits->cd();
0267 TCanvas * canvas = new TCanvas(temp->GetName()+TString("_canvas"), temp->GetTitle(), 1000, 800);
0268 temp->Draw();
0269 TF1 * gaussian = temp->GetFunction("gaus");
0270 gaussian->SetLineColor(kRed);
0271 gaussian->Draw("same");
0272
0273 temp->Write();
0274 }
0275 }
0276 target->cd();
0277 profile->Write();
0278 h1->Write();
0279 h1RMS->Write();
0280 }
0281 }
0282 }
0283 }
0284 else if ( obj->IsA()->InheritsFrom( "TDirectory" ) ) {
0285
0286
0287
0288
0289
0290 target->cd();
0291 TDirectory *newdir = target->mkdir( obj->GetName(), obj->GetTitle() );
0292
0293
0294
0295
0296 draw( newdir, sourcelist, vecNames, doHalfEta, minEntries, rebinX, rebinY );
0297 }
0298 else {
0299
0300
0301
0302 }
0303
0304
0305
0306
0307
0308 if ( obj ) {
0309 target->cd();
0310
0311 if( obj->IsA()->InheritsFrom( "TH1" ) ) {
0312
0313
0314 }
0315 }
0316
0317 }
0318
0319
0320 std::map<TString, std::vector<TH1D*> >::const_iterator mapIt = fitHistograms.begin();
0321 for( ; mapIt != fitHistograms.end(); ++mapIt ) {
0322 TCanvas * canvas = new TCanvas(mapIt->first, mapIt->first, 1000, 800);
0323
0324 int sizeCheck = mapIt->second.size();
0325 int x = int(sqrt(sizeCheck));
0326 int y = x;
0327 if( x*y < sizeCheck ) y += 1;
0328 if( x*y < sizeCheck ) x += 1;
0329 canvas->Divide(x,y);
0330 std::vector<TH1D*>::const_iterator histoIt = mapIt->second.begin();
0331 int histoNum = 1;
0332 for( ; histoIt != mapIt->second.end(); ++histoIt, ++histoNum ) {
0333 canvas->cd(histoNum);
0334 (*histoIt)->Draw();
0335 }
0336 canvas->Write();
0337 }
0338
0339
0340 target->SaveSelf(kTRUE);
0341 TH1::AddDirectory(status);
0342 }