/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86

Bug Summary

File:	alld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h
Warning:	line 398, column 7 Null pointer passed to 2nd parameter expecting 'nonnull'

Annotated Source Code

Press '?' to see keyboard shortcuts

Show analyzer invocation

clang -cc1 -cc1 -triple x86_64-unknown-linux-gnu -analyze -disable-free -main-file-name DCCALShower_factory.cc -analyzer-store=region -analyzer-opt-analyze-nested-blocks -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=cplusplus -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -fhalf-no-semantic-interposition -mframe-pointer=none -fmath-errno -fno-rounding-math -mconstructor-aliases -munwind-tables -target-cpu x86-64 -tune-cpu generic -fno-split-dwarf-inlining -debugger-tuning=gdb -resource-dir /w/halld-scifs17exp/home/sdobbs/clang/llvm-project/install/lib/clang/12.0.0 -D HAVE_CCDB -D HAVE_RCDB -D HAVE_EVIO -D HAVE_TMVA=1 -D RCDB_MYSQL=1 -D RCDB_SQLITE=1 -D SQLITE_USE_LEGACY_STRUCT=ON -I .Linux_CentOS7.7-x86_64-gcc4.8.5/libraries/CCAL -I libraries/CCAL -I . -I libraries -I libraries/include -I /w/halld-scifs17exp/home/sdobbs/clang/halld_recon/Linux_CentOS7.7-x86_64-gcc4.8.5/include -I external/xstream/include -I /usr/include/tirpc -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/root/root-6.08.06/include -I /w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/ccdb/ccdb_1.06.06/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/rcdb/rcdb_0.06.00/cpp/include -I /usr/include/mysql -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/sqlitecpp/SQLiteCpp-2.2.0^bs130/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/sqlite/sqlite-3.13.0^bs130/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/hdds/hdds-4.9.0/Linux_CentOS7.7-x86_64-gcc4.8.5/src -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/xerces-c/xerces-c-3.1.4/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/evio/evio-4.4.6/Linux-x86_64/include -internal-isystem /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5 -internal-isystem /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5/x86_64-redhat-linux -internal-isystem /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5/backward -internal-isystem /usr/local/include -internal-isystem /w/halld-scifs17exp/home/sdobbs/clang/llvm-project/install/lib/clang/12.0.0/include -internal-externc-isystem /include -internal-externc-isystem /usr/include -O2 -std=c++11 -fdeprecated-macro -fdebug-compilation-dir /home/sdobbs/work/clang/halld_recon/src -ferror-limit 19 -fgnuc-version=4.2.1 -fcxx-exceptions -fexceptions -vectorize-loops -vectorize-slp -analyzer-output=html -faddrsig -o /tmp/scan-build-2021-01-21-110224-160369-1 -x c++ libraries/CCAL/DCCALShower_factory.cc

libraries/CCAL/DCCALShower_factory.cc

→


2/*
*  File: DCCALHit_factory.cc
*
* Created on 11/25/18 by A.S.
* use structure similar to FCAL
*/

9#include "CCAL/DCCALShower_factory.h"


12static mutex CCAL_MUTEX;
13//static bool CCAL_PROFILE_LOADED = false;


16//==========================================================
17//
18//   Constructor
19//
20//==========================================================

22DCCALShower_factory::DCCALShower_factory()
23{
// Set defaults:

  	MIN_CLUSTER_BLOCK_COUNT   =  2;
MAX_HITS_FOR_CLUSTERING   =  80;
  	MIN_CLUSTER_SEED_ENERGY   =  0.035;  /* [GeV] */
MIN_CLUSTER_ENERGY        =  0.05;   /* [GeV] */
MAX_CLUSTER_ENERGY        =  15.9;   /* [GeV] */
TIME_CUT                  =  15.0;   /*  [ns] */

SHOWER_DEBUG              =  0;
DO_NONLINEAR_CORRECTION   =  1;

CCAL_RADIATION_LENGTH     =  0.86;
CCAL_CRITICAL_ENERGY      =  1.1e-3;

LOG_POS_CONST             =  4.2;


gPARMS->SetDefaultParameter( "CCAL:SHOWER_DEBUG", SHOWER_DEBUG );
gPARMS->SetDefaultParameter( "CCAL:MIN_CLUSTER_BLOCK_COUNT", MIN_CLUSTER_BLOCK_COUNT, 
	"minimum number of blocks to form a cluster" );
gPARMS->SetDefaultParameter( "CCAL:MIN_CLUSTER_SEED_ENERGY", MIN_CLUSTER_SEED_ENERGY, 
	"minimum energy for a block to be considered as a seed for a cluster" );
gPARMS->SetDefaultParameter( "CCAL:MIN_CLUSTER_ENERGY", MIN_CLUSTER_ENERGY, 
	"minimum allowed cluster energy" );
gPARMS->SetDefaultParameter( "CCAL:MAX_CLUSTER_ENERGY", MAX_CLUSTER_ENERGY, 
	"maximum allowed cluster energy" );
gPARMS->SetDefaultParameter( "CCAL:MAX_HITS_FOR_CLUSTERING", MAX_HITS_FOR_CLUSTERING, 
	"maximum hits allowed to call clustering algorithm" );
gPARMS->SetDefaultParameter( "CCAL:TIME_CUT", TIME_CUT, 
	"time cut for associating CCAL hits together into a cluster" );
gPARMS->SetDefaultParameter( "CCAL:DO_NONLINEAR_CORRECTION", DO_NONLINEAR_CORRECTION, 
	"set this to zero when no nonlinear correction is desired" );
gPARMS->SetDefaultParameter( "CCAL:CCAL_RADIATION_LENGTH", CCAL_RADIATION_LENGTH );
gPARMS->SetDefaultParameter( "CCAL:CCAL_CRITICAL_ENERGY", CCAL_CRITICAL_ENERGY );
gPARMS->SetDefaultParameter( "CCAL:LOG_POS_CONST", LOG_POS_CONST );

VERBOSE = 0;              ///< >0 once off info ; >2 event by event ; >3 everything
gPARMS->SetDefaultParameter("DFCALShower:VERBOSE", VERBOSE, "Verbosity level for DFCALShower objects and factories");
63}




68//==========================================================
69//
70//   brun
71//
72//==========================================================

74jerror_t DCCALShower_factory::brun(JEventLoop *locEventLoop, int32_t runnumber)
75{
  // Only print messages for one thread whenever run number change
  static pthread_mutex_t print_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } };
  static set<int> runs_announced;
  pthread_mutex_lock(&print_mutex);
  bool print_messages = false;
  if(runs_announced.find(runnumber) == runs_announced.end()){
      print_messages = true;
      runs_announced.insert(runnumber);
  }
  pthread_mutex_unlock(&print_mutex);


DApplication *dapp = dynamic_cast<DApplication*>(eventLoop->GetJApplication());
  	const DGeometry *geom = dapp->GetDGeometry(runnumber);

if (geom) {
    	  geom->GetTargetZ(m_zTarget);
    	  geom->GetCCALPosition(m_CCALdX,m_CCALdY,m_CCALfront);
  	}
  	else{
    	  jerr << "No geometry accessible." << endl;
    	  return RESOURCE_UNAVAILABLE;
  	}


//------------------------------------------------------//
//-----------   Read in shower profile data  -----------//

std::unique_lock<std::mutex> lck(CCAL_MUTEX);

string ccal_profile_file;
gPARMS->SetDefaultParameter("CCAL_PROFILE_FILE", ccal_profile_file, 
"CCAL profile data file name");

// follow similar procedure as other resources (DMagneticFieldMapFineMesh)

map< string,string > profile_file_name;
JCalibration *jcalib = dapp->GetJCalibration(runnumber);

if( jcalib->GetCalib("/CCAL/profile_data/profile_data_map", profile_file_name) ) 
{
 jout << "Can't find requested /CCAL/profile_data/profile_data_map in CCDB for this run!"
 	<< endl;

} else if( profile_file_name.find("map_name") != profile_file_name.end() 
&& profile_file_name["map_name"] != "None" ) 
{
 JResourceManager *jresman = dapp->GetJResourceManager(runnumber);
 ccal_profile_file = jresman->GetResource(profile_file_name["map_name"]);
}

if(print_messages)
jout<<"Reading CCAL profile data from "<<ccal_profile_file<<" ..."<<endl;

// check to see if we actually have a file
if(ccal_profile_file.empty()) 
{
 if(print_messages)
   jerr << "Empty file..." << endl;
 return RESOURCE_UNAVAILABLE;
}

ifstream ccal_profile(ccal_profile_file.c_str());
for(int i=0; i<=500; i++) {
    for(int j=0; j<=i; j++) {
	int id1, id2;
	double fcell_hyc, fd2c;
	
	ccal_profile >> id1 >> id2 >> fcell_hyc >> fd2c;

	acell[id1][id2] = fcell_hyc;
	acell[id2][id1] = fcell_hyc;
	ad2c[id1][id2] = fd2c;
	ad2c[id2][id1] = fd2c;
}
}
ccal_profile.close();

lck.unlock();




//------------------------------------------------------//
//----------  Initialize channel status array ----------//

for(int icol = 0; icol < MCOL12; ++icol) {
 for(int irow = 0; irow < MROW12; ++irow) {
   if(icol>=5 && icol<=6 && irow>=5 && irow<=6) { stat_ch[irow][icol] = -1; }
   else { stat_ch[irow][icol] = 0; }
 }
}




//------------------------------------------------------//
//----------  Read shower timewalk parameters ----------//

vector< vector<double> > timewalk_params;
if( eventLoop->GetCalib("/CCAL/shower_timewalk_correction",timewalk_params) )
 jout << "Error loading /CCAL/shower_timewalk_correction !" << endl;
else {
 if( (int)timewalk_params.size() != 2 ) {
   cout << "DCCALShower_factory: Wrong number of entries to timewalk correction table (should be 144)." << endl;
   for( int ii = 0; ii < 2; ++ii ) {
     timewalk_p0.push_back(0.0);
     timewalk_p1.push_back(0.0);
     timewalk_p2.push_back(0.0);
     timewalk_p3.push_back(0.0);
   }
 } else {
   for( vector< vector<double> >::const_iterator iter = timewalk_params.begin(); 
   	iter != timewalk_params.end(); ++iter ) {
     if( iter->size() != 4 ) {
       cout << "DCCALShower_factory: Wrong number of values in timewalk correction table (should be 4)" << endl;
              continue;
     }
   
     timewalk_p0.push_back( (*iter)[0] );
     timewalk_p1.push_back( (*iter)[1] );
     timewalk_p2.push_back( (*iter)[2] );
     timewalk_p3.push_back( (*iter)[3] );
   
   }
 }
}




//------------------------------------------------------//
//---------  Read in non-linearity parameters ----------//

vector< vector<double> > nonlin_params;
if( eventLoop->GetCalib("/CCAL/nonlinear_energy_correction",nonlin_params) )
 jout << "Error loading /CCAL/nonlinear_energy_correction !" << endl;
else {
 if( (int)nonlin_params.size() != CCAL_CHANS ) {
   cout << "DCCALShower_factory: Wrong number of entries to nonlinear energy correction table (should be 144)." << endl;
   for( int ii = 0; ii < CCAL_CHANS; ++ii ) {
     Nonlin_p0.push_back(0.0);
     Nonlin_p1.push_back(0.0);
     Nonlin_p2.push_back(0.0);
     Nonlin_p3.push_back(0.0);
   }
 } else {
   for( vector< vector<double> >::const_iterator iter = nonlin_params.begin(); 
   	iter != nonlin_params.end(); ++iter ) {
     if( iter->size() != 4 ) {
       cout << "DCCALShower_factory: Wrong number of values in nonlinear energy correction table (should be 4)" << endl;
              continue;
     }
   
     Nonlin_p0.push_back( (*iter)[0] );
     Nonlin_p1.push_back( (*iter)[1] );
     Nonlin_p2.push_back( (*iter)[2] );
     Nonlin_p3.push_back( (*iter)[3] );
   
   }
 }
}

if( SHOWER_DEBUG ) {
 cout << "\n\nNONLIN_P0  NONLIN_P1  NONLIN_P2  NONLIN_P3" << endl;
 for(int ii = 0; ii < (int)Nonlin_p0.size(); ii++) {
   cout << Nonlin_p0[ii] << " " << Nonlin_p1[ii] << " " << Nonlin_p2[ii] << " " << 
   	Nonlin_p3[ii] << endl;
 }
 cout << "\n\n";
}




return NOERROR;
252}




257//==========================================================
258//
259//   evnt
260//
261//==========================================================

263jerror_t DCCALShower_factory::evnt(JEventLoop *locEventLoop, uint64_t eventnumber)
264{


vector< const DCCALGeometry* > ccalGeomVect;
  	locEventLoop->Get( ccalGeomVect );
1
Calling 'JEventLoop::Get'→
if (ccalGeomVect.size() < 1)
    	  return OBJECT_NOT_AVAILABLE;
  	const DCCALGeometry& ccalGeom = *(ccalGeomVect[0]);




//------- Get the CCALHits and organize hit pattern -------//


vector< const DCCALHit* > ccalhits;
locEventLoop->Get( ccalhits );

int n_hits = static_cast<int>( ccalhits.size() );
if( n_hits < 1 || n_hits > MAX_HITS_FOR_CLUSTERING ) return NOERROR;


vector< vector< const DCCALHit* > > hitPatterns;
getHitPatterns( ccalhits, hitPatterns );

int n_patterns = static_cast<int>( hitPatterns.size() );
if( n_patterns < 1 ) return NOERROR;

vector< ccalcluster_t > ccalClusters; // will hold all clusters
vector< cluster_t > clusterStorage;   // will hold the constituents of every cluster




//-------  Call island clusterizer on each pattern  -------//


for( int ipat = 0; ipat < n_patterns; ipat++ ) 
{

 
 /*-----  prepare data in dimensionless format  -----*/
 
 vector< const DCCALHit* > rawHitPattern = hitPatterns[ipat];
 vector< const DCCALHit* > locHitPattern;
 cleanHitPattern( rawHitPattern, locHitPattern );
 
 int n_hits = static_cast<int>( locHitPattern.size() );
 
 vector< int > ia;
 vector< int > id;
 	  
 for( int ih = 0; ih < n_hits; ih++ ) {
   const DCCALHit *ccalhit = locHitPattern[ih];
   int row = 12 - ccalhit->row;
   int col = 12 - ccalhit->column;
   int ie  = static_cast<int>( ccalhit->E*10. + 0.5 );
   if( ie > 0 ) {
     int address = 100*col + row;
     ia.push_back( address );
     id.push_back( ie );
   }
 }
 
 
 /*-------------     call to island     -------------*/
 	  
 vector< gamma_t > gammas; gammas.clear();// Output of main_island (holds reconstructed photons)
 main_island( ia, id, gammas );


 /*------------  post-island processing  ------------*/
 	
 int init_clusters = static_cast<int>( gammas.size() );
 if( !init_clusters ) continue;
 
 processShowers( gammas, ccalGeom, locHitPattern, eventnumber, 
 	ccalClusters, clusterStorage );
 
 
} // end looping over hit patterns




//--------------   Fill DCCALShower Object   --------------//

int n_clusters = static_cast<int>( ccalClusters.size() );

for( int k = 0; k < n_clusters; k++ ) {

 DCCALShower *shower = new DCCALShower;
 
 shower->E        =   ccalClusters[k].E;
 shower->Esum     =   ccalClusters[k].Esum;
 
 shower->x        =   ccalClusters[k].x+m_CCALdX;
 shower->y        =   ccalClusters[k].y+m_CCALdY;
 shower->x1       =   ccalClusters[k].x1+m_CCALdX;
 shower->y1       =   ccalClusters[k].y1+m_CCALdY;
 shower->z        =   m_CCALfront;
 
 shower->chi2     =   ccalClusters[k].chi2;
 shower->sigma_E  =   ccalClusters[k].sigma_E;
 shower->Emax     =   ccalClusters[k].emax;
 shower->time     =   ccalClusters[k].time;
 
 shower->dime     =   ccalClusters[k].nhits;
 shower->idmax    =   ccalClusters[k].idmax;
 shower->id       =   ccalClusters[k].id;
 
 shower->type     =   ccalClusters[k].type;
 
 int showTypeVal     = (ccalClusters[k].type)/10;	  
 shower->ClusterType = static_cast<ClusterType_t>( showTypeVal/2 );
 shower->PeakType    = static_cast<PeakType_t>( showTypeVal%2 );
 //shower->ExyztCovariance(i, j)
 float xerr = 0.1016;// / sqrt(ccalClusters[k].E) + 0.2219; // HARD CODED VALUE!!!!   
 float yerr = 0.1016;// / sqrt(ccalClusters[k].E) + 0.2219; // HARD CODED VALUE!!!!   
 float zerr = 2.5; // HARD CODED VALUE!!!!   
 float terr = 0.2; // HARD CODED VALUE!!!!   
 float eerr = 0.5;//pow(ccalClusters[k].E, 2) * (0.01586 / ccalClusters[k].E + 0.0002342 / pow(ccalClusters[k].E, 2) + 1.696e-6); // HARD CODED VALUE!!!!   
 //copy xyz errors into covariance matrix
 shower->ExyztCovariance.ResizeTo(5,5);
 for (int col = 0; col < 5; col ++) {
   for (int row = 0; row < 5; row ++) {
     if (col != row) shower->ExyztCovariance[col][row] = 0;
   }
 }
 shower->ExyztCovariance[0][0] = pow(xerr, 2);
 shower->ExyztCovariance[1][1] = pow(yerr, 2);
 shower->ExyztCovariance[2][2] = pow(zerr, 2);
 shower->ExyztCovariance[3][3] = pow(terr, 2);
 shower->ExyztCovariance[4][4] = pow(eerr, 2);
 
 if (VERBOSE>2) {printf("(E,x,y,z,t)    "); shower->ExyztCovariance.Print(); }

 for( int icell = 0; icell < ccalClusters[k].nhits; icell++ ) {
   
   int hitID   = clusterStorage[k].id[icell];
   float hitX = ccalGeom.positionOnFace( hitID ).X();
   float hitY = ccalGeom.positionOnFace( hitID ).Y();
   int hitROW  = ccalGeom.row( hitY );
   int hitCOL  = ccalGeom.column( hitX );
   	    
   DCCALHit *clusHit = new DCCALHit;
   clusHit->row    = hitROW;
   clusHit->column = hitCOL;
   clusHit->x      = hitX;
   clusHit->y      = hitY;
   clusHit->E      = static_cast<float>( 1000.*clusterStorage[k].E[icell] );
   clusHit->t      = static_cast<float>( clusterStorage[k].t[icell] );
   
   shower->AddAssociatedObject( clusHit );
   
 }

 _data.push_back( shower );
}


return NOERROR;

427}




