StEEmcClusterMaker.cxx

#include "StEEmcClusterMaker.h"

#include <algorithm>
#include <iostream>

#include "TMath.h"

/* StEvent stuff */
#include "StEvent/StEvent.h"
#include "StEvent/StEmcCollection.h"
#include "StEvent/StEmcDetector.h"
#include "StEvent/StEmcModule.h"
#include "StEvent/StEmcClusterCollection.h"
#include "StEvent/StEmcCluster.h"

#define DEBUG 0

ClassImp(StEEmcClusterMaker);

// ----------------------------------------------------------------------------
StEEmcClusterMaker::StEEmcClusterMaker(const Char_t *name):StMaker(name)
{

  mFillStEvent = 0;

  mSuppress = 0;

  mLoose=false;
  mSkip=true;

  mClusterId = 0;

  const Float_t eseeds[] = { 0.6, 1.0/1000, 1.0/1000, 1.0/1000., 0.3/1000., 0.3/1000. };
  for ( Int_t i = 0; i < 6; i++ ) seedEnergy( eseeds[i], i );

  mMaxExtent = 6; 
  mSeedFloor = 2.0; 


  StEEmcClusterVec_t t;
  std::vector< StEEmcClusterVec_t > layers;
  for ( Int_t i = 0; i < 4; i++ ) layers.push_back(t);
  for ( Int_t i = 0; i < 12; i++ ) mTowerClusters.push_back(layers);

  StEEmcSmdClusterVec_t s;
  std::vector< StEEmcSmdClusterVec_t > planes;
  planes.push_back(s);
  planes.push_back(s);
  for ( Int_t i = 0; i < 12; i++ ) mSmdClusters.push_back(planes);

  StEEmcTowerVec_t tow;
  std::vector< StEEmcTowerVec_t > lay;
  for ( Int_t i = 0; i < 4; i++ ) lay.push_back(tow);
  for ( Int_t i = 0; i < 12; i++ ) mSeedTowers.push_back(lay);

  mEEtow=new EEmcGeomSimple();
  mEEsmd=EEmcSmdGeom::instance();
  mEEmap=EEmcSmdMap::instance();

  mMinStrips=3;
  
}

// ----------------------------------------------------------------------------
Int_t StEEmcClusterMaker::Init()
{
  mEEanalysis=(const StEEmcA2EMaker*)GetMaker(mAnalysisName);
  assert(mEEanalysis);

  return StMaker::Init();
}

// ----------------------------------------------------------------------------
Int_t StEEmcClusterMaker::Make()
{
  StMaker::Make();
  if ( !buildTowerClusters() ) return kStWarn;

  if ( !buildSmdClusters() ) return kStWarn; 

  if ( mFillStEvent ) fillStEvent();

  if ( mFillStEvent ) 
    if ( !verifyStEvent() ) Warning("Make","StEvent not properly copied");

  return kStOK;
}

// ----------------------------------------------------------------------------
void StEEmcClusterMaker::Clear( Option_t *opts )
{
  StMaker::Clear();
  for ( Int_t sector=0; sector<12; sector++ ) {
    for ( Int_t layer=0; layer<4; layer++ )
      mTowerClusters[sector][layer].clear();    
    for ( Int_t plane=0; plane<2; plane++ )
      mSmdClusters[sector][plane].clear();
  }

  for ( Int_t i=0;i<6;i++ ) mNumberOfClusters[i]=0;

  mClusterId = 0;

  return;
}


