#include <StdHepInfoAna.h>
Inheritance diagram for StdHepInfoAna:
Public Member Functions | |
| StdHepInfoAna (StdHepInfo &sv) | |
| virtual | ~StdHepInfoAna () |
| void | Analyze (int evtn, RecRecordImp< RecCandHeader > *srobj) |
| void | Analyze (int evtn, NtpStRecord *srobj) |
Private Attributes | |
| StdHepInfo & | fStdHepInfo |
|
|
Definition at line 17 of file StdHepInfoAna.cxx. 00017 : 00018 fStdHepInfo(shi) 00019 {}
|
|
|
Definition at line 21 of file StdHepInfoAna.cxx. 00022 {}
|
|
||||||||||||
|
Definition at line 38 of file StdHepInfoAna.cxx. References abs(), StdHepInfo::baryonxf, calbaryonpt(), checkParent(), countpi0(), StdHepInfo::egeant, StdHepInfo::ekaon, StdHepInfo::elep, StdHepInfo::emcount, StdHepInfo::emenergy, StdHepInfo::emfrac, StdHepInfo::eneut, StdHepInfo::epi0, StdHepInfo::epim, StdHepInfo::epip, StdHepInfo::eprot, StdHepInfo::etot, fStdHepInfo, ANtpRecoNtpManipulator::GetMCManipulator(), ANtpMCManipulator::GetNtpMCTruth(), ANtpMCManipulator::GetNtpTHEvent(), HughCode(), NtpMCStdHep::IdHEP, NtpMCStdHep::index, NtpMCTruth::index, NtpMCTruth::iresonance, isLeptonic(), NtpMCStdHep::IstHEP, StdHepInfo::lep2leptype, StdHepInfo::lepegeant, StdHepInfo::lepekaon, StdHepInfo::lepelep, StdHepInfo::lepeneut, StdHepInfo::lepepi0, StdHepInfo::lepepim, StdHepInfo::lepepip, StdHepInfo::lepeprot, StdHepInfo::lepetot, StdHepInfo::lepnfs, StdHepInfo::lepnfs_he, StdHepInfo::lepngeant, StdHepInfo::lepngeant_he, StdHepInfo::lepnkaon, StdHepInfo::lepnkaon_he, StdHepInfo::lepnlep, StdHepInfo::lepnlep_he, StdHepInfo::lepnneut, StdHepInfo::lepnneut_he, StdHepInfo::lepnpi0, StdHepInfo::lepnpi0_he, StdHepInfo::lepnpim, StdHepInfo::lepnpim_he, StdHepInfo::lepnpip, StdHepInfo::lepnpip_he, StdHepInfo::lepnprot, StdHepInfo::lepnprot_he, StdHepInfo::lepptgeant, StdHepInfo::lepptkaon, StdHepInfo::lepptlep, StdHepInfo::lepptneut, StdHepInfo::lepptpi0, StdHepInfo::lepptpim, StdHepInfo::lepptpip, StdHepInfo::lepptprot, StdHepInfo::leppttot, StdHepInfo::leppzgeant, StdHepInfo::leppzkaon, StdHepInfo::leppzlep, StdHepInfo::leppzneut, StdHepInfo::leppzpi0, StdHepInfo::leppzpim, StdHepInfo::leppzpip, StdHepInfo::leppzprot, StdHepInfo::leppztot, StdHepInfo::lepwfs, StdHepInfo::lepwgeant, StdHepInfo::lepwkaon, StdHepInfo::lepwlep, StdHepInfo::lepwneut, StdHepInfo::lepwpi0, StdHepInfo::lepwpim, StdHepInfo::lepwpip, StdHepInfo::lepwprot, NtpMCStdHep::mass, NtpMCStdHep::mc, NtpTHEvent::neumc, StdHepInfo::nfs, StdHepInfo::nfs_he, StdHepInfo::ngeant, StdHepInfo::ngeant_he, StdHepInfo::nkaon, StdHepInfo::nkaon_he, StdHepInfo::nlep, StdHepInfo::nlep_he, StdHepInfo::nmult, StdHepInfo::nneut, StdHepInfo::nneut_he, StdHepInfo::npi0, StdHepInfo::npi0_he, StdHepInfo::npim, StdHepInfo::npim_he, StdHepInfo::npip, StdHepInfo::npip_he, StdHepInfo::nprot, StdHepInfo::nprot_he, NtpMCStdHep::p4, NtpMCStdHep::parent, StdHepInfo::ptgeant, StdHepInfo::ptkaon, StdHepInfo::ptlep, StdHepInfo::ptneut, StdHepInfo::ptpi0, StdHepInfo::ptpim, StdHepInfo::ptpip, StdHepInfo::ptprot, StdHepInfo::pttot, StdHepInfo::pzgeant, StdHepInfo::pzkaon, StdHepInfo::pzlep, StdHepInfo::pzneut, StdHepInfo::pzpi0, StdHepInfo::pzpim, StdHepInfo::pzpip, StdHepInfo::pzprot, StdHepInfo::pztot, StdHepInfo::Reset(), NtpMCTruth::stdhep, NtpStRecord::stdhep, NtpMCTruth::w2, StdHepInfo::wfs, StdHepInfo::wgeant, StdHepInfo::wkaon, StdHepInfo::wlep, StdHepInfo::wneut, StdHepInfo::wpi0, StdHepInfo::wpim, StdHepInfo::wpip, StdHepInfo::wprot, and StdHepInfo::Zero(). 00039 {
00040 fStdHepInfo.Reset();
00041 if(srobj==0){
00042 return;
00043 }
00044 NtpMCTruth *mctruth = 0;
00045 NtpMCStdHep *ntpStdHep = 0;
00046 NtpTHEvent *thevent = 0;
00047 ANtpRecoNtpManipulator ntpManipulator(srobj);
00048 thevent = ntpManipulator.GetMCManipulator()->GetNtpTHEvent(evtn);
00049 if(thevent){
00050 mctruth = ntpManipulator.GetMCManipulator()->GetNtpMCTruth(thevent->neumc);
00051 if(mctruth){
00052 fStdHepInfo.Zero();
00053 TClonesArray& stdhepArray = *(srobj->stdhep);
00054 Int_t nStdHep = stdhepArray.GetEntries();
00055 HughCode(mctruth, stdhepArray, fStdHepInfo);
00056 Double_t maxmom = 0;
00057 Double_t maxke = 0;
00058 Double_t lepmaxmom = 0;
00059 Double_t lepmaxke = 0;
00060 Double_t initialNuEnergy=0;
00061 //count pi0's
00062 countpi0(mctruth->index,stdhepArray,fStdHepInfo);
00063 //calculate baryon pt
00064 calbaryonpt(mctruth->index,stdhepArray,fStdHepInfo);
00065 if (mctruth->w2>0) fStdHepInfo.baryonxf *= 2./sqrt(mctruth->w2);
00066
00067 for(int i=mctruth->stdhep[0];i<=mctruth->stdhep[1]&&i<nStdHep;i++){
00068 ntpStdHep = dynamic_cast<NtpMCStdHep *>(stdhepArray[i]);
00069
00070 //get initial neutrino energy:
00071 if(ntpStdHep->mc==mctruth->index &&ntpStdHep->IstHEP==0&& (abs(ntpStdHep->IdHEP)==12 || abs(ntpStdHep->IdHEP)==14 || abs(ntpStdHep->IdHEP)==16)){
00072 initialNuEnergy=TMath::Abs(ntpStdHep->p4[3]);
00073
00074 }
00075
00076
00077 //check that stdhep entry corresponds to this mc event
00078 //also don't include final state neutrinos in counting
00079
00080 if(ntpStdHep->mc==mctruth->index && abs(ntpStdHep->IdHEP)!=12 &&
00081 abs(ntpStdHep->IdHEP)!=14 && abs(ntpStdHep->IdHEP)!=16){
00082 //look for neugen intermediate particle to get multiplicity
00083 if(ntpStdHep->IstHEP==3 &&
00084 (mctruth->iresonance!=1002 ||
00085 TMath::Abs(ntpStdHep->IdHEP)!=15)){
00086 Int_t hfs = ntpStdHep->IdHEP%1000;
00087 fStdHepInfo.nmult = Int_t((hfs-128)/4)-20;
00088 }
00089
00090
00091 // look for pi0, e, and photons to make the emcount, emfrac, emenergy variables
00092 //check final decays, but dont double count energies (dont count a particle if its parents energy is also to be counted
00093 //i am using a recursive method of checking - i do not know how far the mc decay tree extends in IstHEP==205
00094 //it goes at least two decays deep in IstHEP==205, possible more
00095 else if(ntpStdHep->IstHEP==205){
00096
00097 if(abs(ntpStdHep->IdHEP)==11||abs(ntpStdHep->IdHEP)==22||abs(ntpStdHep->IdHEP)==111){
00098 Double_t energy = TMath::Abs(ntpStdHep->p4[3]);
00099 if(!checkParent(ntpStdHep->parent[0],stdhepArray)){
00100 if(ntpStdHep->parent[0]!=ntpStdHep->parent[1]){
00101 if(!checkParent(ntpStdHep->parent[1],stdhepArray)){
00102 fStdHepInfo.emenergy+=energy;
00103 fStdHepInfo.emcount++;
00104 }
00105 }else{
00106 fStdHepInfo.emenergy+=energy;
00107 fStdHepInfo.emcount++;
00108 }
00109 }
00110 }
00111 }
00112
00113
00114
00115
00116
00117 else if((ntpStdHep->IstHEP==1 ) && //get final state particles:
00118 ntpStdHep->IdHEP<1000000000){ //excluding nuclei
00119 Double_t energy = TMath::Abs(ntpStdHep->p4[3]);
00120 Double_t pz = TMath::Abs(ntpStdHep->p4[2]);
00121 Double_t pt = ( ntpStdHep->p4[0]*ntpStdHep->p4[0] +
00122 ntpStdHep->p4[1]*ntpStdHep->p4[1] );
00123 Double_t mom = TMath::Sqrt(pz*pz + pt);
00124 Double_t ke = energy - ntpStdHep->mass;
00125
00126 if(mom>maxmom) maxmom = mom;
00127 if(ke>maxke) maxke = ke;
00128 if(ke>0.5) fStdHepInfo.nfs_he += 1;
00129
00130 fStdHepInfo.nfs += 1;
00131 fStdHepInfo.etot += energy;
00132 fStdHepInfo.pztot += pz;
00133 fStdHepInfo.pttot += pt;
00134 fStdHepInfo.wfs += ke;
00135
00136
00137 // by steve cavanaugh
00138
00139 //electrons, protons, and pi0 contribute to em frac
00140 //decay particle of these which also fit 11,22,111 do not contribute as per code dealing with IstHEP===205 above
00141 if(abs(ntpStdHep->IdHEP)==11||abs(ntpStdHep->IdHEP)==22||abs(ntpStdHep->IdHEP)==111){
00142 fStdHepInfo.emenergy+=energy;
00143 fStdHepInfo.emcount++;
00144 }
00145
00146
00147
00148
00149 //see if the current particle is a decay product of the original neutrino
00150 if (isLeptonic(ntpStdHep->index, stdhepArray)){
00151
00152 fStdHepInfo.lepnfs += 1;
00153 fStdHepInfo.lepetot += energy;
00154 fStdHepInfo.leppztot += pz;
00155 fStdHepInfo.leppttot += pt;
00156 fStdHepInfo.lepwfs += ke;
00157
00158 if(mom>lepmaxmom) lepmaxmom = mom;
00159 if(ke>lepmaxke) lepmaxke = ke;
00160 if(ke>0.5) fStdHepInfo.lepnfs_he += 1;
00161
00162
00163
00164 //if not a tau, check if electron or muon, if so set flags and energies
00165 if(abs(ntpStdHep->IdHEP)==11||abs(ntpStdHep->IdHEP)==13||abs(ntpStdHep->IdHEP)==15){ // should i be checking tau also? probably only appears with IstHEP=2
00166 fStdHepInfo.lep2leptype =ntpStdHep->IdHEP;
00167 }
00168
00169
00170
00171
00172 //otherwise do standard particle contribution ...
00173
00174 if(abs(ntpStdHep->IdHEP)==11 || abs(ntpStdHep->IdHEP)==13 ||
00175 abs(ntpStdHep->IdHEP)==15) { //ive probably already accounted for tau energyies from the daughter particles.... so i probably do not want them here!!
00176 fStdHepInfo.lepnlep += 1;
00177 if(ke>0.5) fStdHepInfo.lepnlep_he += 1;
00178 fStdHepInfo.lepelep += energy;
00179 fStdHepInfo.leppzlep += pz;
00180 fStdHepInfo.lepptlep += pt;
00181 fStdHepInfo.lepwlep += ke;
00182 }
00183 else if(ntpStdHep->IdHEP==28) {
00184 fStdHepInfo.lepngeant += 1;
00185 if(ke>0.5) fStdHepInfo.lepngeant_he += 1;
00186 fStdHepInfo.lepegeant += energy;
00187 fStdHepInfo.leppzgeant += pz;
00188 fStdHepInfo.lepptgeant += pt;
00189 fStdHepInfo.lepwgeant += ke;
00190 }
00191 else if(ntpStdHep->IdHEP==2212) {
00192 fStdHepInfo.lepnprot += 1;
00193 if(ke>0.5) fStdHepInfo.lepnprot_he += 1;
00194 fStdHepInfo.lepeprot += energy;
00195 fStdHepInfo.leppzprot += pz;
00196 fStdHepInfo.lepptprot += pt;
00197 fStdHepInfo.lepwprot += ke;
00198 }
00199 else if(ntpStdHep->IdHEP==2112) {
00200 fStdHepInfo.lepnneut += 1;
00201 if(ke>0.5) fStdHepInfo.lepnneut_he += 1;
00202 fStdHepInfo.lepeneut += energy;
00203 fStdHepInfo.leppzneut += pz;
00204 fStdHepInfo.lepptneut += pt;
00205 fStdHepInfo.lepwneut += ke;
00206 }
00207 else if(ntpStdHep->IdHEP==211) {
00208 fStdHepInfo.lepnpip += 1;
00209 if(ke>0.5) fStdHepInfo.lepnpip_he += 1;
00210 fStdHepInfo.lepepip += energy;
00211 fStdHepInfo.leppzpip += pz;
00212 fStdHepInfo.lepptpip += pt;
00213 fStdHepInfo.lepwpip += ke;
00214 }
00215 else if(ntpStdHep->IdHEP==-211) {
00216 fStdHepInfo.lepnpim += 1;
00217 if(ke>0.5) fStdHepInfo.lepnpim_he += 1;
00218 fStdHepInfo.lepepim += energy;
00219 fStdHepInfo.leppzpim += pz;
00220 fStdHepInfo.lepptpim += pt;
00221 fStdHepInfo.lepwpim += ke;
00222 }
00223 else if(ntpStdHep->IdHEP==111) {
00224 fStdHepInfo.lepnpi0 += 1;
00225 if(ke>0.5) fStdHepInfo.lepnpi0_he += 1;
00226 fStdHepInfo.lepepi0 += energy;
00227 fStdHepInfo.leppzpi0 += pz;
00228 fStdHepInfo.lepptpi0 += pt;
00229 fStdHepInfo.lepwpi0 += ke;
00230 }
00231 else if(TMath::Abs(ntpStdHep->IdHEP)==321 || //K+/-
00232 ntpStdHep->IdHEP==130 || //K0_L
00233 ntpStdHep->IdHEP==310 || //K0_S
00234 ntpStdHep->IdHEP==311 ) { //K0
00235 fStdHepInfo.lepnkaon += 1;
00236 if(ke>0.5) fStdHepInfo.lepnkaon_he += 1;
00237 fStdHepInfo.lepekaon += energy;
00238 fStdHepInfo.leppzkaon += pz;
00239 fStdHepInfo.lepptkaon += pt;
00240 fStdHepInfo.lepwkaon += ke;
00241 }
00242
00243
00244 }
00245
00246
00247
00248
00249 //end edit by steve cavanaugh
00250
00251
00252 if(abs(ntpStdHep->IdHEP)==11 || abs(ntpStdHep->IdHEP)==13 ||
00253 abs(ntpStdHep->IdHEP)==15) {
00254 fStdHepInfo.nlep += 1;
00255 if(ke>0.5) fStdHepInfo.nlep_he += 1;
00256 fStdHepInfo.elep += energy;
00257 fStdHepInfo.pzlep += pz;
00258 fStdHepInfo.ptlep += pt;
00259 fStdHepInfo.wlep += ke;
00260 }
00261 else if(ntpStdHep->IdHEP==28) {
00262 fStdHepInfo.ngeant += 1;
00263 if(ke>0.5) fStdHepInfo.ngeant_he += 1;
00264 fStdHepInfo.egeant += energy;
00265 fStdHepInfo.pzgeant += pz;
00266 fStdHepInfo.ptgeant += pt;
00267 fStdHepInfo.wgeant += ke;
00268 }
00269 else if(ntpStdHep->IdHEP==2212) {
00270 fStdHepInfo.nprot += 1;
00271 if(ke>0.5) fStdHepInfo.nprot_he += 1;
00272 fStdHepInfo.eprot += energy;
00273 fStdHepInfo.pzprot += pz;
00274 fStdHepInfo.ptprot += pt;
00275 fStdHepInfo.wprot += ke;
00276 }
00277 else if(ntpStdHep->IdHEP==2112) {
00278 fStdHepInfo.nneut += 1;
00279 if(ke>0.5) fStdHepInfo.nneut_he += 1;
00280 fStdHepInfo.eneut += energy;
00281 fStdHepInfo.pzneut += pz;
00282 fStdHepInfo.ptneut += pt;
00283 fStdHepInfo.wneut += ke;
00284 }
00285 else if(ntpStdHep->IdHEP==211) {
00286 fStdHepInfo.npip += 1;
00287 if(ke>0.5) fStdHepInfo.npip_he += 1;
00288 fStdHepInfo.epip += energy;
00289 fStdHepInfo.pzpip += pz;
00290 fStdHepInfo.ptpip += pt;
00291 fStdHepInfo.wpip += ke;
00292 }
00293 else if(ntpStdHep->IdHEP==-211) {
00294 fStdHepInfo.npim += 1;
00295 if(ke>0.5) fStdHepInfo.npim_he += 1;
00296 fStdHepInfo.epim += energy;
00297 fStdHepInfo.pzpim += pz;
00298 fStdHepInfo.ptpim += pt;
00299 fStdHepInfo.wpim += ke;
00300 }
00301 else if(ntpStdHep->IdHEP==111) {
00302 fStdHepInfo.npi0 += 1;
00303 if(ke>0.5) fStdHepInfo.npi0_he += 1;
00304 fStdHepInfo.epi0 += energy;
00305 fStdHepInfo.pzpi0 += pz;
00306 fStdHepInfo.ptpi0 += pt;
00307 fStdHepInfo.wpi0 += ke;
00308 }
00309 else if(TMath::Abs(ntpStdHep->IdHEP)==321 || //K+/-
00310 ntpStdHep->IdHEP==130 || //K0_L
00311 ntpStdHep->IdHEP==310 || //K0_S
00312 ntpStdHep->IdHEP==311 ) { //K0
00313 fStdHepInfo.nkaon += 1;
00314 if(ke>0.5) fStdHepInfo.nkaon_he += 1;
00315 fStdHepInfo.ekaon += energy;
00316 fStdHepInfo.pzkaon += pz;
00317 fStdHepInfo.ptkaon += pt;
00318 fStdHepInfo.wkaon += ke;
00319 }
00320 }
00321 }
00322 }
00323 if(maxke>0){
00324 fStdHepInfo.wfs /= maxke;
00325 fStdHepInfo.wgeant /= maxke;
00326 fStdHepInfo.wlep /= maxke;
00327 fStdHepInfo.wprot /= maxke;
00328 fStdHepInfo.wneut /= maxke;
00329 fStdHepInfo.wpip /= maxke;
00330 fStdHepInfo.wpim /= maxke;
00331 fStdHepInfo.wpi0 /= maxke;
00332 fStdHepInfo.wkaon /= maxke;
00333 }
00334
00335 if(lepmaxke>0){
00336 fStdHepInfo.lepwfs /= lepmaxke;
00337 fStdHepInfo.lepwgeant /= lepmaxke;
00338 fStdHepInfo.lepwlep /= lepmaxke;
00339 fStdHepInfo.lepwprot /= lepmaxke;
00340 fStdHepInfo.lepwneut /= lepmaxke;
00341 fStdHepInfo.lepwpip /= lepmaxke;
00342 fStdHepInfo.lepwpim /= lepmaxke;
00343 fStdHepInfo.lepwpi0 /= lepmaxke;
00344 fStdHepInfo.lepwkaon /= lepmaxke;
00345 }
00346
00347
00348 if(initialNuEnergy>0){
00349 fStdHepInfo.emfrac=fStdHepInfo.emenergy /initialNuEnergy;
00350
00351 }else{
00352 fStdHepInfo.emfrac= ANtpDefVal::kDouble;
00353 }
00354
00355
00356
00357 if(fStdHepInfo.pttot>0)
00358 fStdHepInfo.pttot = TMath::Sqrt(fStdHepInfo.pttot);
00359 if(fStdHepInfo.ptgeant>0)
00360 fStdHepInfo.ptgeant = TMath::Sqrt(fStdHepInfo.ptgeant);
00361 if(fStdHepInfo.ptlep>0)
00362 fStdHepInfo.ptlep = TMath::Sqrt(fStdHepInfo.ptlep);
00363 if(fStdHepInfo.ptprot>0)
00364 fStdHepInfo.ptprot = TMath::Sqrt(fStdHepInfo.ptprot);
00365 if(fStdHepInfo.ptneut>0)
00366 fStdHepInfo.ptneut = TMath::Sqrt(fStdHepInfo.ptneut);
00367 if(fStdHepInfo.ptpip>0)
00368 fStdHepInfo.ptpip = TMath::Sqrt(fStdHepInfo.ptpip);
00369 if(fStdHepInfo.ptpim>0)
00370 fStdHepInfo.ptpim = TMath::Sqrt(fStdHepInfo.ptpim);
00371 if(fStdHepInfo.ptpi0>0)
00372 fStdHepInfo.ptpi0 = TMath::Sqrt(fStdHepInfo.ptpi0);
00373 if(fStdHepInfo.ptkaon>0)
00374 fStdHepInfo.ptkaon = TMath::Sqrt(fStdHepInfo.ptkaon);
00375
00376 if(fStdHepInfo.leppttot>0)
00377 fStdHepInfo.leppttot = TMath::Sqrt(fStdHepInfo.leppttot);
00378 if(fStdHepInfo.lepptgeant>0)
00379 fStdHepInfo.lepptgeant = TMath::Sqrt(fStdHepInfo.lepptgeant);
00380 if(fStdHepInfo.lepptlep>0)
00381 fStdHepInfo.lepptlep = TMath::Sqrt(fStdHepInfo.lepptlep);
00382 if(fStdHepInfo.lepptprot>0)
00383 fStdHepInfo.lepptprot = TMath::Sqrt(fStdHepInfo.lepptprot);
00384 if(fStdHepInfo.lepptneut>0)
00385 fStdHepInfo.lepptneut = TMath::Sqrt(fStdHepInfo.lepptneut);
00386 if(fStdHepInfo.lepptpip>0)
00387 fStdHepInfo.lepptpip = TMath::Sqrt(fStdHepInfo.lepptpip);
00388 if(fStdHepInfo.lepptpim>0)
00389 fStdHepInfo.lepptpim = TMath::Sqrt(fStdHepInfo.lepptpim);
00390 if(fStdHepInfo.lepptpi0>0)
00391 fStdHepInfo.lepptpi0 = TMath::Sqrt(fStdHepInfo.lepptpi0);
00392 if(fStdHepInfo.lepptkaon>0)
00393 fStdHepInfo.lepptkaon = TMath::Sqrt(fStdHepInfo.lepptkaon);
00394 }
00395 }
00396 }
|
|
||||||||||||
|
Implements NueAnaBase. Definition at line 29 of file StdHepInfoAna.cxx. Referenced by NueDisplayModule::Analyze(), and NueRecordAna::FillTrue(). 00030 {
00031 if(srobj==0){
00032 return;
00033 }
00034 NtpStRecord *st = dynamic_cast<NtpStRecord *>(srobj);
00035 Analyze(event,st);
00036 }
|
|
|
Definition at line 19 of file StdHepInfoAna.h. Referenced by Analyze(). |
1.3.9.1