432//==========================================================
433//
434//   getHitPatterns
435//
436//==========================================================

438void DCCALShower_factory::getHitPatterns( vector< const DCCALHit* > hitarray, 
vector< vector< const DCCALHit* > > &hitPatterns ) 
440{

/*------------------------------

Method for sorting hit patterns:

1. Find hit with largest energy.
2. Sort hit vector in order of increasing time-difference from this hit
3. Push all hits within 15 ns of max hit to a new vector (first few elements of sorted vec)
4. Erase the elements that were moved from the previous vector
- maybe steps 3&4 can be done simultaneously? 
5. Push the new vector back to the hitPatterns vector
6. Repeat steps 1-5 until no hits remain in original hitarray

--------------------------------*/


int n_hits = static_cast<int>( hitarray.size() );

if( n_hits < 1 ) return;
if( n_hits < 2 ) {
 vector< const DCCALHit* > hitVec;
 hitVec.push_back( hitarray[0] );
 hitPatterns.push_back( hitVec );
 return;
}


vector< const DCCALHit* > clonedHitArray = hitarray;

while( clonedHitArray.size() ) 
{

 vector< const DCCALHit* > locHitVec;

 float maxE  = -1.;
 float maxT  = 1.e6;
 
 for( unsigned int ih = 0; ih < clonedHitArray.size(); ih++ ) {
   float trialE = clonedHitArray[ih]->E;
   if( trialE > maxE ) { maxE = trialE; maxT = clonedHitArray[ih]->t; }
 }
 
 if( maxE < 0. ) break;
 
 
 sortByTime( clonedHitArray, maxT );
 
 
 n_hits = static_cast<int>( clonedHitArray.size() );
 
 int n_good_hits = 0;
 for( int ih = 0; ih < n_hits; ih++ ) {
 
   const DCCALHit *locHit = clonedHitArray[ih];
   float timeDiff = fabs( locHit->t - maxT );
   
   if( timeDiff < TIME_CUT ) {
     locHitVec.push_back( locHit );
     n_good_hits++;
   } else { break; }
 
 }
 
 if( locHitVec.size() ) hitPatterns.push_back( locHitVec );
 
 clonedHitArray.erase( clonedHitArray.begin(), clonedHitArray.begin() + n_good_hits );
 
}


return;

513}




518//==========================================================
519//
520//   sortByTime
521//
522//==========================================================

524void DCCALShower_factory::sortByTime( vector< const DCCALHit* > &hitarray, float hitTime )
525{

int nhits = static_cast<int>( hitarray.size() );

if( nhits < 2 ) return; // nothing to sort

for( int ih = 1; ih < nhits; ih++ ) 
{
 
 float timeDiff     = fabs( hitarray[ih]->t - hitTime );
 float lastTimeDiff = fabs( hitarray[ih-1]->t - hitTime );
 
 if( timeDiff <= lastTimeDiff ) 
 {
 
   const DCCALHit *Hit = hitarray[ih];
   
   for( int ii = ih-1; ii >= -1; ii-- ) { 
     
     if( ii >= 0 ) {
     
       const DCCALHit *locHit = hitarray[ii];
       float locTimeDiff = fabs( locHit->t - hitTime );
     
       if( timeDiff < locTimeDiff ) {
         hitarray[ii+1] = locHit;
       } else {
         hitarray[ii+1] = Hit;
  break;
       }
     } else {
       hitarray[0] = Hit;
     }
   }
   
 } // end if statement

} // end loop over hits

return;
565}




570//==========================================================
571//
572//   cleanHitPattern
573//
574//==========================================================

576void DCCALShower_factory::cleanHitPattern( vector< const DCCALHit* > hitarray, 
vector< const DCCALHit* > &hitarrayClean ) 
578{

for( vector< const DCCALHit* >::const_iterator iHit = hitarray.begin(); 
iHit != hitarray.end(); ++iHit ) {
 
 int id12 = ((*iHit)->row)*12 + (*iHit)->column;
 int findVal = -1;
 for( vector<const DCCALHit*>::size_type ii = 0; ii != hitarrayClean.size(); ++ii ) {
   int id = 12*(hitarrayClean[ii]->row) + hitarrayClean[ii]->column;
   if( id == id12 ) { findVal = (int)ii; break; }
 }
 if( findVal >= 0 ) {
   if( (*iHit)->E > hitarrayClean[findVal]->E ) {
     hitarrayClean.erase( hitarrayClean.begin()+findVal );
     hitarrayClean.push_back( (*iHit) );
   }
 } else {
   hitarrayClean.push_back( (*iHit) );
 }
 
}

return;

602}




607//==========================================================
608//
609//   processShowers
610//
611//==========================================================

613void DCCALShower_factory::processShowers( vector< gamma_t > gammas, DCCALGeometry ccalGeom, 
vector< const DCCALHit* > locHitPattern, int eventnumber, 
vector< ccalcluster_t > &ccalClusters, vector< cluster_t > &clusterStorage )
616{


//-------------   Do some post-island processing   -------------//


int n_clusters = 0;
int n_hits = static_cast<int>( locHitPattern.size() );

int init_clusters = static_cast<int>( gammas.size() );
for( int k = 0; k < init_clusters; k++ )  {
 
 ccalcluster_t locCluster;    // stores cluster parameters
 cluster_t locClusterStorage; // stores hit information of cluster cells
 
 int type     = gammas[k].type;
 int dime     = gammas[k].dime;
 int id       = gammas[k].id;
 double chi2  = gammas[k].chi2;
 double e     = gammas[k].energy;
 double x     = gammas[k].x;
 double y     = gammas[k].y;
 double xc    = gammas[k].xc;
 double yc    = gammas[k].yc;
 
 
 // check that shower is not just a single module and that energy is reasonable:
 
 if( dime < MIN_CLUSTER_BLOCK_COUNT ) { continue; } 
 if( e < MIN_CLUSTER_ENERGY || e > MAX_CLUSTER_ENERGY ) { continue; }
 
 n_clusters++;
 
 
 //------------   Find cell with max energy   ------------//
 
 double ecellmax = -1; int idmax = -1;
 double e1 = 0.0;
 for( int j = 0; j < (dime>MAX_CC60 ? MAX_CC60 : dime); j++ ) {
 
   
   double ecell = 1.e-4*static_cast<double>( gammas[k].icl_en[j] );
   
   int ccal_id = gammas[k].icl_in[j];
   int kx  = 12 - (ccal_id/100), ky = 12 - (ccal_id%100);
   
   ccal_id = ky*12 + kx;
   
   e1 += ecell;
   if( ecell > ecellmax ) {
     ecellmax = ecell;
     idmax = ccal_id;
   }
 
 }
 
 double xmax = ccalGeom.positionOnFace(idmax).X();
 double ymax = ccalGeom.positionOnFace(idmax).Y();
 
 
 //-----------   Loop over constituent hits   ------------//
 
 double sW   = 0.0;
 double xpos = 0.0;
 double ypos = 0.0;
 double W;
 
 for( int j = 0; j < (dime>MAX_CC60 ? MAX_CC60 : dime); j++ ) {
 
   int ccal_id = gammas[k].icl_in[j];
   int kx  = 12 - (ccal_id/100), ky = 12 - (ccal_id%100);
   ccal_id = ky*12 + kx;
 
   double ecell  = 1.e-4*static_cast<double>( gammas[k].icl_en[j] );
   double xcell  = ccalGeom.positionOnFace( ccal_id ).X();
   double ycell  = ccalGeom.positionOnFace( ccal_id ).Y();
   
   if(id%10 == 1 || id%10 == 2) {
     xcell += xc;
     ycell += yc;
   }
   
   double hittime = 0.;
   for( int ihit = 0; ihit < n_hits; ihit++ ) {
     int trialid = 12*( locHitPattern[ihit]->row) + locHitPattern[ihit]->column;
     if( trialid == ccal_id ) {
       hittime = locHitPattern[ihit]->t;
break;
     }
   }
   
   locClusterStorage.id[j] = ccal_id;
   locClusterStorage.E[j]  = ecell;
   locClusterStorage.x[j]  = xcell;
   locClusterStorage.y[j]  = ycell;
   locClusterStorage.t[j]  = hittime;
   
   
   // The shower position is calculated using logarithmic weighting:
   
   if( ecell > 0.009 && fabs(xcell-xmax) < 6. && fabs(ycell-ymax) < 6.) {
     W = LOG_POS_CONST + log(ecell/e);
     if( W > 0 ) {
sW   += W;
xpos += xcell*W;
ypos += ycell*W;
     }
   }
   
 }
 
 for( int j = dime; j < MAX_CC60; j++ )  // zero the rest
   locClusterStorage.id[j] = -1;
 
 double weightedTime = getEnergyWeightedTime( locClusterStorage, dime );
 double showerTime   = getCorrectedTime( weightedTime, e );
 
 
 
 //-------  Get position at surface of Calorimeter -------//
 
 DVector3 vertex(0.0, 0.0, m_zTarget); // for now, use center of target as vertex
 
 double x1, y1;
 double zV = vertex.Z();
 double z0 = m_CCALfront - zV;
 if(sW) {
   double dz = getShowerDepth( e );
   double zk = 1. / (1. + dz/z0);
   x1 = zk*xpos/sW;
   y1 = zk*ypos/sW;
 } else {
   printf("WRN bad cluster log. coord at event %i: center id = %i, energy = %f\n",
   	eventnumber, idmax, e);
   x1 = 0.0;
   y1 = 0.0;
 }
 
 
 
 
 //--------  Calculate nonlinear-corrected energy --------//
 
 
 if( idmax < 0 ) {
   printf("WRN negative idmax recorded at event %i; energy = %f\n", eventnumber, e);
 }
 
 double ecorr = e;
 if( DO_NONLINEAR_CORRECTION ) ecorr = getCorrectedEnergy( e, idmax );
 
 if( SHOWER_DEBUG ) {
   cout << "\n\nShower energy before correction: " << e << " GeV" << endl;
   cout << "Shower energy after  correction: " << ecorr << " GeV\n\n" << endl;
 }
 
 
 
 
 //--------  Get energy resolution (needs updating) --------//
 
 double se = sqrt( 0.9*0.9*ecorr*ecorr + 2.5*2.5*ecorr + 1.0 ); 
 	// from HYCAL reconstruction, need to tune
    	  se /= 100.;
 
    	  if( (type%10)==1 )
          se *= 1.5;
    	  else if( (type%10)==2 )
          se *= 1.25;
 
 
 
 //-----------------  Fill cluster bank  -----------------//
 
 locCluster.type    = type;
 locCluster.nhits   = dime;
 locCluster.id      = id;
 locCluster.idmax   = idmax;
 
 locCluster.E       = ecorr;
 locCluster.Esum    = e1;
 locCluster.x       = x;
 locCluster.y       = y;
 locCluster.x1      = x1;
 locCluster.y1      = y1;
 locCluster.chi2    = chi2;
 locCluster.time    = showerTime;
 locCluster.emax    = ecellmax;
 locCluster.sigma_E = se;
 
 clusterStorage.push_back( locClusterStorage );
 ccalClusters.push_back( locCluster );

}


return;
813}




818//==========================================================
819//
820//   getEnergyWeightedTime
821//
822//==========================================================

824double DCCALShower_factory::getEnergyWeightedTime( cluster_t clusterStorage, int nHits )
825{

double weightedtime = 0.;
double totEn = 0;
for( int j = 0; j < (nHits > MAX_CC60 ? MAX_CC60 : nHits); j++ ) {
 weightedtime += clusterStorage.t[j]*clusterStorage.E[j];
 totEn += clusterStorage.E[j];
}
weightedtime /= totEn;

return weightedtime;

837}




842//==========================================================
843//
844//   getCorrectedTime
845//
846//==========================================================

848double DCCALShower_factory::getCorrectedTime( double time, double energy ) 
849{
// timewalk correction:

int iPar;
if( energy < 1.0 ) iPar = 0;
else iPar = 1;

double dt = timewalk_p0[iPar]*exp( timewalk_p1[iPar] + timewalk_p2[iPar]*energy ) + timewalk_p3[iPar];
double t_cor = time - dt;

return t_cor;

861}




866//==========================================================
867//
868//   getShowerDepth
869//
870//==========================================================

872double DCCALShower_factory::getShowerDepth( double energy ) 
873{

double z0 = CCAL_RADIATION_LENGTH, e0 = CCAL_CRITICAL_ENERGY;
double depth = (energy > 0.) ? z0*log(1.+energy/e0) : 0.;
return depth;

879}




884//==========================================================
885//
886//   getCorrectedEnergy
887//
888//==========================================================

890double DCCALShower_factory::getCorrectedEnergy( double energy, int id ) 
891{

if( id < 0 ) return energy;

if( Nonlin_p1[id] == 0. && Nonlin_p2[id] == 0. && Nonlin_p3[id] == 0.) return energy;
if( Nonlin_p0[id] == 0. ) return energy;
	//if( energy < 0.5 || energy > 12. ) return energy;


double emin = 0., emax = 12.;
	double e0 = (emin+emax)/2.;

	double de1 = energy - emin*nonlin_func( emin, id );
	double de2 = energy - emax*nonlin_func( emax, id );
	double de  = energy - e0*nonlin_func( e0, id );

	while( fabs(emin-emax) > 1.e-5 ) {
  	  if( de1*de > 0. && de2*de < 0.) {
    	    emin = e0;
    	    de1 = energy - emin*nonlin_func( emin, id );
  	  } else {
    	    emax = e0;
    	    de2 = energy - emax*nonlin_func( emax, id );
  	  }
  	  e0 = (emin+emax)/2.;
  	  de  = energy - e0*nonlin_func( e0, id );
	}

	return e0;

921}




926//==========================================================
927//
928//   nonlin_func
929//
930//==========================================================

932double DCCALShower_factory::nonlin_func( double e, int id ) 
933{

	return pow( (e/Nonlin_p0[id]), Nonlin_p1[id] + Nonlin_p2[id]*e + Nonlin_p3[id]*e*e );

937}

























963//==========================================================
964//
965//   main_island()
966//
967//==========================================================

969void DCCALShower_factory::main_island( vector<int> &ia, vector<int> &id, vector<gamma_t> &gammas ) 
970{

//----------------------------------------------
// Initialize some useful variables:

int nhits;          // number of hits in detector
int ncl;	    // number of clusters
vector<int> lencl;  // length of a cluster


nhits = static_cast<int>(ia.size());
if( nhits == 0 ) return;

/*
The vector 'ia' will hold the addresses of modules that were hit.
They're defined as 100*(i+1)+(j+1) where i is column and j is row of the hit module.
Row 0, column 0 is bottom right corner of CCAL (looking upstream).

The 'id' vector holds the energies of the hit modules (in units of 0.1 MeVs).
*/


//------------------ Search for clusters: ------------------//


ncl = clus_hyc( nhits, ia, id, lencl );

/*
At this point, the vectors ia and id are sorted such that they are grouped
together by simply-connected clusters. For example, the first lencl[0] elements of ia give 
the addresses of the elements in the first cluster. The next lencl[1] elements
after that give the addresses of the elements in the second cluster.

The clusters here are just created by grouping everything that is physically 
next to each other into one cluster - all simply-connected cells are joined in a cluster.
*/


//-------------------- Process Cluster: --------------------//

int nadcgam = 0; // number of found clusters

if( ncl <= 0 ) return;

int ipncl = 0;
for( int icl = 0; icl < ncl; icl++ ) {
 
 int ecl = 0;
 for( int ii = 0; ii < lencl[icl]; ii++ ) 
   ecl += id[ipncl+ii];
 if( ecl > MIN_ENERGY0. ) {
   
   vector< int > icl_a; // addresses of current cluster
   vector< int > icl_d; // energies of current cluster
   
   icl_a.insert( icl_a.begin(), ia.begin()+ipncl, ia.begin()+ipncl+lencl[icl] );
   icl_d.insert( icl_d.begin(), id.begin()+ipncl, id.begin()+ipncl+lencl[icl] );
   
   
   if( SHOWER_DEBUG ) {
   
     cout << "\n\n======================" << endl;
     cout << "Processing Cluster " << icl << ":" << endl;
     for( unsigned int ih = 0; ih < icl_a.size(); ih++ ) {
       cout << icl_a[ih] << " " << icl_d[ih] << endl;
     }
   
   }
   
   int before = nadcgam;
   gams_hyc( lencl[icl], icl_a, icl_d, nadcgam, gammas );
   int after  = nadcgam;
   
   
   if( SHOWER_DEBUG ) {
     cout << "Reconstructed " << after-before << " gammas. \n\n" << endl;
   }
   
   
   if( nadcgam > MADCGAM50 ) {
     nadcgam = MADCGAM50;
     break;
   }
 }
 ipncl = ipncl + lencl[icl];
}


//-------------------- Prepare gammas for final processing --------------------//


// convert position to units of cm:

for( int ig = 0; ig < nadcgam; ig++ ) {
 gammas[ig].energy /= 10000.;
 gammas[ig].x  = (static_cast<double>(MCOL12+1)-2.*gammas[ig].x)*xsize2.09/2.;
 gammas[ig].y  = (static_cast<double>(MROW12+1)-2.*gammas[ig].y)*ysize2.09/2.;
 gammas[ig].xc = -gammas[ig].xc*xsize2.09;
 gammas[ig].yc = -gammas[ig].yc*ysize2.09;
}


if( nadcgam < 1 ) return;	


// sort gammas in order of decreasing energy:

for( int ig = 1; ig < nadcgam; ig++ ) 
{ 
 if( gammas[ig].energy > gammas[ig-1].energy ) { 
   gamma_t ref_gam = gammas[ig]; 
   
   for( int ii = ig-1; ii >= -1; ii-- ) {	      
     if( ii >= 0 ) {
       if( ref_gam.energy > gammas[ii].energy ) {
         gammas[ii+1] = gammas[ii];
       } else {
         gammas[ii+1] = ref_gam;
  break;
       }
     } else {
       gammas[0] = ref_gam;
     }
   }
 }
}
 


return;
1100}





1106//==========================================================
1107//
1108//   clus_hyc
1109//
1110//==========================================================

