RecoMuon/test/remakeComposites.C

0001
0002 #include <iostream>
0003 #include <string>
0004 #include "TFile.h"
0005 #include "THashList.h"
0006 #include "TH1.h"
0007 #include "TKey.h"
0008 #include "TClass.h"
0009 #include "TSystem.h"
0010
0011 int remakeComposites( TString fname )
0012 {
0013   int a = remakeValidTrack(fname);
0014   int b = remakeValidMuon(fname);
0015   return a+b;
0016 }
0017
0018 int remakeValidTrack( TString fname)
0019 {
0020
0021   TFile* source = TFile::Open( fname , "UPDATE");
0022   if( source==0 ){
0023     return 1;
0024   }
0025   bool multiTrack = source->cd("DQMData/Track");
0026   if(!multiTrack) return -1;
0027
0028   TString path( (char*)strstr( gDirectory->GetPath(), ":" ) );
0029   path.Remove( 0, 2 );
0030
0031   TDirectory *current_sourcedir = source->GetDirectory(path);
0032   if (!current_sourcedir) {
0033     continue;
0034   }
0035
0036   // loop over all keys in this directory
0037   TChain *globChain = 0;
0038   TIter nextkey( current_sourcedir->GetListOfKeys() );
0039   TKey *key, *oldkey=0;
0040   //gain time, do not add the objects in the list in memory
0041   //TH1::AddDirectory(kFALSE);
0042   THashList allNames;
0043   while ( (key = (TKey*)nextkey())) {
0044     if (oldkey && !strcmp(oldkey->GetName(),key->GetName())) continue;
0045
0046     // read object from first source file
0047     current_sourcedir->cd();
0048     TObject *obj = key->ReadObj();
0049     if ( obj->IsA()->InheritsFrom( "TDirectory" ) ) {
0050       // it's a subdirectory
0051       allNames.Add(new TObjString(key->GetName()));
0052       cout << "Found subdirectory " << obj->GetName() << endl;
0053       ((TDirectory*)obj)->cd();
0054
0055       TH1F * sim2reco_eta, *reco2sim_eta, *reco_eta, *sim_eta;
0056       TH1F * sim2reco_pt, *reco2sim_pt, *reco_pt, *sim_pt;
0057       TH1F *effic, *effic_pt, *fakerate;
0058
0059       gDirectory->GetObject("num_assoc(simToReco)_eta",sim2reco_eta);
0060       gDirectory->GetObject("num_assoc(recoToSim)_eta",reco2sim_eta);
0061       gDirectory->GetObject("num_assoc(simToReco)_pT",sim2reco_pt);
0062       //gDirectory->GetObject("num_assoc(recoToSim)_pT",reco2sim_pt);
0063       gDirectory->GetObject("num_reco_eta",reco_eta);
0064       gDirectory->GetObject("num_simul_eta",sim_eta);
0065       //gDirectory->GetObject("num_reco_pT",reco_pt);
0066       gDirectory->GetObject("num_simul_pT",sim_pt);
0067
0068       gDirectory->GetObject("effic",effic); effic->Reset();
0069       gDirectory->GetObject("efficPt",effic_pt); effic_pt->Reset();
0070       gDirectory->GetObject("fakerate",fakerate); fakerate->Reset();
0071
0072       //refill efficiency plot vs eta
0073
0074       makeEffTH1(sim2reco_eta,sim_eta,effic);
0075       makeEffTH1(sim2reco_pt,sim_pt,effic_pt);
0076       makeFakeTH1(reco2sim_eta,reco_eta,fakerate);
0077
0078       effic->Write("",TObject::kOverwrite);
0079       effic_pt->Write("",TObject::kOverwrite);
0080       fakerate->Write("",TObject::kOverwrite);
0081
0082       //refill slices
0083       computeSlice("dxyres_vs_eta","sigmadxy");
0084       computeSlice("dxyres_vs_pt","sigmadxyPt");
0085       computeSlice("ptres_vs_eta","sigmapt");
0086       computeSlice("ptres_vs_pt","sigmaptPt");
0087       computeSlice("dzres_vs_eta","sigmadz");
0088       computeSlice("dzres_vs_pt","sigmadzPt");
0089       computeSlice("phires_vs_eta","sigmaphi");
0090       computeSlice("phires_vs_pt","sigmaphiPt");
0091       computeSlice("cotThetares_vs_eta","sigmacotTheta");
0092       computeSlice("cotThetares_vs_pt","sigmacotThetaPt");
0093       //TH1F* ph1 = (TH1F*) chi2_vs_eta->ProfileX();
0094       //TH1F* ph2 = (TH1F*) nhits_vs_eta->ProfileX();
0095       computeSlice("dxypull_vs_eta","h_dxypulleta");
0096       computeSlice("ptpull_vs_eta","h_ptpulleta","h_ptshifteta");
0097       computeSlice("dzpull_vs_eta","h_dzpulleta");
0098       computeSlice("phipull_vs_eta","h_phipulleta");
0099       computeSlice("thetapull_vs_eta","h_thetapulleta");
0100
0101     }
0102   }
0103
0104   source->Close();
0105
0106   return 0;
0107 }
0108
0109 int remakeValidMuon( TString fname)
0110 {
0111
0112
0113   TFile* source = TFile::Open( fname , "UPDATE");
0114   if( source==0 ){
0115     return 1;
0116   }
0117   bool singleMuon = source->cd("DQMData/RecoMuonTask");
0118   if(!singleMuon) return -1;
0119
0120
0121   TString path( (char*)strstr( gDirectory->GetPath(), ":" ) );
0122   path.Remove( 0, 2 );
0123
0124   TDirectory *current_sourcedir = source->GetDirectory(path);
0125   if (!current_sourcedir) {
0126     continue;
0127   }
0128
0129   // loop over all keys in this directory
0130   TChain *globChain = 0;
0131   TIter nextkey( current_sourcedir->GetListOfKeys() );
0132   TKey *key, *oldkey=0;
0133   //gain time, do not add the objects in the list in memory
0134   //TH1::AddDirectory(kFALSE);
0135   THashList allNames;
0136   while ( (key = (TKey*)nextkey())) {
0137     //keep only the highest cycle number for each key
0138     if (oldkey && !strcmp(oldkey->GetName(),key->GetName())) continue;
0139
0140     // read object from first source file
0141     current_sourcedir->cd();
0142     TObject *obj = key->ReadObj();
0143     if ( obj->IsA()->InheritsFrom( "TDirectory" ) ) {
0144       // it's a subdirectory
0145       allNames.Add(new TObjString(key->GetName()));
0146       cout << "Found subdirectory " << obj->GetName() << endl;
0147       ((TDirectory*)obj)->cd();
0148
0149       computeMuonEff("Glb","Sim");
0150       computeMuonEff("Glb","Seed");
0151       computeMuonEff("Glb","Sta");
0152       computeMuonEff("Glb","Tk");
0153
0154       computeMuonEff("Sta","Sim");
0155       computeMuonEff("Sta","Seed");
0156
0157       computeMuonEff("Seed","Sim");
0158
0159       computeSlice("GlbEtaVsErrQPt","GlbPtResSigma","GlbPtResMean");
0160       computeSlice("StaEtaVsErrQPt","StaPtResSigma","StaPtResMean");
0161       computeSlice("SeedEtaVsErrQPt","SeedPtResSigma","SeedPtResMean");
0162
0163       computeSlice("GlbEtaVsPullPt","GlbPullPtSigma","GlbPullPtMean");
0164       computeSlice("StaEtaVsPullPt","StaPullPtSigma","StaPullPtMean");
0165       computeSlice("SeedEtaVsPullPt","SeedPullPtSigma","SeedPullPtMean");
0166
0167     }
0168   }
0169
0170   source->Close();
0171
0172   return 0;
0173 }
0174
0175 void computeSlice(TString name1, TString name2="", TString name3="")
0176 {
0177   TH2F *theTh2;
0178   cout << "computeSlice " << name1.Data() << endl;
0179   gDirectory->GetObject(name1,theTh2);
0180
0181
0182   theTh2->FitSlicesY();
0183   TH1* h1=(TH1*)gDirectory->Get(TString(theTh2->GetName())+"_1")->Clone();
0184   TH1* h2=(TH1*)gDirectory->Get(TString(theTh2->GetName())+"_2")->Clone();
0185   TH1* h3=(TH1*)gDirectory->Get(TString(theTh2->GetName())+"_chi2")->Clone();
0186
0187   h1->SetTitle(TString(theTh2->GetTitle())+" Gaussian Mean");
0188   h2->SetTitle(TString(theTh2->GetTitle())+" Gaussian Width");
0189   h3->SetTitle(TString(theTh2->GetTitle())+" Gaussian fit #chi^{2}");
0190
0191   TH1* h1a = h1->Clone();
0192   TH1* h2a = h2->Clone();
0193   h1a->Write(name3,TObject::kOverwrite);
0194   h2a->Write(name2,TObject::kOverwrite);
0195
0196   h1->Write("",TObject::kOverwrite);
0197   h2->Write("",TObject::kOverwrite);
0198   h3->Write("",TObject::kOverwrite);
0199 }
0200
0201
0202 void computeSlice(TH2F* theTh2)
0203 {
0204   theTh2->FitSlicesY();
0205   TH1* h1=(TH1*)gDirectory->Get(TString(theTh2->GetName())+"_1")->Clone();
0206   TH1* h2=(TH1*)gDirectory->Get(TString(theTh2->GetName())+"_2")->Clone();
0207   TH1* h3=(TH1*)gDirectory->Get(TString(theTh2->GetName())+"_chi2")->Clone();
0208
0209   h1->SetTitle(TString(theTh2->GetTitle())+" Gaussian Mean");
0210   h2->SetTitle(TString(theTh2->GetTitle())+" Gaussian Width");
0211   h3->SetTitle(TString(theTh2->GetTitle())+" Gaussian fit #chi^{2}");
0212
0213   h1->Write("",TObject::kOverwrite);
0214   h2->Write("",TObject::kOverwrite);
0215   h3->Write("",TObject::kOverwrite);
0216 }
0217
0218 void computeMuonEff(TString first, TString second)
0219 {
0220
0221   TH2F *numTH2, *denTH2;
0222   TH1D *num, *den;
0223   TH1F *effic;
0224
0225   gDirectory->GetObject(first+"EtaVsPhi",numTH2);
0226   gDirectory->GetObject(second+"EtaVsPhi",denTH2);
0227   gDirectory->GetObject(first+second+"_effEta",effic);
0228   num  = numTH2 ->ProjectionX();
0229   den  = denTH2 ->ProjectionX();
0230
0231   makeEffTH1(num,den,effic);
0232   effic->Write("",TObject::kOverwrite);
0233 }
0234
0235 void makeEffTH1(TH1* num, TH1* den, TH1* effic)
0236 {
0237
0238      //fill efficiency plot vs eta
0239       double eff,err;
0240       int nBins = num->GetNbinsX();
0241       for(int bin = 0; bin <= nBins + 1 ; bin++) {
0242         if (den->GetBinContent(bin) != 0 ){
0243           eff = ((double) num->GetBinContent(bin))/((double) den->GetBinContent(bin));
0244       err = sqrt(eff*(1-eff)/((double) den->GetBinContent(bin)));
0245           effic->SetBinContent(bin, eff);
0246           effic->SetBinError(bin,err);
0247         }
0248         else {
0249           effic->SetBinContent(bin, 0);
0250         }
0251       }
0252 }
0253
0254 void makeFakeTH1(TH1F* num, TH1F* den, TH1F* fakerate)
0255 {
0256       //fill fakerate plot
0257       double frate,ferr;
0258       int nBins = num->GetNbinsX();
0259       for (int bin=0; bin <= nBins + 1; bin++){
0260         if (den->GetBinContent(bin) != 0 ){
0261           frate = 1-((double) num->GetBinContent(bin))/((double) den->GetBinContent(bin));
0262       ferr = sqrt( frate*(1-frate)/(double) den->GetBinContent(bin) );
0263           fakerate->SetBinContent(bin, frate);
0264       fakerate->SetBinError(bin,ferr);
0265         }
0266         else {
0267           fakerate->SetBinContent(bin, 0.);
0268         }
0269       }
0270 }