loon: fill_stdhep.h File Reference

00051                                                                                         {
00052   
00053   // add p to stdhep, returning 
00054   // position in list 
00055   // -1 if the particle was ignored.
00056   // -2 in case of trouble
00057   
00058   if(p.IstHEP>200) return -1; // geant made these, or whacky 999 entry
00059   // see if this particle was created from one with IstHEP>200
00060   // this is why
00061   if( (p.parent[0]==p.parent[1]) && p.parent[0]>0 && p.parent[1]>0){
00062     bool dont_add=false;  
00063     for(int j=0; j<2; j++){          
00064       const int m=p.parent[j];
00065       if(m>=0 && m<v.GetEntriesFast()){
00066         const NtpMCStdHep& mommy = dynamic_cast<const NtpMCStdHep&>( *(v[m]));
00067         if(mommy.IstHEP>200) dont_add=true;
00068       }
00069     }
00070     if(dont_add) return -1;
00071   }
00072   const int i=hepevt_.nhep; // postion of new entry
00073   hepevt_.isthep[i]=p.IstHEP;
00074   hepevt_.idhep[i]=p.IdHEP;
00075   // remember, C convention for 2D array indexing is opposite of FORTRAN
00076 
00077   for(int j=0; j<2; j++){        
00078     hepevt_.jmohep[i][j]=p.parent[j];
00079     hepevt_.jdahep[i][j]=p.child[j];
00080     // make indices start at 1
00081     if(p.parent[j]!=-1) hepevt_.jmohep[i][j] -=(offset-1);
00082     if(p.child[j]!=-1) hepevt_.jdahep[i][j]-=(offset-1);
00083     // stdhep package checks for stable particles
00084     // by looking for jdahep = 0 (not -1 and not a test on isthep!!)
00085     // translate -1 to 0 here
00086     if(hepevt_.jmohep[i][j]== -1) hepevt_.jmohep[i][j]=0;
00087     if(hepevt_.jdahep[i][j]== -1) hepevt_.jdahep[i][j]=0;
00088   }
00089   for(int j=0; j<4; j++){
00090     hepevt_.phep[i][j] = p.p4[j];
00091     hepevt_.vhep[i][j] = p.vtx[j];
00092   }
00093   double psq=0.0;
00094   for(int j=0; j<3; j++) psq+=pow(p.p4[j],2);
00095   hepevt_.phep[i][4] =  pow(p.p4[3],2) - psq; 
00096   // precision errors a problem here
00097   // we care about > 1 MeV
00098   const double epsilon = pow(1e-3,2);
00099   if(fabs(hepevt_.phep[i][4])<epsilon) hepevt_.phep[i][4]=0;
00100   // geantinos not on shell
00101   if(hepevt_.idhep[i]==28) hepevt_.phep[i][4]=0;  
00102   // mass zero particles
00103   if(abs(hepevt_.idhep[i])==14 || abs(hepevt_.idhep[i])==12 
00104      || abs(hepevt_.idhep[i])==16 
00105      || abs(hepevt_.idhep[i])==22) hepevt_.phep[i][4]=0;
00106   // by hand for electrons (rounding problems)
00107   if(abs(hepevt_.idhep[i])==11)    hepevt_.phep[i][4] = 0.000502;
00108 
00109   if( hepevt_.phep[i][4] < 0.0) {
00110       std::cerr<<"phep["<<i<<"][4] = "<<hepevt_.phep[i][4]<<" < 0!!"<<std::endl;
00111       p.Print(std::cerr);
00112       inuke_reweight::print_stdhep(2);
00113       return -2;
00114       
00115   }
00116   else hepevt_.phep[i][4] = sqrt(hepevt_.phep[i][4]);
00117   hepevt_.nhep++;
00118   
00119   return i;
00120 
00121   
00122 }

00124                                                {
00125   if(!hardclear) { hepevt_.nhep=0; return; }
00126   // hardclear erases the existing block
00127   for(int i=0; i<hepevt_.nhep; i++){
00128     hepevt_.isthep[i]=0;
00129     hepevt_.idhep[i]=0;
00130     for(int j=0; j<2; j++){
00131       hepevt_.jmohep[i][j]= -1;
00132       hepevt_.jdahep[i][j]= -1;
00133     }
00134     for(int j=0; j<4; j++){
00135       hepevt_.phep[i][j] = 0;
00136       hepevt_.vhep[i][j] = 0;      
00137     }
00138       hepevt_.phep[i][4] = 0;
00139   }
00140   hepevt_.nhep=0;
00141 }

00024                                                                          {
00025   // return number of entries filled
00026   // -1 in case of error
00027   // ilow is the starting index
00028   // ihigh is the ending index (not 1 past the end)
00029   if(ilow>ihigh || ilow<0 || (ihigh-1)>v.GetEntriesFast()){
00030     std::cerr<<"Index problem!! ilow ihigh entries : "
00031              <<ilow<<" "<<ihigh<<" "<<v.GetEntriesFast()<<std::endl;
00032     return -1;
00033   }
00034   int nfill=0;
00035   for(int i=ilow; i<=ihigh; i++){  
00036     const NtpMCStdHep* pp = dynamic_cast<const NtpMCStdHep* >(v[i]);
00037     if(!pp) return -1;  
00038     const NtpMCStdHep& p = *pp;
00039     int r=add_to_stdhep(p, ilow, v);
00040     if(r==-2) return -1;
00041     
00042     else if(r>-1) nfill++;
00043   }
00044   stdhep_to_nuclear_coordinates();
00045   truncate_geant_daughters();
00046   rotate_vhep();
00047 
00048   return nfill;
00049 }