1112int DCCALShower_factory::clus_hyc( int nw, vector<int> &ia, vector<int> &id, vector<int> &lencl ) 
1113{

//---------------- Local Declarations ---------------//

int maxcl = 200;    // maximumum number of clusters allowed
int ncl;            // number of clusters
int next, iak;      // 
int ib, ie;         // 
int ias, iaf;       // 
int last, lastcl;   // 
int leng;           // 
                   // 
int loclencl[200];  // stores the lengths of clusters locally


//---------------- Event Analysis Code --------------//

ncl = 0;
if( nw < 1 ) return ncl;
if( nw < 2 ) { // if only one hit
 ncl = 1;
 lencl.push_back(1);
 return ncl;
}

order_hyc( nw, ia, id ); // sort the addresses (ia) in increasing order

ncl  = 0;
next = 0;

for( int k = 1; k < (nw+1); k++ ) {
 
 if( k < nw ) iak = ia[k];
 if( (iak-ia[k-1] <= 1) && (k < nw) ) continue;
 
 ib   = next;
 ie   = k-1;
 next = k;
 
 if( ncl >= maxcl ) return ncl; 
 ncl++;
 
 loclencl[ncl-1] = next-ib;
 if(ncl == 1) continue;
 
 ias    = ia[ib];
 iaf    = ia[ie];
 last   = ib-1;
 lastcl = ncl-1;
 
 for( int icl = lastcl; icl > 0; icl-- ) {
   
   leng = loclencl[icl-1];
   if( (ias-ia[last]) > 100 ) break; // no subclusters to be glued
   
   for( int ii = last; ii >= last-leng+1; ii-- ) {
     if( ( ias-ia[ii] )  > 100 ) break;
     if( ( iaf-ia[ii] ) >= 100 ) {
     
       if( (icl < (ncl-1)) && (leng <= 10800) ) {  

  vector< int > iawork;
  ucopy1( ia, iawork, last-leng+1, leng );
  ucopy2( ia, last+1, last+1-leng, ib-last-1 );
  ucopy3( iawork, ia, ib-leng, leng );
  
  vector< int > idwork;
  ucopy1( id, idwork, last-leng+1, leng );
  ucopy2( id, last+1, last+1-leng, ib-last-1 );
  ucopy3( idwork, id, ib-leng, leng );
  
  for( int jj = icl; jj < ncl-1; jj++ ) {
                  loclencl[jj-1] = loclencl[jj];
                }

}

ib = ib-leng;
              loclencl[ncl-2] = loclencl[ncl-1]+leng;
              ncl = ncl-1;
              break;
     }
   }
   
   last = last-leng;
   
 } // end loop over previous subclusters
} // end loop over all hits

for( int icl = 0; icl < ncl; icl++ ) {
 lencl.push_back( loclencl[icl] );
}


return ncl;
1208}





1214//==========================================================
1215//
1216//   order_hyc
1217//
1218//==========================================================

1220void DCCALShower_factory::order_hyc( int nw, vector<int> &ia, vector<int> &id ) 
1221{
// sort ia and id in order of increasing address

if( nw < 2 ) return; // only one hit

for( int k = 1; k < nw; k++ ) { // loop over hits
 
 if( ia[k] <= ia[k-1] ) { // check if address is less than previous entry
   int iak = ia[k];
   int idk = id[k];
   
   for( int ii = k-1; ii >= -1; ii = ii-1 ) { // loop over the previous entries
     
     if( ii >= 0 ) {
       if( iak < ia[ii] ) {
         ia[ii+1] = ia[ii];
  id[ii+1] = id[ii];
       } else {
         ia[ii+1] = iak;
  id[ii+1] = idk;
  break;
       }
     } else {
       ia[0] = iak;
       id[0] = idk;
     }
     
   } // end loop over previous entries
 } // endif
} // end loop over hits


return;
1254}





1260//==========================================================
1261//
1262//   gams_hyc
1263//
1264//==========================================================

1266void DCCALShower_factory::gams_hyc( int nadc, vector<int> &ia, vector<int> &id, 
int &nadcgam, vector<gamma_t> &gammas )
1268{

//-------------Local Declarations------------//

int ngam0;                // number of gammas found before
int niter;                // max number of iterations (6)
int npk;                  // number of peaks in the cluster
int ipnpk[10];            // counter number with max energy in a peak
int igmpk[2][10];         // 
int minpk;                // min energy of a counter in a cluster
int idelta;               // 
int itype;                // type of peak
int leng;                 // 
int ixypk, ixpk, iypk;    // 
int ic, idc, in;          // 
int ixy, ixymax, ixymin;  // 
int ix, iy, iyc;          // 
int iwk;                  // 
int iia;                  // 
int ide;                  // 
int idecorr;              // 

double fw[13];	
double chisq;             // current value of chi2
double chisq1;            // value of chi2 for preliminary gammas separation
double chisq2;            // value of chi2 for final gammas separation
double ratio;             // 
double eg;                // 
double epk[10];           // 
double xpk[10];           // 
double ypk[10];           // 
double a, dx, dy;         // 
double e1, x1, y1;        // 
double e2, x2, y2;        // 
double fe, fia;           // 

int idsum;	

idelta = 0;
chisq1 = 90.;
chisq2 = 50.;
niter  = 6;

ngam0 = nadcgam;



vector<int> iwrk[13];     // working array for resolved peaks
vector<int> idp[13];      // energy of each cell of the island belonging to each peak
vector<double> fwrk[13];  // working array for resolved peaks

// allocate memory for each vector:

for( int ii=0; ii<13; ii++ ) {
 iwrk[ii].reserve(nadc);
 fwrk[ii].reserve(nadc);
 idp[ii].reserve(nadc);
 for( int ih=0; ih<nadc; ih++ ) {
   iwrk[ii].push_back(0);
   fwrk[ii].push_back(0.);
   idp[ii].push_back(0);
 }
}




//------------------------------------------
// peaks search:

order_hyc( nadc, ia, id );

idsum = 0;
for( int ic = 0; ic < nadc; ic++ )
 idsum += id[ic];
 
if( nadc < 3 )
 minpk = 1;
else {
 int trial = 7.*log(1. + 0.0001*static_cast<double>(idsum)) + 0.5;
 if( trial > 1 ) minpk = trial;
 else minpk = 1;
}
minpk *= 100;


npk = 0;

for( ic = 0; ic < nadc; ic++ ) {
 idc = id[ic];
 if( idc < minpk ) continue;
 ixy = ia[ic];
 ixymax = ixy + 100 + 1;
 ixymin = ixy - 100 - 1;
 iyc = ixy - (ixy/100)*100;
 
 int peakVal = 1;
 
 in = ic+1;
 while( in < nadc && ia[in] <= ixymax ) {
   iy = ia[in] - (ia[in]/100)*100;
   if( abs(iy-iyc) <= 1 && id[in] >= idc ) peakVal = 0;
   in++;
 }
 
 in = ic-1;
 while( in >= 0 && ia[in] >= ixymin ) {
   iy = ia[in] - (ia[in]/100)*100;
   if( abs(iy-iyc) <= 1 && id[in] > idc ) peakVal = 0;
   in -= 1;
 }
 
 if( !peakVal ) continue;
 
 npk += 1;
 ipnpk[npk-1] = ic;
 if( npk == 10 || npk >= 10000/nadc-3 ) break;
   	
}

if( npk <= 0 ) return;

if( SHOWER_DEBUG ) cout << "Found " << npk << " peaks. Now processing..." << endl;

//------------------------------------------
// gammas search for one peak:

if( npk == 1 ) {

 if( nadcgam >= MADCGAM50-1 ) return;
 nadcgam = nadcgam+1;
 chisq = chisq2;
 
 ic = ipnpk[0];
 ix = ia[ic]/100;
 iy = ia[ic] - ix*100;

 itype = peak_type( ix, iy );
 
 e2 = 0.;
 gamma_hyc( nadc, ia, id, chisq,
 		e1, x1, y1, 
	e2, x2, y2 );
	
 gamma_t gam1;
 gamma_t gam2;
 
 gam1.type   = itype;
 gam1.dime   = nadc;
 gam1.id     = 0;
   
 gam1.chi2   = chisq;
 gam1.energy = e1;
 gam1.x      = x1;
 gam1.y      = y1;
 gam1.xc     = 0.;
 gam1.yc     = 0.;
 
 if( e2 > 0. && nadcgam <= MADCGAM50-1 ) {
   nadcgam = nadcgam+1;
   
   gam2.type   = itype+10;
   gam2.dime   = nadc;
   gam2.id     = 2;
   
   gam2.chi2   = chisq;
   gam2.energy = e2;
   gam2.x      = x2;
   gam2.y      = y2;
   gam2.xc     = 0.5*(x2-x1);
   gam2.yc     = 0.5*(y2-y1);
   
   gam1.type   = itype+10;
   gam1.id     = 1;
   gam1.xc     = 0.5*(x1-x2);
   gam1.yc     = 0.5*(y1-y2);
   
   for( int jj = 0; jj < nadc; jj++ ) {
     if( jj < MAX_CC60 ) {
       gam1.icl_in[jj] = ia[jj];
gam2.icl_in[jj] = ia[jj];
gam1.icl_en[jj] = static_cast<int>(static_cast<double>(id[jj])*e1/(e1+e2) + 0.5);
gam2.icl_en[jj] = static_cast<int>(static_cast<double>(id[jj])*e2/(e1+e2) + 0.5);
     }
   }
   
   gammas.push_back( gam1 );
   gammas.push_back( gam2 );
   
 } else {
   for( int jj = 0; jj < nadc; jj++ ) {
     if( jj < MAX_CC60 ) {
       gam1.icl_in[jj] = ia[jj];
gam1.icl_en[jj] = id[jj];
     }
   }
   
   gammas.push_back( gam1 );
   
 }


} else { // cluster with more than one peak

//------------------------------------------
/*
First step - 1 gamma in each peak.
Do a preliminary estimation of (E,x,y) of each peak, and split each peak into two hits 
only if it is badly needed (chi2 improvement is too high).
If this split occurs, it is only for better (E,x,y) estimation, as it will be 
rejoined and reanalyzed in the second step.
*/
//------------------------------------------


 if( nadcgam >= MADCGAM50-1 ) return;
 
 ratio = 1.;
 for( int iter = 0; iter < niter; iter++ ) {
   for( int ii = 0; ii < nadc; ii++ ) {
     iwrk[0].push_back(0);
     fwrk[0].push_back(0.);
   }
   
   for( int ipk = 0; ipk < npk; ipk++ ) {
   
     ic = ipnpk[ipk];
     if( iter > 0 ) ratio = fwrk[ipk+1][ic]/fwrk[npk+1][ic];
     eg = ratio*static_cast<double>(id[ic]);
     ixypk = ia[ic];
     ixpk  = ixypk/100;
     iypk  = ixypk - ixpk*100;
     epk[ipk] = eg;
     xpk[ipk] = eg*static_cast<double>(ixpk);
     ypk[ipk] = eg*static_cast<double>(iypk);
     
     if( ic < nadc-1 ) {
       for( int ii = ic+1; ii < nadc; ii++ ) {
  ixy = ia[ii];
  ix  = ixy/100;
  iy  = ixy - ix*100;
  if( ixy-ixypk > 100+1 ) break;
  if( abs(iy-iypk) <= 1 ) {
    if( iter != 0 ) ratio = fwrk[ipk+1][ii]/fwrk[npk+1][ii];
    eg = ratio*static_cast<double>(id[ii]);
    epk[ipk] = epk[ipk] + eg;
    xpk[ipk] = xpk[ipk] + eg*static_cast<double>(ix);
    ypk[ipk] = ypk[ipk] + eg*static_cast<double>(iy);
  }
}
     }
     
     if( ic > 0 ) {
       for( int ii = ic-1; ii >= 0; ii--) {
  ixy = ia[ii];
  ix  = ixy/100;
  iy  = ixy - ix*100;
  if( ixypk-ixy > 100+1 ) break;
  if( abs(iy-iypk) <= 1 ) {
    if( iter != 0 ) ratio = fwrk[ipk+1][ii]/fwrk[npk+1][ii];
    eg = ratio*static_cast<double>(id[ii]);
    epk[ipk] = epk[ipk] + eg;
    xpk[ipk] = xpk[ipk] + eg*static_cast<double>(ix);
    ypk[ipk] = ypk[ipk] + eg*static_cast<double>(iy);
  }
}
     }
     
     if( epk[ipk] > 0. ) {
       xpk[ipk] = xpk[ipk]/epk[ipk];
ypk[ipk] = ypk[ipk]/epk[ipk];
     }
     
     for( int ii = 0; ii < nadc; ii++ ) {
       ixy = ia[ii];
ix  = ixy/100;
iy  = ixy - ix*100;
dx  = fabs( static_cast<double>(ix) - xpk[ipk] );
dy  = fabs( static_cast<double>(iy) - ypk[ipk] );

a = epk[ipk]*cell_hyc( dx, dy );

fwrk[ipk+1][ii] = a;
fwrk[0][ii]     = fwrk[0][ii] + fwrk[ipk+1][ii];
iwrk[ipk+1][ii] = static_cast<int>(a + 0.5);
iwrk[0][ii]     = iwrk[0][ii] + iwrk[ipk+1][ii];
     }
     
     
   } // end loop over peaks
   
   for( int ii = 0; ii < nadc; ii++ ) {
     iwk = iwrk[0][ii];
     if( iwk < 1 ) iwk = 1;
     iwrk[npk+1][ii] = iwk;
     
     if( fwrk[0][ii] > 1.e-2 ) 
       fwrk[npk+1][ii] = fwrk[0][ii];
     else 
       fwrk[npk+1][ii] = 1.e-2;

   }
 } // end of iterations to separate peaks in a cluster
 
 
 
 if( SHOWER_DEBUG ) {
 
 
   cout << "\n\n\nAfter 6 iterations: " << endl;
   for( int ipk = 0; ipk < npk; ipk++ ) {
     cout << "peak " << ipk+1 << ": " << endl;
     for( int jj = 0; jj < nadc; jj++ ) {
       cout << ia[jj] <<"; "<< id[jj] <<"; "<< fwrk[ipk+1][jj] << endl;
     }
   }
 
 }
 
 
 
 
 for( int ipk = 0; ipk < npk; ipk++ ) {
   
   vector<int> iwrk_a;
   vector<int> iwrk_d;
   
   leng = 0;
   for( int ii = 0; ii < nadc; ii++ ) {
     if( fwrk[0][ii] > 1.e-2 ) {
       ixy = ia[ii];
fe  = static_cast<double>(id[ii])*fwrk[ipk+1][ii]/fwrk[0][ii];

if( fe > idelta ) {
  leng = leng+1;
  iwrk_a.push_back( ixy );
  iwrk_d.push_back( static_cast<int>(fe+0.5) );
}
     
     }
   }
   
   if( nadcgam >= MADCGAM50-1 ) return;
   
   igmpk[1][ipk] = 0;
   if( leng == 0 ) continue;
   nadcgam = nadcgam + 1;
   chisq = chisq1;
   
   ic = ipnpk[ipk];
   ix = ia[ic]/100;
   iy = ia[ic] - ix*100;
   
   itype = peak_type( ix, iy );
   
   e2 = 0.;
   gamma_hyc( leng, iwrk_a, iwrk_d, chisq,
   		e1, x1, y1, e2, x2, y2 );
	
   gamma_t gam1;
   gamma_t gam2;
 
   gam1.chi2   = chisq;
   gam1.type   = itype+40;
   gam1.energy = e1;
   gam1.x      = x1;
   gam1.y      = y1;
   gam1.dime   = leng;
   
   gam1.id     = 90;
   
   igmpk[0][ipk] = nadcgam;
   igmpk[1][ipk] = nadcgam;
 
   if( e2 > 0. && nadcgam <= MADCGAM50-1 ) {
     nadcgam = nadcgam+1;
   
     gam2.chi2   = chisq;
     gam2.type   = itype+50;
     gam2.energy = e2;
     gam2.x      = x2;
     gam2.y      = y2;
     gam2.xc     = 0.5*(x2-x1);
     gam2.yc     = 0.5*(y2-y1);
     gam1.xc     = 0.5*(x1-x2);
     gam1.yc     = 0.5*(y1-y2);
     
     gam2.id     = 92;
     gam1.id     = 91;
     
     gam2.dime   = leng;
     igmpk[1][ipk] = nadcgam;
     
     gammas.push_back( gam1 );
     gammas.push_back( gam2 );
     
   } else { gammas.push_back( gam1 ); }
   
 } // end loop over peaks
 
 
 
 /*
 This is the second step: ( 1 or 2 gamma in each peak )
 (e,x,y) of hits were preliminarily estimated in the first step.	  
 */
 
 for( int ii = 0; ii < nadc; ii++ ) {
   iwrk[0][ii] = 0;
   fwrk[0][ii] = 0.;
   idp[0][ii]  = 0;
 }
 
 for( int ipk = 0; ipk < npk; ipk++ ) {
   for( int ii = 0; ii < nadc; ii++ ) {
     iwrk[ipk+1][ii] = 0;
     fwrk[ipk+1][ii] = 0.;
     idp[ipk+1][ii]  = 0;
     
     if( igmpk[1][ipk] == 0 ) continue;
     
     for( int ig = igmpk[0][ipk]; ig <= igmpk[1][ipk]; ig++ ) {
       ixy = ia[ii];
ix  = ixy/100;
iy  = ixy-(ix*100);
dx  = static_cast<double>(ix) - gammas[ig-1].x;
dy  = static_cast<double>(iy) - gammas[ig-1].y;

fia = gammas[ig-1].energy*cell_hyc( dx, dy );
iia = static_cast<int>(fia+0.5);

// part of gamma 'ig' energy belonging to cell 'i' from peak 'ipk':
iwrk[ipk+1][ii] += iia; 
fwrk[ipk+1][ii] += fia;
idp[ipk+1][ii]  += iia;

iwrk[0][ii] += iia;
fwrk[0][ii] += fia;

     }
     
   } // end loop over hits
 } // end loop over peaks
 
 
 
 
 // Recover working array and renormalize total sum to the original cell energy:
 
 for( int ii = 0; ii < nadc; ii++ ) {
   
   idp[0][ii] = 0;
   for( int ipk = 0; ipk < npk; ipk++ ) {
     idp[0][ii] += idp[ipk+1][ii];
   }
   
   ide = id[ii] - idp[0][ii];
   if( ide == 0 ) continue;
   if( fwrk[0][ii] == 0. ) continue;
   
   for( int ipk = 0; ipk < npk; ipk++ ) {
     fw[ipk+1] = fwrk[ipk+1][ii]/fwrk[0][ii];
   }
   
   idecorr = 0;
   for( int ipk = 0; ipk < npk; ipk++ ) {
     fia = ide*fw[ipk+1];
     if( (fwrk[ipk+1][ii] + fia) > 0. ) {
       fwrk[ipk+1][ii] += fia;
fwrk[0][ii] += fia;
     }
     iia = static_cast<int>(fia+0.5);
     if( (iwrk[ipk+1][ii] + iia) > 0 ) {
       iwrk[ipk+1][ii] += iia;
iwrk[0][ii] += iia;
idecorr += iia;
     } else if( (iwrk[ipk+1][ii] + iia) < 0 ) {
       //cout << "WARNING NEGATIVE CORR: ia = " << ia[ii] << "; id = " << id[ii] << endl;
     }
   } // end loop over peaks
   
 } // end loop over hits
 
 
 
 // erase the gammas found in the previous step:
 
 nadcgam = ngam0;
 gammas.erase( gammas.begin()+nadcgam, gammas.end() );
 
 
 
 
 
 
 for( int ipk = 0; ipk < npk; ipk++ ) {
   leng = 0;
   
   vector<int> iwrk_a;
   vector<int> iwrk_d;
   
   for( int ii = 0; ii < nadc; ii++ ){
     if( iwrk[0][ii] > 0 ) {
       fe = id[ii]*fwrk[ipk+1][ii]/fwrk[0][ii];
if( fe > idelta ) {
  leng++;
  iwrk_a.push_back( ia[ii] );
  iwrk_d.push_back( static_cast<int>(fe+0.5) );
}
     }
   }
 
   if( nadcgam >= MADCGAM50-1 ) return;
   
   
   if( leng == 0 ) continue; 
   
   nadcgam++;
   
   chisq = chisq2;
   
   ic = ipnpk[ipk];
   ix = ia[ic]/100;
   iy = ia[ic] - ix*100;
   
   itype = peak_type( ix, iy );
   
   e2 = 0.;
   gamma_hyc( leng, iwrk_a, iwrk_d, chisq,
   		e1, x1, y1, e2, x2, y2 );
	
   gamma_t gam1;
   gamma_t gam2;
   
   gam1.type   = itype+20;
   gam1.dime   = leng;
   gam1.id     = 10;
   
   gam1.chi2   = chisq;
   gam1.energy = e1;
   gam1.x      = x1;
   gam1.y      = y1;
   gam1.xc     = 0.;
   gam1.yc     = 0.;
   
   if( e2 > 0. && nadcgam <= MADCGAM50-1 ) {
     nadcgam = nadcgam+1;
     
     gam2.type   = itype+30;
     gam2.dime   = leng;
     gam2.id     = 12;
   
     gam2.chi2   = chisq;
     gam2.energy = e2;
     gam2.x      = x2;
     gam2.y      = y2;
     gam2.xc     = 0.5*(x2-x1);
     gam2.yc     = 0.5*(y2-y1);
     
     gam1.type   = itype+30;
     gam1.id     = 11;
     gam1.xc     = 0.5*(x1-x2);
     gam1.yc     = 0.5*(x1-x2);
     
     for( int jj = 0; jj < leng; jj++ ) {
       if( jj < MAX_CC60 ) {
  gam1.icl_in[jj] = iwrk_a[jj];
  gam2.icl_in[jj] = iwrk_a[jj];
  gam1.icl_en[jj] = static_cast<int>( 
  	static_cast<double>(iwrk_d[jj])*e1/(e1+e2) + 0.5 );
  gam2.icl_en[jj] = static_cast<int>( 
  	static_cast<double>(iwrk_d[jj])*e2/(e1+e2) + 0.5 );
}
     }
     
     gammas.push_back( gam1 );
     gammas.push_back( gam2 );
   
   } else {
     for( int jj = 0; jj < leng; jj++ ) {
       if( jj < MAX_CC60 ) {
         gam1.icl_in[jj] = iwrk_a[jj];
  gam1.icl_en[jj] = iwrk_d[jj];
       }
     }
     
     gammas.push_back( gam1 );
     
   }
   	    
   
 } // end loop over peaks
 
 
 	
} // end looping over multi-peak cluster





return;
1870}