// ----------------------------------------------------------------------------
Bool_t StEEmcClusterMaker::buildTowerClusters()
{

  //static const Int_t layer=0;


  for ( Int_t layer=0;layer<4;layer++)
  {

    Float_t weights[720]; for ( Int_t i=0;i<720;i++ ) weights[i]=0.;

    StEEmcClusterVec_t myClusters;


    StEEmcTowerVec_t towers = mEEanalysis -> towers( layer );



    std::sort(towers.begin(),towers.end());
    std::reverse(towers.begin(),towers.end());
   

    StEEmcTowerVec_t::iterator last = towers.begin();
    while ( last != towers.end() ) {

      if ( (*last).energy() < mSeedEnergy[layer] ) break;
      //$$$      mSeedTowers[layer][ (*last).sector() ].push_back((*last));
#if DEBUG
      std::cout << "-- Seed tower ----------------------" << std::endl;
      (*last).print();
#endif
      last++;
    }
    
    StEEmcTowerVec_t::iterator iter = towers.begin();
    while ( iter != towers.end() ) {

      if ( iter==last ) break;
      

      if ( weights[ (*iter).index() ] > 0. ) {
        iter++;
        continue;
      }

      for ( Int_t in=0; in < (*iter).numberOfNeighbors(); in++ ) {
        StEEmcTower t=(*iter).neighbor(in);
        weights[ t.index() ] += (*iter).energy();
      }

     
      iter++;

    }// loop over towers to init weights


    iter=towers.begin();
    while ( iter != towers.end() ) {

      if ( iter==last ) break;

      if ( weights[ (*iter).index() ] > 0. ) {
        iter++;
        continue;
      }

      StEEmcCluster cluster;
      StEEmcTower seed=(*iter);
#if DEBUG
      std::cout << "--- Clustering ----------------" << std::endl;
      seed.print();
#endif

      TVector3 momentum;
      cluster.add(seed,1.0);
      UInt_t sec,sub,eta;
      sec=(UInt_t)seed.sector(); 
      sub=(UInt_t)seed.subsector(); 
      eta=(UInt_t)seed.etabin();
      TVector3 d=mEEtow->getTowerCenter(sec,sub,eta).Unit();
      momentum += ( seed.energy() * d );

      for ( Int_t in=0; in<seed.numberOfNeighbors(); in++ ) {
        StEEmcTower t=seed.neighbor(in);
        sec=(UInt_t)t.sector(); 
        sub=(UInt_t)t.subsector(); 
        eta=(UInt_t)t.etabin();
        d=mEEtow->getTowerCenter(sec,sub,eta).Unit();
        Float_t weight = seed.energy() / weights[ t.index() ];
        momentum += ( t.energy() * d * weight );
        cluster.add(t, weight);
#if DEBUG
        std::cout << "adding " << t.name() << " E=" << t.energy() << " W=" << weight << std::endl;
#endif
      }


      cluster.momentum( momentum );
      cluster.key( mClusterId++ );
#if DEBUG
      cluster.print();
#endif

      mTowerClusters[ seed.sector() ][ layer ].push_back( cluster );
      mNumberOfClusters[layer]++;
     
      iter++;

    }
          
  }// end of eventual loop over layers
  

  return true;

}

