1106 if (
sqrts<13) assert(0);
1110 TString TREE_NAME =
"ZZTree/candTree";
1111 TString COUNTERS_NAME =
"ZZTree/Counters";
1112 const bool writeFinalTree =
false && debug;
1114 short Z1Flav, Z2Flav, njets;
1115 float varTrue, varReco, ZZMass, weight;
1118 std::vector<T>* LepPt=0;
1119 std::vector<T>* LepEta=0;
1120 std::vector<T>* LepPhi=0;
1121 std::vector<short>* LepLepId=0;
1122 std::vector<T>* JetPt=0;
1123 std::vector<T>* JetEta=0;
1124 std::vector<T>* JetPhi=0;
1125 std::vector<T>* JetMass=0;
1127 std::vector<T>* LHEMotherPz=0;
1128 std::vector<T>* LHEMotherE=0;
1129 std::vector<short>* LHEMotherId=0;
1130 std::vector<T>* LHEDaughterPt=0;
1131 std::vector<T>* LHEDaughterEta=0;
1132 std::vector<T>* LHEDaughterPhi=0;
1133 std::vector<T>* LHEDaughterMass=0;
1134 std::vector<short>* LHEDaughterId=0;
1135 std::vector<T>* LHEAssociatedParticlePt=0;
1136 std::vector<T>* LHEAssociatedParticleEta=0;
1137 std::vector<T>* LHEAssociatedParticlePhi=0;
1138 std::vector<T>* LHEAssociatedParticleMass=0;
1139 std::vector<short>* LHEAssociatedParticleId=0;
1141 float xsec=1, overallEventWeight=1;
1143 vector<SimpleEntry> index;
1144 TString strchannel[3]={
"4mu",
"4e",
"2mu2e" };
1146 vector<TString> dumappend;
1147 vector<TString> strSamples;
1150 appendVector<TString>(strSamples, dumappend);
1153 else if (strprod==
"VBF"){
1155 appendVector<TString>(strSamples, dumappend);
1157 else if (strprod==
"ZH"){
1159 appendVector<TString>(strSamples, dumappend);
1161 else if (strprod==
"WH"){
1163 appendVector<TString>(strSamples, dumappend);
1165 else if (strprod==
"JJEW"){
1167 appendVector<TString>(strSamples, dumappend);
1169 appendVector<TString>(strSamples, dumappend);
1171 appendVector<TString>(strSamples, dumappend);
1174 cerr <<
"Production " << strprod <<
" is unknown." << endl;
1178 TFile* foutput = TFile::Open(Form(
"resolution_mJJ_recoVStrue_%s_%.0fTeV.root", strprod.Data(),
sqrts),
"recreate");
1180 vector<TFile*> finputList;
1181 vector<TTree*> treeList;
1183 TString cinput_main;
1185 unordered_map<TTree*, pair<float, float>> nGenMap;
1186 unordered_map<TTree*, float> mass_map;
1187 unordered_map<float, vector<TTree*>> mass_sample_map;
1190 TREE_NAME =
"SelectedTree";
1194 for (
int is=0; is<(
int)strSamples.size(); is++){
1195 for (
unsigned int ic=0; ic<3; ic++){
1196 TString cinput = Form(
"%s/%s/%s", cinput_main.Data(), strchannel[ic].Data(), (strSamples[is]).Data());
1197 TFile* finput = TFile::Open(cinput,
"read");
1198 cout <<
"Opening file " << cinput <<
"..." << endl;
1201 if (finput->IsOpen() && !finput->IsZombie()){
1202 cout << cinput <<
" opened. Extracting tree " << TREE_NAME <<
"..." << endl;
1203 tree = (TTree*)finput->Get(TREE_NAME);
1205 cout << TREE_NAME <<
" is found." << endl;
1206 tree->SetBranchStatus(
"*", 0);
1207 tree->SetBranchStatus(
"MC_weight", 1); tree->SetBranchAddress(
"MC_weight", &overallEventWeight);
1208 tree->SetBranchStatus(
"ZZMass", 1); tree->SetBranchAddress(
"ZZMass", &ZZMass);
1209 tree->SetBranchStatus(
"NJets30", 1); tree->SetBranchAddress(
"nCleanedJetsPt30", &njets);
1210 tree->SetBranchStatus(
"JetPt", 1); tree->SetBranchAddress(
"JetPt", &JetPt);
1211 tree->SetBranchStatus(
"JetEta", 1); tree->SetBranchAddress(
"JetEta", &JetEta);
1212 tree->SetBranchStatus(
"JetPhi", 1); tree->SetBranchAddress(
"JetPhi", &JetPhi);
1213 tree->SetBranchStatus(
"JetMass", 1); tree->SetBranchAddress(
"JetMass", &JetMass);
1215 TH1F* htmp = (TH1F*)finput->Get(COUNTERS_NAME);
1216 pair<float, float> nsum(1, 1);
1219 nEntries += tree->GetEntries();
1220 treeList.push_back(tree);
1221 finputList.push_back(finput);
1224 cout <<
"Pole mass = " << polemass << endl;
1225 mass_map[tree]=polemass;
1226 if (mass_sample_map.find(polemass)==mass_sample_map.end()){
1227 cout <<
"Inserting new pole mass sample array" << endl;
1228 vector<TTree*> dumarr;
1229 mass_sample_map[polemass] = dumarr;
1231 mass_sample_map[polemass].push_back(tree);
1233 else finput->Close();
1235 else if (finput->IsOpen()) finput->Close();
1242 TREE_NAME =
"ZZTree/candTree";
1245 for (
int is=0; is<(
int)strSamples.size(); is++){
1246 TString cinput = Form(
"%s/%s/ZZ4lAnalysis.root", cinput_main.Data(), (strSamples[is]).Data());
1247 cout <<
"Opening file " << cinput <<
"..." << endl;
1248 TFile* finput = TFile::Open(cinput,
"read");
1249 cout <<
"File open attempt passed." << endl;
1252 if (finput->IsOpen() && !finput->IsZombie()){
1253 cout << cinput <<
" opened. Extracting tree " << TREE_NAME <<
"..." << endl;
1254 tree = (TTree*)finput->Get(TREE_NAME);
1256 cout << TREE_NAME <<
" is found." << endl;
1257 tree->SetBranchStatus(
"*", 0);
1258 tree->SetBranchStatus(
"xsec", 1); tree->SetBranchAddress(
"xsec", &xsec);
1259 tree->SetBranchStatus(
"overallEventWeight", 1); tree->SetBranchAddress(
"overallEventWeight", &overallEventWeight);
1260 tree->SetBranchStatus(
"ZZMass", 1); tree->SetBranchAddress(
"ZZMass", &ZZMass);
1262 tree->SetBranchStatus(
"LepPt", 1); tree->SetBranchAddress(
"LepPt", &LepPt);
1263 tree->SetBranchStatus(
"LepEta", 1); tree->SetBranchAddress(
"LepEta", &LepEta);
1264 tree->SetBranchStatus(
"LepPhi", 1); tree->SetBranchAddress(
"LepPhi", &LepPhi);
1265 tree->SetBranchStatus(
"LepLepId", 1); tree->SetBranchAddress(
"LepLepId", &LepLepId);
1267 tree->SetBranchStatus(
"nCleanedJetsPt30", 1); tree->SetBranchAddress(
"nCleanedJetsPt30", &njets);
1268 tree->SetBranchStatus(
"JetPt", 1); tree->SetBranchAddress(
"JetPt", &JetPt);
1269 tree->SetBranchStatus(
"JetEta", 1); tree->SetBranchAddress(
"JetEta", &JetEta);
1270 tree->SetBranchStatus(
"JetPhi", 1); tree->SetBranchAddress(
"JetPhi", &JetPhi);
1271 tree->SetBranchStatus(
"JetMass", 1); tree->SetBranchAddress(
"JetMass", &JetMass);
1273 tree->SetBranchStatus(
"LHEMotherPz", 1); tree->SetBranchAddress(
"LHEMotherPz", &LHEMotherPz);
1274 tree->SetBranchStatus(
"LHEMotherE", 1); tree->SetBranchAddress(
"LHEMotherE", &LHEMotherE);
1275 tree->SetBranchStatus(
"LHEMotherId", 1); tree->SetBranchAddress(
"LHEMotherId", &LHEMotherId);
1277 tree->SetBranchStatus(
"LHEDaughterPt", 1); tree->SetBranchAddress(
"LHEDaughterPt", &LHEDaughterPt);
1278 tree->SetBranchStatus(
"LHEDaughterEta", 1); tree->SetBranchAddress(
"LHEDaughterEta", &LHEDaughterEta);
1279 tree->SetBranchStatus(
"LHEDaughterPhi", 1); tree->SetBranchAddress(
"LHEDaughterPhi", &LHEDaughterPhi);
1280 tree->SetBranchStatus(
"LHEDaughterMass", 1); tree->SetBranchAddress(
"LHEDaughterMass", &LHEDaughterMass);
1281 tree->SetBranchStatus(
"LHEDaughterId", 1); tree->SetBranchAddress(
"LHEDaughterId", &LHEDaughterId);
1283 tree->SetBranchStatus(
"LHEAssociatedParticlePt", 1); tree->SetBranchAddress(
"LHEAssociatedParticlePt", &LHEAssociatedParticlePt);
1284 tree->SetBranchStatus(
"LHEAssociatedParticleEta", 1); tree->SetBranchAddress(
"LHEAssociatedParticleEta", &LHEAssociatedParticleEta);
1285 tree->SetBranchStatus(
"LHEAssociatedParticlePhi", 1); tree->SetBranchAddress(
"LHEAssociatedParticlePhi", &LHEAssociatedParticlePhi);
1286 tree->SetBranchStatus(
"LHEAssociatedParticleMass", 1); tree->SetBranchAddress(
"LHEAssociatedParticleMass", &LHEAssociatedParticleMass);
1287 tree->SetBranchStatus(
"LHEAssociatedParticleId", 1); tree->SetBranchAddress(
"LHEAssociatedParticleId", &LHEAssociatedParticleId);
1290 cout <<
"Cross section = " << xsec << endl;
1292 TH1F* htmp = (TH1F*)finput->Get(COUNTERS_NAME);
1293 pair<float, float> nsum(htmp->GetBinContent(1), htmp->GetBinContent(40));
1296 nEntries += tree->GetEntries();
1297 treeList.push_back(tree);
1298 finputList.push_back(finput);
1301 cout <<
"Pole mass = " << polemass << endl;
1302 mass_map[tree]=polemass;
1303 if (mass_sample_map.find(polemass)==mass_sample_map.end()){
1304 cout <<
"Inserting new pole mass sample array" << endl;
1305 vector<TTree*> dumarr;
1306 mass_sample_map[polemass] = dumarr;
1308 mass_sample_map[polemass].push_back(tree);
1310 else finput->Close();
1312 else if (finput->IsOpen()) finput->Close();
1317 for (
auto it = mass_sample_map.begin(); it != mass_sample_map.end(); ++it){
1320 for (
unsigned int ix=0; ix<it->second.size(); ix++){
1321 it->second.at(ix)->GetEntry(0);
1322 sum_ngen += nGenMap[it->second.at(ix)].first;
1325 float overallWeight = sum_ngen/sum_xsec;
1326 for (
unsigned int ix=0; ix<it->second.size(); ix++){
1327 cout <<
"Sum Hep MC weights in tree " << ix <<
" / " << it->second.size() <<
" was " << nGenMap[it->second.at(ix)].second <<
" over " << nGenMap[it->second.at(ix)].first <<
" total gen. events." << endl;
1328 nGenMap[it->second.at(ix)].first = overallWeight/nGenMap[it->second.at(ix)].second;
1329 cout <<
"Event scale for tree " << ix <<
" / " << it->second.size() <<
" at pole mass " << it->first <<
" = " << nGenMap[it->second.at(ix)].first << endl;
1333 cout <<
"NEntries = " << nEntries <<
" over " << treeList.size() <<
" trees." << endl;
1337 if (writeFinalTree){
1338 newtree =
new TTree(
"FinalTree",
"");
1339 newtree->Branch(
"varTrue", &varTrue);
1340 newtree->Branch(
"varReco", &varReco);
1341 newtree->Branch(
"ZZMass", &ZZMass);
1342 newtree->Branch(
"weight", &weight);
1346 for (
int ev=0;
ev<nEntries;
ev++){
1349 njets>=2 && LHEAssociatedParticlePt->size()>=2 && LepPt->size()==4
1351 if (!doProcess)
continue;
1353 weight = nGenMap[tree].first*xsec*overallEventWeight;
1355 vector<MELAParticle*> partList;
1358 vector<MELAParticle*> ttot, tg, tq, tothers;
1359 for (
unsigned int i=0; i<LHEAssociatedParticlePt->size(); i++){
1360 int id = LHEAssociatedParticleId->at(i);
1361 TLorentzVector tmp(0, 0, 0, 0);
1363 LHEAssociatedParticlePt->at(i),
1364 LHEAssociatedParticleEta->at(i),
1365 LHEAssociatedParticlePhi->at(i),
1366 LHEAssociatedParticleMass->at(i)
1372 else tothers.push_back(part);
1373 partList.push_back(part);
1377 vector<MELAParticle*> daus;
1378 TLorentzVector pSumDaus(0, 0, 0, 0);
1380 for (
unsigned int i=0; i<LepPt->size(); i++){
1381 int id = LepLepId->at(i);
1386 TLorentzVector tmp(0, 0, 0, 0);
1394 daus.push_back(part);
1395 pSumDaus = pSumDaus + tmp;
1397 partList.push_back(part);
1401 vector<MELACandidate*> candList;
1404 vector<MELAParticle*> tdaus;
1405 TLorentzVector pSumTrueDaus(0, 0, 0, 0);
1407 vector<MELAParticle*> tdaus_split[2][2];
1408 for (
unsigned int i=0; i<LHEDaughterPt->size(); i++){
1409 int id = LHEDaughterId->at(i);
1410 TLorentzVector tmp(0, 0, 0, 0);
1412 LHEDaughterPt->at(i),
1413 LHEDaughterEta->at(i),
1414 LHEDaughterPhi->at(i),
1415 LHEDaughterMass->at(i)
1418 tdaus.push_back(part);
1421 if (part->
id>0) ipart=0;
1423 if (std::abs(
id)==11) jdau=0;
1424 else if (std::abs(
id)==13) jdau=1;
1425 if (ipart>=0 && jdau>=0) tdaus_split[jdau][ipart].push_back(part);
1427 partList.push_back(part);
1430 for (
auto it=tothers.begin(); it<tothers.end(); it++){
1434 if (part->
id>0) ipart=0;
1436 if (std::abs(part->
id)==11) jdau=0;
1437 else if (std::abs(part->
id)==13) jdau=1;
1438 if (ipart>=0 && jdau>=0) tdaus_split[jdau][ipart].push_back(part);
1442 std::vector<MELAParticle*> tmpVhandle;
1443 for (
int c=0; c<2; c++){
1444 for (
unsigned int i=0; i<tdaus_split[c][0].size(); i++){
1445 for (
unsigned int j=0; j<tdaus_split[c][1].size(); j++){
1446 TLorentzVector pV = tdaus_split[c][0].at(i)->
p4+tdaus_split[c][1].at(j)->p4;
1450 tmpVhandle.push_back(V);
1455 for (
unsigned int i=0; i<tmpVhandle.size(); i++){
1456 for (
unsigned int j=i+1; j<tmpVhandle.size(); j++){
1461 if (Vi1==Vj1 || (Vi2==Vj2 && Vi2 != 0))
continue;
1463 TLorentzVector pH(0, 0, 0, 0);
1464 if (Vi1!=0) pH = pH + Vi1->
p4;
1465 if (Vi2!=0) pH = pH + Vi2->
p4;
1466 if (Vj1!=0) pH = pH + Vj1->
p4;
1467 if (Vj2!=0) pH = pH + Vj2->
p4;
1479 candList.push_back(cand);
1483 for (
unsigned int i=0; i<tmpVhandle.size(); i++)
delete tmpVhandle.at(i);
1485 for (
unsigned int icand=0; icand<candList.size(); icand++){
1487 float DM = fabs(cand->
m()-ZZMass);
1502 if (tq.size()+tg.size()<2 || Z1Z2id!=tZ1Z2id || theCand==0){
1503 for (
unsigned int icand=0; icand<candList.size(); icand++)
delete candList.at(icand);
1504 for (
unsigned int i=0; i<partList.size(); i++)
delete partList.at(i);
1509 TVector3 boostTrueDaus = pSumTrueDaus.BoostVector();
1510 TVector3 boostDaus = pSumDaus.BoostVector();
1513 vector<MELAParticle*> tmothers;
1514 for (
unsigned int i=0; i<LHEMotherPz->size(); i++){
1515 int id = LHEMotherId->at(i);
1516 TLorentzVector tmp(0, 0, LHEMotherPz->at(i), LHEMotherE->at(i));
1518 tmothers.push_back(part);
1519 partList.push_back(part);
1522 vector<MELAParticle*> jets;
1523 for (
unsigned int i=0; i<(
unsigned int)min(2, (
int)JetPt->size()); i++){
1525 TLorentzVector tmp(0, 0, 0, 0);
1533 jets.push_back(part);
1534 partList.push_back(part);
1537 for (
unsigned int i=0; i<daus.size(); i++){ daus.at(i)->p4.Boost(-boostDaus); daus.at(i)->p4.Boost(boostTrueDaus); }
1538 for (
unsigned int i=0; i<jets.size(); i++){ jets.at(i)->p4.Boost(-boostDaus); jets.at(i)->p4.Boost(boostTrueDaus); }
1539 varReco = (jets.at(0)->p4 + jets.at(1)->p4).M();
1541 if (ev_acc%10000==0)
cout <<
"Pre-processing event " <<
ev << endl;
1543 for (
unsigned int i=0; i<tg.size(); i++){
1544 int jfound=-1;
float dRmin=9e10;
1547 for (
unsigned int j=0; j<tq.size(); j++){
1549 float dR = 1./fabs(quark->
dot(gluon));
1555 for (
unsigned int j=0; j<tmothers.size(); j++){
1558 float dR = 1./fabs(quark->
dot(gluon));
1566 int&
id = tq.at(jfound)->id;
1567 TLorentzVector mom = tq.at(jfound)->p4 + gluon->
p4;
1571 ttot.push_back(composite);
1572 partList.push_back(composite);
1576 TLorentzVector sum(0, 0, 0, 0);
1577 vector<MELAParticle*> matchedpartons;
1578 for (
unsigned int ijet=0; ijet<2; ijet++){
1583 for (
unsigned int i=0; i<ttot.size(); i++){
1585 bool unmatched=
true;
1586 for (
unsigned int im=0; im<matchedpartons.size(); im++){
1588 matchedpartons.at(im)==part
1594 unmatched=
false;
break;
1597 if (!unmatched)
continue;
1598 float dR = part->
deltaR(jet);
1605 for (
unsigned int i=0; i<tq.size(); i++){
1607 bool unmatched=
true;
1608 for (
unsigned int im=0; im<matchedpartons.size(); im++){
1609 if (matchedpartons.at(im)==part) unmatched=
false;
1610 for (
int idau=0; idau<matchedpartons.at(im)->getNDaughters(); idau++){
1611 if (matchedpartons.at(im)->getDaughter(idau)==part){ unmatched=
false;
break; }
1613 if (!unmatched)
break;
1615 if (!unmatched)
continue;
1616 float dR = part->
deltaR(jet);
1623 for (
unsigned int i=0; i<tg.size(); i++){
1625 bool unmatched=
true;
1626 for (
unsigned int im=0; im<matchedpartons.size(); im++){
1627 if (matchedpartons.at(im)==part) unmatched=
false;
1628 for (
int idau=0; idau<matchedpartons.at(im)->getNDaughters(); idau++){
1629 if (matchedpartons.at(im)->getDaughter(idau)==part){ unmatched=
false;
break; }
1631 if (!unmatched)
break;
1633 if (!unmatched)
continue;
1634 float dR = part->
deltaR(jet);
1642 sum = sum + matched->
p4;
1643 matchedpartons.push_back(matched);
1649 if (varTrue<0. || ev_acc%100000==0){
1650 cout <<
"True mJJ = " << varTrue <<
", reco mJJ = " << varReco <<
". Event info:" << endl;
1651 cout <<
"- Best cand mass = " << pSumTrueDaus.M() <<
" vs. reco mass = " << ZZMass << endl;
1652 cout <<
" - Reco jets:" << endl;
1653 for (
unsigned int i=0; i<jets.size(); i++)
cout <<
" - " << jets.at(i)->p4.X() <<
" " << jets.at(i)->p4.Y() <<
" " << jets.at(i)->p4.Z() <<
" " << jets.at(i)->p4.T() << endl;
1654 cout <<
" - " <<
" - Matched truth:" << endl;
1655 for (
unsigned int i=0; i<matchedpartons.size(); i++)
cout <<
" - " << matchedpartons.at(i)->id <<
" " << matchedpartons.at(i)->p4.X() <<
" " << matchedpartons.at(i)->p4.Y() <<
" " << matchedpartons.at(i)->p4.Z() <<
" " << matchedpartons.at(i)->p4.T() << endl;
1656 cout <<
" - " <<
" - True partons:" << endl;
1657 for (
unsigned int i=0; i<tg.size(); i++)
cout <<
" - " << tg.at(i)->id <<
" " << tg.at(i)->p4.X() <<
" " << tg.at(i)->p4.Y() <<
" " << tg.at(i)->p4.Z() <<
" " << tg.at(i)->p4.T() << endl;
1658 for (
unsigned int i=0; i<tq.size(); i++)
cout <<
" - " << tq.at(i)->id <<
" " << tq.at(i)->p4.X() <<
" " << tq.at(i)->p4.Y() <<
" " << tq.at(i)->p4.Z() <<
" " << tq.at(i)->p4.T() << endl;
1659 cout <<
" - " <<
" - Merged partons:" << endl;
1660 for (
unsigned int i=0; i<ttot.size(); i++){
1661 cout <<
" - " << ttot.at(i)->id <<
" " << ttot.at(i)->p4.X() <<
" " << ttot.at(i)->p4.Y() <<
" " << ttot.at(i)->p4.Z() <<
" " << ttot.at(i)->p4.T() << endl;
1662 for (
int j=0; j<ttot.at(i)->getNDaughters(); j++)
cout <<
" ---Dau" << j+1 <<
": " << ttot.at(i)->getDaughter(j)->id <<
" " << ttot.at(i)->getDaughter(j)->p4.X() <<
" " << ttot.at(i)->getDaughter(j)->p4.Y() <<
" " << ttot.at(i)->getDaughter(j)->p4.Z() <<
" " << ttot.at(i)->getDaughter(j)->p4.T() << endl;
1664 cout <<
" - " <<
" - Other true asssociated:" << endl;
1665 for (
unsigned int i=0; i<tothers.size(); i++)
cout <<
" - " << tothers.at(i)->id <<
" " << tothers.at(i)->p4.X() <<
" " << tothers.at(i)->p4.Y() <<
" " << tothers.at(i)->p4.Z() <<
" " << tothers.at(i)->p4.T() << endl;
1666 cout <<
"True candidate summary:" << endl;
1670 for (
unsigned int icand=0; icand<candList.size(); icand++)
delete candList.at(icand);
1671 for (
unsigned int i=0; i<partList.size(); i++)
delete partList.at(i);
1678 cout <<
"Number of valid entries: " << ev_acc << endl;
1680 float firstVal=index.at(0).trueval;
1681 float lastVal=index.at(index.size()-1).trueval;
1682 float infimum = std::floor(firstVal);
1683 float supremum = std::ceil(lastVal);
1684 cout <<
"Nentries = " << nEntries <<
" | truth = " << firstVal <<
" - " << lastVal <<
"(" << infimum <<
", " << supremum <<
")" << endl;
1686 float divisor=20000;
1687 int nbins = index.size()/divisor;
1688 const int nbins_th=10;
1689 while (nbins<nbins_th){
1690 if (divisor>1000) divisor -= 1000;
1691 else if (divisor>100) divisor -= 100;
1693 nbins=index.size()/divisor;
1695 cout <<
"nbins=" << nbins << endl;
1696 if (nbins<3) cerr <<
"Not enough bins!" << endl;
1697 vector<ExtBin> binList;
1698 float* binning =
new float[nbins+1];
1700 binning[nbins]=supremum;
1701 int ev_stepsize = index.size()/nbins;
1702 cout <<
"Event step size: " << ev_stepsize << endl;
1703 cout <<
"Boundary (" << 0 <<
") = " << binning[0] << endl;
1704 for (
int ix=1; ix<nbins; ix++){
1705 binning[ix]=(index[ix*ev_stepsize-1].trueval+index[ix*ev_stepsize].trueval)*0.5;
1707 tmpbin.
binlow = binning[ix-1]; tmpbin.
binhigh = binning[ix];
1708 for (
int bin=0; bin<ev_stepsize; bin++){ tmpbin.
addEvent(index[(ix-1)*ev_stepsize+bin]); }
1709 binList.push_back(tmpbin);
1710 cout <<
"Boundary (" << ix <<
")= " << binning[ix] <<
" [event " << index[ix*ev_stepsize].id <<
", step " << ix*ev_stepsize <<
"]" << endl;
1713 tmpbin.
binlow = binning[nbins-1]; tmpbin.
binhigh = binning[nbins];
1714 for (
unsigned int bin=(nbins-1)*ev_stepsize; bin<index.size(); bin++){ tmpbin.
addEvent(index[bin]); }
1715 binList.push_back(tmpbin);
1716 cout <<
"Boundary (" << nbins <<
") = " << binning[nbins] << endl;
1717 cout <<
"Bin list has the following bins:" << endl;
1718 for (
unsigned int ib=0; ib<binList.size(); ib++){
1719 cout << ib <<
" / " << binList.size() <<
": [" << binList.at(ib).binlow <<
"," << binList.at(ib).binhigh <<
"]" << endl;
1724 TProfile* hvar =
new TProfile(
"varTrue",
"", nbins, binning); hvar->Sumw2();
1725 TProfile* havg_varReco =
new TProfile(
"avg_varReco",
"", nbins, binning); havg_varReco->Sumw2();
1726 TH2F* hvarTrue_varReco =
new TH2F(
"varTrue_varReco",
"", nbins, binning, nbins, binning); hvarTrue_varReco->Sumw2();
1728 vector<TGraphErrors*> tglist;
1731 for (
unsigned int bin=0; bin<binList.size(); bin++){
1732 cout <<
"Bin " << bin <<
" is now being scrutinized..." << endl;
1734 binList.at(bin).sift(0.99,
true);
1735 binList.at(bin).adjustWeights();
1737 for (
unsigned int ev=0;
ev<binList.at(bin).collection.size();
ev++){
1738 varTrue = binList.at(bin).collection.at(
ev).trueval;
1739 varReco = binList.at(bin).collection.at(
ev).recoval;
1740 ZZMass = binList.at(bin).collection.at(
ev).recotrackval;
1741 weight = binList.at(bin).collection.at(
ev).weight;
1742 havg_varReco->Fill(varTrue, varReco, weight);
1743 hvar->Fill(varTrue, varTrue, weight);
1744 hvarTrue_varReco->Fill(varTrue, varReco, weight);
1745 if (writeFinalTree) newtree->Fill();
1748 float avgbin_varTrue = hvar->GetBinContent(bin+1);
1750 vector<pair<float, int>> recominustruevals_entries;
1751 binList.at(bin).sortByRecoMinusTrueVals(recominustruevals_entries);
1752 vector<float> recominustruevalbounds;
1753 unsigned int nrecototal = recominustruevals_entries.size();
1754 for (
unsigned int ev=0;
ev<nrecototal;
ev++) recominustruevals_entries.at(
ev).first = (recominustruevals_entries.at(
ev).first+1.)*avgbin_varTrue;
1756 float recodivisor=1600;
1757 unsigned int nrecobins = nrecototal/recodivisor;
1758 const unsigned int nrecobins_th=12;
1759 while (nrecobins<nrecobins_th){
1760 if (recodivisor>100) recodivisor -= 100;
1761 else if (recodivisor>10) recodivisor -= 10;
1763 nrecobins = nrecototal/recodivisor;
1765 unsigned int recoinc = nrecototal / nrecobins;
1766 float minval = min(
float(0), std::floor(recominustruevals_entries.at(0).first));
1767 float maxval = max(std::ceil(recominustruevals_entries.at(nrecototal-1).first),
float(
sqrts*1000.));
1768 recominustruevalbounds.push_back(minval);
1769 for (
unsigned int i=1; i<nrecobins; i++){
1770 float a = recominustruevals_entries.at(i*recoinc).first;
1771 float b = recominustruevals_entries.at(i*recoinc-1).first;
1772 recominustruevalbounds.push_back((a+b)*0.5);
1774 recominustruevalbounds.push_back(maxval);
1775 cout <<
"Bin " << bin <<
" has the following reco-true boundaries:";
1776 for (
unsigned int j=0; j<recominustruevalbounds.size(); j++)
cout <<
" " << recominustruevalbounds.at(j);
1781 TH1F* h_RecoMinusTrueVals =
new TH1F(Form(
"RecoMinusTrueVal_Bin%i", bin),
"", nrecobins, recominustruevalbounds.data()); h_RecoMinusTrueVals->Sumw2();
1782 TProfile* p_RecoMinusTrueVals =
new TProfile(Form(
"avg_RecoMinusTrueVal_Bin%i", bin),
"", nrecobins, recominustruevalbounds.data()); p_RecoMinusTrueVals->Sumw2();
1783 for (
unsigned int i=0; i<nrecototal; i++){
1784 float& val = recominustruevals_entries.at(i).first;
1785 int&
ev = recominustruevals_entries.at(i).second;
1786 float& weight = binList.at(bin).collection.at(
ev).weight;
1787 if (weight==0.) cerr <<
"Weight is 0!" << endl;
1788 else if (
isnan(weight) || isinf(weight)) cerr <<
"Weight " << weight <<
" is nan or inf" << endl;
1789 if (
isnan(val) || isinf(val)) cerr <<
"Value " << val <<
" is nan or inf" << endl;
1790 h_RecoMinusTrueVals->Fill(val, weight);
1791 p_RecoMinusTrueVals->Fill(val, val, weight);
1793 const TAxis* axis = h_RecoMinusTrueVals->GetXaxis();
1794 for (
int bin=1; bin<=h_RecoMinusTrueVals->GetNbinsX(); bin++){
1795 float width = axis->GetBinWidth(bin);
1796 h_RecoMinusTrueVals->SetBinContent(bin, h_RecoMinusTrueVals->GetBinContent(bin)/width);
1797 h_RecoMinusTrueVals->SetBinError(bin, h_RecoMinusTrueVals->GetBinError(bin)/width);
1800 foutput->WriteTObject(h_RecoMinusTrueVals);
1801 foutput->WriteTObject(p_RecoMinusTrueVals);
1804 TGraphErrors* tgbin =
makeGraphFromTH1(p_RecoMinusTrueVals, h_RecoMinusTrueVals, Form(
"tg_%s", h_RecoMinusTrueVals->GetName()));
1805 tgbin->SetTitle(Form(
"<m^{true}_{JJ}> = %.3f #pm %.3f GeV, <m^{reco}_{JJ}> = %.3f #pm %.3f GeV", hvar->GetBinContent(bin+1), hvar->GetBinError(bin+1), havg_varReco->GetBinContent(bin+1), havg_varReco->GetBinError(bin+1)));
1806 tgbin->GetYaxis()->SetTitle(
"Sum of weights");
1807 tgbin->GetXaxis()->SetTitle(
"m^{reco}_{JJ} / m^{true}_{JJ} - 1");
1808 cout <<
"Adding point (0,1)" << endl;
1810 for (
int ip=0; ip<tgbin->GetN(); ip++) tgbin->GetY()[ip] = log(tgbin->GetY()[ip]);
1812 double dfirst, dlast;
1813 dfirst = (tgbin->GetY()[1]-tgbin->GetY()[0])/(tgbin->GetX()[1]-tgbin->GetX()[0]);
1814 dlast = (0-tgbin->GetY()[tgbin->GetN()-1])/(maxval-tgbin->GetX()[tgbin->GetN()-1]);
1815 TSpline3* sp_tgbin =
new TSpline3(
"sp_tgbin", tgbin,
"b1e1", dfirst, dlast);
1816 unsigned int ninsert;
1817 vector<double> xcoord;
1818 for (
int ip=0; ip<tgbin->GetN(); ip++) xcoord.push_back(tgbin->GetX()[ip]);
1829 unsigned int binzero=0;
1830 for (
unsigned int ip=0; ip<xcoord.size()-1; ip++){
1831 if (xcoord[ip]<avgbin_varTrue && xcoord[ip+1]>=avgbin_varTrue){
1836 for (
unsigned int ip=0; ip<xcoord.size()-1; ip++){
1837 if (ip==binzero || ip==binzero+1 || ip==(xcoord.size()-2)) ninsert=9;
1838 else if (ip<binzero) ninsert=5;
1841 sp_tgbin =
new TSpline3(
"sp_tgbin", tgbin,
"b1e1", dfirst, dlast);
1842 const double xwidth = (xcoord[ip+1] - xcoord[ip])/((
double)(ninsert+1));
1843 for (
unsigned int j=1; j<=ninsert; j++){
1844 double xtmp = xwidth*j+xcoord[ip];
1845 double ytmp = sp_tgbin->Eval(xtmp);
1846 cout <<
"Inserting (x,y) = ( " << xtmp <<
" , " << ytmp <<
" )" << endl;
1851 for (
int ip=0; ip<tgbin->GetN(); ip++) tgbin->GetY()[ip] = exp(tgbin->GetY()[ip]);
1854 dlast = (tgbin->GetY()[tgbin->GetN()-1]-tgbin->GetY()[tgbin->GetN()-2])/(tgbin->GetX()[tgbin->GetN()-1]-tgbin->GetX()[tgbin->GetN()-2]);
1855 float tgxmax = tgbin->GetX()[tgbin->GetN()-1];
1856 float tgyatxmax = tgbin->GetY()[tgbin->GetN()-1];
1857 float tgyatinf = tgyatxmax*0.001;
1858 float tgZeroPoint = (tgyatinf-tgyatxmax)/dlast + tgxmax;
1859 float tgLastPoint = min(tgZeroPoint, (
float)(
sqrts*1000.));
1861 for (
unsigned int j=1; j<=ninsert; j++){
1862 const double xwidth = (tgLastPoint-tgxmax)/((
double)ninsert);
1863 double xtmp = tgxmax+xwidth*j;
1864 double ytmp = tgyatxmax + dlast*(xtmp-tgxmax);
1867 xtmp = min(xtmp, (
double)(
sqrts*1000.));
1868 cout <<
"Inserting (x,y) = ( " << xtmp <<
" , " << ytmp <<
" )" << endl;
1870 for (
unsigned int k=0; k<1; k++){
1871 if (xtmp+xwidth<=
sqrts*1000.){
1873 cout <<
"Inserting (x,y) = ( " << xtmp <<
" , " << ytmp <<
" )" << endl;
1880 cout <<
"Inserting (x,y) = ( " << xtmp <<
" , " << ytmp <<
" )" << endl;
1885 for (
unsigned int j=1; j<=ninsert; j++){
1886 const double xwidth = (
sqrts*1000.-tgLastPoint)/((
double)ninsert);
1887 double xtmp = tgLastPoint+xwidth*j;
1888 double ytmp = tgyatinf;
1889 cout <<
"Inserting (x,y) = ( " << xtmp <<
" , " << ytmp <<
" )" << endl;
1892 tgbin->SetMarkerStyle(20);
1893 tgbin->SetMarkerSize(1.2);
1894 tgbin->SetLineColor(kBlack);
1895 tgbin->SetLineWidth(2);
1897 RooRealVar* reco_minus_true_var =
new RooRealVar(
"plotVar",
"", minval, maxval);
1898 reco_minus_true_var->setBins((
int)(maxval-minval));
1902 double ymin=9e9, ymax=-9e9;
1907 RooRealIntegral* pdf_int =
new RooRealIntegral(
"pdf_int",
"", *spPDF_bin, RooArgSet(*reco_minus_true_var));
1909 double pdfintval = pdf_int->getVal();
1910 cout <<
"PDF integral: " << pdfintval << endl;
1911 cout <<
"Sum of weights: " << h_RecoMinusTrueVals->Integral(
"width") << endl;
1912 double riemannintval=0;
1913 for (
int ip=0; ip<tgbin->GetN(); ip++){
1914 tgbin->GetY()[ip] /= pdfintval;
1915 tgbin->GetEY()[ip] /= pdfintval;
1916 ymin = min(ymin, tgbin->GetY()[ip]);
1917 ymax = max(ymax, tgbin->GetY()[ip]);
1919 if (ip!=tgbin->GetN()-1) riemannintval += (tgbin->GetY()[ip]+tgbin->GetY()[ip-1])*(tgbin->GetX()[ip]-tgbin->GetX()[ip-1])*0.5;
1921 cout <<
"Post-normalization Riemann integral: " << riemannintval << endl;
1924 spPDF_bin = fac_sp->
getFunc();
1927 pdf_int =
new RooRealIntegral(
"pdf_int",
"", *spPDF_bin, RooArgSet(*reco_minus_true_var));
1929 spPDF_bin->forceNumInt(
true);
1931 pdf_int =
new RooRealIntegral(
"pdf_int",
"", *spPDF_bin, RooArgSet(*reco_minus_true_var));
1935 RooPlot* plot_bin = reco_minus_true_var->frame();
1936 cout <<
"Frame created for [" << reco_minus_true_var->getMin() <<
" , " << reco_minus_true_var->getMax() <<
" , " << reco_minus_true_var->getBins() <<
"]" << endl;
1937 plot_bin->GetXaxis()->CenterTitle();
1938 plot_bin->GetYaxis()->SetTitleOffset(1.2);
1939 plot_bin->GetYaxis()->CenterTitle();
1940 plot_bin->GetXaxis()->SetTitle(tgbin->GetXaxis()->GetTitle());
1941 plot_bin->GetYaxis()->SetTitle(tgbin->GetYaxis()->GetTitle());
1942 plot_bin->SetTitle(tgbin->GetTitle()); tgbin->SetTitle(
"");
1944 spPDF_bin->plotOn(plot_bin, LineColor(kRed), LineWidth(2), LineStyle(1));
1945 plot_bin->GetYaxis()->SetRangeUser(ymin, ymax*1.2);
1946 plot_bin->GetXaxis()->SetNdivisions(-505);
1948 TCanvas* canvas =
new TCanvas(Form(
"c_%s", h_RecoMinusTrueVals->GetName()),
"", 800, 800);
1950 tgbin->Draw(
"e1psame");
1953 if (debug) foutput->WriteTObject(canvas);
1958 reco_minus_true_var->setRange(-1, reco_minus_true_var->getMax()/avgbin_varTrue-1.);
1959 for (
int ip=0; ip<tgbin->GetN(); ip++){
1960 tgbin->GetY()[ip] *= avgbin_varTrue;
1961 tgbin->GetEY()[ip] *= avgbin_varTrue;
1962 tgbin->GetX()[ip] = tgbin->GetX()[ip]/avgbin_varTrue-1.;
1963 tgbin->GetEX()[ip] /= avgbin_varTrue;
1965 tglist.push_back(tgbin);
1968 delete reco_minus_true_var;
1969 delete p_RecoMinusTrueVals;
1970 delete h_RecoMinusTrueVals;
1973 RooArgList splineList;
1974 vector<MELALinearInterpFunc::T> truthList;
1975 vector<MELANCSplineFactory_1D*> facList;
1977 RooRealVar finalX(
"true",
"", 0.1,
sqrts*1000.);
1978 finalX.setBins((
int)(finalX.getMax()-finalX.getMin())/50);
1979 RooRealVar finalY(
"reco",
"", finalX.getMin(), finalX.getMax());
1980 finalY.setBins(finalX.getBins());
1981 RooFormulaVar finalYoverX(
"RecoOverTrueVar",
"@1/max(0.1, @0)-1.", RooArgList(finalX, finalY));
1983 for (
int is=-1; is<=
int(tglist.size()); is++){
1986 if (is>=0 && is<
int(tglist.size())){
1987 tval = hvar->GetBinContent(is+1);
1988 tg = (TGraphErrors*)tglist[is]->Clone(Form(
"tg_bin%i", is+1));
1991 tval=finalX.getMin();
1992 tg = (TGraphErrors*)tglist[is+1]->Clone(Form(
"tg_bin%i", is+1));
1995 tval = finalX.getMax();
1996 tg = (TGraphErrors*)tglist[is-1]->Clone(Form(
"tg_bin%i", is+1));
1998 TSpline3 sptmp(
"sptmp", tg,
"b2e2", 0, 0);
2001 for (
int ip=0; ip<tg->GetN(); ip++){
if (tg->GetX()[ip]<0.) n++; }
2002 for (
unsigned int i=0; i<2; i++) xy[i]=
new double[n];
2003 for (
int ip=0; ip<n; ip++){
2004 if (ip==n-1){ xy[0][ip]=0; xy[1][ip]=sptmp.Eval(xy[0][ip]); }
2005 else{ xy[0][ip]=tg->GetX()[ip]; xy[1][ip]=tg->GetY()[ip]; }
2008 tg =
new TGraph(n, xy[0], xy[1]);
2009 tg->SetName(Form(
"tg_bin%i", is+1));
2011 for (
unsigned int i=0; i<2; i++)
delete[] xy[i];
2013 tg->SetName(Form(
"tg_bin%i", is+1));
2014 tg->SetTitle(Form(
"Truth at %.2f", tval));
2015 tg->GetXaxis()->SetTitle(
"Reco/Truth-1");
2016 tg->GetYaxis()->SetTitle(
"Probability");
2023 finalX.setVal(tval);
2024 RooRealIntegral sp_int = RooRealIntegral(
"pdftmp_int",
"", *functmp, RooArgSet(finalY));
2025 double pdftmpintval = sp_int.getVal();
2026 cout <<
"Bin " << is+1 <<
" at truth = " << tval <<
" renormalized by 1/" << pdftmpintval << endl;
2027 for (
int ip=0; ip<tg->GetN(); ip++){ tg->GetY()[ip] /= pdftmpintval; }
2033 splineList.add(*pdftmp);
2034 truthList.push_back(tval);
2035 facList.push_back(fac_tmp);
2037 foutput->WriteTObject(tg);
2046 finalX.setRange(35, 535);
2047 finalY.setRange(35, 535);
2048 finalX.setBins(finalX.getMax()-finalX.getMin());
2049 finalY.setBins(finalY.getMax()-finalY.getMin());
2050 TH2F* h_frm =
new TH2F(
"hFinalResolution",
"", finalX.getBins(), finalX.getMin(), finalX.getMax(), finalY.getBins(), finalY.getMin(), finalY.getMax());
2051 for (
int ix=1; ix<=h_frm->GetNbinsX(); ix++){
2052 double xval = h_frm->GetXaxis()->GetBinCenter(ix);
2053 finalX.setVal(xval);
2054 for (
int iy=1; iy<=h_frm->GetNbinsY(); iy++){
2055 double yval = h_frm->GetYaxis()->GetBinCenter(iy);
2056 finalY.setVal(yval);
2057 h_frm->SetBinContent(ix, iy, thePdf->getVal());
2061 h_frm->GetXaxis()->SetTitle(
"True m_{JJ} (GeV)");
2062 h_frm->GetYaxis()->SetTitle(
"Reco. m_{JJ} (GeV)");
2063 h_frm->SetTitle(TString(
"m_{JJ} resolution for ")+strprod);
2064 canvas =
new TCanvas(TString(
"cFinalResolution")+strprod,
"", 800, 800);
2065 h_frm->Draw(
"colz");
2068 canvas->SaveAs(Form(
"%s%s", canvas->GetName(),
".pdf"));
2070 if (debug) foutput->WriteTObject(h_frm);
2074 if (strprod==
"WH" || strprod==
"ZH"){
2075 RooRealVar* mV =
new RooRealVar(
"mV",
"", 0, 1000.*
sqrts);
2076 RooRealVar* GaV =
new RooRealVar(
"GaV",
"", 0, 1000.*
sqrts);
2078 mV->setVal(
mela.getPrimaryMass(23));
2079 GaV->setVal(
mela.getPrimaryWidth(23));
2082 mV->setVal(
mela.getPrimaryMass(24));
2083 GaV->setVal(
mela.getPrimaryWidth(24));
2085 mV->setConstant(
true);
2086 GaV->setConstant(
true);
2087 cout <<
"Constructing a BW on top of the resolution. mV = " << mV->getVal() <<
", GammaV = " << GaV->getVal() << endl;
2090 RooGenericPdf* trueBW =
new RooGenericPdf(
"trueBW",
"2.*@0/(pow(pow(@0,2)-pow(@1,2),2) + pow(@1*@2,2))", RooArgList(finalX, *mV, *GaV));
2091 vector<float> samplePointsList;
2093 vector<pair<double, double>> generateRanges;
2094 vector<pair<int, int>> generateReqStep;
2095 generateRanges.push_back(pair<double, double>(0.2, 64.9)); generateReqStep.push_back(pair<int, int>(100, 20));
2096 generateRanges.push_back(pair<double, double>(65., mV->getVal()-0.05)); generateReqStep.push_back(pair<int, int>(100, 10));
2097 generateRanges.push_back(pair<double, double>(mV->getVal()+0.05, 120.)); generateReqStep.push_back(pair<int, int>(100, 10));
2098 generateRanges.push_back(pair<double, double>(120.1, 200.)); generateReqStep.push_back(pair<int, int>(25, 40));
2099 generateRanges.push_back(pair<double, double>(200.1, 600.)); generateReqStep.push_back(pair<int, int>(25, 40));
2100 generateRanges.push_back(pair<double, double>(600.1, 1500.)); generateReqStep.push_back(pair<int, int>(25, 40));
2101 generateRanges.push_back(pair<double, double>(1500.1, 4000.)); generateReqStep.push_back(pair<int, int>(25, 40));
2102 generateRanges.push_back(pair<double, double>(4000.1, 1000.*
sqrts-0.1)); generateReqStep.push_back(pair<int, int>(50, 20));
2103 for (
unsigned int ig=0; ig<generateRanges.size(); ig++){
2104 int& nrequested = generateReqStep[ig].first;
2105 int& nstep = generateReqStep[ig].second;
2106 double& rmin = generateRanges[ig].first;
2107 double&
rmax = generateRanges[ig].second;
2108 finalX.setRange(rmin,
rmax);
2109 RooDataSet* dsSamplePoints = trueBW->generate(finalX, nrequested*nstep);
2110 vector<float> tmpList;
2111 for (
int ev=0;
ev<dsSamplePoints->numEntries();
ev++){
2112 float val =
dynamic_cast<RooAbsReal*
>(dsSamplePoints->get(
ev)->find(finalX))->getVal();
2115 delete dsSamplePoints;
2116 for (
int ev=0;
ev<(
int)tmpList.size();
ev+=nstep) samplePointsList.push_back(tmpList.at(
ev));
2118 finalX.setRange(0.1, 1000.*
sqrts);
2119 finalY.setRange(0.1, 1000.*
sqrts);
2120 addByLowest(samplePointsList,
float(finalX.getMin()),
true);
2121 addByLowest(samplePointsList,
float(finalX.getMax()),
true);
2123 vector<pair<double, double>> idealBWpoints;
2124 for (
unsigned int ev=0;
ev<samplePointsList.size();
ev++){
2125 finalX.setVal(samplePointsList[
ev]);
2126 idealBWpoints.push_back(pair<double, double>(finalX.getVal(), trueBW->getVal()));
2127 cout <<
"Sampled point " << idealBWpoints.at(
ev).first <<
" , " << idealBWpoints.at(
ev).second <<
" from ideal BW." << endl;
2129 vector<pair<double, double>> recoBWpoints;
2131 for (
unsigned int ev=0;
ev<samplePointsList.size();
ev++){
2132 double recoval = samplePointsList[
ev];
2133 finalY.setVal(recoval);
2135 vector<pair<stype, stype>> idealBWconvpoints;
2136 for (
unsigned int ip=0; ip<idealBWpoints.size(); ip++){
2137 double& trueval = idealBWpoints.at(ip).first;
2138 finalX.setVal(trueval);
2139 double resoval = thePdf->getVal();
2140 idealBWconvpoints.push_back(pair<stype, stype>((stype)trueval, (stype)(idealBWpoints.at(ip).second*resoval)));
2143 truespline.setPoints(idealBWconvpoints);
2144 double recobwval = truespline.getPDF()->analyticalIntegral(2);
2145 recoBWpoints.push_back(pair<double, double>(recoval, recobwval));
2147 cout <<
"Conv. ideal BW integral at reco = " << recoval <<
" : " << recobwval << endl;
2161 tg_recoBW->SetMarkerStyle(20);
2162 tg_recoBW->SetMarkerSize(1.2);
2163 tg_recoBW->SetMarkerColor(kRed);
2164 tg_recoBW->SetLineColor(kRed);
2165 tg_recoBW->GetXaxis()->SetRangeUser(0.1, 250.);
2166 tg_recoBW->GetXaxis()->SetTitle(
"m_{JJ} (GeV)");
2167 tg_recoBW->GetYaxis()->SetTitle(
"BW (GeV^{-3})");
2168 tg_recoBW->SetTitle(TString(
"Reco. BW for ")+strprod);
2169 tg_idealBW->SetMarkerStyle(21);
2170 tg_idealBW->SetMarkerSize(0.8);
2171 tg_idealBW->SetMarkerColor(kBlack);
2172 tg_idealBW->SetLineColor(kBlack);
2173 tg_idealBW->GetXaxis()->SetRangeUser(0.1, 250.);
2174 tg_idealBW->GetXaxis()->SetTitle(
"m_{JJ} (GeV)");
2175 tg_idealBW->GetYaxis()->SetTitle(
"BW (GeV^{-3})");
2176 tg_idealBW->SetTitle(TString(
"True BW for ")+strprod);
2177 foutput->WriteTObject(tg_idealBW);
2178 foutput->WriteTObject(tg_recoBW);
2180 canvas =
new TCanvas(TString(
"cFinalBWComparison_")+strprod,
"", 800, 800);
2181 tg_idealBW->SetTitle(TString(
"BW for ")+strprod);
2182 tg_recoBW->SetTitle(
"");
2183 tg_idealBW->Draw(
"acp");
2184 tg_recoBW->Draw(
"cpsame");
2187 canvas->SaveAs(Form(
"%s%s", canvas->GetName(),
".pdf"));
2200 for (
unsigned int is=0; is<facList.size(); is++)
delete facList.at(is);
2201 for (
unsigned int is=0; is<tglist.size(); is++)
delete tglist.at(is);
2203 TGraphErrors* tg_avg_varReco =
makeGraphFromTH1(hvar, havg_varReco, Form(
"tg_%s", havg_varReco->GetName()));
2204 foutput->WriteTObject(tg_avg_varReco);
2205 delete tg_avg_varReco;
2207 for (
int ix=0; ix<=hvarTrue_varReco->GetNbinsX()+1; ix++){
2208 double integral = hvarTrue_varReco->Integral(ix, ix, 0, hvarTrue_varReco->GetNbinsY()+1);
2210 for (
int iy=0; iy<=hvarTrue_varReco->GetNbinsY()+1; iy++){
2211 double bincontent = hvarTrue_varReco->GetBinContent(ix, iy);
2212 double binerror = hvarTrue_varReco->GetBinError(ix, iy);
2213 bincontent /= integral;
2214 binerror /= integral;
2215 hvarTrue_varReco->SetBinContent(ix, iy, bincontent);
2216 hvarTrue_varReco->SetBinError(ix, iy, binerror);
2221 foutput->WriteTObject(hvar);
2222 foutput->WriteTObject(havg_varReco);
2223 foutput->WriteTObject(hvarTrue_varReco);
2224 if (writeFinalTree) foutput->WriteTObject(newtree);
2226 if (writeFinalTree)
delete newtree;
2227 delete hvarTrue_varReco;
2228 delete havg_varReco;
2232 for (
unsigned int f=0;
f<finputList.size();
f++) finputList.at(
f)->Close();