1876//==========================================================
1877//
1878//   peak_type
1879//
1880//==========================================================

1882int DCCALShower_factory::peak_type( int ix, int iy )
1883{
/*
itype = 2 : if the peak is in the most outer layer
itype = 1 : if the peak is in the most inner layer
itype = 0 : if the peak is anywhere else
*/

int itype = 0;
if( (ix == MCOL12/2-1 || ix == MCOL12/2+2) && (iy >= MROW12/2-1 && iy <= MROW12/2+2) ) itype = 1;
if( (iy == MROW12/2-1 || iy == MROW12/2+2) && (ix >= MCOL12/2-1 && ix <= MCOL12/2+2) ) itype = 1;
if( ix == 1 || ix == MCOL12 || iy == 1 || iy == MROW12 ) itype = 2;

return itype;
1896}





1902//==========================================================
1903//
1904//   gamma_hyc
1905//
1906//==========================================================

1908void DCCALShower_factory::gamma_hyc( int nadc, vector<int> ia, vector<int> id, double &chisq, 
double &e1, double &x1, double &y1, 
double &e2, double &x2, double &y2 )
1911{
//-------------Local Declarations------------//

int dof;

double dxy;                    // initial step for iteration
double stepmin;                // minimum step for iteration
double stepx, stepy;           // current steps
double parx, pary;             // 
double chimem, chisq0, chi0;   //
double chiold;                 // 
double chi00;                  // 
double x0, y0;                 // 
double ee, xx, yy;             // 
double d2, xm2, xm2cut;        // 
double chir, chil, chiu, chid; //

dxy     = 0.05;
stepmin = 0.002;
xm2cut  = 1.7;

int nzero;
vector<int> iaz;

//-------------Event Analysis Code------------//

e2 = 0.;
x2 = 0.;
y2 = 0.;	


fill_zeros( nadc, ia, nzero, iaz );
mom1_pht( nadc, ia, id, nzero, iaz, e1, x1, y1 ); // calculate initial values of (E,x,y)

if( nadc <= 0 ) return;

chimem = chisq;
chisq1_hyc( nadc, ia, id, nzero, iaz, e1, x1, y1, chi0 ); // initial value of chi2


chisq0 = chi0;
dof = nzero + nadc - 2;
if( dof < 1 ) dof = 1;
chisq = chi0/dof;
x0 = x1;
y0 = y1;


// start of iterations

int rounds = 0;
while( 1 ) {

 chisq1_hyc( nadc, ia, id, nzero, iaz, e1, x0+dxy, y0, chir );
 chisq1_hyc( nadc, ia, id, nzero, iaz, e1, x0-dxy, y0, chil );
 chisq1_hyc( nadc, ia, id, nzero, iaz, e1, x0, y0+dxy, chiu );
 chisq1_hyc( nadc, ia, id, nzero, iaz, e1, x0, y0-dxy, chid );
 
 if( chi0 > chir || chi0 > chil ) {
   stepx = dxy;
   if( chir > chil ) stepx = -stepx;
 } else {
   stepx = 0.;
   parx = chir + chil - 2.*chi0;
   if( parx > 0. ) stepx = -dxy*(chir-chil)/(2.*parx);
 }
 if( chi0 > chiu || chi0 > chid ) {
   stepy = dxy;
   if( chiu > chid ) stepy = -stepy;
 } else {
   stepy = 0.;
   pary = chiu + chid - 2.*chi0;
   if( pary > 0. ) stepy = -dxy*(chiu-chid)/(2.*pary);
 }


 // if steps at minimum, end iterations

 if( fabs(stepx) < stepmin && fabs(stepy) < stepmin ) break;
   
 chisq1_hyc( nadc, ia, id, nzero, iaz, e1, x0+stepx, y0+stepy, chi00 );
 
 // if chi2 at minimum, end iterations
 
 if( chi00 >= chi0 ) break;
 
 chi0 = chi00;
 x0 = x0+stepx;
 y0 = y0+stepy;
 
 rounds++;
 if( rounds > 10000 ) { cout << "max rounds" << endl; break; }
   
}

if( chi0 < chisq0 ) { // if chi2 improved, then fix the improved values
 x1 = x0;
 y1 = y0;
 chisq = chi0/dof;
}

// if chi2 is less than maximum allowed for one gamma in a peak, return.
// otherwise, try separating the peak into two gammas:

if( chisq <= chimem ) return; 

chiold = chisq;
tgamma_hyc( nadc, ia, id, nzero, iaz, chisq, ee, xx, yy, e2, x2, y2 );

if( e2 > 0. ) {   // if chi2 improved, decide if the separation
          // has physical meaning by calculating the 
	  // effective mass of the two gammas
	  
 d2 = pow((xx-x2)*xsize2.09, 2.0) + pow((yy-y2)*ysize2.09, 2.0);
 xm2 = ee*e2*d2;
 
 if( xm2 > 0. ) xm2 = sqrt(xm2)/1270.*0.1; // mass in MeV; 1270 = zccal
 
 if( xm2 > xm2cut*xsize2.09 ) { // if the separation has physical meaning
   e1 = ee;		     // fix the parameters of the first gamma
   x1 = xx;
   y1 = yy;
 } else {		     // if no physical meaning e2=0 
   e2 = 0.;		     // (second gamma is killed)
   chisq = chiold;
 }
}

return;
2040}





2046//==========================================================
2047//
2048//   fill_zeros
2049//
2050//==========================================================

2052void DCCALShower_factory::fill_zeros( int nadc, vector<int> ia, int &nneib, vector<int> &iaz )
2053{

int ix, iy, nneibnew;
vector<int> ian;

nneib = 0;
for( int ii = 0; ii < nadc; ii++ ) {
 ix = ia[ii]/100;
 iy = ia[ii] - ix*100;
 
 if( ix > 1 ) { // fill left neib
   nneib = nneib+1;
   ian.push_back(iy + (ix-1)*100);
 
   if( iy > 1 ) { // fill bottom left neib
     nneib = nneib+1;
     ian.push_back(iy-1 + (ix-1)*100);
   }
   if( iy < MROW12 ) { // fill top left neib
     nneib = nneib+1;
     ian.push_back(iy+1 + (ix-1)*100);
   }
 }
 if( ix < MCOL12 ) {
   nneib = nneib+1;
   ian.push_back(iy + (ix+1)*100);
   
   if( iy > 1 ) { // fill bottom right neib
     nneib = nneib+1;
     ian.push_back(iy-1 + (ix+1)*100);
   }
   if( iy < MROW12 ) { // fill top right neib
     nneib = nneib+1;
     ian.push_back(iy+1 + (ix+1)*100);
   }
 }
 
 if( iy > 1 ) { // fill bottom neib
   nneib = nneib+1;
   ian.push_back(iy-1 + ix*100);
 }
 if( iy < MROW12 ) { // fill top neib
   nneib = nneib+1;
   ian.push_back(iy+1 + ix*100);
 }
}


for( int ii = 0; ii < nneib; ii++ ) {
 for( int jj = 0; jj < nadc; jj++ ) {
   if( ia[jj] == ian[ii] ) ian[ii] = -1;
 }
}

for( int ii = 0; ii < nneib; ii++ ) {
 if( ian[ii] == -1 ) continue;
 for( int jj = ii+1; jj < nneib; jj++ ) {
   if( ian[jj] == ian[ii] ) ian[jj] = -1;
 }
}

nneibnew = 0;
for( int ii = 0; ii < nneib; ii++ ) {
 ix = ian[ii]/100;
 iy = ian[ii] - ix*100;
 if( ian[ii] != -1 ) {
   if( stat_ch[iy-1][ix-1] == 0 ) {
     nneibnew = nneibnew+1;
     iaz.push_back( ian[ii] );
   }
 }
}
nneib = nneibnew;


return;
2129}





2135//==========================================================
2136//
2137//   mom1_pht
2138//
2139//==========================================================

2141void DCCALShower_factory::mom1_pht( int nadc, vector<int> ia, vector<int> id, 
int nzero, vector<int> iaz, double &e1, double &x1, double &y1 )
2143{
//-------------Local Declarations------------//

int ix, iy;
double dx, dy;

double a;          // 
double corr;       // correction to energy

//--------------Event Analysis Code-------------//

e1 = 0.;
x1 = 0.;
y1 = 0.;

if( nadc <= 0 ) return;
for( int ii = 0; ii < nadc; ii++ ) {
 a  = static_cast<double>(id[ii]);
 ix = ia[ii]/100;
 iy = ia[ii] - ix*100;
 e1 = e1 + a;
 x1 = x1 + a*static_cast<double>(ix);
 y1 = y1 + a*static_cast<double>(iy);
}
if( e1 <= 0. ) return;
x1 = x1/e1;
y1 = y1/e1;
corr = 0.;
for( int ii = 0; ii < nadc; ii++ ) {
 ix = ia[ii]/100;
 iy = ia[ii] - ix*100;
 dx = static_cast<double>(ix) - x1;
 dy = static_cast<double>(iy) - y1;
 corr = corr + cell_hyc( dx, dy );
}

if( COUNT_ZERO1 ) {
 for( int ii = 0; ii < nzero; ii++ ) {
   ix = iaz[ii]/100;
   iy = iaz[ii] - ix*100;
   dx = static_cast<double>(ix) - x1;
   dy = static_cast<double>(iy) - y1;
   corr = corr + cell_hyc( dx, dy );
 }
}

corr = corr/1.006;
if( SHOWER_DEBUG ) cout << "e, corr = " << e1 << ", " << corr << endl;


if( corr < 0.8 ) {
 if( SHOWER_DEBUG ) {
   cout << "corr = " << corr << ", " << e1 << ", " << x1 << ", " << y1;
   cout << "! Too many around central hole!" << endl;
 }
 corr = 0.8;
} else if( corr > 1.0 ) {
 corr = 1.;
}
e1 = e1/corr;


return;
2206}





2212//==========================================================
2213//
2214//    chisq1_hyc
2215//
2216//==========================================================

2218void DCCALShower_factory::chisq1_hyc( int nadc, vector<int> ia, vector<int> id, 
int nneib, vector<int> iaz, double e1, double x1, double y1, double &chisq )
2220{
//-------------Local Declarations------------//

int ix, iy;
double dx, dy;
double fa, fcell;


//--------------Event Analysis Code-------------//

chisq = 0.;

for( int ii = 0; ii < nadc; ii++ ) {
 fa = static_cast<double>(id[ii]);
 ix = ia[ii]/100;
 iy = ia[ii] - ix*100;
 dx = x1 - static_cast<double>(ix);
 dy = y1 - static_cast<double>(iy);
 if( e1 != 0. ) {
   if( fabs(dx) <= 6.0 && fabs(dy) <= 6.0 ) {
     fcell = cell_hyc( dx, dy );
     chisq = chisq + e1*pow((fcell-(fa/e1)), 2.)/sigma2(dx, dy, fcell, e1);
   }
 } else {
   chisq = chisq + fa*fa/9.;
   cout << "case 0 ch" << endl;
   //if( SHOWER_DEBUG ) cout << "case 0 ch" << endl;
 }
}

for(int ii = 0; ii < nneib; ii++ ) {
 ix = iaz[ii]/100;
 iy = iaz[ii] - ix*100;
 dx = x1 - static_cast<double>(ix);
 dy = y1 - static_cast<double>(iy);
 if( e1 != 0. ) {
   if( fabs(dx) < 6.0 && fabs(dy) < 6.0 ) {
     fcell = cell_hyc( dx, dy );
     chisq = chisq + e1*fcell*fcell/sigma2(dx, dy, fcell, e1);
   }
 } else {
   chisq = chisq + id[ii]*id[ii]/9.;
   //if( SHOWER_DEBUG ) cout << "case 0 ch" << endl;
 }
}


return;
2268}





2274//==========================================================
2275//
2276//   sigma2, d2c, & cell_hyc
2277//
2278//==========================================================

2280double DCCALShower_factory::sigma2( double dx, double dy, double fc, double e ) 
2281{
double sig2;
double alp = 0.816;
double bet1 = 32.1;
double bet2 = 1.72;

double r = dx*dx + dy*dy;
if( r > 25. ) {
 sig2 = 100.;
 return sig2;
}

sig2 = alp*fc + (bet1 + bet2*sqrt(e/100.))*d2c( dx, dy ) + 0.2/(e/100.); 
if( TEST_P1 ) sig2 = sig2/pow(0.0001*e, 0.166);
sig2 *= 100.;

return sig2;
2298}


2301double DCCALShower_factory::d2c( double dx, double dy ) 
2302{
int i, j;
double ax, ay, wx, wy;
double d2c;

ax = fabs(dx*100.);
ay = fabs(dy*100.);
i = int(ax);
j = int(ay);

if( i < 499 && j < 499 && i >= 0 && j >= 0 ) {
//if( i < 500. && j < 500. ) {

 wx = ax-static_cast<double>(i);
 wy = ay-static_cast<double>(j);
 
 d2c = ad2c[i][j]     * (1.-wx) * (1.-wy) + 
 	ad2c[i+1][j]   *     wx  * (1.-wy) + 
ad2c[i][j+1]   * (1.-wx) *     wy  +
ad2c[i+1][j+1] *     wx  *     wy; 

} else d2c = 1.;

return d2c;
2326}


2329double DCCALShower_factory::cell_hyc( double dx, double dy )
2330{
int i, j;
double ax, ay, wx, wy;
double cell_hyc;

ax = fabs(dx*100.);
ay = fabs(dy*100.);
i = static_cast<int>(ax);
j = static_cast<int>(ay);

if( i < 499 && j < 499 && i >= 0 && j >= 0 ) {

 wx = ax-static_cast<double>(i);
 wy = ay-static_cast<double>(j);
 //std::cout << "i " << i << " j " << j << " wx " << wx << " wy " << wy << std::endl;
 cell_hyc = acell[i][j]     * (1.-wx) * (1.-wy) + 
 	     acell[i+1][j]   *     wx  * (1.-wy) + 
     acell[i][j+1]   * (1.-wx) *     wy  +
     acell[i+1][j+1] *     wx  *     wy; 

} else cell_hyc = 0.;

return cell_hyc;
2353}





2359//==========================================================
2360//
2361//   tgamma_hyc( int, vector<int>, vector<int>, int, vector<int>, double,
2362//		double, double, double, double, double, double )
2363//
2364//==========================================================

