| 1 |
#include <fstream>
|
| 2 |
#include <iostream>
|
| 3 |
#include <math.h>
|
| 4 |
#include <sstream>
|
| 5 |
#include <stdlib.h>
|
| 6 |
#include "TreeProducer.h"
|
| 7 |
#include "TwoBodyDecaySpy.h"
|
| 8 |
#include "histograms.h"
|
| 9 |
#include <TLorentzVector.h>
|
| 10 |
|
| 11 |
#define FITMODE 0
|
| 12 |
|
| 13 |
|
| 14 |
double cos(double a1, double b1, double c1, double a2, double b2, double c2)
|
| 15 |
{
|
| 16 |
double n1 = sqrt(a1*a1+b1*b1+c1*c1);
|
| 17 |
double n2 = sqrt(a2*a2+b2*b2+c2*c2);
|
| 18 |
double product = a1*a2+b1*b2+c1*c2;
|
| 19 |
return product/(n1*n2);
|
| 20 |
}
|
| 21 |
|
| 22 |
|
| 23 |
double stringTodouble(const std::string& word)
|
| 24 |
{
|
| 25 |
std::stringstream str;
|
| 26 |
str<<word;
|
| 27 |
double a;
|
| 28 |
str >> a;
|
| 29 |
return a;
|
| 30 |
}
|
| 31 |
|
| 32 |
double massInvariant(const double& px1, const double& py1, const double& pz1,
|
| 33 |
const double& px2, const double& py2, const double& pz2)
|
| 34 |
{
|
| 35 |
double p1=sqrt(px1*px1+py1*py1+pz1*pz1);
|
| 36 |
double p2=sqrt(px2*px2+py2*py2+pz2*pz2);
|
| 37 |
|
| 38 |
|
| 39 |
double m_muon=105.658e-3;
|
| 40 |
|
| 41 |
double e1 = sqrt(pow(m_muon,2)+pow(p1,2));
|
| 42 |
double e2 = sqrt(pow(m_muon,2)+pow(p2,2));
|
| 43 |
double E2 = pow(e1+e2,2);
|
| 44 |
double P2 = pow(px1+px2,2)
|
| 45 |
+ pow(py1+py2,2)
|
| 46 |
+ pow(pz1+pz2,2);
|
| 47 |
return sqrt(E2-P2);
|
| 48 |
}
|
| 49 |
|
| 50 |
double pT(const double& px1, const double& py1)
|
| 51 |
{
|
| 52 |
double PT = sqrt( pow(px1,2)
|
| 53 |
+ pow(py1,2));
|
| 54 |
return PT;
|
| 55 |
}
|
| 56 |
|
| 57 |
double p(const double& px1, const double& py1, const double& pz1,
|
| 58 |
const double& px2, const double& py2, const double& pz2)
|
| 59 |
{
|
| 60 |
double P = sqrt( pow(px1+px2,2)
|
| 61 |
+ pow(py1+py2,2)
|
| 62 |
+ pow(pz1+pz2,2));
|
| 63 |
return P;
|
| 64 |
}
|
| 65 |
|
| 66 |
|
| 67 |
//################################ Eric and Jeremy ############################################
|
| 68 |
|
| 69 |
/*
|
| 70 |
void LoopRef(TreeProducer* treeRef, HistoCollection& histos)
|
| 71 |
{
|
| 72 |
Long64_t nentries = treeRef->GetEntriesFast();
|
| 73 |
|
| 74 |
Long64_t nbytes1 = 0, nb1 = 0;
|
| 75 |
Long64_t nbytes2 = 0, nb2 = 0;
|
| 76 |
for (Long64_t jentry=0; jentry<nentries;jentry++)
|
| 77 |
{
|
| 78 |
Long64_t ientry1 = treeRef->LoadTree(jentry);
|
| 79 |
if (ientry1 < 0) break;
|
| 80 |
|
| 81 |
nb1 = treeRef->GetEntry(jentry); nbytes1 += nb1;
|
| 82 |
|
| 83 |
double m=massInvariant(treeRef->Track_px[0], treeRef->Track_py[0], treeRef->Track_pz[0], treeRef->Track_px[1], treeRef->Track_py[1], treeRef->Track_pz[1]);
|
| 84 |
double mgen=massInvariant(treeRef->TrackGen_px[0], treeRef->TrackGen_py[0], treeRef->TrackGen_pz[0],treeRef->TrackGen_px[1], treeRef->TrackGen_py[1], treeRef->TrackGen_pz[1]);
|
| 85 |
histos.DeltaM_None->Fill((m-mgen)/mgen);
|
| 86 |
|
| 87 |
histos.mref_histo->Fill( massInvariant(treeRef->Track_px[0], treeRef->Track_py[0], treeRef->Track_pz[0], treeRef->Track_px[1], treeRef->Track_py[1], treeRef->Track_pz[1]) );
|
| 88 |
|
| 89 |
|
| 90 |
if (treeRef->Track_charge[0]==treeRef->Track_charge[1]) std::cout << "aie" << std::endl;
|
| 91 |
if (treeRef->Track_charge[0]<0)
|
| 92 |
{
|
| 93 |
histos.mVSphi_m->Fill(treeRef->Track_phi[0],massInvariant(treeRef->Track_px[0], treeRef->Track_py[0], treeRef->Track_pz[0], treeRef->Track_px[1], treeRef->Track_py[1], treeRef->Track_pz[1]) );
|
| 94 |
}else
|
| 95 |
{
|
| 96 |
histos.mVSphi_p->Fill(treeRef->Track_phi[0],massInvariant(treeRef->Track_px[0], treeRef->Track_py[0], treeRef->Track_pz[0], treeRef->Track_px[1], treeRef->Track_py[1], treeRef->Track_pz[1]) );
|
| 97 |
}
|
| 98 |
|
| 99 |
if (treeRef->Track_charge[1]<0)
|
| 100 |
{
|
| 101 |
histos.mVSphi_m->Fill(treeRef->Track_phi[1],massInvariant(treeRef->Track_px[0], treeRef->Track_py[0], treeRef->Track_pz[0], treeRef->Track_px[1], treeRef->Track_py[1], treeRef->Track_pz[1]) );
|
| 102 |
}else
|
| 103 |
{
|
| 104 |
histos.mVSphi_p->Fill(treeRef->Track_phi[1],massInvariant(treeRef->Track_px[0], treeRef->Track_py[0], treeRef->Track_pz[0], treeRef->Track_px[1], treeRef->Track_py[1], treeRef->Track_pz[1]) );
|
| 105 |
}
|
| 106 |
|
| 107 |
|
| 108 |
histos.mgen_histo->Fill( massInvariant(treeRef->TrackGen_px[0], treeRef->TrackGen_py[0], treeRef->TrackGen_pz[0],treeRef->TrackGen_px[1], treeRef->TrackGen_py[1], treeRef->TrackGen_pz[1]) );
|
| 109 |
|
| 110 |
|
| 111 |
for (unsigned int i=0;i<2;i++)
|
| 112 |
{
|
| 113 |
unsigned int j=0;
|
| 114 |
double theMax1 = cos(treeRef->Track_px[i],treeRef->Track_py[i],treeRef->Track_pz[i],treeRef->TrackGen_px[0],treeRef->TrackGen_py[0],treeRef->TrackGen_pz[0]);
|
| 115 |
double theMax2 = cos(treeRef->Track_px[i],treeRef->Track_py[i],treeRef->Track_pz[i],treeRef->TrackGen_px[1],treeRef->TrackGen_py[1],treeRef->TrackGen_pz[1]);
|
| 116 |
if (theMax1>theMax2) j=0; else j=1;
|
| 117 |
|
| 118 |
double pt = sqrt(pow(treeRef->Track_px[i],2)+ pow(treeRef->Track_py[i],2));
|
| 119 |
double ptgen = sqrt(pow(treeRef->TrackGen_px[j],2) +pow(treeRef->TrackGen_py[j],2));
|
| 120 |
|
| 121 |
double p = sqrt(pow(treeRef->Track_px[i],2)+ pow(treeRef->Track_py[i],2) + pow(treeRef->Track_pz[i],2));
|
| 122 |
double pgen = sqrt(pow(treeRef->TrackGen_px[j],2) +pow(treeRef->TrackGen_py[j],2) +pow(treeRef->TrackGen_pz[j],2));
|
| 123 |
|
| 124 |
histos.rptref_histo->Fill((ptgen-pt)/ptgen);
|
| 125 |
histos.rpref_histo->Fill((pgen-p)/pgen);
|
| 126 |
}
|
| 127 |
|
| 128 |
|
| 129 |
}
|
| 130 |
}
|
| 131 |
|
| 132 |
|
| 133 |
|
| 134 |
|
| 135 |
void Loop(TreeProducer* tree, HistoCollection& histos)
|
| 136 |
{
|
| 137 |
Long64_t nentries = tree->GetEntriesFast();
|
| 138 |
|
| 139 |
Long64_t nbytes1 = 0, nb1 = 0;
|
| 140 |
Long64_t nbytes2 = 0, nb2 = 0;
|
| 141 |
for (Long64_t jentry=0; jentry<nentries;jentry++)
|
| 142 |
{
|
| 143 |
Long64_t ientry2 = tree->LoadTree(jentry);
|
| 144 |
if (ientry2 < 0) break;
|
| 145 |
|
| 146 |
nb2 = tree->GetEntry(jentry); nbytes2 += nb2;
|
| 147 |
|
| 148 |
|
| 149 |
|
| 150 |
double m=massInvariant(tree->Track_px[0], tree->Track_py[0], tree->Track_pz[0], tree->Track_px[1], tree->Track_py[1], tree->Track_pz[1]);
|
| 151 |
double mgen=massInvariant(tree->TrackGen_px[0], tree->TrackGen_py[0], tree->TrackGen_pz[0],tree->TrackGen_px[1], tree->TrackGen_py[1], tree->TrackGen_pz[1]);
|
| 152 |
histos.DeltaM_Refit->Fill((m-mgen)/mgen);
|
| 153 |
|
| 154 |
histos.m_histo->Fill( massInvariant(tree->Track_px[0], tree->Track_py[0], tree->Track_pz[0],
|
| 155 |
tree->Track_px[1], tree->Track_py[1], tree->Track_pz[1]) );
|
| 156 |
|
| 157 |
|
| 158 |
|
| 159 |
for (unsigned int i=0;i<2;i++)
|
| 160 |
{
|
| 161 |
unsigned int j=0;
|
| 162 |
double theMax1 = cos(tree->Track_px[i],tree->Track_py[i],tree->Track_pz[i],tree->TrackGen_px[0],tree->TrackGen_py[0],tree->TrackGen_pz[0]);
|
| 163 |
double theMax2 = cos(tree->Track_px[i],tree->Track_py[i],tree->Track_pz[i],tree->TrackGen_px[1],tree->TrackGen_py[1],tree->TrackGen_pz[1]);
|
| 164 |
if (theMax1>theMax2) j=0; else j=1;
|
| 165 |
|
| 166 |
double pt = sqrt(pow(tree->Track_px[i],2)+ pow(tree->Track_py[i],2));
|
| 167 |
double ptgen = sqrt(pow(tree->TrackGen_px[j],2) +pow(tree->TrackGen_py[j],2));
|
| 168 |
double p = sqrt(pow(tree->Track_px[i],2)+ pow(tree->Track_py[i],2) + pow(tree->Track_pz[i],2));
|
| 169 |
double pgen = sqrt(pow(tree->TrackGen_px[j],2) +pow(tree->TrackGen_py[j],2) +pow(tree->TrackGen_pz[j],2));
|
| 170 |
|
| 171 |
histos.rpt_histo->Fill((ptgen-pt)/ptgen);
|
| 172 |
histos.rp_histo->Fill((pgen-p)/pgen);
|
| 173 |
}
|
| 174 |
|
| 175 |
|
| 176 |
}
|
| 177 |
|
| 178 |
}
|
| 179 |
|
| 180 |
|
| 181 |
|
| 182 |
|
| 183 |
void LoopTBD(TwoBodyDecaySpy* tree, HistoCollection& histos)
|
| 184 |
{
|
| 185 |
Long64_t nentries = tree->GetEntriesFast();
|
| 186 |
|
| 187 |
Long64_t nbytes1 = 0, nb1 = 0;
|
| 188 |
Long64_t nbytes2 = 0, nb2 = 0;
|
| 189 |
for (Long64_t jentry=0; jentry<nentries;jentry++)
|
| 190 |
{
|
| 191 |
Long64_t ientry2 = tree->LoadTree(jentry);
|
| 192 |
if (ientry2 < 0) break;
|
| 193 |
|
| 194 |
nb2 = tree->GetEntry(jentry); nbytes2 += nb2;
|
| 195 |
|
| 196 |
histos.pp_histo->Fill(tree->NewTrack_Eric_p[0]);
|
| 197 |
histos.pm_histo->Fill(tree->NewTrack_Eric_p[1]);
|
| 198 |
histos.pall_histo->Fill(tree->NewTrack_Eric_p[0]);
|
| 199 |
histos.pall_histo->Fill(tree->NewTrack_Eric_p[1]);
|
| 200 |
|
| 201 |
histos.ep_histo->Fill(tree->NewTrack_Eric_pErr[0]/tree->NewTrack_Eric_p[0]);
|
| 202 |
histos.em_histo->Fill(tree->NewTrack_Eric_pErr[1]/tree->NewTrack_Eric_p[1]);
|
| 203 |
histos.eall_histo->Fill(tree->NewTrack_Eric_pErr[0]/tree->NewTrack_Eric_p[0]);
|
| 204 |
histos.eall_histo->Fill(tree->NewTrack_Eric_pErr[1]/tree->NewTrack_Eric_p[1]);
|
| 205 |
histos.mass_histo->Fill(tree->Zmass);
|
| 206 |
|
| 207 |
histos.massrec_histo->Fill(massInvariant(tree->NewTrack_px[0],tree->NewTrack_py[0],tree->NewTrack_pz[0],tree->NewTrack_px[1],tree->NewTrack_py[1],tree->NewTrack_pz[1]));
|
| 208 |
histos.mbefore_histo->Fill(massInvariant(tree->OldTrack_px[0],tree->OldTrack_py[0],tree->OldTrack_pz[0],tree->OldTrack_px[1],tree->OldTrack_py[1],tree->OldTrack_pz[1]));
|
| 209 |
|
| 210 |
}
|
| 211 |
|
| 212 |
}
|
| 213 |
|
| 214 |
|
| 215 |
void LoopCommon(TreeProducer* tree, TwoBodyDecaySpy * treeTBD, HistoCollection& histos)
|
| 216 |
{
|
| 217 |
Long64_t nentries = tree->GetEntriesFast();
|
| 218 |
|
| 219 |
Long64_t nbytes1 = 0, nb1 = 0;
|
| 220 |
Long64_t nbytes2 = 0, nb2 = 0;
|
| 221 |
for (Long64_t jentry=0; jentry<nentries;jentry++)
|
| 222 |
{
|
| 223 |
Long64_t ientry1 = treeTBD->LoadTree(jentry);
|
| 224 |
Long64_t ientry2 = tree->LoadTree(jentry);
|
| 225 |
if (ientry2 < 0) break;
|
| 226 |
if (ientry1 < 0) break;
|
| 227 |
|
| 228 |
nb1 = treeTBD->GetEntry(jentry); nbytes1 += nb1;
|
| 229 |
nb2 = tree->GetEntry(jentry); nbytes2 += nb2;
|
| 230 |
|
| 231 |
|
| 232 |
double m=massInvariant(treeTBD->NewTrack_px[0], treeTBD->NewTrack_py[0], treeTBD->NewTrack_pz[0], treeTBD->NewTrack_px[1], treeTBD->NewTrack_py[1], treeTBD->NewTrack_pz[1]);
|
| 233 |
double mgen=massInvariant(tree->TrackGen_px[0], tree->TrackGen_py[0], tree->TrackGen_pz[0],tree->TrackGen_px[1], tree->TrackGen_py[1], tree->TrackGen_pz[1]);
|
| 234 |
histos.DeltaM_TBD->Fill((m-mgen)/mgen);
|
| 235 |
|
| 236 |
|
| 237 |
std::vector<int> sign(2,0);
|
| 238 |
for (unsigned int i=0;i<2;i++)
|
| 239 |
{
|
| 240 |
double Max1 = cos(treeTBD->NewTrack_px[i],treeTBD->NewTrack_py[i],treeTBD->NewTrack_pz[i],tree->Track_px[0],tree->Track_py[0],tree->Track_pz[0]);
|
| 241 |
double Max2 = cos(treeTBD->NewTrack_px[i],treeTBD->NewTrack_py[i],treeTBD->NewTrack_pz[i],tree->Track_px[1],tree->Track_py[1],tree->Track_pz[1]);
|
| 242 |
if (Max1>Max2) sign[i]=tree->Track_charge[0]; else sign[i]=tree->Track_charge[1];
|
| 243 |
}
|
| 244 |
|
| 245 |
TVector3 vec1(treeTBD->NewTrack_px[0], treeTBD->NewTrack_py[0], treeTBD->NewTrack_pz[0]);
|
| 246 |
TVector3 vec2(treeTBD->NewTrack_px[1], treeTBD->NewTrack_py[1], treeTBD->NewTrack_pz[1]);
|
| 247 |
|
| 248 |
if (sign[0]<0)
|
| 249 |
{
|
| 250 |
histos.mVSphi_m_TBD->Fill(vec1.Phi(), massInvariant(treeTBD->NewTrack_px[0], treeTBD->NewTrack_py[0], treeTBD->NewTrack_pz[0], treeTBD->NewTrack_px[1], treeTBD->NewTrack_py[1], treeTBD->NewTrack_pz[1]) );
|
| 251 |
}
|
| 252 |
else
|
| 253 |
{
|
| 254 |
histos.mVSphi_p_TBD->Fill(vec1.Phi(), massInvariant(treeTBD->NewTrack_px[0], treeTBD->NewTrack_py[0], treeTBD->NewTrack_pz[0], treeTBD->NewTrack_px[1], treeTBD->NewTrack_py[1], treeTBD->NewTrack_pz[1]));
|
| 255 |
}
|
| 256 |
|
| 257 |
if (sign[1]<0)
|
| 258 |
{
|
| 259 |
histos.mVSphi_m_TBD->Fill(vec2.Phi(),massInvariant(treeTBD->NewTrack_px[0], treeTBD->NewTrack_py[0], treeTBD->NewTrack_pz[0], treeTBD->NewTrack_px[1], treeTBD->NewTrack_py[1], treeTBD->NewTrack_pz[1]) );
|
| 260 |
}else
|
| 261 |
{
|
| 262 |
histos.mVSphi_p_TBD->Fill(vec2.Phi(), massInvariant(treeTBD->NewTrack_px[0], treeTBD->NewTrack_py[0], treeTBD->NewTrack_pz[0], treeTBD->NewTrack_px[1], treeTBD->NewTrack_py[1], treeTBD->NewTrack_pz[1]) );
|
| 263 |
}
|
| 264 |
|
| 265 |
for (unsigned int i=0; i<2; i++)
|
| 266 |
{
|
| 267 |
unsigned int j=0;
|
| 268 |
double theMax1 = cos(treeTBD->NewTrack_px[i],treeTBD->NewTrack_py[i],treeTBD->NewTrack_pz[i],tree->TrackGen_px[0],tree->TrackGen_py[0],tree->TrackGen_pz[0]);
|
| 269 |
double theMax2 = cos(treeTBD->NewTrack_px[i],treeTBD->NewTrack_py[i],treeTBD->NewTrack_pz[i],tree->TrackGen_px[1],tree->TrackGen_py[1],tree->TrackGen_pz[1]);
|
| 270 |
if (theMax1>theMax2) j=0; else j=1;
|
| 271 |
|
| 272 |
double pt = sqrt(pow(treeTBD->NewTrack_px[i],2)+ pow(treeTBD->NewTrack_py[i],2));
|
| 273 |
double ptgen = sqrt(pow(tree->TrackGen_px[j],2) +pow(tree->TrackGen_py[j],2));
|
| 274 |
double p = sqrt(pow(treeTBD->NewTrack_px[i],2)+ pow(treeTBD->NewTrack_py[i],2) + pow(treeTBD->NewTrack_pz[i],2));
|
| 275 |
double pgen = sqrt(pow(tree->TrackGen_px[j],2) +pow(tree->TrackGen_py[j],2) +pow(tree->TrackGen_pz[j],2));
|
| 276 |
|
| 277 |
histos.rptTBD_histo->Fill((ptgen-pt)/ptgen);
|
| 278 |
histos.rpTBD_histo->Fill((pgen-p)/pgen);
|
| 279 |
}
|
| 280 |
}
|
| 281 |
|
| 282 |
}
|
| 283 |
*/
|
| 284 |
|
| 285 |
//############################# Christophe #################################################
|
| 286 |
|
| 287 |
|
| 288 |
void LoopRef(TreeProducer* treeRef, HistoCollection& histos)
|
| 289 |
{
|
| 290 |
Long64_t nentries = treeRef->GetEntriesFast();
|
| 291 |
|
| 292 |
Long64_t nbytes1 = 0, nb1 = 0;
|
| 293 |
Long64_t nbytes2 = 0, nb2 = 0;
|
| 294 |
for (Long64_t jentry=0; jentry<nentries;jentry++)
|
| 295 |
{
|
| 296 |
Long64_t ientry1 = treeRef->LoadTree(jentry);
|
| 297 |
if (ientry1 < 0) break;
|
| 298 |
|
| 299 |
nb1 = treeRef->GetEntry(jentry); nbytes1 += nb1;
|
| 300 |
|
| 301 |
histos.mRef_histo->Fill( massInvariant(treeRef->Track_px[0], treeRef->Track_py[0], treeRef->Track_pz[0], treeRef->Track_px[1], treeRef->Track_py[1], treeRef->Track_pz[1]) );
|
| 302 |
|
| 303 |
histos.mGen_histo->Fill( massInvariant(treeRef->TrackGen_px[0], treeRef->TrackGen_py[0], treeRef->TrackGen_pz[0],treeRef->TrackGen_px[1], treeRef->TrackGen_py[1], treeRef->TrackGen_pz[1]) );
|
| 304 |
|
| 305 |
|
| 306 |
for (unsigned int i=0;i<2;i++)
|
| 307 |
{
|
| 308 |
unsigned int j=0;
|
| 309 |
double theMax1 = cos(treeRef->Track_px[i],treeRef->Track_py[i],treeRef->Track_pz[i],treeRef->TrackGen_px[0],treeRef->TrackGen_py[0],treeRef->TrackGen_pz[0]);
|
| 310 |
double theMax2 = cos(treeRef->Track_px[i],treeRef->Track_py[i],treeRef->Track_pz[i],treeRef->TrackGen_px[1],treeRef->TrackGen_py[1],treeRef->TrackGen_pz[1]);
|
| 311 |
if (theMax1>theMax2) j=0; else j=1;
|
| 312 |
|
| 313 |
double pt = pT(treeRef->Track_px[i], treeRef->Track_py[i]);
|
| 314 |
double ptgen = pT(treeRef->TrackGen_px[j], treeRef->TrackGen_py[j]);
|
| 315 |
|
| 316 |
double p = sqrt(pow(treeRef->Track_px[i],2)+ pow(treeRef->Track_py[i],2) + pow(treeRef->Track_pz[i],2));
|
| 317 |
double pgen = sqrt(pow(treeRef->TrackGen_px[j],2) +pow(treeRef->TrackGen_py[j],2) +pow(treeRef->TrackGen_pz[j],2));
|
| 318 |
|
| 319 |
histos.rptRef_histo->Fill((ptgen-pt)/ptgen);
|
| 320 |
histos.rpRef_histo->Fill((pgen-p)/pgen);
|
| 321 |
}
|
| 322 |
|
| 323 |
|
| 324 |
}
|
| 325 |
}
|
| 326 |
|
| 327 |
|
| 328 |
void LoopVtx(TreeProducer* treeVtx, HistoCollection& histos)
|
| 329 |
{
|
| 330 |
Long64_t nentries = treeVtx->GetEntriesFast();
|
| 331 |
|
| 332 |
Long64_t nbytes1 = 0, nb1 = 0;
|
| 333 |
Long64_t nbytes2 = 0, nb2 = 0;
|
| 334 |
for (Long64_t jentry=0; jentry<nentries;jentry++)
|
| 335 |
{
|
| 336 |
Long64_t ientry1 = treeVtx->LoadTree(jentry);
|
| 337 |
if (ientry1 < 0) break;
|
| 338 |
|
| 339 |
nb1 = treeVtx->GetEntry(jentry); nbytes1 += nb1;
|
| 340 |
|
| 341 |
histos.mVtx_histo->Fill( massInvariant(treeVtx->Track_px[0], treeVtx->Track_py[0], treeVtx->Track_pz[0], treeVtx->Track_px[1], treeVtx->Track_py[1], treeVtx->Track_pz[1]) );
|
| 342 |
|
| 343 |
//histos.mGen_histo->Fill( massInvariant(treeVtx->TrackGen_px[0], treeVtx->TrackGen_py[0], treeVtx->TrackGen_pz[0],treeVtx->TrackGen_px[1], treeVtx->TrackGen_py[1], treeVtx->TrackGen_pz[1]) );
|
| 344 |
|
| 345 |
|
| 346 |
for (unsigned int i=0;i<2;i++)
|
| 347 |
{
|
| 348 |
unsigned int j=0;
|
| 349 |
double theMax1 = cos(treeVtx->Track_px[i],treeVtx->Track_py[i],treeVtx->Track_pz[i],treeVtx->TrackGen_px[0],treeVtx->TrackGen_py[0],treeVtx->TrackGen_pz[0]);
|
| 350 |
double theMax2 = cos(treeVtx->Track_px[i],treeVtx->Track_py[i],treeVtx->Track_pz[i],treeVtx->TrackGen_px[1],treeVtx->TrackGen_py[1],treeVtx->TrackGen_pz[1]);
|
| 351 |
if (theMax1>theMax2) j=0; else j=1;
|
| 352 |
|
| 353 |
double pt = pT(treeVtx->Track_px[i], treeVtx->Track_py[i]);
|
| 354 |
double ptgen = pT(treeVtx->TrackGen_px[j], treeVtx->TrackGen_py[j]);
|
| 355 |
|
| 356 |
double p = sqrt(pow(treeVtx->Track_px[i],2)+ pow(treeVtx->Track_py[i],2) + pow(treeVtx->Track_pz[i],2));
|
| 357 |
double pgen = sqrt(pow(treeVtx->TrackGen_px[j],2) +pow(treeVtx->TrackGen_py[j],2) +pow(treeVtx->TrackGen_pz[j],2));
|
| 358 |
|
| 359 |
histos.rptVtx_histo->Fill((ptgen-pt)/ptgen);
|
| 360 |
histos.rpVtx_histo->Fill((pgen-p)/pgen);
|
| 361 |
}
|
| 362 |
|
| 363 |
|
| 364 |
}
|
| 365 |
}
|
| 366 |
|
| 367 |
|
| 368 |
void LoopMom(TreeProducer* treeMom, HistoCollection& histos)
|
| 369 |
{
|
| 370 |
Long64_t nentries = treeMom->GetEntriesFast();
|
| 371 |
|
| 372 |
Long64_t nbytes1 = 0, nb1 = 0;
|
| 373 |
Long64_t nbytes2 = 0, nb2 = 0;
|
| 374 |
for (Long64_t jentry=0; jentry<nentries;jentry++)
|
| 375 |
{
|
| 376 |
Long64_t ientry1 = treeMom->LoadTree(jentry);
|
| 377 |
if (ientry1 < 0) break;
|
| 378 |
|
| 379 |
nb1 = treeMom->GetEntry(jentry); nbytes1 += nb1;
|
| 380 |
|
| 381 |
histos.mMom_histo->Fill( massInvariant(treeMom->Track_px[0], treeMom->Track_py[0], treeMom->Track_pz[0], treeMom->Track_px[1], treeMom->Track_py[1], treeMom->Track_pz[1]) );
|
| 382 |
|
| 383 |
//histos.mGen_histo->Fill( massInvariant(treeMom->TrackGen_px[0], treeMom->TrackGen_py[0], treeMom->TrackGen_pz[0],treeMom->TrackGen_px[1], treeMom->TrackGen_py[1], treeMom->TrackGen_pz[1]) );
|
| 384 |
|
| 385 |
|
| 386 |
for (unsigned int i=0;i<2;i++)
|
| 387 |
{
|
| 388 |
unsigned int j=0;
|
| 389 |
double theMax1 = cos(treeMom->Track_px[i],treeMom->Track_py[i],treeMom->Track_pz[i],treeMom->TrackGen_px[0],treeMom->TrackGen_py[0],treeMom->TrackGen_pz[0]);
|
| 390 |
double theMax2 = cos(treeMom->Track_px[i],treeMom->Track_py[i],treeMom->Track_pz[i],treeMom->TrackGen_px[1],treeMom->TrackGen_py[1],treeMom->TrackGen_pz[1]);
|
| 391 |
if (theMax1>theMax2) j=0; else j=1;
|
| 392 |
|
| 393 |
double pt = pT(treeMom->Track_px[i], treeMom->Track_py[i]);
|
| 394 |
double ptgen = pT(treeMom->TrackGen_px[j], treeMom->TrackGen_py[j]);
|
| 395 |
|
| 396 |
double p = sqrt(pow(treeMom->Track_px[i],2)+ pow(treeMom->Track_py[i],2) + pow(treeMom->Track_pz[i],2));
|
| 397 |
double pgen = sqrt(pow(treeMom->TrackGen_px[j],2) +pow(treeMom->TrackGen_py[j],2) +pow(treeMom->TrackGen_pz[j],2));
|
| 398 |
|
| 399 |
histos.rptMom_histo->Fill((ptgen-pt)/ptgen);
|
| 400 |
histos.rpMom_histo->Fill((pgen-p)/pgen);
|
| 401 |
}
|
| 402 |
|
| 403 |
|
| 404 |
}
|
| 405 |
}
|
| 406 |
|
| 407 |
|
| 408 |
void LoopFull(TreeProducer* treeFull, HistoCollection& histos)
|
| 409 |
{
|
| 410 |
Long64_t nentries = treeFull->GetEntriesFast();
|
| 411 |
|
| 412 |
Long64_t nbytes1 = 0, nb1 = 0;
|
| 413 |
Long64_t nbytes2 = 0, nb2 = 0;
|
| 414 |
for (Long64_t jentry=0; jentry<nentries;jentry++)
|
| 415 |
{
|
| 416 |
Long64_t ientry1 = treeFull->LoadTree(jentry);
|
| 417 |
if (ientry1 < 0) break;
|
| 418 |
|
| 419 |
nb1 = treeFull->GetEntry(jentry); nbytes1 += nb1;
|
| 420 |
|
| 421 |
histos.mFull_histo->Fill( massInvariant(treeFull->Track_px[0], treeFull->Track_py[0], treeFull->Track_pz[0], treeFull->Track_px[1], treeFull->Track_py[1], treeFull->Track_pz[1]) );
|
| 422 |
|
| 423 |
//histos.mGen_histo->Fill( massInvariant(treeFull->TrackGen_px[0], treeFull->TrackGen_py[0], treeFull->TrackGen_pz[0],treeFull->TrackGen_px[1], treeFull->TrackGen_py[1], treeFull->TrackGen_pz[1]) );
|
| 424 |
|
| 425 |
|
| 426 |
for (unsigned int i=0;i<2;i++)
|
| 427 |
{
|
| 428 |
unsigned int j=0;
|
| 429 |
double theMax1 = cos(treeFull->Track_px[i],treeFull->Track_py[i],treeFull->Track_pz[i],treeFull->TrackGen_px[0],treeFull->TrackGen_py[0],treeFull->TrackGen_pz[0]);
|
| 430 |
double theMax2 = cos(treeFull->Track_px[i],treeFull->Track_py[i],treeFull->Track_pz[i],treeFull->TrackGen_px[1],treeFull->TrackGen_py[1],treeFull->TrackGen_pz[1]);
|
| 431 |
if (theMax1>theMax2) j=0; else j=1;
|
| 432 |
|
| 433 |
double pt = pT(treeFull->Track_px[i], treeFull->Track_py[i]);
|
| 434 |
double ptgen = pT(treeFull->TrackGen_px[j], treeFull->TrackGen_py[j]);
|
| 435 |
|
| 436 |
double p = sqrt(pow(treeFull->Track_px[i],2)+ pow(treeFull->Track_py[i],2) + pow(treeFull->Track_pz[i],2));
|
| 437 |
double pgen = sqrt(pow(treeFull->TrackGen_px[j],2) +pow(treeFull->TrackGen_py[j],2) +pow(treeFull->TrackGen_pz[j],2));
|
| 438 |
|
| 439 |
histos.rptFull_histo->Fill((ptgen-pt)/ptgen);
|
| 440 |
histos.rpFull_histo->Fill((pgen-p)/pgen);
|
| 441 |
}
|
| 442 |
|
| 443 |
|
| 444 |
}
|
| 445 |
}
|
| 446 |
|
| 447 |
|
| 448 |
void LoopKinFit(TreeProducer* treeKinFit, HistoCollection& histos)
|
| 449 |
{
|
| 450 |
Long64_t nentries = treeKinFit->GetEntriesFast();
|
| 451 |
|
| 452 |
Long64_t nbytes1 = 0, nb1 = 0;
|
| 453 |
Long64_t nbytes2 = 0, nb2 = 0;
|
| 454 |
for (Long64_t jentry=0; jentry<nentries;jentry++)
|
| 455 |
{
|
| 456 |
Long64_t ientry1 = treeKinFit->LoadTree(jentry);
|
| 457 |
if (ientry1 < 0) break;
|
| 458 |
|
| 459 |
nb1 = treeKinFit->GetEntry(jentry); nbytes1 += nb1;
|
| 460 |
|
| 461 |
histos.mKinFit_histo->Fill( massInvariant(treeKinFit->Track_px[0], treeKinFit->Track_py[0], treeKinFit->Track_pz[0], treeKinFit->Track_px[1], treeKinFit->Track_py[1], treeKinFit->Track_pz[1]) );
|
| 462 |
|
| 463 |
//histos.mGen_histo->Fill( massInvariant(treeKinFit->TrackGen_px[0], treeKinFit->TrackGen_py[0], treeKinFit->TrackGen_pz[0],treeKinFit->TrackGen_px[1], treeKinFit->TrackGen_py[1], treeKinFit->TrackGen_pz[1]) );
|
| 464 |
|
| 465 |
|
| 466 |
for (unsigned int i=0;i<2;i++)
|
| 467 |
{
|
| 468 |
unsigned int j=0;
|
| 469 |
double theMax1 = cos(treeKinFit->Track_px[i],treeKinFit->Track_py[i],treeKinFit->Track_pz[i],treeKinFit->TrackGen_px[0],treeKinFit->TrackGen_py[0],treeKinFit->TrackGen_pz[0]);
|
| 470 |
double theMax2 = cos(treeKinFit->Track_px[i],treeKinFit->Track_py[i],treeKinFit->Track_pz[i],treeKinFit->TrackGen_px[1],treeKinFit->TrackGen_py[1],treeKinFit->TrackGen_pz[1]);
|
| 471 |
if (theMax1>theMax2) j=0; else j=1;
|
| 472 |
|
| 473 |
double pt = pT(treeKinFit->Track_px[i], treeKinFit->Track_py[i]);
|
| 474 |
double ptgen = pT(treeKinFit->TrackGen_px[j], treeKinFit->TrackGen_py[j]);
|
| 475 |
|
| 476 |
double p = sqrt(pow(treeKinFit->Track_px[i],2)+ pow(treeKinFit->Track_py[i],2) + pow(treeKinFit->Track_pz[i],2));
|
| 477 |
double pgen = sqrt(pow(treeKinFit->TrackGen_px[j],2) +pow(treeKinFit->TrackGen_py[j],2) +pow(treeKinFit->TrackGen_pz[j],2));
|
| 478 |
|
| 479 |
histos.rptKinFit_histo->Fill((ptgen-pt)/ptgen);
|
| 480 |
histos.rpKinFit_histo->Fill((pgen-p)/pgen);
|
| 481 |
}
|
| 482 |
|
| 483 |
|
| 484 |
}
|
| 485 |
}
|
| 486 |
|
| 487 |
|
| 488 |
int init()
|
| 489 |
{
|
| 490 |
|
| 491 |
return 0;
|
| 492 |
}
|
| 493 |
|
| 494 |
|
| 495 |
int main(int argc, char *argv[])
|
| 496 |
{
|
| 497 |
std::cout << "Analysis begins ...." << std::endl;
|
| 498 |
std::cout << "Check argument ..." << std::endl;
|
| 499 |
|
| 500 |
|
| 501 |
if (argc!=2)
|
| 502 |
{
|
| 503 |
std::cout << "ERROR : wrong number of arguments" << std::endl;
|
| 504 |
std::cout << "CORRECT SYNTAX : ./analysis input.root" << std::endl;
|
| 505 |
return 1;
|
| 506 |
}
|
| 507 |
|
| 508 |
std::string filename = argv[1];
|
| 509 |
std::cout << "Input file : " << filename << std::endl;
|
| 510 |
|
| 511 |
/// First Fit (without constraint)
|
| 512 |
TreeProducer treeRef(filename,"treeproducer1/ttree_NoAlgo");
|
| 513 |
std::cout << "N_entries for treeRef :" << treeRef.GetEntriesFast() << std::endl;
|
| 514 |
|
| 515 |
|
| 516 |
/// For Momentum constraint
|
| 517 |
// TwoBodyDecaySpy treeTBDMom(filename,"TwoBodyDecayMomConstraint/TwoBodyDecaySpyMom");
|
| 518 |
// std::cout << "N_entries for treeTBDMom :" << treeTBDMom.GetEntriesFast() << std::endl;
|
| 519 |
|
| 520 |
TreeProducer treeMom(filename,"treeproducer3/ttree_MomTBD");
|
| 521 |
std::cout << "N_entries for treeMom :" << treeMom.GetEntriesFast() << std::endl;
|
| 522 |
|
| 523 |
|
| 524 |
|
| 525 |
/// For Vertex constraint
|
| 526 |
|
| 527 |
// TwoBodyDecaySpy treeTBDVtx(filename,"TwoBodyDecayVertexConstraint/TwoBodyDecaySpyVtx");
|
| 528 |
// std::cout << "N_entries for treeTBDVtx :" << treeTBDVtx.GetEntriesFast() << std::endl;
|
| 529 |
|
| 530 |
TreeProducer treeVtx(filename,"treeproducer5/ttree_VtxTBD");
|
| 531 |
std::cout << "N_entries for treeVtx :" << treeVtx.GetEntriesFast() << std::endl;
|
| 532 |
|
| 533 |
/// For Full constraint
|
| 534 |
TreeProducer treeFull(filename,"treeproducer2/ttree_FullTBD");
|
| 535 |
std::cout << "N_entries for treeFull :" << treeFull.GetEntriesFast() << std::endl;
|
| 536 |
|
| 537 |
/// For KinFitter (mass)constraint
|
| 538 |
|
| 539 |
// TwoBodyDecaySpy treeKinFitSpy(filename,"doConstraintKinFit/KinFitterSpy");
|
| 540 |
// std::cout << "N_entries for treeKinFitSpy :" << treeTBDVtx.GetEntriesFast() << std::endl;
|
| 541 |
|
| 542 |
TreeProducer treeKinFit(filename,"treeproducer4/ttree_KinFit");
|
| 543 |
std::cout << "N_entries for treeKinFit :" << treeKinFit.GetEntriesFast() << std::endl;
|
| 544 |
|
| 545 |
|
| 546 |
std::cout << "Initializing histogram" << std::endl;
|
| 547 |
HistoCollection histos;
|
| 548 |
histos.Initialize();
|
| 549 |
|
| 550 |
std::cout << "Filling"<<std::endl;
|
| 551 |
|
| 552 |
LoopRef(&treeRef,histos); // Fill histos related to measurements at First Fit and at Generator Level
|
| 553 |
LoopVtx(&treeVtx,histos);
|
| 554 |
LoopMom(&treeMom,histos);
|
| 555 |
LoopFull(&treeFull,histos);
|
| 556 |
LoopKinFit(&treeKinFit,histos);
|
| 557 |
|
| 558 |
//LoopTBD(&treeTBD,histos);
|
| 559 |
//LoopCommon(&tree,&treeTBD,histos);
|
| 560 |
std::cout << std::endl;
|
| 561 |
std::cout << "Drawing plots ... "<< std::endl;
|
| 562 |
histos.DrawInit();
|
| 563 |
histos.Finalize();
|
| 564 |
|
| 565 |
std::cout << "End"<<std::endl;
|
| 566 |
}
|
| 567 |
|