00018                                                                          {
00019   const int ilow=t.stdhep[0];
00020   const int ihigh=t.stdhep[1];
00021   return fill_stdhep(v,ilow,ihigh);
00022 }

00147                                           {
00148   heplst_(&iopt);
00149 }

00218                                 {
00219   // neugen generates neutrinos oriented along z
00220   // neugen_kine.F then rotates phep[][] along the direction
00221   // indicated by gnumi but neglects to rotate vhep[][]
00222   // this isn't a problem for isthep=1 since those are assigned
00223   // global coordinates
00224   //
00225   // so: rotate vhep to agree with neutrino direction
00226   // 
00227   // this is a copy of geant's gdrot
00228   //
00229   
00230   const double costh = hepevt_.phep[0][2]/hepevt_.phep[0][3];
00231   const double sinth = sqrt(1.0-costh*costh);
00232   const double phi=atan2(hepevt_.phep[0][1],hepevt_.phep[0][0]);
00233   const double cosph = cos(phi);
00234   const double sinph = sin(phi);
00235   
00236   for(int i=1; i<hepevt_.nhep; i++){
00237     double v[3]={hepevt_.vhep[i][0],hepevt_.vhep[i][1],hepevt_.vhep[i][2]};
00238     hepevt_.vhep[i][0]=v[0]*costh*cosph - v[1]*sinph + v[2]*sinth*cosph;
00239     hepevt_.vhep[i][1]=v[0]*costh*sinph + v[1]*cosph + v[2]*sinth*sinph;
00240     hepevt_.vhep[i][2]=-v[0]*sinth + v[2]*costh;
00241 
00242   }
00243   
00244   
00245 }

00180                                                   {
00181 
00182   for(int i=1; i<hepevt_.nhep; i++){
00183     if( hepevt_.isthep[i]!=14 
00184         && hepevt_.isthep[i]!=11 && hepevt_.isthep[i]!=3){
00185       // remove external vertex position
00186       // preserve any relative changes
00187       for(int j=0; j<4; j++){
00188         hepevt_.vhep[i][j] -= hepevt_.vhep[0][j];
00189       }
00190     }
00191     else{
00192       for(int j=0; j<3; j++){
00193         // convert to fm
00194         hepevt_.vhep[i][j] *= 1e15;
00195       }
00196       hepevt_.vhep[i][3]=0.0;
00197     }
00198   }
00199   // zero first entry
00200   for(int j=0; j<4; j++){
00201     hepevt_.vhep[0][j] = 0;
00202   }
00203 
00204 
00205 }

00152                                                                     {
00153 
00154   TFile f(filename);
00155   TTree* t = (TTree*) f.Get("NtpSt");
00156   NtpStRecord* rec=0;
00157   t->SetBranchAddress("NtpStRecord",&rec);
00158   t->GetEntry(entry);
00159 
00160   const TClonesArray& mc_array= *(rec->mc);
00161   const TClonesArray& stdheps = *(rec->stdhep);
00162 
00163   for(int i=0; i<stdheps.GetEntries(); i++){
00164     const NtpMCStdHep* s = dynamic_cast<const NtpMCStdHep*>(stdheps[i]);
00165     s->Print();
00166   }
00167   
00168   int nmc = rec->mchdr.nmc;
00169   for(int i=0; i<nmc; i++){
00170     const NtpMCTruth* mct = dynamic_cast<const NtpMCTruth*>(mc_array[i]);
00171     mct->Print();
00172     int nfill = inuke_reweight::fill_stdhep(*mct, stdheps);
00173     std::cout<<"MC Event: "<<i<<"  filled "<<nfill<<" entries"<<std::endl;
00174     inuke_reweight::print_stdhep(2);
00175     inuke_reweight::clear_stdhep();
00176   }
00177   return;
00178 }

00208                                              {
00209 
00210   for(int i=0; i<hepevt_.nhep; i++){
00211 
00212     if(hepevt_.jdahep[i][0] > hepevt_.nhep) hepevt_.jdahep[i][0]=0;
00213     if(hepevt_.jdahep[i][1] > hepevt_.nhep) hepevt_.jdahep[i][1]=0;
00214   }
00215 
00216 }


Namespaces
namespace	inuke_reweight
Functions
int	fill_stdhep (const NtpMCTruth &t, const TClonesArray &stdheps)
int	fill_stdhep (const TClonesArray &stdheps, int ilow, int ihigh)
int	add_to_stdhep (const NtpMCStdHep &p, int offset, const TClonesArray &stdheps)
void	clear_stdhep (bool hardclear=false)
void	print_stdhep (int iopt)
void	test_fill_stdhep (const char *filename, int entry=0)
void	stdhep_to_nuclear_coordinates ()
void	truncate_geant_daughters ()
void	rotate_vhep ()

fill_stdhep.h File Reference

Namespaces

Functions

Function Documentation