2366void DCCALShower_factory::tgamma_hyc( int nadc, vector<int> ia, vector<int> id, 
int nzero, vector<int> iaz, double &chisq, double &e1, double &x1, double &y1, 
double &e2, double &x2, double &y2 )
2369{
//-------------- Local Declarations -------------//

int ix, iy;
int dof;

double dx, dy;
double dxy;
double dxc, dyc;
double dx0, dy0;
double dx1, dy1;
double dx2, dy2;
double u, r, rsq, rsq2;
double epsc, eps0, eps1, eps2;
double stepmin, epsmax;
double delch;
double step;
double cosi, scal;
double chisq2, chisqc;
double dchi, dchida, dchi0;
double a1, a2;
double ex;
double e1c, x1c, y1c;
double e2c, x2c, y2c;
double gr, gre, grx, gry;
double grc;
double grec, grxc, gryc;
double gx1, gy1;
double gx2, gy2;
double e0, x0, y0;
double xx, yy, yx;

double f1c, f2c, f1x, f2x, f1y, f2y;
double chisqt, chisqt0, chisqt1, chisqt2;
double chisqtx1, chisqtx2, chisqty1, chisqty2;
double dchidax1, dchidax2, dchiday1, dchiday2;

stepmin = 0.5;
epsmax  = 0.9999;
delch   = 10.;

//-------------- Event Analysis Code -------------//

yx   = 0.;
grx  = 0.; gry = 0.;
gre  = 0.; gr  = 0.;
dx0  = 0.; dy0 = 0.;
eps0 = 0.;

mom2_pht( nadc, ia, id, nzero, iaz, e0, x0, y0, xx, yy, yx );

e2 = 0.;
x2 = 0.;
y2 = 0.;

if( nadc <= 0 ) return;

// choosing of the starting point

dxy  = xx - yy;
rsq2 = dxy*dxy + 4.*yx*yx;
if( rsq2 < 1.e-20 ) rsq2 = 1.e-20;
rsq = sqrt(rsq2);
dxc = -sqrt((rsq+dxy)*2.);
dyc =  sqrt((rsq-dxy)*2.);
if( yx >= 0. ) dyc = -dyc;
r = sqrt(dxc*dxc + dyc*dyc);
epsc = 0.;
for( int ii = 0; ii < nadc; ii++ ) {
 ix   = ia[ii]/100;
 iy   = ia[ii] - ix*100;
 dx   = static_cast<double>(ix) - x0;
 dy   = static_cast<double>(iy) - y0;
 u    = dx*dxc/r + dy*dyc/r;
 epsc = epsc - 0.01*id[ii]*u*fabs(u);
}
epsc = epsc/(0.01*e0*rsq);
if(  epsc > 0.8 ) epsc =  0.8;
if( epsc < -0.8 ) epsc = -0.8;
dxc = dxc/sqrt(1.-(epsc*epsc));
dyc = dyc/sqrt(1.-(epsc*epsc));



// start of iterations:

step   = 0.1;
cosi   = 0.0;
chisq2 = 1.e35;

for( int iter = 0; iter < 100; iter++ ) {
 
 c3to5_pht( e0, x0, y0, epsc, dxc, dyc, e1c, x1c, y1c, e2c, x2c, y2c );
 
 eps1   = (1. + epsc)/2.;
 eps2   = (1. - epsc)/2.;
 chisqc =  0.;
 
 for( int ii = 0; ii < nadc+nzero; ii++ ) { // chi2 calculation
   
   if( ii < nadc ) {
     ex = static_cast<double>(id[ii]);
     ix = ia[ii]/100;
     iy = ia[ii] - ix*100;
   } else {
     ex = 0.;
     ix = iaz[ii-nadc]/100;
     iy = iaz[ii-nadc] - ix*100;
   }
   
   dx1 = x1c - static_cast<double>(ix);
   dy1 = y1c - static_cast<double>(iy);
   dx2 = x2c - static_cast<double>(ix);
   dy2 = y2c - static_cast<double>(iy);
   
   f1c = cell_hyc( dx1, dy1 );
   f2c = cell_hyc( dx2, dy2 );
   
   chisq2t_hyc( ex, e1c, dx1, dy1, e2c, dx2, dy2, f1c, f2c, chisqt );
   chisqc += chisqt;
   
 }
 
 if( chisqc >= chisq2 ) { // new step if no improvement
   
   if( iter > 0 ) {
     dchi  = chisqc - chisq2;
     dchi0 = gr*step;
     step  = 0.5*step/sqrt(1. + dchi/dchi0); 
   }
   step = 0.5*step;
 
 } else { // calculate gradient
   
   grec = 0.;
   grxc = 0.;
   gryc = 0.;
   
   for( int ii = 0; ii < nadc+nzero; ii++ ) {
     
     if( ii < nadc ) {
       ex = static_cast<double>(id[ii]);
ix = ia[ii]/100;
iy = ia[ii] - ix*100;
     } else {
       ex = 0.;
       ix = iaz[ii-nadc]/100;
       iy = iaz[ii-nadc] - ix*100;
     }
     
     dx1 = x1c - static_cast<double>(ix);
     dy1 = y1c - static_cast<double>(iy);
     dx2 = x2c - static_cast<double>(ix);
     dy2 = y2c - static_cast<double>(iy);
     
     f1c = cell_hyc( dx1, dy1 );
     f2c = cell_hyc( dx2, dy2 );
   
     a1 = e1c*f1c;
     a2 = e2c*f2c;
     //a0 = a1 + a2;
     
     chisq2t_hyc( ex, e1c, dx1, dy1, e2c, dx2, dy2, f1c, f2c, chisqt0 );
     chisq2t_hyc( ex, e1c+1., dx1, dy1, e2c, dx2, dy2, f1c, f2c, chisqt1 );
     chisq2t_hyc( ex, e1c, dx1, dy1, e2c+1., dx2, dy2, f1c, f2c, chisqt2 );
     
     f1x = cell_hyc( dx1+0.05, dy1 );
     f2x = cell_hyc( dx2+0.05, dy2 );
     f1y = cell_hyc( dx1, dy1+0.05 );
     f2y = cell_hyc( dx2, dy2+0.05 );
     
     chisq2t_hyc( ex, e1c, dx1+0.05, dy1, e2c, dx2, dy2, f1x, f2c, chisqtx1 );
     chisq2t_hyc( ex, e1c, dx1, dy1, e2c, dx2+0.05, dy2, f1c, f2x, chisqtx2 );
     chisq2t_hyc( ex, e1c, dx1, dy1+0.05, e2c, dx2, dy2, f1y, f2c, chisqty1 );
     chisq2t_hyc( ex, e1c, dx1, dy1, e2c, dx2, dy2+0.05, f1c, f2y, chisqty2 );
     
     dchidax1 = 20.*(chisqtx1 - chisqt0);
     dchidax2 = 20.*(chisqtx2 - chisqt0);
     dchiday1 = 20.*(chisqty1 - chisqt0);
     dchiday2 = 20.*(chisqty2 - chisqt0);
     dchida   = 0.5*(chisqt1 + chisqt2 - chisqt0);
     
     gx1 = (e1c*f1x - a1)*dchidax1;
     gx2 = (e2c*f2x - a2)*dchidax2;
     gy1 = (e1c*f1y - a1)*dchiday1;
     gy2 = (e2c*f2y - a2)*dchiday2;
     
     grec += dchida*(f1c-f2c)*e0 - ( (gx1+gx2)*dxc + (gy1+gy2)*dyc );
     grxc += gx1*eps2 - gx2*eps1;
     gryc += gy1*eps2 - gy2*eps1;
     
   }
   
   grc = sqrt( grec*grec + grxc*grxc + gryc*gryc );
   if( grc < 1.e-6 ) grc = 1.e-6;
   if( iter > 0 ) {
     cosi = ( gre*grec + grx*grxc + gry*gryc )/(gr*grc);
     scal = fabs((gr/grc) - cosi);
     if( scal < 0.1 ) scal = 0.1;
     step = step/scal;
   }
   chisq2 = chisqc;
   eps0   = epsc;
   dx0    = dxc;
   dy0    = dyc;
   gre    = grec;
   grx    = grxc;
   gry    = gryc;
   gr     = grc;
 }
 epsc = eps0 - step*gre/gr;
 while( fabs(epsc) > epsmax ) {
   step = step/2.;
   epsc = eps0 - step*gre/gr;
 }
 dxc = dx0 - step*grx/gr;
 dyc = dy0 - step*gry/gr;
 if( step*gr < stepmin ) break;
}

if( (chisq*(nadc+nzero-2) - chisq2) < delch ) return;
dof = nzero+nadc-5;
if( dof < 1 ) dof = 1;
chisq = chisq2/dof;

c3to5_pht( e0, x0, y0, eps0, dx0, dy0, e1, x1, y1, e2, x2, y2 );

return;
2597}





2603//==========================================================
2604//
2605//   mom2_pht( int, vector<int>, vector<int>, int, vector<int>,
2606//		double, double, double, double, double, double )
2607//
2608//==========================================================

2610void DCCALShower_factory::mom2_pht( int nadc, vector<int> ia, vector<int> id, 
int nzero, vector<int> iaz, double &a0, double &x0, double &y0, 
double &xx, double &yy, double &yx)
2613{
//-------------- Local Declarations --------------//

int ix, iy;
double dx, dy;
double a;
double corr;

//-------------- Event Analysis Code -------------//

a0 = 0.;
x0 = 0.;
y0 = 0.;
xx = 0.;
yy = 0.;
yx = 0.;

if( nadc <= 0 ) return;
for( int ii = 0; ii < nadc; ii++ ) {
 a  = static_cast<double>(id[ii]);
 ix = ia[ii]/100;
 iy = ia[ii] - ix*100;
 a0 = a0 + a;
 x0 = x0 + a*static_cast<double>(ix);
 y0 = y0 + a*static_cast<double>(iy);
}
if( a0 <= 0. ) return;
x0 = x0/a0;
y0 = y0/a0;

for( int ii = 0; ii < nadc; ii++ ) {
 a  = static_cast<double>(id[ii])/a0;
 ix = ia[ii]/100;
 iy = ia[ii] - ix*100;
 dx = static_cast<double>(ix) - x0;
 dy = static_cast<double>(iy) - y0;
 xx = xx + a*dx*dx;
 yy = yy + a*dy*dy;
 yx = yx + a*dx*dy;
}

corr = 0.;
for( int ii = 0; ii < nadc; ii++ ) {
 ix = ia[ii]/100;
 iy = ia[ii] - ix*100;
 dx = static_cast<double>(ix) - x0;
 dy = static_cast<double>(iy) - y0;
 corr = corr + cell_hyc( dx, dy );
}
if( COUNT_ZERO1 ) {
 for( int ii = 0; ii < nzero; ii++ ) {
   ix = iaz[ii]/100;
   iy = iaz[ii] - ix*100;
   dx = static_cast<double>(ix) - x0;
   dy = static_cast<double>(iy) - y0;
   corr = corr + cell_hyc( dx, dy );
 }
}

corr = corr/1.006;
if( corr < 0.8 ) {
 corr = 0.8;
} else if( corr > 1. ) {
 corr = 1.0;
}
a0 = a0/corr;


return;
2682}





2688//==========================================================
2689//
2690//   c3to5_pht( double, double, double, double, double, double,
2691//		double, double, double, double, double, double )
2692//
2693//==========================================================

2695void DCCALShower_factory::c3to5_pht( double e0, double x0, double y0, double eps, 
double dx, double dy, double &e1, double &x1, double &y1, double &e2, 
double &x2, double &y2 )
2698{
e1 = e0*(1.+eps)/2.;
e2 = e0 - e1;
x1 = x0 + 0.5*dx*(1.-eps);
y1 = y0 + 0.5*dy*(1.-eps);
x2 = x0 - 0.5*dx*(1.+eps);
y2 = y0 - 0.5*dy*(1.+eps);

return;
2707}





2713//==========================================================
2714//
2715//   chisq2t_hyc( double, double, double, double, double, double,
2716//		double, double, double, double )
2717//
2718//==========================================================

2720void DCCALShower_factory::chisq2t_hyc( double ecell, double e1, double dx1, double dy1, 
double e2, double dx2, double dy2, double f1, double f2, double &chisqt )
2722{
double s;
double p1, p2;

if( TEST_P1 ) {
 p1 = pow(0.0001*e1, 0.166);
 p2 = pow(0.0001*e2, 0.166);
} else {
 p1 = 1.;
 p2 = 1.;
}

if( e1 != 0. && e2 != 0. )
 s = e1*sigma2(dx1, dy1, f1, e1)/p1 + e2*sigma2(dx2, dy2, f2, e2)/p2;
else if( e1 == 0. && e2 == 0. )
 s = 90000.;
else if( e1 == 0. )
 s = e2*sigma2(dx2, dy2, f2, e2)/p2;
else 
 s = e1*sigma2(dx1, dy1, f1, e1)/p1;

chisqt = pow((e1*f1 + e2*f2 - ecell), 2.)/s;

return;
2746}






2753//==========================================================
2754//
2755//   ucopy1( vector<int>, vector<int>, int, int )
2756//
2757//==========================================================

2759void DCCALShower_factory::ucopy1( vector<int> &ia, vector<int> &iwork, int start, int length )
2760{

for( int ii = 0; ii < length; ii++ ) {
 iwork.push_back( ia[start+ii] );
}


return;
2768}



2772//==========================================================
2773//
2774//   ucopy2( vector<int>, int, int, int )
2775//
2776//==========================================================

2778void DCCALShower_factory::ucopy2( vector<int> &ia, int start1, int start2, int length )
2779{

vector<int> work;
for( int ii = 0; ii < length; ii++ ) {
 work.push_back( ia[start1+ii] );
}

for( int ii = 0; ii < length; ii++ ) {
 ia[start2+ii] = work[ii];
}


return;
2792}




2797//==========================================================
2798//
2799//   ucopy3( vector<int>, vector<int>, int, int )
2800//
2801//==========================================================

2803void DCCALShower_factory::ucopy3( vector<int> &iwork, vector<int> &ia, int start, int length )
2804{

for( int ii = 0; ii < length; ii++ ) {
 ia[start+ii] = iwork[ii];	
}


return;
2812}























←

/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h

→

1// $Id: JEventLoop.h 1763 2006-05-10 14:29:25Z davidl $
2//
3//    File: JEventLoop.h
4// Created: Wed Jun  8 12:30:51 EDT 2005
5// Creator: davidl (on Darwin wire129.jlab.org 7.8.0 powerpc)
6//

8#ifndef _JEventLoop_
9#define _JEventLoop_

11#include <sys/time.h>

13#include <vector>
14#include <list>
15#include <string>
16#include <utility>
17#include <typeinfo>
18#include <string.h>
19#include <map>
20#include <utility>
21using std::vector;
22using std::list;
23using std::string;
24using std::type_info;

26#include <JANA/jerror.h>
27#include <JANA/JObject.h>
28#include <JANA/JException.h>
29#include <JANA/JEvent.h>
30#include <JANA/JThread.h>
31#include <JANA/JFactory_base.h>
32#include <JANA/JCalibration.h>
33#include <JANA/JGeometry.h>
34#include <JANA/JResourceManager.h>
35#include <JANA/JStreamLog.h>

37// The following is here just so we can use ROOT's THtml class to generate documentation.
38#include "cint.h"


