Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

EnergyModule - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: EnergyModule.cc,v 1.74 2005/11/14 17:02:49 naumov Exp $
// Naumov Dmitry created Jul,  4 2004

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  EnergyModule                                                            //
//                                                                          //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////
//
/* */ //

//
/* */ //


#include "EnergyModule.hh"
#include "RecoEvent.hh"
#include "RecoPixelData.hh"
#include "Config.hh"
#include "EOpticsResponse.hh"
#include "EConst.hh"
#include "RecoRootEvent.hh"
#include "RecoShowerTrack.hh"

// ---> Begin of Energy Module helpers
#include "EnergyViewer.hh"
#include "EnergyFitUtils.hh"
// <--- End of Energy Module helpers

// ---> Begin of Root includes
#include <TCutG.h>
#include <TMinuit.h>
#include <TSpectrum.h>
#include <TProfile.h>
// <--- End of Root includes
#include <map>
#include <vector>

using namespace TMath;
using namespace EConst;
using namespace sou;

ClassImp(EnergyModule)

//_____________________________________________________________________________
 EnergyModule::EnergyModule() : RecoModule("Energy") {
    // ctor
    const string& ORdbname = Conf()->GetStr("EnergyModule.fORname");
    fOpticsResponse = new EOpticsResponse(ORdbname.c_str());
    Constructor();
}

//_____________________________________________________________________________
 void EnergyModule::Constructor() {
    //
    // Default ctor initialization
    fEvent             = NULL;
    fShowerTrack       = NULL;
    fGraphTheta        = NULL;
    fGraphPhi          = NULL;
    fGraphXY           = NULL;
    fGraphX            = NULL;
    fGraphY            = NULL;
    fHistoOverlap      = NULL;
    fHistoIdealOverlap = NULL;
    fUseTrueAngles     = kFALSE;
}
//_____________________________________________________________________________
 EnergyModule::~EnergyModule() {

    // dtor
    SafeDelete(fShowerTrack);
    ClearGraphs();
}

