StEEmcPointFitMaker.cxx

#include "StEEmcPointFitMaker.h"
#include "EEmcSectorFit.h" 
#include "StEEmcPool/StEEmcA2EMaker/StEEmcA2EMaker.h"
#include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
#include <algorithm>
#include "TString.h"
#include "StEEmcPool/StEEmcClusterMaker/StEEmcSmdCluster.h" 

ClassImp(StEEmcPointFitMaker);

// ----------------------------------------------------------------------------
StEEmcPointFitMaker::StEEmcPointFitMaker(const Char_t *n):StEEmcPointMaker(n)
{
  for ( Int_t sec=0;sec<12;sec++ )
    {
      mSectorFit[sec]=0;
    }

  mPermutations=true;

  mLimitFits=10;

}

// ----------------------------------------------------------------------------
Int_t StEEmcPointFitMaker::Init()
{
    return StEEmcPointMaker::Init();
}

// ----------------------------------------------------------------------------
Int_t StEEmcPointFitMaker::Make()
{
    mClusterId=10000; 
  for ( Int_t sec=0;sec<12;sec++ )
    {
      mSectorFit[sec]=new EEmcSectorFit(10); 
      mSectorFit[sec]->doPermutations=mPermutations;
      //mSectorFit[sec]->SetPrintLevel(-1); 
    }

  for ( Int_t i=0;i<12;i++ ) 
    {
      FitSector(i); 
    }

  if ( mEnergyMode == 0 )
    shareEnergySmd();
  else if ( mEnergyMode == 1 )
    shareEnergySimple();


  return kStOK;
}
// ----------------------------------------------------------------------------
Int_t StEEmcPointFitMaker::FitSector(Int_t sector)
{

  EEmcSectorFit *fitter = mSectorFit[sector];

  Float_t etowers = mEEanalysis->energy(sector,0);
  Float_t eproject = etowers * 0.007 * 1000.0; // expected smd E [MeV]
  Float_t nproject = eproject / 1.3; // expected smd E [nmips]


  TString uname="h";uname+=(sector+1<10)?"0":"";uname+=sector+1;uname+="U";
  TString vname="h";vname+=(sector+1<10)?"0":"";vname+=sector+1;vname+="V";
  TH1F *hU=new TH1F(uname,"SMD-u energy distribution [nmips]",288,0.,288.);
  TH1F *hV=new TH1F(vname,"SMD-v energy distribution [nmips]",288,0.,288.);

  for ( Int_t hit=0;hit<mEEanalysis->numberOfHitStrips(sector,0);hit++ ) {
    StEEmcStrip strip=mEEanalysis->hitstrip(sector,0,hit);
    if ( strip.fail() ) continue;
    hU->SetBinContent( strip.index()+1, strip.energy() * 1000.0 / 1.3 );
  }
  for ( Int_t hit=0;hit<mEEanalysis->numberOfHitStrips(sector,1);hit++ ) {
    StEEmcStrip strip=mEEanalysis->hitstrip(sector,1,hit);
    if ( strip.fail() ) continue;
    hV->SetBinContent( strip.index()+1, strip.energy() * 1000.0 / 1.3 );
  }

  mSectorFit[sector]->SetHistograms(hU,hV);


  if ( etowers < 1.0 ) return kStOK;
  

  Float_t nfit = 0.;
  Float_t nlast = 5.0;
  Int_t count = 0;

  while ( nfit < 0.90 * nproject && nlast >= 4.0 && count < mLimitFits ) {
    
    Double_t res5umax=0.;
    Double_t res5vmax=0.;
    Double_t res5u=0.;
    Double_t res5v=0.;
    Int_t maxu=-1;
    Int_t maxv=-1;
    for ( Int_t i=0;i<288;i++ )
      {
        res5u=fitter->Residual(i,0,2);
        res5v=fitter->Residual(i,1,2);
        // resu=fitter->Residual(i,0);
        // resv=fitter->Residual(i,1);
        if ( res5u > res5umax ) {
          res5umax=res5u;
          maxu=i;
        }
        if ( res5v > res5vmax ) {
          res5vmax=res5v;
          maxv=i;
        }
      }

    if ( maxu<0 || maxv<0 ) {
      //      std::cout << "maxu || maxv < 0, abort" << std::endl;
      break;
    }
    if ( res5umax < 4.0 && res5vmax < 4.0 ) {
      //      std::cout << "Not enough signal, abort" << std::endl;
      break;
    }

    //  std::cout << "Adding candidate, now I expect some output!" << std::endl;
    fitter -> AddCandidate( (res5umax+res5vmax)/2., 0.85, (Double_t)maxu+0.5, (Double_t)maxv+0.5 );

    fitter -> Migrad();

    Double_t e,s,u,v;
    fitter->GetLastCandidate( e, s, u, v );
    nlast = e;

    if ( count && count < 5 ) fitter -> TryPermutations();

    nfit = 0.;
    for ( Int_t i=0;i<fitter->numberOfCandidates();i++ ){

      fitter->GetCandidate(i,e,s,u,v);
      nfit += e;

    }
    
    count++;
    if ( nlast < 4.0 ) break;

  }


  //$$$  std::cout << "found ncand=" << fitter->numberOfCandidates() << std::endl;


  std::vector<Bool_t> drop( fitter->numberOfCandidates(), false );
  for ( Int_t ifit=0;ifit<fitter->numberOfCandidates()-1;ifit++ )
    {
      Double_t e1,s1,u1,v1;
      fitter->GetCandidate(ifit,e1,s1,u1,v1);
      for ( Int_t jfit=ifit+1;jfit<fitter->numberOfCandidates();jfit++ )
        {
          Double_t e2,s2,u2,v2;
          fitter->GetCandidate(jfit,e2,s2,u2,v2);

          Float_t dist = TMath::Sqrt( (u1-u2)*(u1-u2) + (v1-v2)*(v1-v2) );

          if ( e1 < 0.15 * e2 ) {           
            Float_t min_dist = TMath::Abs(s1) * 3.5;
            if ( dist < min_dist ) drop[ifit]=true;
          }

          if ( e2 < 0.15 * e1 ) {
            Float_t min_dist = TMath::Abs(s2) * 3.5;
            if ( dist < min_dist ) drop[jfit]=true;
          }

        }
    }




  for ( Int_t ifit=0;ifit<fitter->numberOfCandidates();ifit++ )
    {

      if ( drop[ifit] ) continue;

      Double_t nmip,sigma,u,v;
      fitter->GetCandidate( ifit, nmip, sigma, u, v );
      StEEmcPoint point;
      point.sector( sector );
      point.sigma( sigma );
      point.u(u);
      point.v(v); 
      point.energy( (Float_t)(nmip * 1.3 / 0.007 / 1000.0) );
      TVector3 position = mEEsmd->getIntersection( sector, (Int_t)u, (Int_t)v );
      point.position( position );
      StEEmcTower *tower = mEEanalysis->tower(position,0);
      if ( !tower ) continue; 
      point.tower( *tower );
      point.residueU( (Float_t)fitter -> Residual( (Int_t)u, 0, 60 ) * 1.3 / 1000.0 );
      point.residueV( (Float_t)fitter -> Residual( (Int_t)v, 1, 60 ) * 1.3 / 1000.0 );
                      

      StEEmcSmdCluster uc;
      uc.energy( (float)(nmip * 1.3 / 1000.0) );
      uc.mean( (float)u );
      uc.sigma( (float)sigma ); 
      uc.sector ( sector );
      uc.plane( 0 ); 
      uc.key( mClusterId++ ); 
      StEEmcSmdCluster vc;
      vc.energy( (float)(nmip * 1.3 / 1000.0) );
      vc.mean ( (float)v );
      vc.sigma( (float)sigma );
      vc.plane( 1 ); 
      vc.key( mClusterId++ ); 
      point.cluster( uc, 0 );
      point.cluster( vc, 1 ); 
      point.energy( (float)(nmip * 1.3 / 1000.0) );
      point.sector( sector ); 
     
      mPoints.push_back(point);

    }

  return kStOK;

}
// ----------------------------------------------------------------------------
void StEEmcPointFitMaker::Clear(Option_t *opts)
{
  for ( Int_t i=0;i<12;i++ )
    if ( mSectorFit[i] ) delete mSectorFit[i];

  StEEmcPointMaker::Clear(opts);
}

// ----------------------------------------------------------------------------
void StEEmcPointFitMaker::print() const
{
    for ( Int_t i=0;i<12;i++ )
    {
        mSectorFit[i]->print(); 
    } 
    for ( UInt_t i=0;i<mPoints.size();i++ )
    {
        mPoints[i].print(); 
    } 
}