41// Place everything in JANA namespace
42namespace jana{


45template<class T> class JFactory;
46class JApplication;
47class JEventProcessor;


50class JEventLoop{
public:

friend class JApplication;

enum data_source_t{
	DATA_NOT_AVAILABLE = 1,
	DATA_FROM_CACHE,
	DATA_FROM_SOURCE,
	DATA_FROM_FACTORY
};

typedef struct{
	string caller_name;
	string caller_tag;
	string callee_name;
	string callee_tag;
	double start_time;
	double end_time;
	data_source_t data_source;
}call_stack_t;

typedef struct{
	const char* factory_name;
	string tag;
	const char* filename;
	int line;
}error_call_stack_t;

                                  JEventLoop(JApplication *app); ///< Constructor
					   virtual ~JEventLoop(); ////< Destructor
		   virtual const char* className(void){return static_className();}
			static const char* static_className(void){return "JEventLoop";}

				 JApplication* GetJApplication(void) const {return app;} ///< Get pointer to the JApplication object
                                void RefreshProcessorListFromJApplication(void); ///< Re-copy the list of JEventProcessors from JApplication
                    virtual jerror_t AddFactory(JFactory_base* factory); ///< Add a factory
					  jerror_t RemoveFactory(JFactory_base* factory); ///< Remove a factory
                      JFactory_base* GetFactory(const string data_name, const char *tag="", bool allow_deftag=true); ///< Get a specific factory pointer
              vector<JFactory_base*> GetFactories(void) const {return factories;} ///< Get all factory pointers
                                void GetFactoryNames(vector<string> &factorynames); ///< Get names of all factories in name:tag format
                                void GetFactoryNames(map<string,string> &factorynames); ///< Get names of all factories in map with key=name, value=tag
                  map<string,string> GetDefaultTags(void) const {return default_tags;}
                            jerror_t ClearFactories(void); ///< Reset all factories in preparation for next event.
					  jerror_t PrintFactories(int sparsify=0); ///< Print a list of all factories.
                            jerror_t Print(const string data_name, const char *tag=""); ///< Print the data of the given type

				 JCalibration* GetJCalibration();
              template<class T> bool GetCalib(string namepath, map<string,T> &vals);
              template<class T> bool GetCalib(string namepath, vector<T> &vals);
              template<class T> bool GetCalib(string namepath, T &val);

                          JGeometry* GetJGeometry();
              template<class T> bool GetGeom(string namepath, map<string,T> &vals);
	    template<class T> bool GetGeom(string namepath, T &val);

                   JResourceManager* GetJResourceManager(void);
                              string GetResource(string namepath);
              template<class T> bool GetResource(string namepath, T vals, int event_number=0);

                                void Initialize(void); ///< Do initializations just before event processing starts
                            jerror_t Loop(void); ///< Loop over events
                            jerror_t OneEvent(uint64_t event_number); ///< Process a specific single event (if source supports it)
                            jerror_t OneEvent(void); ///< Process a single event
                         inline void Pause(void){pause = 1;} ///< Pause event processing
                         inline void Resume(void){pause = 0;} ///< Resume event processing
                         inline void Quit(void){quit = 1;} ///< Clean up and exit the event loop
                         inline bool GetQuit(void) const {return quit;}
                                void QuitProgram(void);

                               // Support for random access of events
                          bool HasRandomAccess(void);
                          void AddToEventQueue(uint64_t event_number){ next_events_to_process.push_back(event_number); }
                          void AddToEventQueue(list<uint64_t> &event_numbers) { next_events_to_process.insert(next_events_to_process.end(), event_numbers.begin(), event_numbers.end()); }
                list<uint64_t> GetEventQueue(void){ return next_events_to_process; }
                          void ClearEventQueue(void){ next_events_to_process.clear(); }

      template<class T> JFactory<T>* GetSingle(const T* &t, const char *tag="", bool exception_if_not_one=true); ///< Get pointer to first data object from (source or factory).
      template<class T> JFactory<T>* Get(vector<const T*> &t, const char *tag="", bool allow_deftag=true); ///< Get data object pointers from (source or factory)
      template<class T> JFactory<T>* GetFromFactory(vector<const T*> &t, const char *tag="", data_source_t &data_source=null_data_source, bool allow_deftag=true); ///< Get data object pointers from factory
          template<class T> jerror_t GetFromSource(vector<const T*> &t, JFactory_base *factory=NULL__null); ///< Get data object pointers from source.
                      inline JEvent& GetJEvent(void){return event;} ///< Get pointer to the current JEvent object.
                         inline void SetJEvent(JEvent *event){this->event = *event;} ///< Set the JEvent pointer.
                         inline void SetAutoFree(int auto_free){this->auto_free = auto_free;} ///< Set the Auto-Free flag on/off
                    inline pthread_t GetPThreadID(void) const {return pthread_id;} ///< Get the pthread of the thread to which this JEventLoop belongs
                              double GetInstantaneousRate(void) const {return rate_instantaneous;} ///< Get the current event processing rate
                              double GetIntegratedRate(void) const {return rate_integrated;} ///< Get the current event processing rate
                              double GetLastEventProcessingTime(void) const {return delta_time_single;}
                        unsigned int GetNevents(void) const {return Nevents;}

                         inline bool CheckEventBoundary(uint64_t event_numberA, uint64_t event_numberB);

                         inline bool GetCallStackRecordingStatus(void){ return record_call_stack; }
                         inline void DisableCallStackRecording(void){ record_call_stack = false; }
                         inline void EnableCallStackRecording(void){ record_call_stack = true; }
                     inline void CallStackStart(JEventLoop::call_stack_t &cs, const string &caller_name, const string &caller_tag, const string callee_name, const string callee_tag);
                     inline void CallStackEnd(JEventLoop::call_stack_t &cs);
         inline vector<call_stack_t> GetCallStack(void){return call_stack;} ///< Get the current factory call stack
                         inline void AddToCallStack(call_stack_t &cs){if(record_call_stack) call_stack.push_back(cs);} ///< Add specified item to call stack record but only if record_call_stack is true
                         inline void AddToErrorCallStack(error_call_stack_t &cs){error_call_stack.push_back(cs);} ///< Add layer to the factory call stack
   inline vector<error_call_stack_t> GetErrorCallStack(void){return error_call_stack;} ///< Get the current factory error call stack
                                void PrintErrorCallStack(void); ///< Print the current factory call stack

                      const JObject* FindByID(JObject::oid_t id); ///< Find a data object by its identifier.
          template<class T> const T* FindByID(JObject::oid_t id); ///< Find a data object by its type and identifier
                      JFactory_base* FindOwner(const JObject *t); ///< Find the factory that owns a data object by pointer
                      JFactory_base* FindOwner(JObject::oid_t id); ///< Find a factory that owns a data object by identifier

                                     // User defined references
              template<class T> void SetRef(T *t);    ///< Add a user reference to this JEventLoop (must be a pointer)
                template<class T> T* GetRef(void);   ///< Get a user-defined reference of a specific type
        template<class T> vector<T*> GetRefsT(void); ///< Get all user-defined refrences of a specicif type
   vector<pair<const char*, void*> > GetRefs(void){ return user_refs; }  ///< Get copy of full list of user-defined references
              template<class T> void RemoveRef(T *t); ///< Remove user reference from list

							   // Convenience methods wrapping JEvent methods of same name
                      uint64_t GetStatus(void){return event.GetStatus();}
                          bool GetStatusBit(uint32_t bit){return event.GetStatusBit(bit);}
                          bool SetStatusBit(uint32_t bit, bool val=true){return event.SetStatusBit(bit, val);}
                          bool ClearStatusBit(uint32_t bit){return event.ClearStatusBit(bit);}
                          void ClearStatus(void){event.ClearStatus();}
                          void SetStatusBitDescription(uint32_t bit, string description){event.SetStatusBitDescription(bit, description);}
                        string GetStatusBitDescription(uint32_t bit){return event.GetStatusBitDescription(bit);}
                          void GetStatusBitDescriptions(map<uint32_t, string> &status_bit_descriptions){return event.GetStatusBitDescriptions(status_bit_descriptions);}
                              string StatusWordToString(void);

private:
JEvent event;
vector<JFactory_base*> factories;
vector<JEventProcessor*> processors;
vector<error_call_stack_t> error_call_stack;
vector<call_stack_t> call_stack;
JApplication *app;
JThread *jthread;
bool initialized;
bool print_parameters_called;
int pause;
int quit;
int auto_free;
pthread_t pthread_id;
map<string, string> default_tags;
vector<pair<string,string> > auto_activated_factories;
bool record_call_stack;
string caller_name;
string caller_tag;
vector<uint64_t> event_boundaries;
int32_t event_boundaries_run; ///< Run number boundaries were retrieved from (possbily 0)
list<uint64_t> next_events_to_process;

uint64_t Nevents;			      ///< Total events processed (this thread)
uint64_t Nevents_rate;		   ///< Num. events accumulated for "instantaneous" rate
double delta_time_single;		///< Time spent processing last event
double delta_time_rate;			///< Integrated time accumulated "instantaneous" rate (partial number of events)
double delta_time;				///< Total time spent processing events (this thread)
double rate_instantaneous;		///< Latest instantaneous rate
double rate_integrated;			///< Rate integrated over all events
 
static data_source_t null_data_source;

vector<pair<const char*, void*> > user_refs;
210};


213// The following is here just so we can use ROOT's THtml class to generate documentation.
214#ifdef G__DICTIONARY
215typedef JEventLoop::call_stack_t call_stack_t;
216typedef JEventLoop::error_call_stack_t error_call_stack_t;
217#endif

219#if !defined(__CINT__) && !defined(__CLING__)

221//-------------
222// GetSingle
223//-------------
224template<class T>
225JFactory<T>* JEventLoop::GetSingle(const T* &t, const char *tag, bool exception_if_not_one)
226{
/// This is a convenience method that can be used to get a pointer to the single
/// object of type T from the specified factory. It simply calls the Get(vector<...>) method
/// and copies the first pointer into "t" (or NULL if something other than 1 object is returned).
/// 
/// This is intended to address the common situation in which there is an interest
/// in the event if and only if there is exactly 1 object of type T. If the event
/// has no objects of that type or more than 1 object of that type (for the specified
/// factory) then an exception of type "unsigned long" is thrown with the value
/// being the number of objects of type T. You can supress the exception by setting
/// exception_if_not_one to false. In that case, you will have to check if t==NULL to
/// know if the call succeeded.
vector<const T*> v;
JFactory<T> *fac = Get(v, tag);

if(v.size()!=1){
t = NULL__null;
if(exception_if_not_one) throw v.size();
}

t = v[0];

return fac;
249}

251//-------------
252// Get
253//-------------
254template<class T> 
255JFactory<T>* JEventLoop::Get(vector<const T*> &t, const char *tag, bool allow_deftag)
256{
/// Retrieve or generate the array of objects of
/// type T for the curent event being processed
/// by this thread.
///
/// By default, preference is given to reading the
/// objects from the data source(e.g. file) before generating
/// them in the factory. A flag exists in the factory
/// however to change this so that the factory is
/// given preference.
///
/// Note that regardless of the setting of this flag,
/// the data are only either read in or generated once.
/// Ownership of the objects will always be with the
/// factory so subsequent calls will always return pointers to
/// the same data.
///
/// If the factory is called on to generate the data,
/// it is done by calling the factory's Get() method
/// which, in turn, calls the evnt() method. 
/// 
/// First, we just call the GetFromFactory() method.
/// It will make the initial decision as to whether
/// it should look in the source first or not. If
/// it returns NULL, then the factory couldn't be
/// found so we automatically try the file.
///
/// Note that if no factory exists to hold the objects
/// from the file, one can be created automatically
/// providing the <i>JANA:AUTOFACTORYCREATE</i>
/// configuration parameter is set.

// Check if a tag was specified for this data type to use for the
// default.
const char *mytag = tag1.1
'tag' is not equal to NULL
1
'tag' is not equal to NULL
1
'tag' is not equal to NULL
==NULL__null ? "":tag; // protection against NULL tags
2
←
'?' condition is false→
if(strlen(mytag)==0 && allow_deftag2.1
'allow_deftag' is true
1
'allow_deftag' is true
1
'allow_deftag' is true
){
3
←
Taking true branch→
map<string, string>::const_iterator iter = default_tags.find(T::static_className());
if(iter!=default_tags.end())tag = iter->second.c_str();
4
←
Assuming the condition is true→
5
←
Taking true branch→
6
←
Value assigned to 'tag'→
}


// If we are trying to keep track of the call stack then we
// need to add a new call_stack_t object to the the list
// and initialize it with the start time and caller/callee
// info.
call_stack_t cs;

// Optionally record starting info of call stack entry
if(record_call_stack) CallStackStart(cs, caller_name, caller_tag, T::static_className(), tag);
7
←
Assuming field 'record_call_stack' is false→
8
←
Taking false branch→

// Get the data (or at least try to)
JFactory<T>* factory=NULL__null;
try{
factory = GetFromFactory(t, tag, cs.data_source, allow_deftag);
9
←
Passing value via 2nd parameter 'tag'→
10
←
Calling 'JEventLoop::GetFromFactory'→
if(!factory){
	// No factory exists for this type and tag. It's possible
	// that the source may be able to provide the objects
	// but it will need a place to put them. We can create a
	// dumb JFactory just to hold the data in case the source
	// can provide the objects. Before we do though, make sure
	// the user condones this via the presence of the
	// "JANA:AUTOFACTORYCREATE" config parameter.
	string p;
	try{
		gPARMS->GetParameter("JANA:AUTOFACTORYCREATE", p);
	}catch(...){}
	if(p.size()==0){
		jout<<std::endl;
		_DBG__std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<324<<std::endl;
		jout<<"No factory of type \""<<T::static_className()<<"\" with tag \""<<tag<<"\" exists."<<std::endl;
		jout<<"If you are reading objects from a file, I can auto-create a factory"<<std::endl;
		jout<<"of the appropriate type to hold the objects, but this feature is turned"<<std::endl;
		jout<<"off by default. To turn it on, set the \"JANA:AUTOFACTORYCREATE\""<<std::endl;
		jout<<"configuration parameter. This can usually be done by passing the"<<std::endl;
		jout<<"following argument to the program from the command line:"<<std::endl;
		jout<<std::endl;
		jout<<"   -PJANA:AUTOFACTORYCREATE=1"<<std::endl;
		jout<<std::endl;
		jout<<"Note that since the most commonly expected occurance of this situation."<<std::endl;
		jout<<"is an error, the program will now throw an exception so that the factory."<<std::endl;
		jout<<"call stack can be printed."<<std::endl;
		jout<<std::endl;
		throw exception();
	}else{
		AddFactory(new JFactory<T>(tag));
		jout<<__FILE__"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"<<":"<<__LINE__341<<" Auto-created "<<T::static_className()<<":"<<tag<<" factory"<<std::endl;
	
		// Now try once more. The GetFromFactory method will call
		// GetFromSource since it's empty.
		factory = GetFromFactory(t, tag, cs.data_source, allow_deftag);
	}
}
}catch(exception &e){
// Uh-oh, an exception was thrown. Add us to the call stack
// and re-throw the exception
error_call_stack_t ecs;
ecs.factory_name = T::static_className();
ecs.tag = tag;
ecs.filename = NULL__null;
error_call_stack.push_back(ecs);
throw e;
}

// If recording the call stack, update the end_time field and add to stack
if(record_call_stack) CallStackEnd(cs);

return factory;
363}

365//-------------
366// GetFromFactory
367//-------------
368template<class T> 
369JFactory<T>* JEventLoop::GetFromFactory(vector<const T*> &t, const char *tag, data_source_t &data_source, bool allow_deftag)
370{
// We need to find the factory providing data type T with
// tag given by "tag". 
vector<JFactory_base*>::iterator iter=factories.begin();
JFactory<T> *factory = NULL__null;
string className(T::static_className());

// Check if a tag was specified for this data type to use for the
// default.
const char *mytag = tag==NULL__null ? "":tag; // protection against NULL tags
11
←
Assuming 'tag' is equal to NULL→
12
←
Assuming pointer value is null→
13
←
'?' condition is true→
if(strlen(mytag)==0 && allow_deftag13.1
'allow_deftag' is true
1
'allow_deftag' is true
1
'allow_deftag' is true
){
14
←
Taking true branch→
map<string, string>::const_iterator iter = default_tags.find(className);
if(iter!=default_tags.end())tag = iter->second.c_str();
15
←
Assuming the condition is false→
16
←
Taking false branch→
}

for(; iter!=factories.end(); iter++){
17
←
Calling 'operator!=<jana::JFactory_base **, std::vector<jana::JFactory_base *>>'→
20
←
Returning from 'operator!=<jana::JFactory_base **, std::vector<jana::JFactory_base *>>'→
21
←
Loop condition is true.  Entering loop body→
// It turns out a long standing bug in g++ makes dynamic_cast return
// zero improperly when used on objects created on one side of
// a dynamically shared object (DSO) and the cast occurs on the 
// other side. I saw bug reports ranging from 2001 to 2004. I saw
// saw it first-hand on LinuxEL4 using g++ 3.4.5. This is too bad
// since it is much more elegant (and safe) to use dynamic_cast.
// To avoid this problem which can occur with plugins, we check
// the name of the data classes are the same. (sigh)
//factory = dynamic_cast<JFactory<T> *>(*iter);
if(className == (*iter)->GetDataClassName())factory = (JFactory<T>*)(*iter);
22
←
Taking true branch→
if(factory == NULL__null)continue;
23
←
Assuming 'factory' is not equal to NULL→
24
←
Taking false branch→
const char *factag = factory->Tag()==NULL__null ? "":factory->Tag();
25
←
Assuming the condition is true→
26
←
'?' condition is true→
if(!strcmp(factag, tag)){
27
←
Null pointer passed to 2nd parameter expecting 'nonnull'
	break;
}else{
	factory=NULL__null;
}
}

// If factory not found, just return now
if(!factory){
data_source = DATA_NOT_AVAILABLE;
return NULL__null;
}

// OK, we found the factory. If the evnt() routine has already
// been called, then just call the factory's Get() routine
// to return a copy of the existing data
if(factory->evnt_was_called()){
factory->CopyFrom(t);
data_source = DATA_FROM_CACHE;
return factory;
}

// Next option is to get the objects from the data source
if(factory->GetCheckSourceFirst()){
// If the object type/tag is found in the source, it
// will return NOERROR, even if there are zero instances
// of it. If it is not available in the source then it
// will return OBJECT_NOT_AVAILABLE.

jerror_t err = GetFromSource(t, factory);
if(err == NOERROR){
	// A return value of NOERROR means the source had the objects
	// even if there were zero of them.(If the source had no
	// information about the objects OBJECT_NOT_AVAILABLE would 
	// have been returned.)
	// The GetFromSource() call will eventually lead to a call to
	// the GetObjects() method of the concrete class derived
	// from JEventSource. That routine should copy the object
	// pointers into the factory using the factory's CopyTo()
	// method which also sets the evnt_called flag for the factory.
	// Note also that the "t" vector is then filled with a call
	// to the factory's CopyFrom() method in JEvent::GetObjects().
	// All we need to do now is just set the factory pointers in
	// the newly generated JObjects and return the factory pointer.

	factory->SetFactoryPointers();
	data_source = DATA_FROM_SOURCE;

	return factory;
}
}

// OK. It looks like we have to have the factory make this.
// Get pointers to data from the factory.
factory->Get(t);
factory->SetFactoryPointers();
data_source = DATA_FROM_FACTORY;

return factory;
457}

459//-------------
460// GetFromSource
461//-------------
462template<class T> 
463jerror_t JEventLoop::GetFromSource(vector<const T*> &t, JFactory_base *factory)
464{
/// This tries to get objects from the event source.
/// "factory" must be a valid pointer to a JFactory
/// object since that will take ownership of the objects
/// created by the source.
/// This should usually be called from JEventLoop::GetFromFactory
/// which is called from JEventLoop::Get. The latter will
/// create a dummy JFactory of the proper flavor and tag if
/// one does not already exist so if objects exist in the
/// file without a corresponding factory to create them, they
/// can still be used.
if(!factory)throw OBJECT_NOT_AVAILABLE;

return event.GetObjects(t, factory);
478}

480//-------------
481// CallStackStart
482//-------------
483inline void JEventLoop::CallStackStart(JEventLoop::call_stack_t &cs, const string &caller_name, const string &caller_tag, const string callee_name, const string callee_tag)
484{
/// This is used to fill initial info into a call_stack_t stucture
/// for recording the call stack. It should be matched with a call
/// to CallStackEnd. It is normally called from the Get() method
/// above, but may also be used by external actors to manipulate
/// the call stack (presumably for good and not evil).

struct itimerval tmr;
getitimer(ITIMER_PROFITIMER_PROF, &tmr);

cs.caller_name    = this->caller_name;
cs.caller_tag     = this->caller_tag;
this->caller_name = cs.callee_name = callee_name;
this->caller_tag  = cs.callee_tag  = callee_tag;
cs.start_time = tmr.it_value.tv_sec + tmr.it_value.tv_usec/1.0E6;
499}

501//-------------
502// CallStackEnd
503//-------------
504inline void JEventLoop::CallStackEnd(JEventLoop::call_stack_t &cs)
505{
/// Complete a call stack entry. This should be matched
/// with a previous call to CallStackStart which was
/// used to fill the cs structure.

struct itimerval tmr;
getitimer(ITIMER_PROFITIMER_PROF, &tmr);
cs.end_time = tmr.it_value.tv_sec + tmr.it_value.tv_usec/1.0E6;
caller_name = cs.caller_name;
caller_tag  = cs.caller_tag;
call_stack.push_back(cs);
516}

518//-------------
519// CheckEventBoundary
520//-------------
521inline bool JEventLoop::CheckEventBoundary(uint64_t event_numberA, uint64_t event_numberB)
522{
/// Check whether the two event numbers span one or more boundaries
/// in the calibration/conditions database for the current run number.
/// Return true if they do and false if they don't. The first parameter
/// "event_numberA" is also checked if it lands on a boundary in which
/// case true is also returned. If event_numberB lands on a boundary,
/// but event_numberA does not, then false is returned.
///
/// This method is not expected to be called by a user. It is, however called,
/// everytime a JFactory's Get() method is called.

// Make sure our copy of the boundaries is up to date
if(event.GetRunNumber()!=event_boundaries_run){
event_boundaries.clear(); // in case we can't get the JCalibration pointer
JCalibration *jcalib = GetJCalibration();
if(jcalib)jcalib->GetEventBoundaries(event_boundaries);
event_boundaries_run = event.GetRunNumber();
}

// Loop over boundaries
for(unsigned int i=0; i<event_boundaries.size(); i++){
uint64_t eb = event_boundaries[i];
if((eb - event_numberA)*(eb - event_numberB) < 0.0 || eb==event_numberA){ // think about it ....
	// events span a boundary or is on a boundary. Return true
	return true;
}
}

return false;
551}

553//-------------
554// FindByID
555//-------------
556template<class T> 
557const T* JEventLoop::FindByID(JObject::oid_t id)
558{
/// This is a templated method that can be used in place
/// of the non-templated FindByID(oid_t) method if one knows
/// the class of the object with the specified id.
/// This method is faster than calling the non-templated
/// FindByID and dynamic_cast-ing the JObject since
/// this will only search the objects of factories that
/// produce the desired data type.
/// This method will cast the JObject pointer to one
/// of the specified type. To use this method,
/// a type is specified in the call as follows:
///
/// const DMyType *t = loop->FindByID<DMyType>(id);

// Loop over factories looking for ones that provide
// specified data type.
for(unsigned int i=0; i<factories.size(); i++){
if(factories[i]->GetDataClassName() != T::static_className())continue;

// This factory provides data of type T. Search it for
// the object with the specified id.
const JObject *my_obj = factories[i]->GetByID(id);
if(my_obj)return dynamic_cast<const T*>(my_obj);
}

return NULL__null;
584}

586//-------------
587// GetCalib (map)
588//-------------
589template<class T>
590bool JEventLoop::GetCalib(string namepath, map<string,T> &vals)
591{
/// Get the JCalibration object from JApplication for the run number of
/// the current event and call its Get() method to get the constants.

// Note that we could do this by making "vals" a generic type T thus, combining
// this with the vector version below. However, doing this explicitly will make
// it easier for the user to understand how to call us.

vals.clear();

JCalibration *calib = GetJCalibration();
if(!calib){
_DBG_std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<603<<" "<<"Unable to get JCalibration object for run "<<event.GetRunNumber()<<std::endl;
return true;
}

return calib->Get(namepath, vals, event.GetEventNumber());
608}

610//-------------
611// GetCalib (vector)
612//-------------
613template<class T> bool JEventLoop::GetCalib(string namepath, vector<T> &vals)
614{
/// Get the JCalibration object from JApplication for the run number of
/// the current event and call its Get() method to get the constants.

vals.clear();

JCalibration *calib = GetJCalibration();
if(!calib){
_DBG_std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<622<<" "<<"Unable to get JCalibration object for run "<<event.GetRunNumber()<<std::endl;
return true;
}

return calib->Get(namepath, vals, event.GetEventNumber());
627}

629//-------------
630// GetCalib (single)
631//-------------
632template<class T> bool JEventLoop::GetCalib(string namepath, T &val)
633{
/// This is a convenience method for getting a single entry. It
/// simply calls the vector version and returns the first entry.
/// It returns true if the vector version returns true AND there
/// is at least one entry in the vector. No check is made for there
/// there being more than one entry in the vector.

vector<T> vals;
bool ret = GetCalib(namepath, vals);
if(vals.empty()) return true;
val = vals[0];

return ret;
646}

648//-------------
649// GetGeom (map)
650//-------------
651template<class T>
652bool JEventLoop::GetGeom(string namepath, map<string,T> &vals)
653{
/// Get the JGeometry object from JApplication for the run number of
/// the current event and call its Get() method to get the constants.

// Note that we could do this by making "vals" a generic type T thus, combining
// this with the vector version below. However, doing this explicitly will make
// it easier for the user to understand how to call us.

vals.clear();

JGeometry *geom = GetJGeometry();
if(!geom){
_DBG_std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<665<<" "<<"Unable to get JGeometry object for run "<<event.GetRunNumber()<<std::endl;
return true;
}

return geom->Get(namepath, vals);
670}

672//-------------
673// GetGeom (atomic)
674//-------------
675template<class T> bool JEventLoop::GetGeom(string namepath, T &val)
676{
/// Get the JCalibration object from JApplication for the run number of
/// the current event and call its Get() method to get the constants.

JGeometry *geom = GetJGeometry();
if(!geom){
_DBG_std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<682<<" "<<"Unable to get JGeometry object for run "<<event.GetRunNumber()<<std::endl;
return true;
}

return geom->Get(namepath, val);
687}

689//-------------
690// SetRef
691//-------------
692template<class T>
693void JEventLoop::SetRef(T *t)
694{
pair<const char*, void*> p(typeid(T).name(), (void*)t);
user_refs.push_back(p);
697}

699//-------------
700// GetResource
701//-------------
702template<class T> bool JEventLoop::GetResource(string namepath, T vals, int event_number)
703{
JResourceManager *resource_manager = GetJResourceManager();
if(!resource_manager){
string mess = string("Unable to get the JResourceManager object (namepath=\"")+namepath+"\")";
throw JException(mess);
}

return resource_manager->Get(namepath, vals, event_number);
711}

713//-------------
714// GetRef
715//-------------
716template<class T>
717T* JEventLoop::GetRef(void)
718{
/// Get a user-defined reference (a pointer)
for(unsigned int i=0; i<user_refs.size(); i++){
if(user_refs[i].first == typeid(T).name()) return (T*)user_refs[i].second;
}

return NULL__null;
725}

727//-------------
728// GetRefsT
729//-------------
730template<class T>
731vector<T*> JEventLoop::GetRefsT(void)
732{
vector<T*> refs;
for(unsigned int i=0; i<user_refs.size(); i++){
if(user_refs[i].first == typeid(T).name()){
	refs.push_back((T*)user_refs[i].second);
}
}

return refs;
741}

743//-------------
744// RemoveRef
745//-------------
746template<class T>
747void JEventLoop::RemoveRef(T *t)
748{
vector<pair<const char*, void*> >::iterator iter;
for(iter=user_refs.begin(); iter!= user_refs.end(); iter++){
if(iter->second == (void*)t){
	user_refs.erase(iter);
	return;
}
}
_DBG_std::cerr<<"/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JEventLoop.h"
<<":"<<756<<" "<<" Attempt to remove user reference not in event loop!" << std::endl;
757}


760#endif //__CINT__  __CLING__

762} // Close JANA namespace