// ----------------------------------------------------------------------------
Bool_t StEEmcClusterMaker::buildSmdClusters()
{

  Int_t max_extent = mMaxExtent;

  for ( Int_t sector=0; sector<12; sector++ ) 
    for ( Int_t plane=0; plane<2; plane++ ) {

        Float_t floor[288]; for ( Int_t i=0; i<288; i++ ) floor[i]=0.;

        StEEmcStripVec_t strips=mEEanalysis->strips(sector,plane);
        StEEmcStripVec_t seeds;
        std::sort(strips.begin(),strips.end());
        std::reverse(strips.begin(),strips.end());
        
        std::vector<Bool_t> seed_flags( strips.size(), false );


        Int_t nstrip=0;
        Int_t nseeds=0;
        StEEmcStripVec_t::iterator istrip=strips.begin();
        while ( istrip!=strips.end() ) {

          Int_t index=(*istrip).index();
          Float_t eseed=(*istrip).energy();

          if ( index <= 3 || index >= 283 ) {
            istrip++;
            nstrip++;
            continue;
          }

          if ( mSkip )
            if ( (*istrip).fail() ) {
              istrip++;
              nstrip++;
              continue;
            }


          if ( eseed < mSeedFloor*floor[ index ] + mSeedEnergy[4+plane] ) {
            istrip++;
            nstrip++;
            continue;
          }

          seeds.push_back( (*istrip) );
          nstrip++;
          nseeds++;


#ifndef LOOSE_CUTS        
          for ( Int_t i=0; i < 288; i++ ) {
            Int_t dx=TMath::Abs(index-i);
            if ( dx<3 ) 
              if ( eseed > floor[i] ) floor[i]=eseed;       
            if ( dx<5 ) 
              if ( 0.20 * eseed > floor[i] ) floor[i] = 0.2 * eseed;
            if ( dx<11 )
              if ( 0.10 * eseed > floor[i] ) floor[i] = 0.1 * eseed;
            if ( dx<21 )
              if ( 0.20*eseed > floor[i] ) floor[i] = 0.05 * eseed;                         
          }
#else
          for ( Int_t i=0; i < 288; i++ ) {

            Int_t dx=TMath::Abs(index-i);

            if ( dx<7 ) 
              if ( 0.05 * eseed > floor[i] ) floor[i] = 0.05 * eseed;

          }
#endif
          istrip++;

        }// found all seed strips

        Bool_t owned[288]; for (Int_t i=0;i<288;i++) owned[i]=false;
        Bool_t xseed[288]; for (Int_t i=0;i<288;i++) xseed[i]=false;

        StEEmcStripVec_t::iterator iseed=seeds.begin();
        while ( iseed != seeds.end() ) {
          
          Int_t index=(*iseed).index();
          if ( owned[index] || (mSuppress&&xseed[index]) ) {
            iseed++;
            continue;
          }

          owned[index]=true;
          xseed[index]=true;
          StEEmcSmdCluster cluster;
          cluster.add( (*iseed) );

          Int_t ind_max = TMath::Min(287,index+max_extent); 
          Int_t ind_min = TMath::Max(0,  index-max_extent); 

          for ( Int_t i=index+1; i<=ind_max; i++ ) {
            StEEmcStrip strip=mEEanalysis->strip(sector,plane,i);

            if ( mSkip ) {
              if (strip.energy()<=0.&&!strip.fail()) break; 
            } else {
              if ( strip.energy()<=0. ) break;
            }

            owned[ strip.index() ] = true;
            xseed[ strip.index() ] = true;
            cluster.add(strip);             
          }
          for ( Int_t i=index-1; i>=ind_min; i-- ) {
            StEEmcStrip strip=mEEanalysis->strip(sector,plane,i);

            if ( mSkip ) {
              if (strip.energy()<=0.&&!strip.fail()) break;
            } else {
              if (strip.energy()<=0.) break;
            }

            owned[ strip.index() ] = true;
            xseed[ strip.index() ] = true;
            cluster.add(strip);             
          }
          
          if ( cluster.size() >= mMinStrips ) {
            cluster.key( mClusterId++ );
            mSmdClusters[ sector ][ plane ].push_back(cluster);
            mNumberOfClusters[4+plane]++;

            // disallow strips on either side of the cluster
            // from forming seeds
            Int_t ns=cluster.numberOfStrips();
            Int_t left=999,right=-999;
            for ( Int_t is=0;is<ns;is++ )
            {
                StEEmcStrip s=cluster.strip(is);
                if ( s.index()<left ) left=s.index();
                if ( s.index()>right ) right=s.index();
            }
            /*
            xseed[left-1]=true;
            xseed[left-2]=true;
            xseed[right+1]=true;
            xseed[right+2]=true;
            */
            for ( Int_t ii=0;ii<=mSuppress;ii++ )
              {
                if ( left-ii>=0   ) xseed[left-ii] = true;
                if ( right+ii<288 ) xseed[right+ii]= true;
              }
            
          }
          

          iseed++;
        }                       

    }// loop over planes/sectors

  return true;

}



// ----------------------------------------------------------------------------