//_____________________________________________________________________________
 Bool_t EnergyModule::Init() {    
    // Get fluorescence calculator
    Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
    fDebugPs           = (Int_t)Conf()->GetNum("EnergyModule.fDebugPs");
    fMinuitOutputLevel = (Int_t)Conf()->GetNum("EnergyModule.fMinuitOutputLevel");
    if ( fMinuitOutputLevel < -1 || fMinuitOutputLevel > 3 ) {
      Msg(EsafMsg::Warning) << "Wrong config value EnergyModule.fMinuitOutputLevel " 
			    << fMinuitOutputLevel << " Setted to 1" << MsgDispatch;
      fMinuitOutputLevel = 1;
    }
    fUseTrueAngles     = (Int_t)Conf()->GetNum("EnergyModule.fUseTrueAngles");

    // Initialization of pixels viewer

    if (fDebugPs >= 1) fEnergyViewer = new EnergyViewer(this);

    fShowerTrack = new RecoShowerTrack;

    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t EnergyModule::PreProcess() {

    fTheta          = -1;
    fPhi            = -1;
    fAltitude       = -1; 
    fEnergyChi2     = -1;
    fThetaFoVAtMax  = -1;
    fPhiFoVAtMax    = -1;
    fGtuTimeAtMax   = -1;
    fBackGround     = -1;
    fFCtime[0] = fFCtime[1] = fFCtime[2] = fFCtime[3] = -1;
    fPixelsCluster.clear();
    fAllActivePixels.clear();

    return kTRUE;
}

//_____________________________________________________________________________
 void EnergyModule::CreateAuxilaryObjects() {

    ClearGraphs();
    fGraphTheta       = new TGraph();
    fGraphPhi         = new TGraph();
    fGraphXY          = new TGraph();
    fGraphX           = new TGraph();
    fGraphY           = new TGraph();

    if (fDebugPs >= 1) fEnergyViewer->InitFilePS();
}

//_____________________________________________________________________________
 void EnergyModule::ClearGraphs() {
    SafeDelete(fGraphTheta);
    SafeDelete(fGraphPhi);
    SafeDelete(fGraphXY);
    SafeDelete(fGraphX);
    SafeDelete(fGraphY);
}

//_____________________________________________________________________________
 Bool_t EnergyModule::LoadModules() {

  vector <Int_t> PixelsClusterInt;
  
    // Load information from previous modules essential for the energy reconstruction
    fTDM = fEvent->GetModuleData("TrackDirection");
    
    // try to get TrackDirection2
    if ( fTDM == NULL ) fTDM = fEvent->GetModuleData("TrackDirection2"); 
    if ( fTDM == NULL ) Msg(EsafMsg::Panic) << "None of TrackDirectionModule and TrackDirectionModule2 is not found." << MsgDispatch;
    
    if (fUseTrueAngles) {
        fTheta = fEvent->GetHeader().GetTrueTheta();
        fPhi   = fEvent->GetHeader().GetTruePhi();
    }
    else {
        fTheta = fTDM->GetDouble("Theta");
        fPhi   = fTDM->GetDouble("Phi");
        if(IsNaN(fTheta) || IsNaN(fPhi)) return kFALSE;
    }
    
    fShowerDir.SetXYZ(Sin(fTheta)*Cos(fPhi),Sin(fTheta)*Sin(fPhi),Cos(fTheta));
    fGtuCM = fEvent->GetModuleData("GTUClustering");
    
    if ( fGtuCM == NULL )  fGtuCM = fEvent->GetModuleData("HoughTransform");
    if ( fGtuCM == NULL )  Msg(EsafMsg::Panic) << "None of GTUClusteringModule and HoughModule is not found." << MsgDispatch;
    
    if ( fGtuCM->GetObj("CluPixels") == NULL ) 
      Msg(EsafMsg::Panic) << "Gtu Cluster Module has NULL CluPixels" << MsgDispatch;

    PixelsClusterInt.clear();

    PixelsClusterInt = *(vector<Int_t>*) fGtuCM->GetObj("CluPixels");
    fAllActivePixels =  (Pixels_t)fEvent->GetRecoPixels();

    for (UInt_t i (0); i < PixelsClusterInt.size(); i++) {
      RecoPixelData *pix = fEvent->GetRecoPixelData(PixelsClusterInt[i]);
      fPixelsCluster.push_back(pix);
    }
    
    if (!fPixelsCluster.size()) 
      Msg(EsafMsg::Panic) << "CluPixels has no pixels" << MsgDispatch;

    Msg(EsafMsg::Debug) << "Processing " << fPixelsCluster.size() 
			<< " points" << MsgDispatch;
    
    fHmaxShapeM  = fEvent->GetModuleData("HmaxByShapeMethod");
    fHmaxProtonM = fEvent->GetModuleData("HmaxForProton");
    
    if (fHmaxShapeM == NULL && fHmaxProtonM == NULL) fAltitude = 5*km;
    if (fHmaxShapeM)  fAltitude = fHmaxShapeM->GetDouble("Altitude");
    if (fHmaxProtonM) fAltitude = fHmaxProtonM->GetDouble("Altitude");

    return kTRUE;
}


//_____________________________________________________________________________
 void EnergyModule::SortGtu() {

    Int_t pmtId_previous = -1;
    Float_t OffSet = 10.; // mm

    map <Int_t, Int_t> EC;
    PMTs_t::iterator itCell, itPmt;
    Pixels_t::iterator itPix;
    map <Int_t, Int_t>::iterator itEC;
    
    fGtuMin = fMyGtuMin =  10000;
    fGtuMax = fMyGtuMax = -10000;
    fXmin   = fYmin = fPhiMin = fThetaMin =  kHuge;
    fXmax   = fYmax = fPhiMax = fThetaMax = -kHuge;

    fGtuPixels.clear(); fGtuMyPixels.clear(); fMyPixelsCluster.clear();
    fPixelsBgr.clear(); fECells.clear(); fPmts.clear(); EC.clear();
    fActivePixels.clear();

    // Sort Gtu and define cluster position

    for (itPix = fPixelsCluster.begin(); itPix != fPixelsCluster.end(); itPix++) {
      Int_t   pmtId = fRunpars->Pmt((*itPix)->GetPixelId());
      Int_t   gtu   = (*itPix)->GetGtu();
      Float_t theta = (*itPix)->GetTheta()*RadToDeg();
      Float_t phi   = (*itPix)->GetPhi()*RadToDeg();

      fGtuPixels[gtu].push_back(*itPix);

      if (fGtuMin   > gtu)   fGtuMin   = gtu;
      if (fGtuMax   < gtu)   fGtuMax   = gtu;
      if (fPhiMin   > phi)   fPhiMin   = phi;
      if (fPhiMax   < phi)   fPhiMax   = phi;
      if (fThetaMin > theta) fThetaMin = theta;
      if (fThetaMax < theta) fThetaMax = theta;

      if (pmtId == pmtId_previous) continue; // Avoid double calculation 
      pmtId_previous = pmtId;                // for PMT's

      for(Int_t iCorner(0); iCorner < 4; iCorner++){
	TVector3 vPmt = fRunpars->PmtCorner(pmtId,iCorner);
	if (fXmin > vPmt.X()) fXmin = vPmt.X();
	if (fYmin > vPmt.Y()) fYmin = vPmt.Y();
	if (fXmax < vPmt.X()) fXmax = vPmt.X();
	if (fYmax < vPmt.Y()) fYmax = vPmt.Y();
      }
      EC[fRunpars->GetPmtGeo(pmtId)->GetEC()] = 0; // list of macrocells
    }

    // Make list of macrocells

    for (itEC = EC.begin(); itEC != EC.end(); itEC++) 
      fECells.push_back(itEC->first);

    fXmin -= OffSet; fYmin -= OffSet;
    fXmax += OffSet; fYmax += OffSet;

    // prepeare list of PMTs and pixels

    for(Int_t i = 1; i <= fRunpars->GetNumPmts(); i++) {
      TVector3 vPmt = fRunpars->GetPmtGeo(i)->GetCenter();
      if (vPmt.X() < fXmin || vPmt.X() > fXmax || 
	  vPmt.Y() < fYmin || vPmt.Y() > fYmax ) continue;

      for (itCell = fECells.begin(); itCell != fECells.end();itCell++) {
	if (*itCell == fRunpars->GetPmtGeo(i)->GetEC()) {
	  fPmts.push_back(i);
	  const EPmtGeo* pmt = fRunpars->GetPmtGeo(i);
	  for (Int_t PixUid = pmt->GetStartUniqueId(); PixUid <= pmt->GetLastUniqueId(); PixUid++)
	    fPixelsBgr.push_back(PixUid);
	  break;
	}
      }
    }

    fMyGtuMin        = fGtuMin;
    fMyGtuMax        = fGtuMax;
    fMyPixelsCluster = fPixelsCluster;
    fGtuMyPixels     = fGtuPixels;

    // Make map of active pixels around cluster

    for (itPix = fAllActivePixels.begin(); itPix != fAllActivePixels.end(); itPix++)
      if ((*itPix)->GetGtu() >= fMyGtuMin && (*itPix)->GetGtu() <= fMyGtuMax) 
	for (itPmt = fPixelsBgr.begin(); itPmt != fPixelsBgr.end(); itPmt++)
	  if (*itPmt == (*itPix)->GetPixelId()) {
	    fActivePixels.push_back(*itPix);
	    break;
	  }

    InitHisto();                 // Histogram initialization
    BuildFoVAngleFunctions();    // Build angles and fit distributions
    BuildMyPixels();             // Define my pixels around signal

    if (fDebugPs >= 1) fEnergyViewer->DisplayFSTrack();
    if (fDebugPs >= 1) fEnergyViewer->DisplayPixels();
    if (fDebugPs >= 2) DumpGtuSorted();

    EC.clear();
}


//_____________________________________________________________________________
 void EnergyModule::InitHisto() {

    Float_t OffSet = 15;
    Float_t OffSetDeg = 2;

    if (fPhiMax > 150. && fPhiMin < -150.) {
      fHistoPhi   = new TH2F("GtuPhi","GtuPhi",fMyGtuMax-fMyGtuMin,
			     fMyGtuMin,fMyGtuMax,3600,-360,360);
    }
    else {
      fHistoPhi   = new TH2F("GtuPhi","GtuPhi",fMyGtuMax-fMyGtuMin,
			     fMyGtuMin,fMyGtuMax,
			     5*(Int_t)(fPhiMax-fPhiMin+2*OffSetDeg),
			     fPhiMin-OffSetDeg,fPhiMax+OffSetDeg);
    }
    
    fHistoTheta = new TH2F("GtuTheta","GtuTheta",fMyGtuMax-fMyGtuMin,
			   fMyGtuMin,fMyGtuMax,
			   5*(Int_t)(fThetaMax-fThetaMin+2*OffSetDeg),
			   fThetaMin-OffSetDeg,fThetaMax+OffSetDeg);
    fHistoX     = new TH2F("GtuX","GtuX",fMyGtuMax-fMyGtuMin,
			   fMyGtuMin,fMyGtuMax,
			   5*(Int_t)(fXmax-fXmin+2*OffSet),
			   fXmin-OffSet,fXmax+OffSet);
    fHistoY     = new TH2F("GtuY","GtuY",fMyGtuMax-fMyGtuMin,
			   fMyGtuMin,fMyGtuMax,
			   5*(Int_t)(fYmax-fYmin+2*OffSet),
			   fYmin-OffSet,fYmax+OffSet);
    fHistoXY    = new TH2F("XY","XY",
			   5*(Int_t)(fXmax-fXmin+2*OffSet),
			   fXmin-OffSet,fXmax+OffSet,
			   5*(Int_t)(fYmax-fYmin+2*OffSet),
			   fYmin-OffSet,fYmax+OffSet);

    fHistoOverlap      = new TProfile("Overlap","Overlap",
				      fMyGtuMax-fMyGtuMin,fMyGtuMin,fMyGtuMax);
    fHistoIdealOverlap = new TProfile("IdealOverlap","Overlap",
				      fMyGtuMax-fMyGtuMin,fMyGtuMin,fMyGtuMax);
}


//_____________________________________________________________________________
 void EnergyModule::BuildFoVAngleFunctions() {
    //
    //
    //
    Float_t binx, biny, bintheta, binphi, phip;
    Pixels_t OneGtuPixels;
    Pixels_t::iterator itPix;
    GtuPixels_t::iterator
    Int_t pixel = 0;

    itPix = (fGtuMyPixels.begin()->second).begin();
    phip = (*itPix)->GetPhi()*RadToDeg();

    for (itGtu = fGtuMyPixels.begin(); itGtu != fGtuMyPixels.end(); itGtu++) {
         OneGtuPixels.clear();
         OneGtuPixels = itGtu->second;
	 Int_t Count = 0;
         binx = biny = bintheta = binphi = 0;

         for (itPix = OneGtuPixels.begin(); itPix != OneGtuPixels.end(); itPix++) {
	   Int_t nc = (*itPix)->GetCounts();
	   Float_t phi = (*itPix)->GetPhi()*RadToDeg();
	   if (phi-phip > 300.) {phi -= 360; phip = phi;}
	   else if (phi-phip < -300.) {phi += 360; phip = phi;}

	   fHistoPhi->Fill(itGtu->first,phi,nc);
	   fHistoTheta->Fill(itGtu->first,(*itPix)->GetTheta()*RadToDeg(),nc);
	   TVector3 PixelCenter=fRunpars->PixelCenter((*itPix)->GetPixelId());
	   fHistoX->Fill(itGtu->first,PixelCenter.X(),nc);
	   fHistoY->Fill(itGtu->first,PixelCenter.Y(),nc);
 	   fHistoXY->Fill(PixelCenter.X(),PixelCenter.Y(),nc);

	   bintheta += (*itPix)->GetTheta()*RadToDeg()*nc;
	   binphi   += phi*nc;
	   binx     += PixelCenter.X()*nc;
	   biny     += PixelCenter.Y()*nc;
	   Count    += nc;
         }

	 for (Int_t icount = 0; icount < Count; icount++) {
	   fGraphTheta->SetPoint(pixel,itGtu->first,bintheta/Count);
	   fGraphPhi->SetPoint(pixel,itGtu->first,binphi/Count);
	   fGraphX->SetPoint(pixel,itGtu->first,binx/Count);
	   fGraphY->SetPoint(pixel,itGtu->first,biny/Count);
	   fGraphXY->SetPoint(pixel,binx/Count,biny/Count);
	   pixel++;
	 }
    }

    fGraphTheta->Fit("pol2","QEM");
    fGraphPhi->Fit("pol2","QEM");
    fGraphX->Fit("pol1","QEM");
    fGraphY->Fit("pol1","QEM");
    fGraphXY->Fit("pol1","QEM");
}


//_____________________________________________________________________________
 void EnergyModule::BuildMyPixels() {

    // Define the box around fitted line

    PMTs_t::iterator itPmt;

    Float_t  Xnew[5], Ynew[5];
    TVector2 VBox1,   VBox2;

    Float_t  HistoSize = 15.*mm;  // It is necessary to define it via Histo
    TVector2 beginLine = XYCluster(fMyGtuMin);
    TVector2 endLine   = XYCluster(fMyGtuMax);
    TVector2 Line      = endLine-beginLine;

    fMyPixels.clear();
    VBox1.SetMagPhi(2.*HistoSize,Line.Phi());

    // Coordinate in box corners

    for (Int_t Rot = 45, i = 0; Rot <= 315; Rot += 90, i++) {
      VBox2 = VBox1.Rotate(Rot*DegToRad());
      Xnew[i] = VBox2.X(); Ynew[i] = VBox2.Y();
      if (!i || i == 3) {Xnew[i] += endLine.X(); Ynew[i] += endLine.Y();}
      else {Xnew[i] += beginLine.X(); Ynew[i] += beginLine.Y();}	
    }

    Xnew[4] = Xnew[0]; Ynew[4] = Ynew[0];
    TCutG BoxCut("BoxCut",5,Xnew,Ynew);

    // Define pixels within box

    for (itPmt = fPixelsBgr.begin(); itPmt != fPixelsBgr.end(); itPmt++) {
      TVector3 VCenter = fRunpars->PixelCenter(*itPmt);
      if (BoxCut.IsInside(VCenter.X(),VCenter.Y()))
	fMyPixels.push_back(*itPmt);
    }
}


//_____________________________________________________________________________
 void EnergyModule::DumpGtuSorted() {
  // dump on screen

  Int_t MaxCounts = 0, GtuTimeAtMax = 0;
  Pixels_t OneGtuPixels;
  Pixels_t::iterator itPix;
  GtuPixels_t::iterator


  Msg(EsafMsg::Debug) << "Dumped maps of pixels" << MsgDispatch;

  for (itGtu = fGtuPixels.begin(); itGtu != fGtuPixels.end(); itGtu++) {
    OneGtuPixels.clear();
    OneGtuPixels = itGtu->second;
    Msg(EsafMsg::Debug) << Form("GTU Entry: %d vector of %d pixels:",
				itGtu->first,OneGtuPixels.size())<<MsgDispatch;

    Int_t sum(0);
    for (itPix = OneGtuPixels.begin(); itPix != OneGtuPixels.end(); itPix++) {
      Msg(EsafMsg::Debug) << Form("          GTU %d PixelUid %d Counts %d", 
				  (*itPix)->GetGtu(),(*itPix)->GetPixelId(),
				  (*itPix)->GetCounts()) << MsgDispatch;
      sum += (*itPix)->GetCounts();
    }

    Msg(EsafMsg::Debug) << Form("GTU Entry: Total counts %d",sum)<<MsgDispatch;
    if (MaxCounts < sum) {
      MaxCounts = sum;
      GtuTimeAtMax = itGtu->first;
    }
  }

  Msg(EsafMsg::Debug) << Form("Most numerous GTU %d with %d count",
			      Int_t(GtuTimeAtMax),MaxCounts) << MsgDispatch;
  Msg(EsafMsg::Debug) << Form("Min and Max GTU  %d  %d ",Int_t(fGtuMin),
			      Int_t(fGtuMax)) << MsgDispatch;
  Msg(EsafMsg::Debug) << Form("DumpGtuSorted(): End Dumping") << MsgDispatch;
  Msg(EsafMsg::Debug) << Form("DumpGtuSorted(): Dumping My map <GtuHit_t,Pixels_t>") << MsgDispatch;

  MaxCounts = 0;
  for (itGtu = fGtuMyPixels.begin(); itGtu != fGtuMyPixels.end(); itGtu++) {
    OneGtuPixels.clear();
    OneGtuPixels = itGtu->second;

    Msg(EsafMsg::Debug) << Form("GTU Entry: %d vector of %d pixels:",
				itGtu->first,OneGtuPixels.size())<<MsgDispatch;

    Int_t sum(0);
    for (itPix = OneGtuPixels.begin(); itPix != OneGtuPixels.end(); itPix++) {
      Msg(EsafMsg::Debug) << Form("          GTU %d PixelUid %d Counts %d",
				  (*itPix)->GetGtu(), (*itPix)->GetPixelId(),
				  (*itPix)->GetCounts()) << MsgDispatch;

      sum += (*itPix)->GetCounts();
    }
    
    Msg(EsafMsg::Debug) << Form("GTU Entry: Total counts %d",sum)<<MsgDispatch;

    if (MaxCounts < sum) {
      MaxCounts = sum;
      GtuTimeAtMax = itGtu->first;
    }
  }
  Msg(EsafMsg::Debug) << Form("Most numerous GTU %d with %d count",
			      Int_t(GtuTimeAtMax),MaxCounts) << MsgDispatch;
  Msg(EsafMsg::Debug) << Form("Min and Max GTU  %d  %d ",Int_t(fMyGtuMin),
			      Int_t(fMyGtuMax)) << MsgDispatch;
  Msg(EsafMsg::Debug) << Form("DumpGtuSorted(): End Dumping") << MsgDispatch;
}


//_____________________________________________________________________________
 TVector2 EnergyModule::XYCluster(Float_t time) {
    //
    if (!fGraphX || !fGraphY )  Msg(EsafMsg::Panic) << "fGraphX or fGraphY does not exist! Check your code." << MsgDispatch;
    
    Float_t XP0 = fGraphX->GetFunction("pol1")->GetParameter(0);
    Float_t XP1 = fGraphX->GetFunction("pol1")->GetParameter(1);
    Float_t YP0 = fGraphY->GetFunction("pol1")->GetParameter(0);
    Float_t YP1 = fGraphY->GetFunction("pol1")->GetParameter(1);

    TVector2 PosCluster(XP0 + XP1*time,YP0 + YP1*time);
    return PosCluster;
}


//_____________________________________________________________________________
 Float_t EnergyModule::ThetaFoV(Float_t time) {
    //
    //
    //
    if (!fGraphTheta) Msg(EsafMsg::Panic) << "fGraphTheta does not exist! Check your code." << MsgDispatch;
    Float_t ThetaP0 = fGraphTheta->GetFunction("pol2")->GetParameter(0);
    Float_t ThetaP1 = fGraphTheta->GetFunction("pol2")->GetParameter(1);
    Float_t ThetaP2 = fGraphTheta->GetFunction("pol2")->GetParameter(2);

    return DegToRad()*(ThetaP0 + ThetaP1*time + ThetaP2*time*time);
}


//_____________________________________________________________________________
 Float_t EnergyModule::PhiFoV(Float_t time) {
    //
    if (!fGraphPhi) Msg(EsafMsg::Panic) << "fGraphPhi does not exist! Check your code." << MsgDispatch;
    Float_t PhiP0 = fGraphPhi->GetFunction("pol2")->GetParameter(0);
    Float_t PhiP1 = fGraphPhi->GetFunction("pol2")->GetParameter(1);
    Float_t PhiP2 = fGraphPhi->GetFunction("pol2")->GetParameter(2);

    return DegToRad()*(PhiP0 + PhiP1*time + PhiP2*time*time);
}


//_____________________________________________________________________________
 Bool_t EnergyModule::BuildRecoveryFunction() {    
    // Build the recovery function along the track
    // for a number of steps along the track compute overlap of efficacy with pixels coverage.

    Float_t  Xnew[5], Ynew[5];

    PMTs_t::iterator itPixUid;
    vector <TCutG> Cuts;
    vector <TCutG>::const_iterator itCuts;

    Int_t   progress    = 0;
    Float_t MeanLambda  = 360*nm;
    Float_t GtuOffSet   = 0., TimeStep = 0.25;
    Float_t InitialTime = fMyGtuMin - GtuOffSet + 0.5*TimeStep;
    Float_t FinalTime   = fMyGtuMax + GtuOffSet - 0.4*TimeStep; // Not error
    Float_t TotalTime   = fMyGtuMax - fMyGtuMin + 2.*GtuOffSet - TimeStep;

    // loop along the track

    for (Float_t time = InitialTime; time <= FinalTime; time += TimeStep) {
      Float_t fThetaFoV  = ThetaFoV(time);
      Float_t fPhiFoV    = PhiFoV(time);
      TVector2 XYcluster = XYCluster(time);

      if (fThetaFoV > 30*DegToRad()) continue;

      TH2F* htemp = fOpticsResponse->MakePsfHist(fThetaFoV,fPhiFoV,MeanLambda);
      if (!htemp) continue;

      Float_t xminPsf = htemp->GetXaxis()->GetXmin();
      Float_t xmaxPsf = htemp->GetXaxis()->GetXmax();
      Float_t yminPsf = htemp->GetYaxis()->GetXmin();
      Float_t ymaxPsf = htemp->GetYaxis()->GetXmax();

      Float_t deltaX = XYcluster.X() - 0.5*(xminPsf+xmaxPsf);
      Float_t deltaY = XYcluster.Y() - 0.5*(yminPsf+ymaxPsf);

      xminPsf += deltaX; xmaxPsf += deltaX;
      yminPsf += deltaY; ymaxPsf += deltaY;

      TH2F* Psf = new TH2F("Psf","",htemp->GetNbinsX(),xminPsf,xmaxPsf,
			   htemp->GetNbinsY(),yminPsf,ymaxPsf);

      for (Int_t j = 1; j <= htemp->GetNbinsX(); j++)
	for (Int_t k = 1; k <= htemp->GetNbinsY(); k++) {
	  Float_t x = htemp->GetXaxis()->GetBinCenter(j) + deltaX;
	  Float_t y = htemp->GetYaxis()->GetBinCenter(k) + deltaY;
	  Psf->Fill(x,y,htemp->GetBinContent(j,k));
	}

      SafeDelete(htemp);

      // Define the square for signal integration

      Float_t PixelSide = 2.*mm;
      Float_t histoX[5], histoY[5];

      histoX[1] = histoX[2] = xmaxPsf + PixelSide;
      histoY[2] = histoY[3] = yminPsf - PixelSide;
      histoX[0] = histoX[3] = histoX[4] = xminPsf - PixelSide;
      histoY[0] = histoY[1] = histoY[4] = ymaxPsf + PixelSide;

      TCutG HistoCut("HistoCut",5,histoX,histoY);

      // Define the pixels inside histogram

      Cuts.clear();

      for (itPixUid = fMyPixels.begin(); itPixUid != fMyPixels.end(); itPixUid++) {
	TVector3 VCenter = fRunpars->PixelCenter(*itPixUid);
	if (HistoCut.IsInside(VCenter.X(),VCenter.Y())) {
	  for(Int_t iCorner = 0; iCorner < 4; iCorner++) {
	    TVector3 dummy = fRunpars->PixelCorner(*itPixUid,iCorner);
	    Xnew[iCorner] = dummy.X(); Ynew[iCorner] = dummy.Y();
	  }
	  Xnew[4] = Xnew[0]; Ynew[4] = Ynew[0];
	  TCutG PixelCut("PixelCut",5,Xnew,Ynew);
	  Cuts.push_back(PixelCut);
	}
      }

      Float_t IntegratedEfficacy(0);
      itCuts = Cuts.begin();

      for (Int_t i = 1; i <= Psf->GetNbinsX(); i++) {
	Float_t x = Psf->GetXaxis()->GetBinCenter(i);
	for (Int_t j = 1; j <= Psf->GetNbinsY(); j++) {
	  Float_t y = Psf->GetYaxis()->GetBinCenter(j);
	  if (itCuts != Cuts.end() && (*itCuts).IsInside(x,y)) 
	    IntegratedEfficacy += Psf->GetBinContent(i,j);
	  else
	    for (itCuts = Cuts.begin(); itCuts != Cuts.end(); itCuts++) 
	      if ((*itCuts).IsInside(x,y)) {
		IntegratedEfficacy += Psf->GetBinContent(i,j);
		break;
	      }
	}
      }

      if (fDebugPs >= 2) fEnergyViewer->DisplayRecoveryDebug(Psf, IntegratedEfficacy, time+0.5*TimeStep);
      fHistoOverlap->Fill(time,IntegratedEfficacy);
      fHistoIdealOverlap->Fill(time,Psf->Integral());

      if (100.*(time-InitialTime)/TotalTime >= progress) {         
	Msg(EsafMsg::Info).SetProgress(progress);
	Msg(EsafMsg::Info) << "BuildRecoveryFunction():" << MsgCount;
	progress+=10;
      }

      SafeDelete(Psf);
    }

    if (progress <= 100) Msg(EsafMsg::Info) << MsgDispatch;
    if (fDebugPs >= 1) fEnergyViewer->DisplayRecovery();

    return kTRUE;
}


//_____________________________________________________________________________
 void EnergyModule::FitTimeDistribution() {
  // Fit temporal distribution
  // define the data point and plot them for the check

  Int_t CurrentGtuCounts, CurrentGtuCountsUnfolded;
  Pixels_t OneGtuPixels;
  Pixels_t::iterator itPix;
  GtuPixels_t::iterator

  // Currently is not used

//     for (itGtu = fGtuPixels.begin(); itGtu != fGtuPixels.end(); itGtu++) {

//          CurrentGtuCounts = 0;
//          OneGtuPixels.clear();
//          OneGtuPixels = itGtu->second;

//          for (itPix = OneGtuPixels.begin(); itPix != OneGtuPixels.end(); itPix++) {
//               CurrentGtuCounts += (*itPix)->GetCounts();
//              if (fBackground) 
// 	         CurrentGtuCounts-= fRunpars->GetNightGlowRateByUId((*itPix)->GetPixelId())* 
// 	         fEvent->GetHeader().GetDetector()->GetElectronics()->GetGtuLength()/microsecond;
//          }

//          point++;
//     }

    // make my cluster

    TH1F* htime    = new TH1F("htime","Time distribution",
			      fMyGtuMax-fMyGtuMin,fMyGtuMin,fMyGtuMax);
    TH1F* htimeS   = new TH1F("htimeS","Unsaturated time distribution",
			      fMyGtuMax-fMyGtuMin,fMyGtuMin,fMyGtuMax);
    TH1F* htimeRec = new TH1F("Recovery Time distribution","",
			      fMyGtuMax-fMyGtuMin,fMyGtuMin,fMyGtuMax);

    for (itGtu = fGtuMyPixels.begin(); itGtu != fGtuMyPixels.end(); itGtu++) {
      CurrentGtuCounts = 0;
      CurrentGtuCountsUnfolded = 0;
      OneGtuPixels.clear();
      OneGtuPixels = itGtu->second;

      for (itPix = OneGtuPixels.begin(); itPix != OneGtuPixels.end();itPix++) {
	CurrentGtuCounts += (*itPix)->GetCounts();
        CurrentGtuCountsUnfolded += UnfoldSaturation((*itPix)->GetCounts());
// 	if (fBackground) 
// 	  CurrentGtuCounts-=fRunpars->GetNightGlowRateByUId((*itPix)->GetPixelId()) * 
// 	    fEvent->GetHeader().GetDetector()->GetElectronics()->GetGtuLength()/microsecond;
      }
      htime->Fill(itGtu->first,CurrentGtuCounts);
      htimeS->Fill(itGtu->first,CurrentGtuCountsUnfolded);
    }

    // Find Chrenkov signal
    
    FindCherenkovPeak(htimeS);
    htimeRec->Divide(htimeS,fHistoOverlap);
    FindCherenkovPeak(htimeRec);

    // Generate Shower Track, fill expected histo for flux of photons on EP

    BuildShowerTrack();

    // Fitting histogram

    TF1 *func_energy = new TF1("func_energy",Fit_Energy,fMyGtuMin,fMyGtuMax,2);
    TF1 *func_time   = new TF1("func_time",TimeDistr,fMyGtuMin,fMyGtuMax,2);
    TF1 *func_bg     = new TF1("func_bg",Background,fMyGtuMin,fMyGtuMax,0);

    gMinuit->SetObjectFit(this);
    func_energy->SetParLimits(0,0,100);
    func_energy->SetParameter(0,5.);
    func_energy->SetParameter(1,22.);
    htimeS->Fit("func_energy","QEMN");

    Double_t *par = func_energy->GetParameters();
    Double_t *err = func_energy->GetParErrors();

    fEnergy      = par[0]*1.e+14;
    fEnergyError = err[0]*1.e+14;
    fEnergyChi2  = func_energy->GetChisquare()/(fMyGtuMax-fMyGtuMin-2.);

    if (fDebugPs >= 1) {

      htimeRec->SetMarkerColor(kBlue); 
      htimeRec->SetMarkerStyle(22);
      htime->SetMarkerColor(kBlack); 
      htime->SetMarkerStyle(21);
      htimeS->SetMarkerColor(kBlack); 
      htimeS->SetMarkerStyle(21);

      func_time->SetLineStyle(1);    
      func_time->SetLineColor(1);
      func_time->SetParameters(par);
      func_bg->SetLineColor(kBlue);
      func_energy->SetLineStyle(2);    
      func_energy->SetLineColor(2);
      fEnergyViewer->DisplayTimeFit(htimeRec,func_time,func_bg);
      fEnergyViewer->DisplayTimeFit(htime,func_energy,func_bg);
      fEnergyViewer->DisplayTimeFit(htimeS,func_energy,func_bg);
    }

    Msg(EsafMsg::Info) << Form("Simulated          Reconstructed Information") 
		       << MsgDispatch; 
    Msg(EsafMsg::Info) << Form("Theta  %.2E        %.2E deg",
			       fEvent->GetHeader().GetTrueTheta()*RadToDeg(),
			       fTheta*RadToDeg()) << MsgDispatch; 
    Msg(EsafMsg::Info) << Form("Phi    %.2E        %.2E deg",
			       fEvent->GetHeader().GetTruePhi()*RadToDeg(),
			       fPhi*RadToDeg()) << MsgDispatch; 
    Msg(EsafMsg::Info) << Form("Energy %.2E        %.2E MeV", 
			       fEvent->GetHeader().GetTrueEnergy(),
			       fEnergy) << MsgDispatch; 

    SafeDelete(htime);
    SafeDelete(htimeRec);
    SafeDelete(htimeS);
    SafeDelete(func_energy);
    SafeDelete(func_time);
    SafeDelete(func_bg);
}

//_____________________________________________________________________________
 void EnergyModule::FindCherenkovPeak(TH1F *htime) {

  Int_t ShowerTime, TimePeak;
  Int_t MaxSignalGtu = 0, npeak = 0, ChBin = 0, ChBinMin = 0;
  Int_t ShowerTime2  = fMyGtuMax - fMyGtuMin;
  Int_t ShowerTime1  = 0;
  Float_t *PeakPosition;


  TSpectrum *spectrum = new TSpectrum(10,1);
  Int_t PeaksNumber = spectrum->Search(htime,2,"");

  if (PeaksNumber) {
    
    // All time unit defines via bin number

    // Determine time caracteristics of event

    PeakPosition = spectrum->GetPositionX();

    // Find a reasonable time of shower evolution
 
    // Searching after last peak

    TimePeak = htime->FindBin(PeakPosition[PeaksNumber-1]);
    
    for (Int_t iTime = TimePeak+1; iTime < fMyGtuMax; iTime++) {
      if (!(htime->GetBinContent(iTime))) {
	ShowerTime2 = iTime;
	break;
      }
    }
    
    // Searching before first peak
    
    TimePeak = htime->FindBin(PeakPosition[0]);
    
    for (Int_t iTime = TimePeak; iTime > fMyGtuMin; iTime--) {
      if (!(htime->GetBinContent(iTime))) {
	ShowerTime1 = iTime;
	break;
      }
    }

    ShowerTime = ShowerTime2 - ShowerTime1; // Time of shower

    // Find maximum signal

    for (Int_t iPeak = 0; iPeak < PeaksNumber; iPeak++) {
      Int_t Time = htime->FindBin(PeakPosition[iPeak]);
      if (htime->GetBinContent(MaxSignalGtu) < htime->GetBinContent(Time)) {
	MaxSignalGtu = Time;
	npeak        = iPeak;
      }
    }

    ChBin = htime->FindBin(PeakPosition[PeaksNumber-1]);

    // Check if the maximum signal is Cherenkov light

    if (npeak == PeaksNumber-1) {
      if (1.*(MaxSignalGtu-ShowerTime1)/ShowerTime > 0.8) {
	if (PeaksNumber > 1) {
	  // Cherenkov peak founded, try to find fluorescent maximum 
	  MaxSignalGtu = 0;
	  for (Int_t iPeak = 0; iPeak < PeaksNumber-1; iPeak++) {
	    Int_t Time = htime->FindBin(PeakPosition[iPeak]);
	    if (htime->GetBinContent(MaxSignalGtu) < htime->GetBinContent(Time))
	      MaxSignalGtu = Time;
	  }
	}
	else {MaxSignalGtu = (Int_t)(0.5*(ShowerTime1+ShowerTime2));}
      }
      else {ChBin = 0;}
    }

    // Find a Cherenkov light

    if (ChBin) {

      // Position of maximum signal from Cherenkov light
            
      for (Int_t iTime = ChBin-5; iTime < ShowerTime2; iTime++)
	if (htime->GetBinContent(ChBin) < htime->GetBinContent(iTime)) 
	  ChBin = iTime;

      // Time of Cherenkov light signal on focal plane

      ChBinMin = ChBin;

      for (Int_t iTime = ChBin; iTime > ChBin - 6; iTime--)
	if (htime->GetBinContent(ChBinMin) > htime->GetBinContent(iTime)) 
	  ChBinMin = iTime;
    }
  }
  else {
    // Put the peak to the point with maximum density
    
    Float_t Sumt = 0;

    for (Int_t iTime = fMyGtuMin+2; iTime < fMyGtuMax-3; iTime++) {
      Float_t Sum = 0;
      for (Int_t iDh = iTime - 2; iDh <= iTime + 2; iDh ++)
	Sum += htime->GetBinContent(iDh);
      if (Sum > Sumt) {
	MaxSignalGtu = iTime;
	Sumt = Sum;
      }
    }
  }

  // Theta and Phi at shower maximum

  fPhiFoVAtMax   = PhiFoV(MaxSignalGtu);
  fThetaFoVAtMax = ThetaFoV(MaxSignalGtu);
  fGtuTimeAtMax  = MaxSignalGtu;

  Msg(EsafMsg::Debug) << Form("FindThetaPhiAtMaximum(): Gtu %d, Theta %3.2f deg, Phi %3.2f deg", MaxSignalGtu,fThetaFoVAtMax*RadToDeg(),fPhiFoVAtMax*RadToDeg()) << MsgDispatch;

  if (fDebugPs >= 1) {

    // Draw polymarker for Cherenkov region and shower maximum

    Int_t nPoints = 3;
    Float_t x[6], y[6];

    x[0] = fMyGtuMin + ShowerTime1;
    y[0] = 0;
    x[1] = fMyGtuMin + MaxSignalGtu-0.5;
    y[1] = htime->GetBinContent(MaxSignalGtu);
    x[2] = fMyGtuMin + ShowerTime2;
    y[2] = 0;

    if (ChBin) {
      x[3] = fMyGtuMin + ChBin-0.5;
      y[3] = htime->GetBinContent(ChBin);
      nPoints = 4;
	
      if (ChBinMin) {
	x[4] = fMyGtuMin + ChBinMin-0.5;
	y[4] = htime->GetBinContent(ChBinMin);
	x[5] = fMyGtuMin + 2*ChBin - ChBinMin - 0.5;
	y[5] = htime->GetBinContent(2*ChBin - ChBinMin);
	nPoints = 6;
      }
    }
    fEnergyViewer->DisplayCherenkov(htime,nPoints,x,y);
    }
    SafeDelete(spectrum);
    fFCtime[0] = fMyGtuMin + MaxSignalGtu-0.5;
    fFCtime[1] = fMyGtuMin + ChBinMin-0.5;
    fFCtime[2] = fMyGtuMin + ChBin-0.5;
    fFCtime[3] = fMyGtuMin + 2*ChBin - ChBinMin - 0.5;
}


//_____________________________________________________________________________
 void EnergyModule::EstimateBackground() {

    PMTs_t BgrPix;
    map <Int_t, Float_t> simBgr;
    map <Int_t,Int_t> Npix, simNpix;
    map <Int_t,TH1F*> hBgrCount;
    PMTs_t::iterator itPmt, itMyPmt;
    Pixels_t::iterator itPix;

    Npix.clear(); simBgr.clear(); simNpix.clear(); 
    BgrPix.clear(); hBgrCount.clear();

    // Make a list of pixels not in signal

    for (itPmt = fPixelsBgr.begin(); itPmt != fPixelsBgr.end(); itPmt++) {
      for (itMyPmt = fMyPixels.begin(); itMyPmt != fMyPixels.end(); itMyPmt++)
	if (*itPmt == *itMyPmt) break;
      BgrPix.push_back(*itPmt);
    }

    // create histogram

    for (itPmt = fECells.begin(); itPmt != fECells.end(); itPmt++)
      hBgrCount[*itPmt] = new TH1F(Form("Count%d",*itPmt),"Gtu",
				   fMyGtuMax-fMyGtuMin,fMyGtuMin,fMyGtuMax);

    // fill background histo

    for (itPix = fActivePixels.begin();itPix != fActivePixels.end(); itPix++) {
      for (itPmt = BgrPix.begin(); itPmt != BgrPix.end(); itPmt++)
	if (*itPmt == (*itPix)->GetPixelId()) {
	  Int_t ECuid = fRunpars->ElementaryCell((*itPix)->GetPixelId());
	  hBgrCount[ECuid]->Fill((*itPix)->GetGtu(),(*itPix)->GetCounts());
	  break;
	}
    }

    // Calculate number of pixels in each elementary cells

    for (itPmt = BgrPix.begin(); itPmt != BgrPix.end(); itPmt++)
      Npix[fRunpars->ElementaryCell(*itPmt)]++;

    // Calculate simulated background

    for (itPmt = fPixelsBgr.begin(); itPmt != fPixelsBgr.end(); itPmt++) {
      simBgr[fRunpars->ElementaryCell(*itPmt)] += 
	fRunpars->GetNightGlowRateByUId(*itPmt);
      simNpix[fRunpars->ElementaryCell(*itPmt)]++;
    }

    // calculate time interval excluding zero bin for each elementary cell

    Int_t nec = 0;
    Float_t rec = 0, sim = 0;

    for (itPmt = fECells.begin(); itPmt != fECells.end(); itPmt++) {

      Int_t Ngtu = 0;
      Float_t integral = 0;

      for (Int_t i = 1; i <= hBgrCount[*itPmt]->GetNbinsX(); i++)
	if (hBgrCount[*itPmt]->GetBinContent(i)) {
	  integral += hBgrCount[*itPmt]->GetBinContent(i);
	  Ngtu++;
	}

//     cout << "Cell " << *itPmt << " rec: " 
// 	 << hBgrCount[*itPmt]->Integral()/(1.*Npix[*itPmt]*Ngtu) 
// 	 << " sim: " << 2.5*NsimBgr[*itPmt]/NpixBgr[*itPmt] << endl;

      rec += integral/(1.*Npix[*itPmt]*Ngtu);
      sim += 2.5*simBgr[*itPmt]/simNpix[*itPmt];
      nec++;

//     hBgrCount[*itPmt]->Scale(1./Npix[*itPmt]);
//     fEnergyViewer->DrawHisto(hBgrPoisson[*itPmt]);
//     fEnergyViewer->DrawHisto(hBgrCount[*itPmt]);

      SafeDelete(hBgrCount[*itPmt]);
    }

    hBgrCount.clear(); BgrPix.clear(); simNpix.clear(); simBgr.clear();

//   Int_t np = 0, bgr = 0;
//   for (itPix = fPixelsCluster.begin();itPix != fPixelsCluster.end(); itPix++) {
//     bgr += (*itPix)->GetCountsNoise(); np++;
//   }
//   cout << "Mean value sim:       " << 1.*sim/nec << endl;
//   cout << "Simulated background: " << 1.*bgr/np << endl;
//   cout << "Reconstructed:        " << 1.*rec/nec << endl;

    fBackGround = 1.*rec/nec;
}


//_____________________________________________________________________________
 Float_t EnergyModule::BgrPerPixelPerGtu(Float_t time) {
    //
    //
    //
    return fBackGround;
}


//_____________________________________________________________________________
 void EnergyModule::BuildShowerTrack() {
    //
    //
    //
    fShowerTrack->SetEnergyModule(this);
    fShowerTrack->SetHmax(fAltitude); 
    fShowerTrack->SetThetaPhiFoV(fThetaFoVAtMax,fPhiFoVAtMax);
    fShowerTrack->SetThetaPhi(fTheta,fPhi);
    fShowerTrack->MakeTrack();   
}

//_____________________________________________________________________________
 Bool_t EnergyModule::Process(RecoEvent *ev) {

    // manager of the energy reconstruction chain

    fEvent      = ev;
    fRunpars    = fEvent->GetHeader().GetRunPars();
    gErrorIgnoreLevel = 1;
    CreateAuxilaryObjects();

    // Load information from previous modules
    if (!LoadModules()) return kFALSE;
    // Make map <GtuHit_t,Pixels_t>
    SortGtu();
    // Estimate Background
    EstimateBackground();
    // Build the recovery function along the track
    BuildRecoveryFunction();
    // Fit the time distribution
    FitTimeDistribution();
    // Save data                                    
    MyData()->Add("Energy",fEnergy);

    if (fDebugPs >= 1) fEnergyViewer->DisplayPixels3D("ClusterSurf");
    if (fDebugPs >= 1) fEnergyViewer->DisplayPixels3D("FocalSurf");
    if (fDebugPs >= 1) fEnergyViewer->CloseFilePS();
    if (fDebugPs >= 1) fEnergyViewer->Clean();

    gErrorIgnoreLevel = 0;

    fShowerTrack->Clear();
    
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t EnergyModule::PostProcess() {

     return kTRUE;
}

//_____________________________________________________________________________
 Bool_t EnergyModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
    //
    // Save desired information into the reco root file
    //
    fRecoRootEvent->GetRecoEnergy().SetEnergy(fEnergy,fEnergyError); 
    fRecoRootEvent->GetRecoEnergy().SetQuality(fEnergyChi2);
    Double_t MCError = fRecoRootEvent->GetTruth().GetEnergy() - fEnergy;
    fRecoRootEvent->GetRecoEnergy().SetEnergyMCError(MCError);
    fRecoRootEvent->GetRecoEnergy().SetThetaFoV(fThetaFoVAtMax);
    fRecoRootEvent->GetRecoEnergy().SetPhiFoV(fPhiFoVAtMax);
    
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t EnergyModule::Done() {
    // 
    // Here we inform that the module is ended its work
    // 
  
    if (fDebugPs >= 1) SafeDelete(fEnergyViewer);
    Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
    return kTRUE;
}

//_____________________________________________________________________________
 void EnergyModule::UserMemoryClean()  {
    // 
    // Here we clear any auxilary objects created during one event energy reconstruction
    //
    fPixelsCluster.clear();
    fMyPixelsCluster.clear();
    fAllActivePixels.clear();
    fActivePixels.clear();
    fShowerDir.Clear();
    fGtuPixels.clear();
    fGtuMyPixels.clear();
    fMyPixels.clear();
    fPixelsBgr.clear();
    fECells.clear();
    fPmts.clear();

    SafeDelete(fHistoTheta);
    SafeDelete(fHistoPhi);
    SafeDelete(fHistoX);
    SafeDelete(fHistoY);
    SafeDelete(fHistoXY);
    SafeDelete(fHistoOverlap);
    SafeDelete(fHistoIdealOverlap);
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005