766#endif // _JEventLoop_


←

/usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5/bits/stl_iterator.h

1// Iterators -*- C++ -*-

3// Copyright (C) 2001-2013 Free Software Foundation, Inc.
4//
5// This file is part of the GNU ISO C++ Library.  This library is free
6// software; you can redistribute it and/or modify it under the
7// terms of the GNU General Public License as published by the
8// Free Software Foundation; either version 3, or (at your option)
9// any later version.

11// This library is distributed in the hope that it will be useful,
12// but WITHOUT ANY WARRANTY; without even the implied warranty of
13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14// GNU General Public License for more details.

16// Under Section 7 of GPL version 3, you are granted additional
17// permissions described in the GCC Runtime Library Exception, version
18// 3.1, as published by the Free Software Foundation.

20// You should have received a copy of the GNU General Public License and
21// a copy of the GCC Runtime Library Exception along with this program;
22// see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
23// <http://www.gnu.org/licenses/>.

25/*
*
* Copyright (c) 1994
* Hewlett-Packard Company
*
* Permission to use, copy, modify, distribute and sell this software
* and its documentation for any purpose is hereby granted without fee,
* provided that the above copyright notice appear in all copies and
* that both that copyright notice and this permission notice appear
* in supporting documentation.  Hewlett-Packard Company makes no
* representations about the suitability of this software for any
* purpose.  It is provided "as is" without express or implied warranty.
*
*
* Copyright (c) 1996-1998
* Silicon Graphics Computer Systems, Inc.
*
* Permission to use, copy, modify, distribute and sell this software
* and its documentation for any purpose is hereby granted without fee,
* provided that the above copyright notice appear in all copies and
* that both that copyright notice and this permission notice appear
* in supporting documentation.  Silicon Graphics makes no
* representations about the suitability of this software for any
* purpose.  It is provided "as is" without express or implied warranty.
*/

51/** @file bits/stl_iterator.h
*  This is an internal header file, included by other library headers.
*  Do not attempt to use it directly. @headername{iterator}
*
*  This file implements reverse_iterator, back_insert_iterator,
*  front_insert_iterator, insert_iterator, __normal_iterator, and their
*  supporting functions and overloaded operators.
*/

60#ifndef _STL_ITERATOR_H1
61#define _STL_ITERATOR_H1 1

63#include <bits/cpp_type_traits.h>
64#include <ext/type_traits.h>
65#include <bits/move.h>