void StEEmcClusterMaker::fillStEvent()
{

  StEvent *stevent=(StEvent*)GetInputDS("StEvent");
  if ( !stevent ) {
    Warning("fillStEvent","called, but no StEvent to be found");
    return;
  }

  std::cout << "Adding tower clusters to StEvent at " << stevent << std::endl;

  StEmcDetector *detector=stevent->emcCollection()->detector(kEndcapEmcTowerId);
  if ( !detector )
    {
      Warning("fillStEvent","detector == NULL, MAJOR StEvent problem, continuing");
      return;
    }
  

  
  if ( mNumberOfClusters[0] > 0 )
    {


      StEmcClusterCollection *collect = detector -> cluster();
      if ( !collect ) 
        {
          //Warning("fillStEvent","StEmcClusterCollection (towers) was NULL, so I'm creating one.");
          collect = new StEmcClusterCollection();
          detector->setCluster( collect );
        }

      assert(collect);
      collect->setDetector( kEndcapEmcTowerId );
      collect->setClusterFinderId( 123 );
      collect->setClusterFinderParamVersion( 123 );
  
      for ( Int_t isector=0; isector<12; isector++ )
        {

          for ( UInt_t iclust=0; iclust<mTowerClusters[isector][0].size(); iclust++ ) 
            {

              StEEmcCluster cl=(mTowerClusters[isector][0])[iclust]; 
              
              StEmcCluster *emccluster = new StEmcCluster();
              emccluster->setEta( cl.momentum().Eta() );
              emccluster->setPhi( cl.momentum().Phi() );
              emccluster->setSigmaEta(-1.);
              emccluster->setSigmaPhi(-1.);
              emccluster->setEnergy( cl.energy() );
              emccluster->SetUniqueID( cl.key() );
#if 1
              for ( Int_t i=0; i< cl.numberOfTowers(); i++ ) 
                {
                  StEmcRawHit *hit=cl.tower(i).stemc();
                  assert( hit );         
                  emccluster->addHit( hit );
                }
#endif
            
              collect->addCluster( emccluster );

              
              mEtoEE[ emccluster ] = cl;
              cl.stemc( emccluster );

            }
          
        }

    }

  else // if ( mNumberOfClusters[0] == 0 )
    {

      detector->setCluster( NULL );

    }




  //  detector->setCluster( collect );

  detector=stevent->emcCollection()->detector(kEndcapEmcPreShowerId);
  if ( !detector )
    {
      Warning("fillStEvent","detector == NULL for pre/post, no clusters for you");
    }
  else if ( mNumberOfClusters[1] > 0 ||
            mNumberOfClusters[2] > 0 || 
            mNumberOfClusters[3] > 0 )
    {

      StEmcClusterCollection *pqr = detector -> cluster();
      if ( !pqr )
        {
          //Warning("fillStEvent","StEmcClusterCollection (pre/post) was NULL, so I'm creating one.");
          pqr = new StEmcClusterCollection();
          detector->setCluster( pqr );
        }
      assert(pqr); 
      pqr -> setDetector( kEndcapEmcPreShowerId );
      pqr -> setClusterFinderId( 123 );
      pqr -> setClusterFinderParamVersion( 321 );
  
      for ( Int_t isector=0; isector<12; isector++ )
    
      for ( Int_t ilayer=1; ilayer<4; ilayer++ ) 
    
      for ( UInt_t iclust=0; iclust<mTowerClusters[isector][ilayer].size(); iclust++ ) 
        {

          StEEmcCluster cl=(mTowerClusters[isector][ilayer])[iclust];

          StEmcCluster *emccluster = new StEmcCluster();
          emccluster->setEta( cl.momentum().Eta() );
          emccluster->setPhi( cl.momentum().Phi() );
          emccluster->setSigmaEta(-1.);
          emccluster->setSigmaPhi(-1.);
          emccluster->setEnergy( cl.energy() );
          emccluster->SetUniqueID( cl.key() );
#if 1
          for ( Int_t i=0; i< cl.numberOfTowers(); i++ ) 
            {
              StEmcRawHit *hit=cl.tower(i).stemc();
              assert( hit );         
              emccluster->addHit( hit );
            }
#endif


          pqr->addCluster( emccluster );

          mEtoEE[ emccluster ] = cl;
          cl.stemc( emccluster );

        }

    }
  else // if ( mNumberOfClusters[1, 2 OR 3] == 0 )
    {

      detector->setCluster( NULL );

    }

  
  StDetectorId ids[]={ kEndcapSmdUStripId, kEndcapSmdVStripId };


  for ( Int_t iplane=0; iplane<2; iplane++ ) 
    {

      detector=stevent->emcCollection()->detector(ids[iplane]);
      if ( !detector ) 
        {
          Warning("fillStEvent","detector == NULL for smd plane, no clusters for you");
        }
      else if ( mNumberOfClusters[4+iplane] > 0 )
        {
          
          StEmcClusterCollection *smdc = detector -> cluster();
          if ( !smdc ) 
            {
              //Warning("fillStEvent","StEmcClusterCollection (smd) was NULL, so I'm creating one.");
              smdc = new StEmcClusterCollection();
              detector->setCluster( smdc );
            }
          
          smdc->setDetector( ids[iplane] );
          smdc->setClusterFinderId( 123 );
          smdc->setClusterFinderParamVersion( 321 );

          
          for ( Int_t isector=0; isector<12; isector++ ) 
            {
              
              for ( UInt_t iclust=0; iclust<mSmdClusters[isector][iplane].size(); iclust++ ) 
                {
                  
                  StEEmcSmdCluster cl = (mSmdClusters[isector][iplane])[iclust];
                  

                  StEmcCluster *emccluster = new StEmcCluster();
                  emccluster->setEta( cl.mean() );
                  emccluster->setPhi( (Float_t)iplane );
                  emccluster->setSigmaEta(-1.);
                  emccluster->setSigmaPhi(-1.);
                  emccluster->setEnergy( cl.energy() );
                  emccluster->SetUniqueID( cl.key() );
                  for ( Int_t i=0; i< cl.numberOfStrips(); i++ ) 
                    {
                      StEmcRawHit *hit=cl.strip(i).stemc();
                      assert( hit );         
                      emccluster->addHit( hit );
                    }             
                  smdc->addCluster( emccluster );

                  mEtoEEsmd[ emccluster ] = cl;
                  cl.stemc( emccluster );
                  
                }
          
              
            }
          
        }

      else
        {
          detector -> setCluster( NULL );
        }



    }


}



