loon: StdHepInfoAna.cxx Source File

00001 #include "CandNtupleSR/NtpSREvent.h"
00002 #include "StandardNtuple/NtpStRecord.h"
00003 #include "MCNtuple/NtpMCTruth.h"
00004 #include "MCNtuple/NtpMCStdHep.h"
00005 #include "TruthHelperNtuple/NtpTHEvent.h"
00006 #include "MessageService/MsgService.h"
00007 #include "StdHepInfoAna.h"
00008 #include "AnalysisNtuples/Module/ANtpRecoNtpManipulator.h"
00009 #include "AnalysisNtuples/ANtpDefaultValue.h"
00010 
00011 CVSID("$Id: StdHepInfoAna.cxx,v 1.16 2009/02/12 04:58:55 tjyang Exp $");
00012 
00013 void CalcDaughters(int start, int end, double &E, double &EmE, int &npi0, int &nch, int &ntot, TClonesArray& stdhepArray);
00014 void HughCode(NtpMCTruth* mctruth, TClonesArray& stdhepArray, StdHepInfo &shi);
00015 
00016 
00017 StdHepInfoAna::StdHepInfoAna(StdHepInfo &shi):
00018   fStdHepInfo(shi)
00019 {}
00020 
00021 StdHepInfoAna::~StdHepInfoAna()
00022 {}
00023 
00024 bool isLeptonic(int parent, TClonesArray& stdhepArray);
00025 bool checkParent(int parent, TClonesArray& stdhepArray);
00026 void countpi0(int index, TClonesArray& stdhepArray, StdHepInfo &fStdHepInfo);
00027 void calbaryonpt(int index, TClonesArray& stdhepArray, StdHepInfo &fStdHepInfo);
00028 
00029 void StdHepInfoAna::Analyze(int event, RecRecordImp<RecCandHeader> *srobj)
00030 {
00031   if(srobj==0){
00032     return;
00033   }
00034   NtpStRecord *st = dynamic_cast<NtpStRecord *>(srobj);
00035   Analyze(event,st);
00036 }
00037 
00038 void StdHepInfoAna::Analyze(int evtn, NtpStRecord *srobj)
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 }
00397 
00398 
00399 
00400 
00401 //this class will determine if a given particle has a parent which is from the original neutrino, pass in parent of current particle
00402 bool isLeptonic(int parent, TClonesArray& stdhepArray)  
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 }
00421 
00422 
00423 //this class will determine if a given particle has a parent with final state 1 or 205 of type electron, pi0, or photon
00424 //returns false if no parent is wanted type
00425 bool checkParent(int parent, TClonesArray& stdhepArray)  
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 }
00445 
00446 
00447 
00448 void countpi0(int index, TClonesArray& stdhepArray, StdHepInfo &fStdHepInfo)
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 }
00546 
00547 void calbaryonpt(int index, TClonesArray& stdhepArray, StdHepInfo &fStdHepInfo)
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 }
00632 
00633 void HughCode(NtpMCTruth* mctruth, TClonesArray& stdhepArray, StdHepInfo& shi)
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 }
00672                                                                                                              
00673 void CalcDaughters(int start, int end, double &E, double &EmE, int &npi0, int &nch, int &ntot, TClonesArray&
00674 stdhepArray)
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 }
00699                                                                                                              
00700