67namespace std _GLIBCXX_VISIBILITY(default)__attribute__ ((__visibility__ ("default")))
68{
69_GLIBCXX_BEGIN_NAMESPACE_VERSION

/**
 * @addtogroup iterators
 * @{
 */

// 24.4.1 Reverse iterators
/**
 *  Bidirectional and random access iterators have corresponding reverse
 *  %iterator adaptors that iterate through the data structure in the
 *  opposite direction.  They have the same signatures as the corresponding
 *  iterators.  The fundamental relation between a reverse %iterator and its
 *  corresponding %iterator @c i is established by the identity:
 *  @code
 *      &*(reverse_iterator(i)) == &*(i - 1)
 *  @endcode
 *
 *  <em>This mapping is dictated by the fact that while there is always a
 *  pointer past the end of an array, there might not be a valid pointer
 *  before the beginning of an array.</em> [24.4.1]/1,2
 *
 *  Reverse iterators can be tricky and surprising at first.  Their
 *  semantics make sense, however, and the trickiness is a side effect of
 *  the requirement that the iterators must be safe.
*/
template<typename _Iterator>
  class reverse_iterator
  : public iterator<typename iterator_traits<_Iterator>::iterator_category,
      typename iterator_traits<_Iterator>::value_type,
      typename iterator_traits<_Iterator>::difference_type,
      typename iterator_traits<_Iterator>::pointer,
                    typename iterator_traits<_Iterator>::reference>
  {
  protected:
    _Iterator current;

    typedef iterator_traits<_Iterator>		__traits_type;

  public:
    typedef _Iterator					iterator_type;
    typedef typename __traits_type::difference_type	difference_type;
    typedef typename __traits_type::pointer		pointer;
    typedef typename __traits_type::reference		reference;

    /**
     *  The default constructor value-initializes member @p current.
     *  If it is a pointer, that means it is zero-initialized.
    */
    // _GLIBCXX_RESOLVE_LIB_DEFECTS
    // 235 No specification of default ctor for reverse_iterator
    reverse_iterator() : current() { }

    /**
     *  This %iterator will move in the opposite direction that @p x does.
    */
    explicit
    reverse_iterator(iterator_type __x) : current(__x) { }

    /**
     *  The copy constructor is normal.
    */
    reverse_iterator(const reverse_iterator& __x)
    : current(__x.current) { }

    /**
     *  A %reverse_iterator across other types can be copied if the
     *  underlying %iterator can be converted to the type of @c current.
    */
    template<typename _Iter>
      reverse_iterator(const reverse_iterator<_Iter>& __x)
: current(__x.base()) { }

    /**
     *  @return  @c current, the %iterator used for underlying work.
    */
    iterator_type
    base() const
    { return current; }

    /**
     *  @return  A reference to the value at @c --current
     *
     *  This requires that @c --current is dereferenceable.
     *
     *  @warning This implementation requires that for an iterator of the
     *           underlying iterator type, @c x, a reference obtained by
     *           @c *x remains valid after @c x has been modified or
     *           destroyed. This is a bug: http://gcc.gnu.org/PR51823
    */
    reference
    operator*() const
    {
_Iterator __tmp = current;
return *--__tmp;
    }

    /**
     *  @return  A pointer to the value at @c --current
     *
     *  This requires that @c --current is dereferenceable.
    */
    pointer
    operator->() const
    { return &(operator*()); }

    /**
     *  @return  @c *this
     *
     *  Decrements the underlying iterator.
    */
    reverse_iterator&
    operator++()
    {
--current;
return *this;
    }

    /**
     *  @return  The original value of @c *this
     *
     *  Decrements the underlying iterator.
    */
    reverse_iterator
    operator++(int)
    {
reverse_iterator __tmp = *this;
--current;
return __tmp;
    }

    /**
     *  @return  @c *this
     *
     *  Increments the underlying iterator.
    */
    reverse_iterator&
    operator--()
    {
++current;
return *this;
    }

    /**
     *  @return  A reverse_iterator with the previous value of @c *this
     *
     *  Increments the underlying iterator.
    */
    reverse_iterator
    operator--(int)
    {
reverse_iterator __tmp = *this;
++current;
return __tmp;
    }

    /**
     *  @return  A reverse_iterator that refers to @c current - @a __n
     *
     *  The underlying iterator must be a Random Access Iterator.
    */
    reverse_iterator
    operator+(difference_type __n) const
    { return reverse_iterator(current - __n); }

    /**
     *  @return  *this
     *
     *  Moves the underlying iterator backwards @a __n steps.
     *  The underlying iterator must be a Random Access Iterator.
    */
    reverse_iterator&
    operator+=(difference_type __n)
    {
current -= __n;
return *this;
    }

    /**
     *  @return  A reverse_iterator that refers to @c current - @a __n
     *
     *  The underlying iterator must be a Random Access Iterator.
    */
    reverse_iterator
    operator-(difference_type __n) const
    { return reverse_iterator(current + __n); }

    /**
     *  @return  *this
     *
     *  Moves the underlying iterator forwards @a __n steps.
     *  The underlying iterator must be a Random Access Iterator.
    */
    reverse_iterator&
    operator-=(difference_type __n)
    {
current += __n;
return *this;
    }

    /**
     *  @return  The value at @c current - @a __n - 1
     *
     *  The underlying iterator must be a Random Access Iterator.
    */
    reference
    operator[](difference_type __n) const
    { return *(*this + __n); }
  };

//@{
/**
 *  @param  __x  A %reverse_iterator.
 *  @param  __y  A %reverse_iterator.
 *  @return  A simple bool.
 *
 *  Reverse iterators forward many operations to their underlying base()
 *  iterators.  Others are implemented in terms of one another.
 *
*/
template<typename _Iterator>
  inline bool
  operator==(const reverse_iterator<_Iterator>& __x,
      const reverse_iterator<_Iterator>& __y)
  { return __x.base() == __y.base(); }

template<typename _Iterator>
  inline bool
  operator<(const reverse_iterator<_Iterator>& __x,
     const reverse_iterator<_Iterator>& __y)
  { return __y.base() < __x.base(); }

template<typename _Iterator>
  inline bool
  operator!=(const reverse_iterator<_Iterator>& __x,
      const reverse_iterator<_Iterator>& __y)
  { return !(__x == __y); }

template<typename _Iterator>
  inline bool
  operator>(const reverse_iterator<_Iterator>& __x,
     const reverse_iterator<_Iterator>& __y)
  { return __y < __x; }

template<typename _Iterator>
  inline bool
  operator<=(const reverse_iterator<_Iterator>& __x,
      const reverse_iterator<_Iterator>& __y)
  { return !(__y < __x); }

template<typename _Iterator>
  inline bool
  operator>=(const reverse_iterator<_Iterator>& __x,
      const reverse_iterator<_Iterator>& __y)
  { return !(__x < __y); }

template<typename _Iterator>
  inline typename reverse_iterator<_Iterator>::difference_type
  operator-(const reverse_iterator<_Iterator>& __x,
     const reverse_iterator<_Iterator>& __y)
  { return __y.base() - __x.base(); }

template<typename _Iterator>
  inline reverse_iterator<_Iterator>
  operator+(typename reverse_iterator<_Iterator>::difference_type __n,
     const reverse_iterator<_Iterator>& __x)
  { return reverse_iterator<_Iterator>(__x.base() - __n); }

// _GLIBCXX_RESOLVE_LIB_DEFECTS
// DR 280. Comparison of reverse_iterator to const reverse_iterator.
template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator==(const reverse_iterator<_IteratorL>& __x,
      const reverse_iterator<_IteratorR>& __y)
  { return __x.base() == __y.base(); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator<(const reverse_iterator<_IteratorL>& __x,
     const reverse_iterator<_IteratorR>& __y)
  { return __y.base() < __x.base(); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator!=(const reverse_iterator<_IteratorL>& __x,
      const reverse_iterator<_IteratorR>& __y)
  { return !(__x == __y); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator>(const reverse_iterator<_IteratorL>& __x,
     const reverse_iterator<_IteratorR>& __y)
  { return __y < __x; }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator<=(const reverse_iterator<_IteratorL>& __x,
      const reverse_iterator<_IteratorR>& __y)
  { return !(__y < __x); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator>=(const reverse_iterator<_IteratorL>& __x,
      const reverse_iterator<_IteratorR>& __y)
  { return !(__x < __y); }

template<typename _IteratorL, typename _IteratorR>
376#if __cplusplus201103L >= 201103L
  // DR 685.
  inline auto
  operator-(const reverse_iterator<_IteratorL>& __x,
     const reverse_iterator<_IteratorR>& __y)
  -> decltype(__y.base() - __x.base())
382#else
  inline typename reverse_iterator<_IteratorL>::difference_type
  operator-(const reverse_iterator<_IteratorL>& __x,
     const reverse_iterator<_IteratorR>& __y)
386#endif
  { return __y.base() - __x.base(); }
//@}

// 24.4.2.2.1 back_insert_iterator
/**
 *  @brief  Turns assignment into insertion.
 *
 *  These are output iterators, constructed from a container-of-T.
 *  Assigning a T to the iterator appends it to the container using
 *  push_back.
 *
 *  Tip:  Using the back_inserter function to create these iterators can
 *  save typing.
*/
template<typename _Container>
  class back_insert_iterator
  : public iterator<output_iterator_tag, void, void, void, void>
  {
  protected:
    _Container* container;

  public:
    /// A nested typedef for the type of whatever container you used.
    typedef _Container          container_type;

    /// The only way to create this %iterator is with a container.
    explicit
    back_insert_iterator(_Container& __x) : container(&__x) { }

    /**
     *  @param  __value  An instance of whatever type
     *                 container_type::const_reference is; presumably a
     *                 reference-to-const T for container<T>.
     *  @return  This %iterator, for chained operations.
     *
     *  This kind of %iterator doesn't really have a @a position in the
     *  container (you can think of the position as being permanently at
     *  the end, if you like).  Assigning a value to the %iterator will
     *  always append the value to the end of the container.
    */
427#if __cplusplus201103L < 201103L
    back_insert_iterator&
    operator=(typename _Container::const_reference __value)
    {
container->push_back(__value);
return *this;
    }
434#else
    back_insert_iterator&
    operator=(const typename _Container::value_type& __value)
    {
container->push_back(__value);
return *this;
    }

    back_insert_iterator&
    operator=(typename _Container::value_type&& __value)
    {
container->push_back(std::move(__value));
return *this;
    }
448#endif

    /// Simply returns *this.
    back_insert_iterator&
    operator*()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    back_insert_iterator&
    operator++()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    back_insert_iterator
    operator++(int)
    { return *this; }
  };

/**
 *  @param  __x  A container of arbitrary type.
 *  @return  An instance of back_insert_iterator working on @p __x.
 *
 *  This wrapper function helps in creating back_insert_iterator instances.
 *  Typing the name of the %iterator requires knowing the precise full
 *  type of the container, which can be tedious and impedes generic
 *  programming.  Using this function lets you take advantage of automatic
 *  template parameter deduction, making the compiler match the correct
 *  types for you.
*/
template<typename _Container>
  inline back_insert_iterator<_Container>
  back_inserter(_Container& __x)
  { return back_insert_iterator<_Container>(__x); }

/**
 *  @brief  Turns assignment into insertion.
 *
 *  These are output iterators, constructed from a container-of-T.
 *  Assigning a T to the iterator prepends it to the container using
 *  push_front.
 *
 *  Tip:  Using the front_inserter function to create these iterators can
 *  save typing.
*/
template<typename _Container>
  class front_insert_iterator
  : public iterator<output_iterator_tag, void, void, void, void>
  {
  protected:
    _Container* container;

  public:
    /// A nested typedef for the type of whatever container you used.
    typedef _Container          container_type;

    /// The only way to create this %iterator is with a container.
    explicit front_insert_iterator(_Container& __x) : container(&__x) { }

    /**
     *  @param  __value  An instance of whatever type
     *                 container_type::const_reference is; presumably a
     *                 reference-to-const T for container<T>.
     *  @return  This %iterator, for chained operations.
     *
     *  This kind of %iterator doesn't really have a @a position in the
     *  container (you can think of the position as being permanently at
     *  the front, if you like).  Assigning a value to the %iterator will
     *  always prepend the value to the front of the container.
    */
517#if __cplusplus201103L < 201103L
    front_insert_iterator&
    operator=(typename _Container::const_reference __value)
    {
container->push_front(__value);
return *this;
    }
524#else
    front_insert_iterator&
    operator=(const typename _Container::value_type& __value)
    {
container->push_front(__value);
return *this;
    }

    front_insert_iterator&
    operator=(typename _Container::value_type&& __value)
    {
container->push_front(std::move(__value));
return *this;
    }
538#endif

    /// Simply returns *this.
    front_insert_iterator&
    operator*()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    front_insert_iterator&
    operator++()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    front_insert_iterator
    operator++(int)
    { return *this; }
  };

/**
 *  @param  __x  A container of arbitrary type.
 *  @return  An instance of front_insert_iterator working on @p x.
 *
 *  This wrapper function helps in creating front_insert_iterator instances.
 *  Typing the name of the %iterator requires knowing the precise full
 *  type of the container, which can be tedious and impedes generic
 *  programming.  Using this function lets you take advantage of automatic
 *  template parameter deduction, making the compiler match the correct
 *  types for you.
*/
template<typename _Container>
  inline front_insert_iterator<_Container>
  front_inserter(_Container& __x)
  { return front_insert_iterator<_Container>(__x); }

/**
 *  @brief  Turns assignment into insertion.
 *
 *  These are output iterators, constructed from a container-of-T.
 *  Assigning a T to the iterator inserts it in the container at the
 *  %iterator's position, rather than overwriting the value at that
 *  position.
 *
 *  (Sequences will actually insert a @e copy of the value before the
 *  %iterator's position.)
 *
 *  Tip:  Using the inserter function to create these iterators can
 *  save typing.
*/
template<typename _Container>
  class insert_iterator
  : public iterator<output_iterator_tag, void, void, void, void>
  {
  protected:
    _Container* container;
    typename _Container::iterator iter;

  public:
    /// A nested typedef for the type of whatever container you used.
    typedef _Container          container_type;

    /**
     *  The only way to create this %iterator is with a container and an
     *  initial position (a normal %iterator into the container).
    */
    insert_iterator(_Container& __x, typename _Container::iterator __i)
    : container(&__x), iter(__i) {}

    /**
     *  @param  __value  An instance of whatever type
     *                 container_type::const_reference is; presumably a
     *                 reference-to-const T for container<T>.
     *  @return  This %iterator, for chained operations.
     *
     *  This kind of %iterator maintains its own position in the
     *  container.  Assigning a value to the %iterator will insert the
     *  value into the container at the place before the %iterator.
     *
     *  The position is maintained such that subsequent assignments will
     *  insert values immediately after one another.  For example,
     *  @code
     *     // vector v contains A and Z
     *
     *     insert_iterator i (v, ++v.begin());
     *     i = 1;
     *     i = 2;
     *     i = 3;
     *
     *     // vector v contains A, 1, 2, 3, and Z
     *  @endcode
    */
628#if __cplusplus201103L < 201103L
    insert_iterator&
    operator=(typename _Container::const_reference __value)
    {
iter = container->insert(iter, __value);
++iter;
return *this;
    }
636#else
    insert_iterator&
    operator=(const typename _Container::value_type& __value)
    {
iter = container->insert(iter, __value);
++iter;
return *this;
    }

    insert_iterator&
    operator=(typename _Container::value_type&& __value)
    {
iter = container->insert(iter, std::move(__value));
++iter;
return *this;
    }
652#endif

    /// Simply returns *this.
    insert_iterator&
    operator*()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    insert_iterator&
    operator++()
    { return *this; }

    /// Simply returns *this.  (This %iterator does not @a move.)
    insert_iterator&
    operator++(int)
    { return *this; }
  };

/**
 *  @param __x  A container of arbitrary type.
 *  @return  An instance of insert_iterator working on @p __x.
 *
 *  This wrapper function helps in creating insert_iterator instances.
 *  Typing the name of the %iterator requires knowing the precise full
 *  type of the container, which can be tedious and impedes generic
 *  programming.  Using this function lets you take advantage of automatic
 *  template parameter deduction, making the compiler match the correct
 *  types for you.
*/
template<typename _Container, typename _Iterator>
  inline insert_iterator<_Container>
  inserter(_Container& __x, _Iterator __i)
  {
    return insert_iterator<_Container>(__x,
			 typename _Container::iterator(__i));
  }

// @} group iterators

691_GLIBCXX_END_NAMESPACE_VERSION
692} // namespace

694namespace __gnu_cxx _GLIBCXX_VISIBILITY(default)__attribute__ ((__visibility__ ("default")))
695{
696_GLIBCXX_BEGIN_NAMESPACE_VERSION

// This iterator adapter is @a normal in the sense that it does not
// change the semantics of any of the operators of its iterator
// parameter.  Its primary purpose is to convert an iterator that is
// not a class, e.g. a pointer, into an iterator that is a class.
// The _Container parameter exists solely so that different containers
// using this template can instantiate different types, even if the
// _Iterator parameter is the same.
using std::iterator_traits;
using std::iterator;
template<typename _Iterator, typename _Container>
  class __normal_iterator
  {
  protected:
    _Iterator _M_current;

    typedef iterator_traits<_Iterator>		__traits_type;

  public:
    typedef _Iterator					iterator_type;
    typedef typename __traits_type::iterator_category iterator_category;
    typedef typename __traits_type::value_type  	value_type;
    typedef typename __traits_type::difference_type 	difference_type;
    typedef typename __traits_type::reference 	reference;
    typedef typename __traits_type::pointer   	pointer;

    _GLIBCXX_CONSTEXPRconstexpr __normal_iterator() : _M_current(_Iterator()) { }

    explicit
    __normal_iterator(const _Iterator& __i) : _M_current(__i) { }

    // Allow iterator to const_iterator conversion
    template<typename _Iter>
      __normal_iterator(const __normal_iterator<_Iter,
	  typename __enable_if<
    	       (std::__are_same<_Iter, typename _Container::pointer>::__value),
      _Container>::__type>& __i)
      : _M_current(__i.base()) { }

    // Forward iterator requirements
    reference
    operator*() const
    { return *_M_current; }

    pointer
    operator->() const
    { return _M_current; }

    __normal_iterator&
    operator++()
    {
++_M_current;
return *this;
    }

    __normal_iterator
    operator++(int)
    { return __normal_iterator(_M_current++); }

    // Bidirectional iterator requirements
    __normal_iterator&
    operator--()
    {
--_M_current;
return *this;
    }

    __normal_iterator
    operator--(int)
    { return __normal_iterator(_M_current--); }

    // Random access iterator requirements
    reference
    operator[](const difference_type& __n) const
    { return _M_current[__n]; }

    __normal_iterator&
    operator+=(const difference_type& __n)
    { _M_current += __n; return *this; }

    __normal_iterator
    operator+(const difference_type& __n) const
    { return __normal_iterator(_M_current + __n); }

    __normal_iterator&
    operator-=(const difference_type& __n)
    { _M_current -= __n; return *this; }

    __normal_iterator
    operator-(const difference_type& __n) const
    { return __normal_iterator(_M_current - __n); }

    const _Iterator&
    base() const
    { return _M_current; }
  };

// Note: In what follows, the left- and right-hand-side iterators are
// allowed to vary in types (conceptually in cv-qualification) so that
// comparison between cv-qualified and non-cv-qualified iterators be
// valid.  However, the greedy and unfriendly operators in std::rel_ops
// will make overload resolution ambiguous (when in scope) if we don't
// provide overloads whose operands are of the same type.  Can someone
// remind me what generic programming is about? -- Gaby

// Forward iterator requirements
template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator==(const __normal_iterator<_IteratorL, _Container>& __lhs,
      const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() == __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator==(const __normal_iterator<_Iterator, _Container>& __lhs,
      const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() == __rhs.base(); }

template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator!=(const __normal_iterator<_IteratorL, _Container>& __lhs,
      const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() != __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator!=(const __normal_iterator<_Iterator, _Container>& __lhs,
      const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() != __rhs.base(); }
18
←
Assuming the condition is true→
19
←
Returning the value 1, which participates in a condition later→

// Random access iterator requirements
template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator<(const __normal_iterator<_IteratorL, _Container>& __lhs,
     const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() < __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator<(const __normal_iterator<_Iterator, _Container>& __lhs,
     const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() < __rhs.base(); }

template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator>(const __normal_iterator<_IteratorL, _Container>& __lhs,
     const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() > __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator>(const __normal_iterator<_Iterator, _Container>& __lhs,
     const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() > __rhs.base(); }

template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator<=(const __normal_iterator<_IteratorL, _Container>& __lhs,
      const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() <= __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator<=(const __normal_iterator<_Iterator, _Container>& __lhs,
      const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() <= __rhs.base(); }

template<typename _IteratorL, typename _IteratorR, typename _Container>
  inline bool
  operator>=(const __normal_iterator<_IteratorL, _Container>& __lhs,
      const __normal_iterator<_IteratorR, _Container>& __rhs)
  { return __lhs.base() >= __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline bool
  operator>=(const __normal_iterator<_Iterator, _Container>& __lhs,
      const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() >= __rhs.base(); }

// _GLIBCXX_RESOLVE_LIB_DEFECTS
// According to the resolution of DR179 not only the various comparison
// operators but also operator- must accept mixed iterator/const_iterator
// parameters.
template<typename _IteratorL, typename _IteratorR, typename _Container>
881#if __cplusplus201103L >= 201103L
  // DR 685.
  inline auto
  operator-(const __normal_iterator<_IteratorL, _Container>& __lhs,
     const __normal_iterator<_IteratorR, _Container>& __rhs)
  -> decltype(__lhs.base() - __rhs.base())
887#else
  inline typename __normal_iterator<_IteratorL, _Container>::difference_type
  operator-(const __normal_iterator<_IteratorL, _Container>& __lhs,
     const __normal_iterator<_IteratorR, _Container>& __rhs)
891#endif
  { return __lhs.base() - __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline typename __normal_iterator<_Iterator, _Container>::difference_type
  operator-(const __normal_iterator<_Iterator, _Container>& __lhs,
     const __normal_iterator<_Iterator, _Container>& __rhs)
  { return __lhs.base() - __rhs.base(); }

template<typename _Iterator, typename _Container>
  inline __normal_iterator<_Iterator, _Container>
  operator+(typename __normal_iterator<_Iterator, _Container>::difference_type
     __n, const __normal_iterator<_Iterator, _Container>& __i)
  { return __normal_iterator<_Iterator, _Container>(__i.base() + __n); }

906_GLIBCXX_END_NAMESPACE_VERSION
907} // namespace

909#if __cplusplus201103L >= 201103L

911namespace std _GLIBCXX_VISIBILITY(default)__attribute__ ((__visibility__ ("default")))
912{
913_GLIBCXX_BEGIN_NAMESPACE_VERSION

/**
 * @addtogroup iterators
 * @{
 */

// 24.4.3  Move iterators
/**
 *  Class template move_iterator is an iterator adapter with the same
 *  behavior as the underlying iterator except that its dereference
 *  operator implicitly converts the value returned by the underlying
 *  iterator's dereference operator to an rvalue reference.  Some
 *  generic algorithms can be called with move iterators to replace
 *  copying with moving.
 */
template<typename _Iterator>
  class move_iterator
  {
  protected:
    _Iterator _M_current;

    typedef iterator_traits<_Iterator>		__traits_type;

  public:
    typedef _Iterator					iterator_type;
    typedef typename __traits_type::iterator_category iterator_category;
    typedef typename __traits_type::value_type  	value_type;
    typedef typename __traits_type::difference_type	difference_type;
    // NB: DR 680.
    typedef _Iterator					pointer;
    typedef value_type&&				reference;

    move_iterator()
    : _M_current() { }

    explicit
    move_iterator(iterator_type __i)
    : _M_current(__i) { }

    template<typename _Iter>
move_iterator(const move_iterator<_Iter>& __i)
: _M_current(__i.base()) { }

    iterator_type
    base() const
    { return _M_current; }

    reference
    operator*() const
    { return std::move(*_M_current); }

    pointer
    operator->() const
    { return _M_current; }

    move_iterator&
    operator++()
    {
++_M_current;
return *this;
    }

    move_iterator
    operator++(int)
    {
move_iterator __tmp = *this;
++_M_current;
return __tmp;
    }

    move_iterator&
    operator--()
    {
--_M_current;
return *this;
    }

    move_iterator
    operator--(int)
    {
move_iterator __tmp = *this;
--_M_current;
return __tmp;
    }

    move_iterator
    operator+(difference_type __n) const
    { return move_iterator(_M_current + __n); }

    move_iterator&
    operator+=(difference_type __n)
    {
_M_current += __n;
return *this;
    }

    move_iterator
    operator-(difference_type __n) const
    { return move_iterator(_M_current - __n); }
  
    move_iterator&
    operator-=(difference_type __n)
    { 
_M_current -= __n;
return *this;
    }

    reference
    operator[](difference_type __n) const
    { return std::move(_M_current[__n]); }
  };

// Note: See __normal_iterator operators note from Gaby to understand
// why there are always 2 versions for most of the move_iterator
// operators.
template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator==(const move_iterator<_IteratorL>& __x,
      const move_iterator<_IteratorR>& __y)
  { return __x.base() == __y.base(); }

template<typename _Iterator>
  inline bool
  operator==(const move_iterator<_Iterator>& __x,
      const move_iterator<_Iterator>& __y)
  { return __x.base() == __y.base(); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator!=(const move_iterator<_IteratorL>& __x,
      const move_iterator<_IteratorR>& __y)
  { return !(__x == __y); }

template<typename _Iterator>
  inline bool
  operator!=(const move_iterator<_Iterator>& __x,
      const move_iterator<_Iterator>& __y)
  { return !(__x == __y); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator<(const move_iterator<_IteratorL>& __x,
     const move_iterator<_IteratorR>& __y)
  { return __x.base() < __y.base(); }

template<typename _Iterator>
  inline bool
  operator<(const move_iterator<_Iterator>& __x,
     const move_iterator<_Iterator>& __y)
  { return __x.base() < __y.base(); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator<=(const move_iterator<_IteratorL>& __x,
      const move_iterator<_IteratorR>& __y)
  { return !(__y < __x); }

template<typename _Iterator>
  inline bool
  operator<=(const move_iterator<_Iterator>& __x,
      const move_iterator<_Iterator>& __y)
  { return !(__y < __x); }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator>(const move_iterator<_IteratorL>& __x,
     const move_iterator<_IteratorR>& __y)
  { return __y < __x; }

template<typename _Iterator>
  inline bool
  operator>(const move_iterator<_Iterator>& __x,
     const move_iterator<_Iterator>& __y)
  { return __y < __x; }

template<typename _IteratorL, typename _IteratorR>
  inline bool
  operator>=(const move_iterator<_IteratorL>& __x,
      const move_iterator<_IteratorR>& __y)
  { return !(__x < __y); }

template<typename _Iterator>
  inline bool
  operator>=(const move_iterator<_Iterator>& __x,
      const move_iterator<_Iterator>& __y)
  { return !(__x < __y); }

// DR 685.
template<typename _IteratorL, typename _IteratorR>
  inline auto
  operator-(const move_iterator<_IteratorL>& __x,
     const move_iterator<_IteratorR>& __y)
  -> decltype(__x.base() - __y.base())
  { return __x.base() - __y.base(); }

template<typename _Iterator>
  inline auto
  operator-(const move_iterator<_Iterator>& __x,
     const move_iterator<_Iterator>& __y)
  -> decltype(__x.base() - __y.base())
  { return __x.base() - __y.base(); }

template<typename _Iterator>
  inline move_iterator<_Iterator>
  operator+(typename move_iterator<_Iterator>::difference_type __n,
     const move_iterator<_Iterator>& __x)
  { return __x + __n; }

template<typename _Iterator>
  inline move_iterator<_Iterator>
  make_move_iterator(_Iterator __i)
  { return move_iterator<_Iterator>(__i); }

template<typename _Iterator, typename _ReturnType
  = typename conditional<__move_if_noexcept_cond
    <typename iterator_traits<_Iterator>::value_type>::value,
              _Iterator, move_iterator<_Iterator>>::type>
  inline _ReturnType
  __make_move_if_noexcept_iterator(_Iterator __i)
  { return _ReturnType(__i); }

// @} group iterators

1137_GLIBCXX_END_NAMESPACE_VERSION
1138} // namespace

1140#define _GLIBCXX_MAKE_MOVE_ITERATOR(_Iter)std::make_move_iterator(_Iter) std::make_move_iterator(_Iter)
1141#define _GLIBCXX_MAKE_MOVE_IF_NOEXCEPT_ITERATOR(_Iter)std::__make_move_if_noexcept_iterator(_Iter) \
std::__make_move_if_noexcept_iterator(_Iter)
1143#else
1144#define _GLIBCXX_MAKE_MOVE_ITERATOR(_Iter)std::make_move_iterator(_Iter) (_Iter)
1145#define _GLIBCXX_MAKE_MOVE_IF_NOEXCEPT_ITERATOR(_Iter)std::__make_move_if_noexcept_iterator(_Iter) (_Iter)
1146#endif // C++11

1148#endif