#include "CandNtupleSR/NtpSREvent.h"#include "StandardNtuple/NtpStRecord.h"#include "MCNtuple/NtpMCTruth.h"#include "MCNtuple/NtpMCStdHep.h"#include "TruthHelperNtuple/NtpTHEvent.h"#include "MessageService/MsgService.h"#include "StdHepInfoAna.h"#include "AnalysisNtuples/Module/ANtpRecoNtpManipulator.h"#include "AnalysisNtuples/ANtpDefaultValue.h"Go to the source code of this file.
Functions | |
| CVSID ("$Id: StdHepInfoAna.cxx,v 1.16 2009/02/12 04:58:55 tjyang Exp $") | |
| void | CalcDaughters (int start, int end, double &E, double &EmE, int &npi0, int &nch, int &ntot, TClonesArray &stdhepArray) |
| void | HughCode (NtpMCTruth *mctruth, TClonesArray &stdhepArray, StdHepInfo &shi) |
| bool | isLeptonic (int parent, TClonesArray &stdhepArray) |
| bool | checkParent (int parent, TClonesArray &stdhepArray) |
| void | countpi0 (int index, TClonesArray &stdhepArray, StdHepInfo &fStdHepInfo) |
| void | calbaryonpt (int index, TClonesArray &stdhepArray, StdHepInfo &fStdHepInfo) |
|
||||||||||||||||
|
Definition at line 547 of file StdHepInfoAna.cxx. References StdHepInfo::baryonpt, StdHepInfo::baryonxf, NtpMCStdHep::IdHEP, NtpMCStdHep::IstHEP, NtpMCStdHep::mass, NtpMCStdHep::mc, StdHepInfo::nhad, StdHepInfo::npi0_neugen, NtpMCStdHep::p4, NtpMCStdHep::parent, and StdHepInfo::totpt. Referenced by StdHepInfoAna::Analyze(). 00548 {
00549
00550 vector<int> pid;
00551 vector<double> px;
00552 vector<double> py;
00553 vector<double> pz;
00554 vector<double> eng;
00555 vector<double> mass;
00556
00557 double hspx = 0.;
00558 double hspy = 0.;
00559 double hspz = 0.;
00560 double hse = 0.;
00561 double beta = 0.;
00562 double gamma = 0.;
00563
00564 //cout<<"counting pi0"<<endl;
00565 Int_t nStdHep = stdhepArray.GetEntries();
00566 NtpMCStdHep *ntpStdHep = 0;
00567 for(int i=0; i<nStdHep; i++){
00568 ntpStdHep = dynamic_cast<NtpMCStdHep *>(stdhepArray[i]);
00569 //check that stdhep entry corresponds to this mc event
00570 if(ntpStdHep->mc!=index) continue;
00571
00572 if(ntpStdHep->IstHEP==3){//resonance decay
00573 if (TMath::Abs(ntpStdHep->p4[3])>hse){
00574 hspx = ntpStdHep->p4[0];
00575 hspy = ntpStdHep->p4[1];
00576 hspz = ntpStdHep->p4[2];
00577 hse = TMath::Abs(ntpStdHep->p4[3]);
00578 }
00579 }
00580 else {
00581 //check if this particle comes from a resonance decay or DIS
00582 bool resdecay = true;
00583 for (int ip = ntpStdHep->parent[0]; ip<=ntpStdHep->parent[1]; ip++){
00584 if (ip==-1) {
00585 resdecay = false;
00586 break;
00587 }
00588 NtpMCStdHep *parent = dynamic_cast<NtpMCStdHep *>(stdhepArray[ip]);
00589 resdecay = resdecay && (parent->IstHEP==3);
00590 if (!resdecay) break;
00591 }
00592 if (resdecay&&ntpStdHep->IdHEP) {
00593 pid.push_back(ntpStdHep->IdHEP);
00594 px.push_back(ntpStdHep->p4[0]);
00595 py.push_back(ntpStdHep->p4[1]);
00596 pz.push_back(ntpStdHep->p4[2]);
00597 eng.push_back(ntpStdHep->p4[3]);
00598 mass.push_back(ntpStdHep->mass);
00599 }
00600 }
00601 }
00602
00603 double maxbaryone = 0;
00604 double Phs = sqrt(pow(hspx,2)+pow(hspy,2)+pow(hspz,2));
00605
00606 if (Phs>=0.&&hse>Phs){
00607 beta = Phs/hse;
00608 gamma = hse/sqrt(pow(hse,2)-pow(Phs,2));
00609 }
00610
00611
00612 for (unsigned ipar = 0; ipar<pid.size(); ipar++){//loop through hadrons
00613 double ptot = sqrt(pow(px[ipar],2)+pow(py[ipar],2)+pow(pz[ipar],2));
00614 double Pz = 0;
00615 if (Phs>0) Pz = (px[ipar]*hspx+py[ipar]*hspy+pz[ipar]*hspz)/Phs;
00616 double pt = pow(ptot,2)-pow(Pz,2);
00617 Pz = gamma*(Pz-beta*eng[ipar]);
00618 if (pt>0.00000001) pt = sqrt(pt);
00619 else pt = 0;
00620 if (TMath::Abs(pid[ipar])==2212||TMath::Abs(pid[ipar])==2112){
00621 if (eng[ipar]>maxbaryone){
00622 maxbaryone = eng[ipar];
00623 fStdHepInfo.baryonpt = pt;
00624 fStdHepInfo.baryonxf = Pz;
00625 }
00626 }
00627 if (pid[ipar]==111) fStdHepInfo.npi0_neugen++;
00628 fStdHepInfo.totpt += pt;
00629 }
00630 fStdHepInfo.nhad = int(pid.size());
00631 }
|
|
||||||||||||||||||||||||||||||||||||
|
Definition at line 673 of file StdHepInfoAna.cxx. References NtpMCStdHep::child, NtpMCStdHep::IdHEP, NtpMCStdHep::IstHEP, and NtpMCStdHep::p4. Referenced by HughCode(). 00675 {
00676 NtpMCStdHep* ntpStdHep = 0;
00677
00678 for(int j = start; j <= end; j++){
00679 ntpStdHep = dynamic_cast<NtpMCStdHep *>(stdhepArray[j]);
00680 int id = ntpStdHep->IdHEP;
00681 int ist = ntpStdHep->IstHEP;
00682 if(ist == 3){
00683 int st = ntpStdHep->child[0];
00684 int en = ntpStdHep->child[1];
00685 CalcDaughters(st, en, E, EmE, npi0, nch, ntot, stdhepArray);
00686 }
00687
00688 float Etemp = ntpStdHep->p4[3];
00689
00690 E += Etemp;
00691 if(id == 111 || id == 22){
00692 EmE += Etemp;
00693 npi0++;
00694 }
00695 if(id == 211 || id == -211) nch++;
00696 ntot++;
00697 }
00698 }
|
|
||||||||||||
|
Definition at line 425 of file StdHepInfoAna.cxx. References abs(), NtpMCStdHep::IdHEP, NtpMCStdHep::IstHEP, and NtpMCStdHep::parent. Referenced by StdHepInfoAna::Analyze(). 00426 {
00427
00428 NtpMCStdHep *ntpStdHepTemp = 0;
00429 ntpStdHepTemp = dynamic_cast<NtpMCStdHep *>(stdhepArray[parent]);
00430 int parentID = ntpStdHepTemp->IdHEP;
00431
00432 if(!(ntpStdHepTemp->IstHEP==1||ntpStdHepTemp->IstHEP==205))return false;
00433
00434 if ((abs(parentID)==22||abs(parentID)==111||abs(parentID)==11))return true; //
00435
00436 if (checkParent(ntpStdHepTemp->parent[0],stdhepArray)) return true;
00437 if (ntpStdHepTemp->parent[0]!=ntpStdHepTemp->parent[1]) //check each parent branch if they are different
00438 {
00439 if (checkParent(ntpStdHepTemp->parent[1],stdhepArray)) return true;
00440 }
00441
00442 return false;
00443
00444 }
|
|
||||||||||||||||
|
Definition at line 448 of file StdHepInfoAna.cxx. References NtpMCStdHep::child, StdHepInfo::epi0_abs, StdHepInfo::epi0_intranuke, StdHepInfo::epi0_neugen, StdHepInfo::epi0_total, NtpMCStdHep::IdHEP, NtpMCStdHep::IstHEP, NtpMCStdHep::mc, NtpMCStdHep::p4, and NtpMCStdHep::parent. Referenced by StdHepInfoAna::Analyze(). 00449 {
00450 //cout<<"counting pi0"<<endl;
00451 Int_t nStdHep = stdhepArray.GetEntries();
00452 NtpMCStdHep *ntpStdHep = 0;
00453 for(int i=0; i<nStdHep; i++){
00454 ntpStdHep = dynamic_cast<NtpMCStdHep *>(stdhepArray[i]);
00455 //check that stdhep entry corresponds to this mc event
00456 if(ntpStdHep->mc!=index) continue;
00457
00458 if(ntpStdHep->IdHEP==111){//a pi0 is spotted
00459 //check if this pi0 comes from a resonance decay or DIS
00460 bool resdecay = true;
00461 for (int ip = ntpStdHep->parent[0]; ip<=ntpStdHep->parent[1]; ip++){
00462 if (ip==-1) {
00463 resdecay = false;
00464 break;
00465 }
00466 NtpMCStdHep *parent = dynamic_cast<NtpMCStdHep *>(stdhepArray[ip]);
00467 resdecay = resdecay && (parent->IstHEP==3);
00468 if (!resdecay) break;
00469 }
00470 if (resdecay) fStdHepInfo.epi0_neugen += TMath::Abs(ntpStdHep->p4[3]);
00471
00472 Double_t etmp = 0; // count epi0 in the child list
00473 for (int ic = ntpStdHep->child[0]; ic<=ntpStdHep->child[1]; ic++){
00474 if (ic==-1) break;
00475 NtpMCStdHep *child = dynamic_cast<NtpMCStdHep *>(stdhepArray[ic]);
00476 if (child->IdHEP==111) etmp += TMath::Abs(child->p4[3]);
00477 }
00478 if (ntpStdHep->IstHEP!=205 && //not a pi0 decay
00479 ntpStdHep->IstHEP!=1205 &&
00480 ntpStdHep->child[0]!=-1 &&
00481 ntpStdHep->child[1]!=-1){
00482 fStdHepInfo.epi0_abs += TMath::Abs(ntpStdHep->p4[3]) - etmp;
00483 }
00484
00485 bool pi0pro = true;
00486 for (int ip = ntpStdHep->parent[0]; ip<=ntpStdHep->parent[1]; ip++){
00487 if (ip==-1) {
00488 pi0pro = false;
00489 break;
00490 }
00491 NtpMCStdHep *parent = dynamic_cast<NtpMCStdHep *>(stdhepArray[ip]);
00492 pi0pro = pi0pro&&(parent->IstHEP == 14 && parent->IdHEP != 111);
00493 if (!pi0pro) break;
00494 }
00495 if (pi0pro) {
00496 fStdHepInfo.epi0_intranuke += TMath::Abs(ntpStdHep->p4[3]);
00497 fStdHepInfo.epi0_total += TMath::Abs(ntpStdHep->p4[3]);
00498 }
00499 else if (ntpStdHep->IstHEP == 1){
00500 fStdHepInfo.epi0_total += TMath::Abs(ntpStdHep->p4[3]);
00501 }
00502
00503 bool decay = true;
00504 for (int ip = ntpStdHep->parent[0]; ip<=ntpStdHep->parent[1]; ip++){
00505 if (ip==-1) {
00506 decay = false;
00507 break;
00508 }
00509 NtpMCStdHep *parent = dynamic_cast<NtpMCStdHep *>(stdhepArray[ip]);
00510 decay = decay&&(parent->IstHEP == 1 &&
00511 (ntpStdHep->IstHEP == 205||ntpStdHep->IstHEP == 1205)
00512 && parent->IdHEP != 111);
00513 if (!decay) break;
00514 }
00515 if (decay) {
00516 fStdHepInfo.epi0_intranuke += TMath::Abs(ntpStdHep->p4[3]);
00517 fStdHepInfo.epi0_total += TMath::Abs(ntpStdHep->p4[3]);
00518 }
00519 }//find a pi0
00520 else if(ntpStdHep->IdHEP==22){//a gamma is spotted
00521 bool decay = true;
00522 for (int ip = ntpStdHep->parent[0]; ip<=ntpStdHep->parent[1]; ip++){
00523 if (ip==-1) {
00524 decay = false;
00525 break;
00526 }
00527 NtpMCStdHep *parent = dynamic_cast<NtpMCStdHep *>(stdhepArray[ip]);
00528 decay = decay&&(parent->IstHEP == 1 &&
00529 (ntpStdHep->IstHEP == 205||ntpStdHep->IstHEP == 1205)
00530 &&parent->IdHEP != 111);
00531 if (!decay) break;
00532 }
00533 if (decay) {
00534 fStdHepInfo.epi0_intranuke += TMath::Abs(ntpStdHep->p4[3]);
00535 fStdHepInfo.epi0_total += TMath::Abs(ntpStdHep->p4[3]);
00536 }
00537 }
00538 }//loop through the particle list
00539 // cout<<"pi0 summary "<<endl;
00540 // cout<<"epi0_total "<<fStdHepInfo.epi0_total<<endl;
00541 // cout<<"epi0_neugen "<<fStdHepInfo.epi0_neugen<<endl;
00542 // cout<<"epi0_intranuke "<<fStdHepInfo.epi0_intranuke<<endl;
00543 // cout<<"epi0_abs "<<fStdHepInfo.epi0_abs<<endl;
00544
00545 }
|
|
||||||||||||
|
|
|
||||||||||||||||
|
Definition at line 633 of file StdHepInfoAna.cxx. References CalcDaughters(), NtpMCStdHep::child, StdHepInfo::e_total, StdHepInfo::em_total, NtpMCTruth::index, NtpMCStdHep::IstHEP, NtpMCStdHep::mc, StdHepInfo::nch_jbhg, StdHepInfo::npi0_jbhg, StdHepInfo::ntot_jbhg, and NtpMCTruth::stdhep. Referenced by StdHepInfoAna::Analyze(). 00634 {
00635 //find the mc information for this event
00636 NtpMCStdHep* ntpStdHep = 0;
00637
00638 bool foundFirst = false;
00639 int index = 0;
00640
00641 double Etot = 0.0;
00642 double EMtot = 0.0;
00643 int npi0 = 0;
00644 int nch = 0;
00645 int ntot = 0;
00646
00647 for(int i=mctruth->stdhep[0];i<=mctruth->stdhep[1]&&!foundFirst;i++){
00648 ntpStdHep = dynamic_cast<NtpMCStdHep *>(stdhepArray[i]);
00649
00650 if(ntpStdHep->mc!=mctruth->index) cout<<"Big Problem"<<endl;
00651
00652 if(!foundFirst && ntpStdHep->IstHEP == 3) foundFirst = true;
00653 else continue;
00654 index = i;
00655 }
00656
00657 if(foundFirst){
00658 int start = ntpStdHep->child[0];
00659 int end = ntpStdHep->child[1];
00660
00661 CalcDaughters(start, end, Etot, EMtot, npi0, nch, ntot, stdhepArray);
00662
00663 Etot = Etot - 0.935;
00664 }
00665
00666 shi.e_total = Etot;
00667 shi.em_total = EMtot;
00668 shi.npi0_jbhg = npi0;
00669 shi.nch_jbhg = nch;
00670 shi.ntot_jbhg = ntot;
00671 }
|
|
||||||||||||
|
Definition at line 402 of file StdHepInfoAna.cxx. References abs(), NtpMCStdHep::IdHEP, NtpMCStdHep::IstHEP, and NtpMCStdHep::parent. Referenced by StdHepInfoAna::Analyze(). 00403 {
00404
00405 NtpMCStdHep *ntpStdHepTemp = 0;
00406 ntpStdHepTemp = dynamic_cast<NtpMCStdHep *>(stdhepArray[parent]);
00407 int parentID = ntpStdHepTemp->IdHEP;
00408
00409 if ((abs(parentID)==12||abs(parentID)==14||abs(parentID)==16)&&ntpStdHepTemp->IstHEP==0)return true; //related to the original neutrino!
00410 if (ntpStdHepTemp->parent[0]==-1||ntpStdHepTemp->parent[1]==-1)return false; //at end of chain, not parent so exit
00411
00412 if (isLeptonic(ntpStdHepTemp->parent[0],stdhepArray)) return true;
00413 if (ntpStdHepTemp->parent[0]!=ntpStdHepTemp->parent[1]) //check each parent branch if they are different
00414 {
00415 if (isLeptonic(ntpStdHepTemp->parent[1],stdhepArray)) return true;
00416 }
00417
00418 return false;
00419
00420 }
|
1.3.9.1