test/Macros/ResolDraw.cc

0001 /**
0002  * This macro draws the resolutions for single muons quantities: pt, cotgTheta and phi.
0003  *
0004  * It does a rebin(4) if NbinsX > 50 and in the case of PtGenVsMu_ResoVSPt also a rebin(8) in y. <br>
0005  * It takes a new histogram (a TH1D) equal to the projection in X of the starting histogram (a TH2F) and it empties it
0006  * (so as to have the binning and the axis already set and an empty histogram). <br>
0007  * Takes also a profileX of the TH2F. <br>
0008  * For the fit in eta it takes the events with eta < 0 on those with eta > 0. <br>
0009  * In any case it takes the projection in y (ProjectionY) of the TH2F in each bin (from x to x, that is a single bin). <br>
0010  * It extracts the rms and its error from the gaussian fit and writes them in the TH1D described above.
0011  */
0012
0013 // Needed to use gROOT in a compiled macro
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 // Use this with the ResolutionAnalyzer files
0034 // TString mainNamePt("PtResolutionGenVSMu");
0035 // TString mainNameCotgTheta("CotgThetaResolutionGenVSMu");
0036 // TString mainNamePhi("PhiResolutionGenVSMu");
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   // in an interactive ROOT session, edit the file names
0048   // Target and FileList, then
0049   // root > .L hadd.C
0050   // root > hadd()
0051
0052   TList *FileList;
0053   // vector of names of the histograms to fit
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 //   gROOT->SetStyle("Plain");
0083 //   gStyle->SetCanvasColor(kWhite);
0084 //   gStyle->SetCanvasBorderMode(0);
0085 //   gStyle->SetPadBorderMode(0);
0086 //   gStyle->SetTitleFillColor(kWhite);
0087 //   gStyle->SetTitleColor(kWhite);
0088 //   gStyle->SetOptStat("nemruoi");
0089
0090   Target = TFile::Open( "redrawed_"+numString+".root", "RECREATE" );
0091
0092   FileList = new TList();
0093
0094   // ************************************************************
0095   // List of Files
0096   FileList->Add( TFile::Open(numString+"_MuScleFit.root") );    // 1
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   //  std::cout << "Target path: " << target->GetPath() << std::endl;
0107   TString path( (char*)strstr( target->GetPath(), ":" ) );
0108   path.Remove( 0, 2 );
0109
0110   TFile *first_source = (TFile*)sourcelist->First();
0111
0112   // Stores all the fits
0113   std::map<TString, std::vector<TH1D*> > fitHistograms;
0114
0115   first_source->cd( path );
0116   TDirectory *current_sourcedir = gDirectory;
0117   //gain time, do not add the objects in the list in memory
0118   Bool_t status = TH1::AddDirectoryStatus();
0119   TH1::AddDirectory(kFALSE);
0120
0121   // loop over all keys in this directory
0122   TIter nextkey( current_sourcedir->GetListOfKeys() );
0123   TKey *key, *oldkey=0;
0124   while ( (key = (TKey*)nextkey())) {
0125
0126     //keep only the highest cycle number for each key
0127     if (oldkey && !strcmp(oldkey->GetName(),key->GetName())) continue;
0128
0129     // read object from first source file
0130     first_source->cd( path );
0131     TObject *obj = key->ReadObj();
0132
0133     if ( obj->IsA()->InheritsFrom( "TH1" ) ) {
0134       // descendant of TH1 -> redraw it
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           // Perform different fits for different profiles
0146           TH2F *h2 = (TH2F*)obj;
0147
0148       // Rebin X
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       // Rebin Y
0159       if( rebinY != 0 ) {
0160             h2->RebinY(rebinY);
0161       }
0162
0163           TH1D * h1 = h2->ProjectionX();
0164           TH1D * h1RMS = h2->ProjectionX();
0165           // h1->Clear();
0166           h1->Reset();
0167           h1->SetName(*namesIt+"_resol");
0168           h1RMS->Reset();
0169           h1RMS->SetName(*namesIt+"_resolRMS");
0170           TProfile * profile = h2->ProfileX();
0171
0172           // This is used to fit half eta, needed to get the resolution function by value
0173           // with greater precision assuming symmetry in eta and putting all values in
0174           // -3,0 (as if -fabs(eta) is used).
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                 // std::cout << "sigma = " << rms << std::endl;
0194                 // std::cout << "sigma error = " << rmsError << std::endl;
0195                 // std::cout << "rms = " << rms << std::endl;
0196                 // std::cout << "rms error = " << rmsError << std::endl;
0197                 // Reverse x in the first half to the second half.
0198                 int xToFill = x;
0199                 // Bin 0 corresponds to bin=binNumber(the last bin, which is also considered in the loop).
0200                 if( *namesIt == mainNamePt+"_ResoVSEta" ) {
0201                   std::cout << mainNamePt+"_ResoVSEta" << std::endl;
0202                   if( x<xBins/2+1 ) xToFill = xBins+1 - x;
0203                 }
0204                 // std::cout << "x = " << x << ", xToFill = " << xToFill << std::endl;
0205                 // std::cout << "rms = " << rms << ", rmsError = " << rmsError << std::endl;
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                 // h2->ProjectionY("_px",x,x)->Write();
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                 //canvas->Write();
0223                 temp->Write();
0224               }
0225             }
0226             target->cd();
0227             profile->Write();
0228             // for ( int i=1; i<=h1->GetNbinsX(); ++i ) {
0229             //   std::cout << "bin["<<i<<"] = " << h1->GetBinContent(i) << std::endl;
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                 // double rms = temp->GetRMS();
0246                 // double rmsError = temp->GetRMSError();
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                   // NaN is the only value different from itself. Infact NaN is "not a number"
0259                   // and it is not equal to any value, including itself.
0260                   h1->SetBinContent(x, sigma);
0261                   h1->SetBinError(x, sigmaError);
0262                   h1RMS->SetBinContent(x, rms);
0263                   h1RMS->SetBinError(x, rmsError);
0264                 }
0265                 // h2->ProjectionY("_px",x,x)->Write();
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                 // canvas->Write();
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       // it's a subdirectory
0286
0287       // std::cout << "Found subdirectory " << obj->GetName() << std::endl;
0288
0289       // create a new subdir of same name and title in the target file
0290       target->cd();
0291       TDirectory *newdir = target->mkdir( obj->GetName(), obj->GetTitle() );
0292
0293       // newdir is now the starting point of another round of merging
0294       // newdir still knows its depth within the target file via
0295       // GetPath(), so we can still figure out where we are in the recursion
0296       draw( newdir, sourcelist, vecNames, doHalfEta, minEntries, rebinX, rebinY );
0297     }
0298     else {
0299       // object is of no type that we know or can handle
0300       // std::cout << "Unknown object type, name: "
0301       //      << obj->GetName() << " title: " << obj->GetTitle() << std::endl;
0302     }
0303
0304     // now write the compared histograms (which are "in" obj) to the target file
0305     // note that this will just store obj in the current directory level,
0306     // which is not persistent until the complete directory itself is stored
0307     // by "target->Write()" below
0308     if ( obj ) {
0309       target->cd();
0310
0311       if( obj->IsA()->InheritsFrom( "TH1" ) ) {
0312         // Write the superimposed histograms to file
0313         // obj->Write( key->GetName() );
0314       }
0315     }
0316
0317   } // while ( ( TKey *key = (TKey*)nextkey() ) )
0318
0319   // Save canvases of the fitted histograms
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   // save modifications to target file
0340   target->SaveSelf(kTRUE);
0341   TH1::AddDirectory(status);
0342 }