// ----------------------------------------------------------------------------
Bool_t StEEmcClusterMaker::verifyStEvent() const
  {

  const StEvent *stevent=(const StEvent*)GetInputDS("StEvent");
  if ( !stevent ) {
    Warning("verifyStEvent","No StEvent found.");
    return true;
  }
  //StEmcCollection *emccollection=stevent->emcCollection();

  Bool_t go = true;
 
  const StEmcDetector *detector=stevent->emcCollection()->detector(kEndcapEmcTowerId);
  if ( !detector )
    {
      Warning("verifyStEvent","detector == NULL for towers");
      return false;
    }

  
  const StEmcClusterCollection *cc=detector->cluster();
  if ( cc )
    {


      
      Float_t emc_sum_towers = 0.;
      Float_t eemc_sum_towers = 0.;
  
      const StSPtrVecEmcCluster &emcClusters=cc->clusters();
      for ( UInt_t i=0; i<emcClusters.size(); i++ ) 
        {
          const StEmcCluster *cl=emcClusters[i];
          assert(cl);
          emc_sum_towers += cl->energy();
        }
      
      for ( Int_t sec=0; sec<12; sec++ ) 
        {
          
          for ( Int_t i=0; i<numberOfClusters(sec,0); i++ )
            eemc_sum_towers += cluster(sec,0,i).energy();
        }
      
      std::cout << "StEEmcClusterMaker tower checksum: ";
      if ( emc_sum_towers == eemc_sum_towers ) {
        std::cout << "passed";
      }
      else {
        std::cout << "FAILED"; go=false;
      }
      std::cout << std::endl;

    }
  else {

    std::cout << "StEEmcClusterMaker tower checksum: NULL collection, nclust=" << mNumberOfClusters[0] << std::endl;
    go &= (mNumberOfClusters[0]==0);

  }

  Float_t emc_sum_smdu = 0.;
  Float_t eemc_sum_smdu = 0.;

  detector=stevent->emcCollection()->detector(kEndcapSmdUStripId);
  if ( !detector )
    {
      Warning("verifyStEvent","detector == NULL for smdu");
      return false;
    }

  cc=detector->cluster();
  if ( cc ) 
    {

      const StSPtrVecEmcCluster &smduClusters=cc->clusters();

      for ( UInt_t i=0; i<smduClusters.size(); i++ ) 
        {
          const StEmcCluster *cl=smduClusters[i];
          assert(cl);
          emc_sum_smdu += cl->energy();
        }
      
  
      for ( Int_t sec=0; sec<12; sec++ ) 
        {

          for ( Int_t i=0; i<numberOfSmdClusters(sec,0); i++ )
        
            {
              eemc_sum_smdu += smdcluster(sec,0,i).energy();

            }

        }

      std::cout << "StEEmcClusterMaker smdu checksum: ";
      if ( emc_sum_smdu == eemc_sum_smdu ) {
        std::cout << "passed";
      }
      else {
        std::cout << "FAILED"; go=false;
      }
      std::cout << std::endl;

    }  
  else 
    {
      std::cout << "StEEmcClusterMaker smdu checksum: NULL collection, nclust=" << mNumberOfClusters[4] << std::endl;
      go &= (mNumberOfClusters[4]==0);
    }



  Float_t emc_sum_smdv = 0.;
  Float_t eemc_sum_smdv = 0.;

  detector=stevent->emcCollection()->detector(kEndcapSmdVStripId);
  if (!detector)
    {
      Warning("verifyStEvent","detector == NULL for smdv");
      return false;
    }

  cc=detector->cluster();

  if ( cc )
    {

      const StSPtrVecEmcCluster &smdvClusters=cc->clusters();
      
      for ( UInt_t i=0; i<smdvClusters.size(); i++ ) 
        {
          const StEmcCluster *cl=smdvClusters[i];
          assert(cl);
          emc_sum_smdv += cl->energy();
        }
      
      
      for ( Int_t sec=0; sec<12; sec++ ) 
        {
          
          for ( Int_t i=0; i<numberOfSmdClusters(sec,1); i++ )
            
            {
              eemc_sum_smdv += smdcluster(sec,1,i).energy();
              
            }
          
        }
      
      std::cout << "StEEmcClusterMaker smdv checksum: ";
      if ( emc_sum_smdv == eemc_sum_smdv ) {
        std::cout << "passed";
      }
      else {
        std::cout << "FAILED"; go=false;
      }
      std::cout << std::endl;

    }
  else 
    {
      std::cout << "StEEmcClusterMaker smdv checksum: NULL collection, nclust=" << mNumberOfClusters[5] << std::endl;
      go &= (mNumberOfClusters[5]==0);
    }
  

  
  return go;
}



// ----------------------------------------------------------------------------
void StEEmcClusterMaker::print() const
{

  std::cout << "StEEmcClusterMaker::print()" << std::endl;
  const Char_t *names[] = { "tower", "pre1", "pre2", "post", "smdu", "smdv" };
  for ( Int_t i=0;i<6;i++ )
    {

      std::cout << "Number of " << names[i] 
                << " clusters = " << mNumberOfClusters[i] 
                << std::endl;

    }

  std::cout << "printout of tower clusters follows:" << std::endl;
  for ( Int_t sec=0;sec<12;sec++)
  for ( Int_t i=0;i<numberOfClusters(sec,0);i++ )
    {
      const StEEmcCluster &clust=cluster(sec,0,i);
      clust.print();
      //      std::cout << "cluster.key()=" << clust.key() << std::endl;
      //      std::cout << "cluster.eta()=" << clust.momentum().Eta() << std::endl;
      //      std::cout << "cluster.phi()=" << clust.momentum().Phi() << std::endl;      
      //      std::cout << "cluster.energy()=" << clust.energy() << std::endl;
    }

  

}