CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_2_SLHC2/src/SimDataFormats/SLHC/src/L1TowerJet.cc

Go to the documentation of this file.
00001 #include "SimDataFormats/SLHC/interface/L1TowerJet.h"
00002 #include <stdlib.h>
00003 
00004 namespace l1slhc
00005 {
00006   
00007   L1TowerJet::L1TowerJet( ):
00008   mIeta( 0 ), 
00009   mIphi( 0 ), 
00010 
00011   mE( 0 ), 
00012   mCentral( true ),
00013   mAsymEta(0),
00014   mAsymPhi(0),
00015   mWeightedIeta( 0 ),
00016   mWeightedIphi( 0 ),
00017   mJetSize( 12 ),
00018   mJetShapeType( square ),
00019   mJetArea( 144 )
00020   {
00021   }
00022   
00023   L1TowerJet::L1TowerJet( const int& aJetSize, const L1TowerJet::tJetShape& aJetShapeType , const int& aJetArea ):
00024   mIeta( 0 ), 
00025   mIphi( 0 ), 
00026 
00027   mE( 0 ), 
00028   mCentral( true ),
00029   mAsymEta(0),
00030   mAsymPhi(0),
00031   mWeightedIeta( 0 ),
00032   mWeightedIphi( 0 ),
00033 
00034   mJetSize( aJetSize ),
00035   mJetShapeType( aJetShapeType ),
00036   mJetArea( aJetArea )
00037   {
00038   }
00039   
00040   L1TowerJet::L1TowerJet( const int& aJetSize, const L1TowerJet::tJetShape& aJetShapeType , const int& aJetArea , const int &iEta, const int &iPhi  ):
00041   mIeta( iEta ), 
00042   mIphi( iPhi ), 
00043 
00044   mE( 0 ), 
00045   mCentral( true ),
00046   mAsymEta(0),
00047   mAsymPhi(0),
00048 
00049   mWeightedIeta( 0 ),
00050   mWeightedIphi( 0 ),
00051 
00052   mJetSize( aJetSize ),
00053   mJetShapeType( aJetShapeType ),
00054   mJetArea( aJetArea )
00055   {
00056   }
00057   
00058   L1TowerJet::~L1TowerJet(  )
00059   {
00060   }
00061   
00062   
00063   
00064   
00065   void L1TowerJet::setP4( const math::PtEtaPhiMLorentzVector & p4 )
00066   {
00067           mP4 = p4;
00068   }
00069   
00070   void L1TowerJet::setCentral( const bool & central )
00071   {
00072           mCentral = central;
00073   }
00074 /*
00075   void L1TowerJet::setE( const int &E )
00076   {
00077           mE = E;
00078   }
00079 */
00080 
00081 
00082   const int &L1TowerJet::iEta(  ) const
00083   {
00084           return mIeta;
00085   }
00086 
00087   const int &L1TowerJet::iPhi(  ) const
00088   {
00089           return mIphi;
00090   }
00091   
00092 
00093   
00094 
00095    const double &L1TowerJet::iWeightedEta( ) const 
00096   {
00097 
00098     return mWeightedIeta;
00099 
00100   }
00101 
00102    const double &L1TowerJet::iWeightedPhi( ) const 
00103   {
00104 
00105     return mWeightedIphi;
00106 
00107   }
00108 
00109    const double &L1TowerJet::WeightedEta( ) const 
00110   {
00111 
00112     return mWeightedEta;
00113 
00114   }
00115 
00116    const double &L1TowerJet::WeightedPhi( ) const 
00117   {
00118 
00119     return mWeightedPhi;
00120 
00121   }
00122 
00123 
00124 
00125 
00126   const int &L1TowerJet::E(  ) const
00127   {
00128           return mE;
00129   }
00130   
00131   const int& L1TowerJet::AsymEta(  ) const
00132   {
00133           return mAsymEta;
00134   }
00135   
00136   const int& L1TowerJet::AsymPhi(  ) const
00137   {
00138           return mAsymPhi;
00139   }
00140 
00141 
00142   const bool & L1TowerJet::central(  ) const
00143   {
00144           return mCentral;
00145   }
00146 
00147   const math::PtEtaPhiMLorentzVector & L1TowerJet::p4(  ) const
00148   {
00149           return mP4;
00150   }
00151   
00152   const int& L1TowerJet::JetSize(  ) const
00153   {
00154           return mJetSize;
00155   }
00156   
00157   const L1TowerJet::tJetShape& L1TowerJet::JetShape(  ) const
00158   {
00159           return mJetShapeType;
00160   }
00161 
00162   const int& L1TowerJet::JetArea(  ) const
00163   {
00164           return mJetArea;
00165   }
00166 
00167 /*
00168   double L1TowerJet::EcalVariance(  ) const
00169   {
00170     double lMean(0.0);
00171     double lMeanSq(0.0);
00172   
00173     for ( L1CaloTowerRefVector::const_iterator lConstituentIt = mConstituents.begin() ; lConstituentIt != mConstituents.end(); ++lConstituentIt ){
00174       lMean += (**lConstituentIt).E();
00175       lMeanSq += ((**lConstituentIt).E() * (**lConstituentIt).E());
00176     }
00177   
00178     lMean /= mConstituents.size();      
00179     lMeanSq /= mConstituents.size();    
00180   
00181     return lMeanSq - (lMean*lMean);
00182   }
00183   
00184   
00185   double L1TowerJet::HcalVariance(  ) const
00186   {
00187     double lMean(0.0);
00188     double lMeanSq(0.0);
00189     
00190     for ( L1CaloTowerRefVector::const_iterator lConstituentIt = mConstituents.begin() ; lConstituentIt != mConstituents.end(); ++lConstituentIt ){
00191       lMean += (**lConstituentIt).H();
00192       lMeanSq += ((**lConstituentIt).H() * (**lConstituentIt).H());
00193     }
00194     
00195     lMean /= mConstituents.size();      
00196     lMeanSq /= mConstituents.size();    
00197     
00198     return lMeanSq - (lMean*lMean);
00199   }
00200 
00201 
00202   double L1TowerJet::EnergyVariance(  ) const
00203   {
00204     double lMean( double(mE) / double(mConstituents.size()) );
00205     double lMeanSq(0.0);
00206     
00207     double lTower;
00208     for ( L1CaloTowerRefVector::const_iterator lConstituentIt = mConstituents.begin() ; lConstituentIt != mConstituents.end(); ++lConstituentIt ){
00209       lTower = (**lConstituentIt).E() + (**lConstituentIt).H();
00210       lMeanSq += ( lTower * lTower );
00211     }
00212     
00213     lMeanSq /= mConstituents.size();    
00214   
00215     return lMeanSq - (lMean*lMean);
00216   }
00217 */
00218 
00219 
00220   double L1TowerJet::EcalMAD() const
00221   {
00222     std::deque< int > lEnergy;
00223     for ( L1CaloTowerRefVector::const_iterator lConstituentIt = mConstituents.begin() ; lConstituentIt != mConstituents.end(); ++lConstituentIt ){
00224       lEnergy.push_back( (**lConstituentIt).E() );
00225     }
00226     lEnergy.resize( mJetArea , 0 );
00227     return MAD( lEnergy );
00228   
00229   }
00230   
00231   double L1TowerJet::HcalMAD() const
00232   {
00233     std::deque< int > lEnergy;
00234     for ( L1CaloTowerRefVector::const_iterator lConstituentIt = mConstituents.begin() ; lConstituentIt != mConstituents.end(); ++lConstituentIt ){
00235       lEnergy.push_back( (**lConstituentIt).H() );
00236     }
00237     lEnergy.resize( mJetArea , 0 );
00238     return MAD( lEnergy );
00239   
00240   }
00241   
00242 
00243   void L1TowerJet::CalcWeightediEta() 
00244   {
00245     double etaSumEt(0); 
00246     double sumEt (0);
00247 
00248     for ( L1CaloTowerRefVector::const_iterator lConstituentIt = mConstituents.begin() ; lConstituentIt != mConstituents.end(); ++lConstituentIt ){
00249       etaSumEt += ( (**lConstituentIt).E() + (**lConstituentIt).H() ) * ( (**lConstituentIt).iEta() );
00250       sumEt += ( (**lConstituentIt).E() + (**lConstituentIt).H() ) ;
00251     }
00252 //      std::cout<<" eta* energy = "<<etaSumEt<<" sum energy: "<<sumEt<<std::endl;
00253     
00254     //discretize weighted Ieta: run from |1-28| (no HF)
00255     mWeightedIeta = etaSumEt/sumEt ; 
00256 
00257     //discrete-ize and account for the fact there is no zero
00258     int discrete_iEta(999);
00259     //add 0.5 so that we get rounding up as well as down
00260 
00261     if( mWeightedIeta>=0 ) discrete_iEta=int(mWeightedIeta+0.5);
00262     else  discrete_iEta=int(mWeightedIeta-0.5); 
00263    
00264     //account for the fact there is no 0
00265     if(mWeightedIeta>=0 && mWeightedIeta<1)      discrete_iEta=1;
00266   
00267     if(mWeightedIeta<0 && mWeightedIeta>(-1))    discrete_iEta=-1;
00268 
00269     //std::cout<<"weighted ieta: "<<mWeightedIeta <<" discrete ieta: "<< discrete_iEta<<std::endl;
00270     mWeightedIeta = double(discrete_iEta);
00271 
00272   }
00273 
00274 
00275   // --------------------------------------------------
00276   // modified version to calculate weighted eta by first converting iEta to eta, than weighting
00277   void L1TowerJet::calculateWeightedEta()
00278   {
00279     
00280     const double endcapEta[8] = { 0.09, 0.1, 0.113, 0.129, 0.15, 0.178, 0.15, 0.35 };
00281     
00282     double etaSumEt(0); 
00283     double sumEt(0);
00284 
00285     double tmpEta(9999);
00286     double abs_eta(9999);
00287 
00288     for (L1CaloTowerRefVector::const_iterator lConstituentIt = mConstituents.begin() ; lConstituentIt != mConstituents.end(); ++lConstituentIt ) {
00289 
00290       abs_eta = fabs((**lConstituentIt).iEta());
00291       if (abs_eta < 21) tmpEta = (0.087*abs_eta - 0.0435);
00292       else {
00293         abs_eta -= 21;
00294         tmpEta = 1.74;
00295         
00296         for (int i=0; i!=int(abs_eta); ++i) tmpEta += endcapEta[i];
00297       
00298         if (mJetSize % 2 == 0) tmpEta += endcapEta[int(abs_eta)] / 2.;
00299         else tmpEta += endcapEta[int(abs_eta)];
00300       }
00301       if ((**lConstituentIt).iEta()<0) tmpEta = (-1)*tmpEta;
00302 
00303       etaSumEt += ( (**lConstituentIt).E() + (**lConstituentIt).H() ) * ( tmpEta );
00304       sumEt += ( (**lConstituentIt).E() + (**lConstituentIt).H() ) ;
00305     }
00306     
00307     mWeightedEta = etaSumEt/sumEt ; 
00308 
00309   }
00310 
00311 
00312   /*
00313   // old version of calculating weighted eta using as input weighted iEta (this does not give the right result)
00314   void L1TowerJet::calculateWeightedEta()
00315   {
00316     double WeightedEta(9999);
00317     double abs_eta =fabs(mWeightedIeta);
00318    if ( abs_eta < 21 )
00319   {
00320     //with discrete iEta this first option should always be the case
00321     if( abs_eta >=1) WeightedEta = (0.087*abs_eta - 0.0435);
00322     else WeightedEta = 0.0435*abs_eta ;
00323 
00324   }
00325   else
00326   {
00327     const double endcapEta[8] = { 0.09, 0.1, 0.113, 0.129, 0.15, 0.178, 0.15, 0.35 };
00328     abs_eta -= 21;
00329 
00330     WeightedEta = 1.74;
00331 
00332     for ( int i = 0; i != int(abs_eta); ++i )
00333     {
00334       WeightedEta += endcapEta[i];
00335     }
00336     if( mJetSize % 2 == 0 ){
00337       WeightedEta += endcapEta[int(abs_eta)] / 2.;
00338     }else{
00339       WeightedEta += endcapEta[int(abs_eta)];
00340     }
00341   }
00342  if(mWeightedIeta<0) WeightedEta=-WeightedEta;
00343 
00344   mWeightedEta = WeightedEta;
00345   }
00346   */
00347 
00348 
00349   void L1TowerJet::CalcWeightediPhi() 
00350   {
00351     double phiSumEt(0); 
00352     double sumEt (0);
00353 
00354     for ( L1CaloTowerRefVector::const_iterator lConstituentIt = mConstituents.begin() ; lConstituentIt != mConstituents.end(); ++lConstituentIt ){
00355      double tower_iPhi =(**lConstituentIt).iPhi(); 
00356      if( iPhi() >= (72-JetSize()) ) { //constituents may go over edge, iPhi>66 for 8x8 jet
00357         if( tower_iPhi  < (72 - JetSize()) ){//if constituent tower is over edge, ie. iPhi>1 
00358            tower_iPhi+=72;  //iPhi=1 -> iPhi=73
00359           }
00360       }
00361       //calculate weighted phi using corrected iPhi value 
00362       phiSumEt += ( (**lConstituentIt).E() + (**lConstituentIt).H() ) * (tower_iPhi  );
00363 
00364       sumEt += ( (**lConstituentIt).E() + (**lConstituentIt).H() ) ;
00365     }
00366 //      std::cout<<"phi sum et: "<<phiSumEt<<"sum Et: "<<sumEt<<std::endl;
00367     mWeightedIphi = phiSumEt/sumEt ; 
00368 
00369     // take weighted Iphi back to 1-72 range
00370     if(mWeightedIphi>72) mWeightedIphi-=72;
00371 
00372 
00373   }
00374 
00375 void L1TowerJet::calculateWeightedPhi( )
00376 {
00377 //  double JetSize = double(mJetSize) / 2.0;
00378   double WeightedPhi = ( ( mWeightedIphi-0.5 ) * 0.087 );
00379   //Need this because 72*0.087 != 2pi: else get uneven phi dist
00380   double pi=(72*0.087)/2;
00381   if(WeightedPhi > pi) WeightedPhi -=2*pi;
00382   mWeightedPhi=WeightedPhi;
00383   //  std::cout<<"weighted IPhi: "<<mWeightedIphi<<" weighted phi: "<<WeightedPhi<<std::endl;
00384 }
00385 
00386 
00387 
00388 
00389 
00390 
00391   double L1TowerJet::EnergyMAD() const
00392   {
00393     std::deque< int > lEnergy;
00394     for ( L1CaloTowerRefVector::const_iterator lConstituentIt = mConstituents.begin() ; lConstituentIt != mConstituents.end(); ++lConstituentIt ){
00395       lEnergy.push_back( (**lConstituentIt).E() + (**lConstituentIt).H() );
00396     }
00397     lEnergy.resize( mJetArea , 0 );
00398     return MAD( lEnergy );
00399   }
00400 
00401 
00402 
00403   double L1TowerJet::MAD( std::deque<int>& aDataSet ) const
00404   {
00405     std::sort( aDataSet.begin() , aDataSet.end() );
00406     
00407     std::size_t lDataSetSize( aDataSet.size() );
00408     
00409     double lMedian(0);
00410     if( lDataSetSize%2 == 0 ){
00411       lMedian = double ( aDataSet[ (lDataSetSize/2) - 1 ] + aDataSet[ lDataSetSize/2 ] ) / 2.0 ;
00412     }else{
00413       lMedian = double( aDataSet[ (lDataSetSize-1)/2 ] );
00414     }
00415     
00416     
00417     std::deque< double > lMedianSubtractedDataSet;
00418     for ( std::deque< int >::const_iterator lIt = aDataSet.begin() ; lIt != aDataSet.end(); ++lIt ){
00419       lMedianSubtractedDataSet.push_back( fabs( double(*lIt) - lMedian ) );
00420     }
00421     
00422     std::sort( lMedianSubtractedDataSet.begin() , lMedianSubtractedDataSet.end() );
00423     
00424     if( lDataSetSize%2 == 0 ){
00425       return double ( lMedianSubtractedDataSet[ (lDataSetSize/2) - 1 ] + lMedianSubtractedDataSet[ lDataSetSize/2 ] ) / 2.0 ;
00426     }else{
00427       return double( lMedianSubtractedDataSet[ (lDataSetSize-1)/2 ] );
00428     }
00429   
00430   }
00431 
00432 
00433   void L1TowerJet::addConstituent( const L1CaloTowerRef & Tower )
00434   { 
00435   
00436     int lHalfJetSizeEta( mJetSize >> 1 );
00437     int lHalfJetSizePhi( mJetSize >> 1 );
00438 
00439     int lTowerEnergy( Tower->E(  ) + Tower->H(  ) );
00440 
00441     //slightly different sizes for HF jets
00442     if( abs( iEta() ) > 28 ){
00443       lHalfJetSizeEta =1; //ie mJetSize/4 as in HF jets 2 in eta
00444     }
00445   
00446     //if iPhi at the edge of the calo wrap round in phi
00447     int ToweriPhi = Tower->iPhi(  );
00448     if(iPhi() > (72 - mJetSize)) { 
00449       if(ToweriPhi>=1 && ToweriPhi<mJetSize) ToweriPhi+=72;
00450     }
00451 
00452     mE += lTowerEnergy;
00453     mConstituents.push_back( Tower );
00454    /* if( (Tower->iPhi())<4){
00455      std::cout<<"JET: ("<<iEta()<<","<<iPhi()<<")"<<std::endl;
00456      std::cout<<"TOWER  = ( " <<Tower->iEta()<< " , " <<Tower->iPhi()<<" ), energy = " <<Tower->E()+Tower->H()<<std::endl;
00457   }*/
00458     //when add a tower, also add the asymmetry terms
00459     //positive asym: top RH corner of jet
00460   
00461     if( mJetSize % 2 == 0 ){ //even jet size
00462     
00463     //if(Tower->iEta(  ) == iEta() ) mAsymEta += 0; //do nothing 
00464       if( Tower->iEta(  ) >= (iEta() + lHalfJetSizeEta) ) {
00465         mAsymEta +=  lTowerEnergy; 
00466       }         
00467       else{ /*if( Tower->iEta(  ) <  iEta() + lHalfJetSize )*/ 
00468         mAsymEta -= lTowerEnergy; 
00469   
00470       }
00471   
00472     //if(  Tower->iPhi(  ) == iPhi() ) mAsymEta += 0; //do nothing
00473 
00474       if( ToweriPhi < (iPhi() + lHalfJetSizePhi) ){
00475         mAsymPhi += lTowerEnergy;
00476   
00477       }
00478       else{ /*if( Tower->iPhi(  ) > iPhi() + lHalfJetSize )*/  
00479         mAsymPhi -= lTowerEnergy;
00480   
00481   
00482       }
00483   
00484     }else{ //odd jet size: miss out central towers
00485       
00486       if( Tower->iEta(  ) ==  (lHalfJetSizeEta + iEta()) ) {
00487         mAsymEta += 0; //do nothing
00488       }
00489       else if( Tower->iEta(  ) > (iEta() + lHalfJetSizeEta) ) {
00490         mAsymEta +=  lTowerEnergy;
00491   
00492       }else /*if( Tower->iEta(  ) <  iEta() + lHalfJetSize )*/ {
00493         mAsymEta -= lTowerEnergy;
00494   
00495       }
00496     // else it is in the middle so does not contribute to the asymmetry
00497   
00498       if( ToweriPhi == (lHalfJetSizePhi + iPhi()) ) {
00499         mAsymPhi += 0; //do nothing
00500   
00501       }
00502       else if( ToweriPhi < (iPhi() + lHalfJetSizePhi) ) {
00503         mAsymPhi += lTowerEnergy;
00504   
00505       }else /*if( Tower->iPhi(  ) > iPhi() + lHalfJetSize )*/  {
00506         mAsymPhi -= lTowerEnergy;
00507   
00508       }
00509   // else it is in the middle so does not contribute to the asymmetry
00510 
00511     }
00512   } 
00513 
00514 
00515   void L1TowerJet::removeConstituent( const int &eta, const int &phi )
00516   {
00517     L1CaloTowerRefVector::iterator lConstituent = getConstituent( eta, phi );
00518     if ( lConstituent != mConstituents.end() )
00519     {
00520       int lHalfJetSizeEta( mJetSize >> 1 );
00521       int lHalfJetSizePhi( mJetSize >> 1 );
00522       int lTowerEnergy( (**lConstituent).E(  ) + (**lConstituent).H(  ) );
00523   
00524       mE -= lTowerEnergy;
00525       mConstituents.erase( lConstituent );
00526     
00527       if( abs( iEta() ) > 28 ){
00528         lHalfJetSizeEta =1; //ie mJetSize/4 as in HF jets 2 in eta
00529       }
00530       int ToweriPhi = phi;
00531       //if iPhi at the edge of the calo wrap round in phi
00532       if(iPhi() > (72 - mJetSize)) { 
00533         if(ToweriPhi>=1 && ToweriPhi < mJetSize) ToweriPhi+=72;
00534       }
00535         
00536       
00537       if( mJetSize % 2 == 0 ){ //even jet size
00538   
00539         if( eta >= (iEta() + lHalfJetSizeEta) ) {
00540           mAsymEta -=  lTowerEnergy;
00541         }               
00542         else{ /*if( Tower->iEta(  ) <  iEta() + lHalfJetSize )*/ 
00543           mAsymEta += lTowerEnergy;
00544         }
00545   
00546   
00547         if( ToweriPhi < (iPhi() + lHalfJetSizePhi) ){
00548           mAsymPhi -= lTowerEnergy;
00549   
00550         }else{ /*if( Tower->iPhi(  ) > iPhi() + lHalfJetSize )*/  
00551           mAsymPhi += lTowerEnergy;
00552         }
00553   
00554   
00555       }else{ //odd jet size: miss out central towers
00556     
00557         if( eta ==  (lHalfJetSizeEta + iEta()) ) {
00558           mAsymEta += 0; //do nothing
00559         }
00560         else if( eta > (iEta() + lHalfJetSizeEta) ) {
00561           mAsymEta -=  lTowerEnergy;
00562   
00563         }else /*if( Tower->iEta(  ) <  iEta() + lHalfJetSize )*/ {
00564           mAsymEta += lTowerEnergy;
00565   
00566         }
00567   
00568         if( ToweriPhi == (lHalfJetSizePhi + iPhi()) ) {
00569           mAsymEta -= 0; //do nothing
00570   
00571         }
00572         else if( ToweriPhi < (iPhi() + lHalfJetSizePhi) ) {
00573           mAsymPhi -= lTowerEnergy;
00574   
00575         }else /*if( Tower->iPhi(  ) > iPhi() + lHalfJetSize )*/  {
00576           mAsymPhi += lTowerEnergy;
00577   
00578         }
00579       }
00580     }
00581   }
00582 
00583 
00584   const L1CaloTowerRefVector & L1TowerJet::getConstituents(  ) const
00585   {
00586     return mConstituents;
00587   }
00588   
00589   
00590   L1CaloTowerRefVector::iterator L1TowerJet::getConstituent( const int &eta, const int &phi )
00591   {
00592     for ( L1CaloTowerRefVector::iterator lConstituentIt = mConstituents.begin() ; lConstituentIt != mConstituents.end(); ++lConstituentIt )
00593       if ( (**lConstituentIt).iEta(  ) == eta + mIeta && (**lConstituentIt).iPhi(  ) == phi + mIphi )
00594         return lConstituentIt;
00595   
00596     return mConstituents.end();
00597   }
00598 
00599 }
00600 
00601 
00602 
00603 
00604 namespace std
00605 {
00606         bool operator<( const l1slhc::L1TowerJet & aLeft, const l1slhc::L1TowerJet & aRight )
00607         {
00608                 if ( aLeft.E(  ) == aRight.E(  ) )
00609                 {
00610                         // for two objects with equal energy, favour the more central one
00611                         return ( abs( aLeft.iEta(  ) ) > abs( aRight.iEta(  ) ) );
00612                 }
00613                 else
00614                 {
00615                         return ( aLeft.E(  ) < aRight.E(  ) );
00616                 }
00617         }
00618 }
00619 
00620 
00621 // pretty print
00622 std::ostream & operator<<( std::ostream & aStream, const l1slhc::L1TowerJet & aL1TowerJet )
00623 {
00624         aStream << "L1TowerJet" 
00625                 << " iEta=" << aL1TowerJet.iEta(  ) 
00626                 << " iPhi=" << aL1TowerJet.iPhi(  ) 
00627                 << "\n with constituents:\n";
00628         for ( l1slhc::L1CaloTowerRefVector::const_iterator i = aL1TowerJet.getConstituents(  ).begin(  ); i < aL1TowerJet.getConstituents(  ).end(  ); ++i )
00629                 aStream << "  iEta=" << ( **i ).iEta(  ) 
00630                         << " iPhi=" << ( **i ).iPhi(  ) 
00631                         << " ET=" << ( **i ).E(  ) 
00632                         << "\n";
00633         return aStream;
00634 }