00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073 #include <iostream>
00074 #include <map>
00075 #include <cmath>
00076 #include <TMath.h>
00077 #include <TFile.h>
00078 #include <TTree.h>
00079 #include <TH1.h>
00080 #ifndef READ_BINNING_CINT
00081 #include <TDOMParser.h>
00082 #include <TXMLDocument.h>
00083 #include "TUnfoldBinningXML.h"
00084 #else
00085 #include "TUnfoldBinning.h"
00086 #endif
00087
00088 using namespace std;
00089
00090
00091 void testUnfold5c()
00092 {
00093
00094 TH1::SetDefaultSumw2();
00095
00096
00097
00098
00099 TFile *outputFile=new TFile("testUnfold5_histograms.root","recreate");
00100
00101
00102
00103
00104
00105 #ifdef READ_BINNING_CINT
00106 TFile *binningSchemes=new TFile("testUnfold5_binning.root");
00107 #endif
00108
00109 TUnfoldBinning *detectorBinning,*generatorBinning;
00110
00111 outputFile->cd();
00112
00113
00114 #ifndef READ_BINNING_CINT
00115 TDOMParser parser;
00116 Int_t error=parser.ParseFile("testUnfold5binning.xml");
00117 if(error) cout<<"error="<<error<<" from TDOMParser\n";
00118 TXMLDocument const *XMLdocument=parser.GetXMLDocument();
00119 detectorBinning=
00120 TUnfoldBinningXML::ImportXML(XMLdocument,"detector");
00121 generatorBinning=
00122 TUnfoldBinningXML::ImportXML(XMLdocument,"generator");
00123 #else
00124 binningSchemes->GetObject("detector",detectorBinning);
00125 binningSchemes->GetObject("generator",generatorBinning);
00126
00127 delete binningSchemes;
00128 #endif
00129 detectorBinning->Write();
00130 generatorBinning->Write();
00131
00132 if(detectorBinning) {
00133 detectorBinning->PrintStream(cout);
00134 } else {
00135 cout<<"could not read 'detector' binning\n";
00136 }
00137 if(generatorBinning) {
00138 generatorBinning->PrintStream(cout);
00139 } else {
00140 cout<<"could not read 'generator' binning\n";
00141 }
00142
00143
00144 const TUnfoldBinning *detectordistribution=
00145 detectorBinning->FindNode("detectordistribution");
00146
00147 const TUnfoldBinning *signalBinning=
00148 generatorBinning->FindNode("signal");
00149
00150 const TUnfoldBinning *bgrBinning=
00151 generatorBinning->FindNode("background");
00152
00153
00154
00155
00156
00157
00158 Float_t etaRec,ptRec,discr,etaGen,ptGen;
00159 Int_t istriggered,issignal;
00160
00161 outputFile->cd();
00162
00163 TH1 *histDataReco=detectorBinning->CreateHistogram("histDataReco");
00164 TH1 *histDataTruth=generatorBinning->CreateHistogram("histDataTruth");
00165
00166 TFile *dataFile=new TFile("testUnfold5_data.root");
00167 TTree *dataTree=(TTree *) dataFile->Get("data");
00168
00169 if(!dataTree) {
00170 cout<<"could not read 'data' tree\n";
00171 }
00172
00173 dataTree->ResetBranchAddresses();
00174 dataTree->SetBranchAddress("etarec",&etaRec);
00175 dataTree->SetBranchAddress("ptrec",&ptRec);
00176 dataTree->SetBranchAddress("discr",&discr);
00177
00178 dataTree->SetBranchAddress("istriggered",&istriggered);
00179
00180 dataTree->SetBranchAddress("etagen",&etaGen);
00181 dataTree->SetBranchAddress("ptgen",&ptGen);
00182 dataTree->SetBranchAddress("issignal",&issignal);
00183 dataTree->SetBranchStatus("*",1);
00184
00185
00186 cout<<"loop over data events\n";
00187
00188 for(Int_t ievent=0;ievent<dataTree->GetEntriesFast();ievent++) {
00189 if(dataTree->GetEntry(ievent)<=0) break;
00190
00191 if(istriggered) {
00192 Int_t binNumber=
00193 detectordistribution->GetGlobalBinNumber(ptRec,etaRec,discr);
00194 histDataReco->Fill(binNumber);
00195 }
00196
00197 if(issignal) {
00198
00199 Int_t binNumber=signalBinning->GetGlobalBinNumber(ptGen,etaGen);
00200 histDataTruth->Fill(binNumber);
00201 } else {
00202
00203 Int_t binNumber=bgrBinning->GetGlobalBinNumber(ptRec,etaRec);
00204 histDataTruth->Fill(binNumber);
00205 }
00206 }
00207
00208 delete dataTree;
00209 delete dataFile;
00210
00211
00212
00213
00214
00215 outputFile->cd();
00216
00217 TH2 *histMCGenRec=TUnfoldBinning::CreateHistogramOfMigrations
00218 (generatorBinning,detectorBinning,"histMCGenRec");
00219
00220 TFile *signalFile=new TFile("testUnfold5_signal.root");
00221 TTree *signalTree=(TTree *) signalFile->Get("signal");
00222
00223 if(!signalTree) {
00224 cout<<"could not read 'signal' tree\n";
00225 }
00226
00227 signalTree->ResetBranchAddresses();
00228 signalTree->SetBranchAddress("etarec",&etaRec);
00229 signalTree->SetBranchAddress("ptrec",&ptRec);
00230 signalTree->SetBranchAddress("discr",&discr);
00231 signalTree->SetBranchAddress("istriggered",&istriggered);
00232 signalTree->SetBranchAddress("etagen",&etaGen);
00233 signalTree->SetBranchAddress("ptgen",&ptGen);
00234 signalTree->SetBranchStatus("*",1);
00235
00236 cout<<"loop over MC signal events\n";
00237
00238 for(Int_t ievent=0;ievent<signalTree->GetEntriesFast();ievent++) {
00239 if(signalTree->GetEntry(ievent)<=0) break;
00240
00241
00242 Int_t genBin=signalBinning->GetGlobalBinNumber(ptGen,etaGen);
00243
00244
00245
00246 Int_t recBin=0;
00247 if(istriggered) {
00248 recBin=detectordistribution->GetGlobalBinNumber(ptRec,etaRec,discr);
00249 }
00250 histMCGenRec->Fill(genBin,recBin);
00251 }
00252
00253 delete signalTree;
00254 delete signalFile;
00255
00256 TFile *bgrFile=new TFile("testUnfold5_background.root");
00257 TTree *bgrTree=(TTree *) bgrFile->Get("background");
00258
00259 if(!bgrTree) {
00260 cout<<"could not read 'background' tree\n";
00261 }
00262
00263 bgrTree->ResetBranchAddresses();
00264 bgrTree->SetBranchAddress("etarec",&etaRec);
00265 bgrTree->SetBranchAddress("ptrec",&ptRec);
00266 bgrTree->SetBranchAddress("discr",&discr);
00267 bgrTree->SetBranchAddress("istriggered",&istriggered);
00268 bgrTree->SetBranchStatus("*",1);
00269
00270 cout<<"loop over MC background events\n";
00271
00272 for(Int_t ievent=0;ievent<bgrTree->GetEntriesFast();ievent++) {
00273 if(bgrTree->GetEntry(ievent)<=0) break;
00274
00275
00276
00277 if(istriggered) {
00278
00279 Int_t genBin=bgrBinning->GetGlobalBinNumber(ptRec,etaRec);
00280
00281 Int_t recBin=detectordistribution->GetGlobalBinNumber
00282 (ptRec,etaRec,discr);
00283 histMCGenRec->Fill(genBin,recBin);
00284 }
00285 }
00286
00287 delete bgrTree;
00288 delete bgrFile;
00289
00290 outputFile->Write();
00291 delete outputFile;
00292
00293 }