Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

EnergyModule - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: EnergyModule.cc,v 1.24 2005/04/19 14:27:31 naumov Exp $
// Naumov Dmitry created Jul,  4 2004

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  EnergyModule                                                            //
//                                                                          //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "EnergyModule.hh"
#include "RecoEvent.hh"
#include "RecoPixelData.hh"
#include "FluoCalculator.hh"
#include "Config.hh"
#include "EsafSpectrum.hh"
#include "LightSourceFactory.hh"
#include "SinglePhoton.hh"
#include "RadiativeProcessesCalculator.hh"
#include "O1_ClearSkyPropagator.hh"
#include "RadiativeFactory.hh"
#include "EOpticsResponse.hh"
#include "EORSample.hh"
#include "EORHeader.hh"
#include "ERunParameters.hh"
#include "EDetector.hh"
#include "EConst.hh"
#include "ShowerTrack.hh"
#include "RecoRootEvent.hh"
#include "SlastShowerGenerator.hh"

#include <TH2F.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TView.h>
#include <TPolyLine.h>
#include <TStyle.h>
#include <TGraphErrors.h>
#include <TText.h>
#include <TCutG.h>
#include <TSpline.h>
#include <TMinuit.h>
#include <TFumili.h>
#include <TVirtualFitter.h>
#include <TPaveText.h>
#include <TLegend.h>

#include <map>
#include <vector>

using namespace TMath;
using namespace EConst;

ClassImp(EnergyModule)
ClassImp(EnergyModuleData)

//______________________________________________________________________________
Double_t PhotoElectronsTotal(Double_t *x, Double_t *par) {
    // time (in GTU units)
    // par[0] is fitted energy
    // par[1] is fitted time at the shower maximum
    // par[2] is fitted time sigma
    EnergyModuleData *fData = (EnergyModuleData*)gMinuit->GetObjectFit();
    if (!fData) return 0;
    Float_t time     = x[0];
    Float_t energy   = par[0];
    Float_t timemax  = par[1];
    Float_t sigma    = par[2];
    
    Float_t       fQuantumEff        = fData->fQuantumEff;
    Float_t       fAttenuation       = fData->fAttenuation;
    Float_t       fYield             = fData->fYield;
    Float_t       fShowerLength      = fData->fShowerLength;
    Float_t       fRmax              = fData->fRmax;
    TSpline3     *fPsfIdeal          = (TSpline3*)fData->fPsfIdeal;
    
    // make the energy to scale around 10**20 eV
    Float_t Theory = energy*1.e14/(1450)*fQuantumEff*fAttenuation*fYield*fShowerLength*fPsfIdeal->Eval(time);
    Theory         = Theory/( 4.e0*Pi()*fRmax*fRmax);
    Theory        *= Exp(-0.5*Power(time-timemax,2)/Power(sigma,2)); 
    return Theory;
}

//______________________________________________________________________________
Double_t PhotoElectronsRecovery(Double_t *x, Double_t *par) {
    // time (in GTU units)
    // par[0] is fitted energy
    // par[1] is fitted time at the shower maximum
    // par[2] is fitted time sigma
    EnergyModuleData *fData = (EnergyModuleData*)gMinuit->GetObjectFit();
    if (!fData) return 0;
    Float_t time     = x[0];
    Float_t energy   = par[0];
    Float_t timemax  = par[1];
    Float_t sigma    = par[2];
    
    Float_t       fQuantumEff        = fData->fQuantumEff;
    Float_t       fAttenuation       = fData->fAttenuation;
    Float_t       fYield             = fData->fYield;
    Float_t       fShowerLength      = fData->fShowerLength;
    Float_t       fRmax              = fData->fRmax;
    TSpline3     *fRecovery          = (TSpline3*)fData->fRecovery;
    
    // make the energy to scale around 10**20 eV
    Float_t Theory = energy*1.e14/(1450)*fQuantumEff*fAttenuation*fYield*fShowerLength*fRecovery->Eval(time);
    Theory         = Theory/( 4.e0*Pi()*fRmax*fRmax);
    Theory        *= Exp(-0.5*Power(time-timemax,2)/Power(sigma,2));
    return Theory;
}

//______________________________________________________________________________
Double_t PhotoElectrons2Total(Double_t *x, Double_t *par) {
    // time (in GTU units)
    // par[0] is fitted energy
    // par[1] is fitted time at the shower maximum
    EnergyModuleData *fData = (EnergyModuleData*)gMinuit->GetObjectFit();
    if (!fData) return 0;
    Float_t time       = x[0];
    Float_t energy     = par[0];
    Float_t timeshift  = par[1];
    
    Float_t       fQuantumEff        = fData->fQuantumEff;
    TSpline3     *fPsfIdeal          = (TSpline3*)fData->fPsfIdeal;
    TH1F         *fFlux              = (TH1F*)fData->fFlux;
    
    Int_t bin      = fFlux->GetXaxis()->FindBin(time-timeshift); 
    Float_t Theory = energy*fPsfIdeal->Eval(time)*fQuantumEff;
    Theory        *= fFlux->GetBinContent(bin);
    return Theory;
}

//______________________________________________________________________________
Double_t PhotoElectrons2Recovery(Double_t *x, Double_t *par) {
    // time (in GTU units)
    // par[0] is fitted energy
    // par[1] is fitted time at the shower maximum
    EnergyModuleData *fData = (EnergyModuleData*)gMinuit->GetObjectFit();
    if (!fData) return 0;
    Float_t time       = x[0];
    Float_t energy     = par[0];
    Float_t timeshift  = par[1];
    
    Float_t       fQuantumEff        = fData->fQuantumEff;
    TSpline3     *fRecovery          = (TSpline3*)fData->fRecovery;
    TH1F         *fFlux              = (TH1F*)fData->fFlux;
    
    Int_t bin      = fFlux->GetXaxis()->FindBin(time-timeshift); 
    Float_t Theory = energy*fRecovery->Eval(time)*fQuantumEff;
    Theory        *= fFlux->GetBinContent(bin);
    return Theory;
}

//______________________________________________________________________________
void PhotoElectronsTimeFit(Int_t &npar, Double_t *deriv,Double_t &f, Double_t *par, Int_t flag) {
    // time (in GTU units)
    // par[0] is fitted energy
    // par[1] is fitted time at the shower maximum
    // par[2] is fitted time sigma
    EnergyModuleData *fData = (EnergyModuleData*)gMinuit->GetObjectFit();
    if (!fData) return;
    Float_t energy   = par[0];
    Float_t timemax  = par[1];
    Float_t sigma    = par[2];
    
    Float_t       fQuantumEff        = fData->fQuantumEff;
    Float_t       fAttenuation       = fData->fAttenuation;
    Float_t       fYield             = fData->fYield;
    Float_t       fShowerLength      = fData->fShowerLength;
    Float_t       fRmax              = fData->fRmax;
    
    TGraphErrors *fPhotoElectrons    = (TGraphErrors*)fData->fMyPhotoElectrons;
    TSpline3     *fRecovery          = (TSpline3*)fData->fRecovery;
  
    Double_t Chi2(0);
    
    for (Int_t i(0); i < fPhotoElectrons->GetN(); i++) {
         Double_t time,pe;
         fPhotoElectrons->GetPoint(i,time,pe);    
         // find efficacy integrated over pixels
	 Float_t Theory = energy*1.e14/1450*fRecovery->Eval(time)*fQuantumEff*fAttenuation*fYield*fShowerLength;
         Theory         = Theory/( 4.e0*Pi()*fRmax*fRmax);
         Theory        *= Exp(-0.5*Power(time-timemax,2)/Power(sigma,2)); 
         Chi2          += Power(pe-Theory,2)/pe;
    }
    f = (Double_t)Chi2/(fPhotoElectrons->GetN() - 3);
}
//______________________________________________________________________________
void PhotoElectronsTimeFit2(Int_t &npar, Double_t *deriv,Double_t &f, Double_t *par, Int_t flag) {
    // time (in GTU units)
    // par[0] is fitted energy
    // par[1] is fitted time at the shower maximum
    EnergyModuleData *fData = (EnergyModuleData*)gMinuit->GetObjectFit();
    if (!fData) return;
    Float_t energy    = par[0];
    Float_t timeshift = par[1];
    
    Float_t       fQuantumEff        = fData->fQuantumEff;
    TGraphErrors *fPhotoElectrons    = (TGraphErrors*)fData->fMyPhotoElectrons;
    TSpline3     *fRecovery          = (TSpline3*)fData->fRecovery;
    TH1F         *fFlux              = (TH1F*)fData->fFlux;
     
    Double_t Chi2(0);
    
    for (Int_t i(0); i < fPhotoElectrons->GetN(); i++) {
         Double_t time,pe;
         fPhotoElectrons->GetPoint(i,time,pe);    
	 Int_t bin      = fFlux->GetXaxis()->FindBin(time-timeshift); 
	 Float_t Theory = energy*fRecovery->Eval(time)*fQuantumEff;
         Theory        *= fFlux->GetBinContent(bin);
         Chi2          += Power(pe-Theory,2)/pe;
    }
    f = (Double_t)Chi2/(fPhotoElectrons->GetN() - 2);
}

//_____________________________________________________________________________
void EnergyModuleData::Clear(){
    fQuantumEff        = 0;
    fAttenuation       = 0;
    fYield             = 0;
    fShowerLength      = 0;
    fRmax              = 0;
    SafeDelete(fRecovery);
    SafeDelete(fPsfIdeal);
    SafeDelete(fFlux);
    SafeDelete(fPhotoElectrons);
    SafeDelete(fMyPhotoElectrons);
}

//_____________________________________________________________________________
EnergyModule::EnergyModule() : RecoModule("Energy") {
    // ctor
    const string& ORdbname = Conf()->GetStr("EnergyModule.fORname");
    fOpticsResponse  = new EOpticsResponse(ORdbname.c_str());
    fEsafSpectrum    = new EsafSpectrum();
    fShowerGenerator = new SlastShowerGenerator;
    fEPspectrum      = NULL;
    fEvent           = NULL;
    fEPphotons       = NULL;
    fShowerTrack     = NULL;
    fCanvas          = NULL;
    fPsfOverlap      = NULL;
    fPsfIdealOverlap = NULL;
    fRecovery        = NULL;
    fPhotoElectrons  = NULL;
    fMyPhotoElectrons  = NULL;
    fBgGraph         = NULL;
    //
}

//_____________________________________________________________________________
EnergyModule::~EnergyModule() {

    // dtor
    SafeDelete(fEsafSpectrum);
    SafeDelete(fEPspectrum);
    SafeDelete(fEPphotons);
    SafeDelete(fOpticsResponse);
    SafeDelete(fShowerTrack);
    SafeDelete(fShowerGenerator);
    SafeDelete(fCanvas);
    ClearGraphs();
}

//_____________________________________________________________________________
 Bool_t EnergyModule::Init() {    
    // Get fluorescence calculator
    Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
    string fluoname         = Conf()->GetStr("EnergyModule.fFluoCalculator");
    fFluoCalculator         = LightSourceFactory::Get()->GetFluoCalculator( fluoname );    
    fRadCalculator          = RadiativeFactory::Get()->GetRadiativeProcessesCalculator();
    fEusoAlt                = Conf()->GetNum("EnergyModule.fEusoAltitude")*km;
    fEusoRadius             = Conf()->GetNum("EnergyModule.fEusoRadius")*mm;
    fSaveEps                = Conf()->GetStr("EnergyModule.fSaveEps");
    fEnergyDistributionName = Conf()->GetStr("EnergyModule.fEnergyDistributionName");
    fViewThreshold          = (Int_t)Conf()->GetNum("EnergyModule.fViewThreshold");
    fThreshold              = (Int_t)Conf()->GetNum("EnergyModule.fThreshold");
    fMinuitOutputLevel      = (Int_t)Conf()->GetNum("EnergyModule.fMinuitOutputLevel");
    fFitMode                = (Int_t)Conf()->GetNum("EnergyModule.fFitMode");
    if ( fMinuitOutputLevel < -1 || fMinuitOutputLevel > 3 ) {
         Msg(EsafMsg::Warning) << "Wrong config value EnergyModule.fMinuitOutputLevel " 
	                       << fMinuitOutputLevel << " Setted to 1" << MsgDispatch;
         fMinuitOutputLevel = 1;
    }
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t EnergyModule::PreProcess() {
    fAltitude       = -1; 
    fFluorMaxPe     = -1;
    fDeltaL         = -1;
    fEnergyChi2     = -1;
    fGtuMin         =  1.e21;
    fGtuMax         = -1.21;
    fGtuAtMax       = -1;
    fGtuTimeAtMax   = -1;
    fXmax           = -1.e20;
    fXmin           =  1.e20;
    fYmax           = -1.e20;
    fYmin           =  1.e20;
    fData.Clear();
    fGtuPixels.clear();
    fGtuMyPixels.clear();
    fPixelLines.clear();
    fPmtLines.clear();
    SafeDelete(fEPspectrum);
    return kTRUE;
}

//_____________________________________________________________________________
 void EnergyModule::CreateAuxilaryObjects() {
    ClearGraphs();
    fCanvas          = new TCanvas("EnergyModule");
    fGraphTheta      = new TGraphErrors();
    fGraphPhi        = new TGraphErrors();
    fGraphXY         = new TGraphErrors();
    fPsfOverlap      = new TGraph();
    fPsfIdealOverlap = new TGraph();
    fBgGraph         = new TGraph();
    fRecovery        = new TSpline3();
    fPhotoElectrons  = new TGraphErrors();
    fMyPhotoElectrons  = new TGraphErrors();
    fEPspectrum        = new TH1F("EPspectrum","EP",100,300*nm,400*nm);

}
//_____________________________________________________________________________
 void EnergyModule::ClearGraphs() {
    SafeDelete(fCanvas);
    SafeDelete(fGraphTheta);
    SafeDelete(fGraphPhi);
    SafeDelete(fGraphXY);
    SafeDelete(fPsfOverlap);
    SafeDelete(fPsfIdealOverlap);
    SafeDelete(fRecovery);
    SafeDelete(fPhotoElectrons);
    SafeDelete(fMyPhotoElectrons);
    SafeDelete(fBgGraph);
}

//_____________________________________________________________________________
 void EnergyModule::LoadModules() {
    // Load information from previous modules essential for the energy reconstruction
    fTDM = fEvent->GetModuleData("TrackDirection");
    
    if ( fTDM == NULL ) Msg(EsafMsg::Panic) << "TrackDirectionModule is not found." << MsgDispatch;
    
    fTheta    = fTDM->GetDouble("Theta");
    fPhi      = fTDM->GetDouble("Phi");
    fShowerDir.SetXYZ(Sin(fTheta)*Cos(fPhi), Sin(fTheta)*Sin(fPhi), Cos(fTheta));
    
    fGtuCM    = fEvent->GetModuleData("GTUClustering");
    
    if ( fGtuCM == NULL )  Msg(EsafMsg::Panic) << "GTUClusteringModule is not found." << MsgDispatch;
    
    if ( fGtuCM->GetObj("CluPixels") == NULL ) 
         Msg(EsafMsg::Panic) << "Gtu Cluster Module has NULL CluPixels" << MsgDispatch;	    
    
    fPixelsCluster.clear();
    fAllActivePixels.clear();
    fPixelsCluster = *(vector<Int_t>*) fGtuCM->GetObj("CluPixels");
    fAllActivePixels = (vector<RecoPixelData*>)fEvent->GetRecoPixels();
    
    if ( (fNumPoints = fPixelsCluster.size()) == 0) Msg(EsafMsg::Panic) << "CluPixels has no pixels" 
                                                                        << MsgDispatch;
    Msg(EsafMsg::Debug) << "Processing " << fNumPoints << " 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");

}

//_____________________________________________________________________________
 Float_t EnergyModule::MakeFluorYield() {
    // Estimate fluorescence yield for the shower maximum 
    TF2* EnergyDistribution = (TF2*)ShowerTrack::GetEnergyDistribution(fEnergyDistributionName);
    TF12 EnergyDistributionMax("EDS_name",EnergyDistribution,1.e0,"X");
    fTotalYield = fFluoCalculator->GetFluoYield(fAltitude,&EnergyDistributionMax,fEsafSpectrum);
    
    // Estimate the fluorescence spectrum at the Entrace Pupil
    Int_t NumTrials(0);
    Float_t H1   = fEusoAlt-fAltitude;
    Float_t ct2 = Power(Cos(fThetaFoVAtMax),2);
    Float_t st2 = Power(Sin(fThetaFoVAtMax),2);
    Float_t tmp1 = fEusoAlt*st2 - EarthRadius()*ct2;
    Float_t tmp2 = Power(fEusoAlt,2)*st2 - (2*EarthRadius()*fAltitude + fAltitude*fAltitude);
    Float_t z   = tmp1 + Sqrt(tmp1*tmp1 - tmp2);
    Float_t rho = H1*Tan(fThetaFoVAtMax);
    Float_t x = rho*Cos(fPhiFoVAtMax);
    Float_t y = rho*Sin(fPhiFoVAtMax);
    
    EarthVector MaximumPosition(x,y,z);
    fMaxPos = MaximumPosition;
    EarthVector ISS(0,0,fEusoAlt);
    EarthVector dir = (ISS - MaximumPosition).Unit();
    fRmax = (ISS - MaximumPosition).Mag();
    fDeltaL = EConst::C()*fEvent->GetHeader().GetGtuLength()/(1 + dir*fShowerDir);
    fMeanAttenuation = 0;
    fMeanLambdaOnEntrancePupil = 0;
    
//     O1_ClearSkyPropagator::GoToDetector(EsafSpectrum& spec,const EarthVector& pos)
    while (kTRUE) {	
        // FIXME An optimization can be done if the whole spectrum at once can be transported.
	// FIXME This is what Lowtran does but this is not fully here
	NumTrials++; 
	Float_t lambda = fEsafSpectrum->GetLambda();	    
	SinglePhoton s(0, lambda, MaximumPosition, dir, (PhotonType)1, (PhotonStatus)1);
	fMeanAttenuation += fRadCalculator->Trans(s,ISS);
	fMeanLambdaOnEntrancePupil += lambda; 
	fEPspectrum->Fill(lambda);
	if (NumTrials > 100) break;
     }
     
    fMeanAttenuation           /= NumTrials;
    fMeanLambdaOnEntrancePupil /= NumTrials; 
    Msg(EsafMsg::Debug) <<  Form("MakeFluorYield(): Begin Dump") << MsgDispatch;  
    Msg(EsafMsg::Debug) <<  Form(" Estimated ShowerMax(x,y,z):        (%2.3f, %2.3f, %2.3f) km",x/km,y/km,z/km) 
                        << MsgDispatch;     
    Msg(EsafMsg::Debug) <<  Form(" Estimated Dir To Euso (Theta,Phi): (%2.3f, %2.3f) deg",    
                                   dir.Theta()*RadToDeg(),dir.Phi()*RadToDeg()) << MsgDispatch;     
    Msg(EsafMsg::Debug) <<  Form(" Estimated Mean Attenuation:        (%2.3f)",fMeanAttenuation) << MsgDispatch; 
    Msg(EsafMsg::Debug) <<  Form(" Estimated Mean Lambda on EP:       (%2.3f) nm",
                                   fMeanLambdaOnEntrancePupil/nm) << MsgDispatch; 
    Msg(EsafMsg::Debug) <<  Form(" Estimated Mean Yield:              (%2.3f) ph/m",fTotalYield*m) 
                        << MsgDispatch;				   
    Msg(EsafMsg::Debug) <<  Form("MakeFluorYield(): End Dump") << MsgDispatch;  

    return fTotalYield;
}

//_____________________________________________________________________________
 void EnergyModule::DefinePixelsView() {
    // returnes vector of RecoPixelData which should be considered in both Viewer and Energy reconstruction
    PMTs_t pmts;      
    //PMTs_t::iterator   itPmt;
    //Pixels_t::iterator  itPixel;
    
    // make list of PMTs which containes found pixels
    // find box with found cluster
    for (UInt_t i (0); i < fPixelsCluster.size(); i++) {
         RecoPixelData *pix = fEvent->GetRecoPixelData(fPixelsCluster[i]);
	 Int_t pmtId = fRunpars->Pmt(pix->GetPixelId());
         for(Int_t iCorner(0); iCorner<4; iCorner++){
             TVector3 dummy = fRunpars->PmtCorner(pmtId,iCorner);
             if (fXmin > dummy.X()) fXmin = dummy.X();
	     if (fYmin > dummy.Y()) fYmin = dummy.Y();
	     if (fXmax < dummy.X()) fXmax = dummy.X();
	     if (fYmax < dummy.Y()) fYmax = dummy.Y();
        }
    }    
    // increase size of the box by OffSet
    Float_t OffSet = 10; // mm
    fXmin -= OffSet;
    fYmin -= OffSet;
    fXmax += OffSet;
    fYmax += OffSet;
    // find pmts in this box
    for( Int_t i(1); i <= fRunpars->GetNumPmts(); i++) {
         TVector3 dummy = fRunpars->GetPmtGeo(i)->GetCenter();
         if ( dummy.X() < fXmin || dummy.X() > fXmax || dummy.Y() < fYmin || dummy.Y() > fYmax ) continue;
         // check if this pmt is not yet there
	 UInt_t pmtsize = pmts.size();
	 Bool_t IsPmtThere = kFALSE;
	 if (pmtsize) {
	     for (UInt_t j(0); j<pmtsize; j++) {
	          if (pmts[i] == i) {
		      IsPmtThere = kTRUE;
		      break;
		  }
	     }
         }
	 if (!IsPmtThere) pmts.push_back(i); 

    }


    // make list of lines to draw arond each PMT
    fPmtLines.clear();
    fPixelLines.clear();
    for (UInt_t i(0); i < pmts.size(); i++) {
        Float_t x[5], y[5];
	Int_t pmtUid = pmts[i];
        for(Int_t iCorner(0); iCorner<4; iCorner++){
            TVector3 dummy = fRunpars->PmtCorner(pmtUid,iCorner);
            x[iCorner] = dummy.X();
            y[iCorner] = dummy.Y();
            if (fXmin > x[iCorner]) fXmin = x[iCorner];
	    if (fYmin > y[iCorner]) fYmin = y[iCorner];
	    if (fXmax < x[iCorner]) fXmax = x[iCorner];
	    if (fYmax < y[iCorner]) fYmax = y[iCorner];
	    
        }
        x[4] = x[0];
        y[4] = y[0];

        TPolyLine *p = new TPolyLine(5,x,y); 
	p->SetLineColor(20);
        fPmtLines[pmtUid] = p;
	
	// for each pmt to draw prepare the list of pixels inside the given pmt
	const EPmtGeo* pmt = fRunpars->GetPmtGeo(pmtUid);
	for (Int_t PixUid = pmt->GetStartUniqueId(); PixUid <= pmt->GetLastUniqueId(); PixUid++) {
             Float_t x[5], y[5];
             for(Int_t iCorner(0); iCorner<4; iCorner++){
                 TVector3 dummy = fRunpars->PixelCorner(PixUid,iCorner);
                 x[iCorner] = dummy.X();
                 y[iCorner] = dummy.Y();
	    
              }
             x[4] = x[0];
             y[4] = y[0];

             TPolyLine *p = new TPolyLine(5,x,y); 
	     p->SetLineColor(kWhite);
	     Color_t color = 20;
	     Style_t style = 3490; // special style for defaut
	     // search if this pixel belongs to active pixels
	     for (UInt_t j(0); j<fAllActivePixels.size(); j++){
	          if (PixUid == fAllActivePixels[j]->GetPixelId()) {
		      // make yellow color for active pixel with number of counts above the ViewThreshold
		      if   (fAllActivePixels[j]->GetCounts() > fViewThreshold) color = 20;  
		      // make grey color for active pixel with number of counts below the ViewThreshold
		      else                                                     color = 20; 
		      // define fill style 
                      for (UInt_t k(0); k<fPixelsCluster.size(); k++){ 
		           RecoPixelData* pixel = fEvent->GetRecoPixelData(fPixelsCluster[k]); 
	                   if (pixel->GetPixelId() == PixUid) {
		               style = 1001; // solid fill for the cluster
			       break;
		           }
		      }
		      break;
		  }
		  else
		      // make black color for not active pixel
		      color = 1; 
	     } 
	     p->SetFillColor(color); 
             fPixelLines[PixUid] = p;
	}
	
    }
 
}
//_____________________________________________________________________________
 void EnergyModule::DisplayCluster() {
    DefinePixelsView();
    if (fCanvas) fCanvas->Clear();
    
    if (gPad) {
        Float_t OffSet = 5; 
        gPad->Range(fXmin-OffSet,fYmin-OffSet,fXmax+OffSet,fYmax+OffSet);
    } 
    
    // Draw PMTs 
    map <Int_t, TPolyLine* >::const_iterator itPmt;
    for (itPmt = fPmtLines.begin(); itPmt != fPmtLines.end(); itPmt++) {
	 itPmt->second->Draw("f");
	 itPmt->second->Draw();
    }
    // compute number of p.e. integrated over each pixel and draw the text in the center of each pixel.
    map <Int_t,Int_t> mapPixelsIntegral;
    for (UInt_t i(0) ; i<fAllActivePixels.size(); i++) {
         RecoPixelData *pix = fAllActivePixels[i];
	 mapPixelsIntegral[pix->GetPixelId()] += pix->GetCounts();
    }    
    
    // Draw Pixels 
    for (itPmt = fPixelLines.begin(); itPmt != fPixelLines.end(); itPmt++) {
	 itPmt->second->Draw("f");
	 itPmt->second->Draw();
	 Int_t counts = mapPixelsIntegral[itPmt->first];
	 TVector3 PixelCenter = fRunpars->PixelCenter(itPmt->first);
	 Float_t x = PixelCenter.X() + 0.2*fRunpars->GetPmtData().GetPadSide();
	 Float_t y = PixelCenter.Y() + 0.2*fRunpars->GetPmtData().GetPadSide();
	 TText *text = new TText(x,y,Form("%d",counts));
         Float_t padMin = Min( gPad->GetWh(),gPad->GetWw() );
	 Float_t textsize = 3.e-5*padMin;
	 text->SetTextSize(textsize);
	 text->Draw();
    }
    // Draw fitted line track
    if (fGraphXY)
        fGraphXY->Draw("PX");
	
    if (fSaveEps == "yes")  
        if(fCanvas) fCanvas->Print(Form("output/check_focalplane.run%d.event%d.eps",
	                     fEvent->GetHeader().GetRun(),fEvent->GetHeader().GetNum()));

}
//_____________________________________________________________________________
 void EnergyModule::SortGtu() {
    // Prepare map <GtuHit_t,Pixels_t> 
    fGtuPixels.clear();
    Pixels_t OneGtuPixels;
    
    for (UInt_t i(0); i < fPixelsCluster.size(); i++) {
         RecoPixelData *pix = fEvent->GetRecoPixelData(fPixelsCluster[i]);
	 if (fGtuMin > pix->GetGtu())  fGtuMin = pix->GetGtu();
	 if (fGtuMax < pix->GetGtu())  fGtuMax = pix->GetGtu();
    }
       
    Int_t previousGtu = fEvent->GetRecoPixelData(fPixelsCluster[0])->GetGtu();
    for (UInt_t i(0); i < fPixelsCluster.size(); i++) {
         RecoPixelData *pix = fEvent->GetRecoPixelData(fPixelsCluster[i]);
	 if (pix->GetGtu() == previousGtu) {
	     OneGtuPixels.push_back(i);
	 }
         else {
             fGtuPixels[previousGtu] = OneGtuPixels;
     	     previousGtu = pix->GetGtu();
	     OneGtuPixels.clear();
	     OneGtuPixels.push_back(i);
	 }	 
    }
    // make another map of pixels from all active pixels 
    fGtuMyPixels.clear();
    previousGtu = fEvent->GetRecoPixelData(fPixelsCluster[0])->GetGtu();
    for (UInt_t i(0); i < fAllActivePixels.size(); i++) {
         RecoPixelData *pix = fAllActivePixels[i];
	 if (pix->GetGtu() == previousGtu) {
	     OneGtuPixels.push_back(i);
	 }
         else {
             fGtuMyPixels[previousGtu] = OneGtuPixels;
     	     previousGtu = pix->GetGtu();
	     OneGtuPixels.clear();
	     OneGtuPixels.push_back(i);
	 }	 
    }
    DumpGtuSorted();
}

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

    Msg(EsafMsg::Info)  << Form("DumpGtuSorted(): Dumping maps <GtuHit_t,Pixels_t> into Output File") 
                        << MsgDispatch;
    Msg(EsafMsg::Debug) << Form("DumpGtuSorted(): Dumping maps <GtuHit_t,Pixels_t>") << MsgDispatch;
    Msg(EsafMsg::Debug) << Form("DumpGtuSorted(): Dumping TDM map <GtuHit_t,Pixels_t>") << MsgDispatch;
    GtuPixels_t::iterator it;
    Pixels_t OneGtuPixels;
    Int_t MaxCounts(0);
    for (it = fGtuPixels.begin(); it != fGtuPixels.end(); it++) {    
         OneGtuPixels.clear();
	 OneGtuPixels = it->second;
	 Msg(EsafMsg::Debug) << Form("GTU Entry: %d vector of %d pixels:",it->first,OneGtuPixels.size())
	                     << MsgDispatch;
	 Int_t sum(0);		     
	 for (UInt_t i(0); i < OneGtuPixels.size(); i++) {
 	      RecoPixelData* pix = fEvent->GetRecoPixelData(fPixelsCluster[OneGtuPixels[i]]);
	      Msg(EsafMsg::Debug) 
	      << Form("          GTU %d PixelUid %d Counts %d",pix->GetGtu(),pix->GetPixelId(),pix->GetCounts()) 
	      << MsgDispatch;
	      sum += pix->GetCounts();
         }
	 Msg(EsafMsg::Debug) << Form("GTU Entry: Total counts %d",sum) << MsgDispatch;
	 if (MaxCounts < sum) {
	     MaxCounts = sum;
	     fGtuTimeAtMax = it->first; 
	 }
	 
    }
    Msg(EsafMsg::Debug) << Form("Most numerous GTU %d with %d count",Int_t(fGtuTimeAtMax),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 (it = fGtuMyPixels.begin(); it != fGtuMyPixels.end(); it++) {    
         OneGtuPixels.clear();
	 OneGtuPixels = it->second;
	 Msg(EsafMsg::Debug) << Form("GTU Entry: %d vector of %d pixels:",it->first,OneGtuPixels.size())
	                     << MsgDispatch;
	 Int_t sum(0);		     
	 for (UInt_t i(0); i < OneGtuPixels.size(); i++) {
 	      RecoPixelData* pix = fAllActivePixels[OneGtuPixels[i]];
	      Msg(EsafMsg::Debug) 
	      << Form("          GTU %d PixelUid %d Counts %d",pix->GetGtu(),pix->GetPixelId(),pix->GetCounts()) 
	      << MsgDispatch;
	      sum += pix->GetCounts();
         }
	 Msg(EsafMsg::Debug) << Form("GTU Entry: Total counts %d",sum) << MsgDispatch;
	 if (MaxCounts < sum) {
	     MaxCounts = sum;
	     fGtuTimeAtMax = it->first; 
	 }
	 
     }
     Msg(EsafMsg::Debug) << Form("Most numerous GTU %d with %d count",Int_t(fGtuTimeAtMax),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;
     
   
}

//_____________________________________________________________________________
 void EnergyModule::BuildFoVAngleFunctions() {
    // 
    // make two maps 
    Int_t pixels(0);
    GtuPixels_t::iterator it;
    Pixels_t OneGtuPixels;
    for (it = fGtuPixels.begin(); it != fGtuPixels.end(); it++) {    
         OneGtuPixels.clear();
	 OneGtuPixels = it->second;
	 for (UInt_t i(0); i < OneGtuPixels.size(); i++) {
 	      RecoPixelData* pix = fEvent->GetRecoPixelData(fPixelsCluster[OneGtuPixels[i]]);
	      fGraphTheta->SetPoint(pixels,it->first,pix->GetTheta()*RadToDeg());
	      fGraphTheta->SetPointError(pixels,0,Sqrt(1.e0/pix->GetCounts()));
	      fGraphPhi->SetPoint(pixels,it->first,pix->GetPhi()*RadToDeg());
	      fGraphPhi->SetPointError(pixels,0,Sqrt(1.e0/pix->GetCounts()));
	      TVector3 PixelCenter = fRunpars->PixelCenter(pix->GetPixelId());
	      fGraphXY->SetPoint(pixels,PixelCenter.X(), PixelCenter.Y());
	      fGraphXY->SetPointError(pixels,0,Sqrt(1.e0/pix->GetCounts()));
	      pixels++;
         }
    }    
    fCanvas->Clear();
    fCanvas->Divide(1,2);
    fCanvas->cd(1);
    fGraphTheta->Fit("pol1","Q");
    fGraphTheta->Draw("A*X");
    fGraphTheta->GetYaxis()->SetTitle("#Theta, [deg]");
    fCanvas->cd(2);
    fGraphPhi->Fit("pol1","Q");
    fGraphPhi->Draw("A*X");
    fGraphPhi->GetYaxis()->SetTitle("#Phi, [deg]");
    fGraphPhi->GetXaxis()->SetTitle("Gtu Hit");
     
    fGraphXY->Fit("pol1","Q");
    fGraphXY->SetMarkerColor(kRed);
    fGraphXY->SetMarkerStyle(kMultiply);
    fGraphXY->SetMarkerSize(1.e-3*Min( gPad->GetWh(),gPad->GetWw() ));
    fGraphXY->GetFunction("pol1")->SetLineStyle(2);
    fGraphXY->GetFunction("pol1")->SetLineWidth(2);
    fGraphXY->GetFunction("pol1")->SetLineColor(kBlue);
    
    if (fSaveEps == "yes")  
        if(fCanvas) fCanvas->Print(Form("output/check_fits.run%d.event%d.eps",
	                     fEvent->GetHeader().GetRun(),fEvent->GetHeader().GetNum()));
     
}
//_____________________________________________________________________________
 void EnergyModule::FindThetaPhiAtMaximum() {
    Int_t MaxCounts(0);
    GtuPixels_t::iterator it;
    Pixels_t OneGtuPixels;
    for (it = fGtuMyPixels.begin(); it != fGtuMyPixels.end(); it++) {    
         OneGtuPixels.clear();
	 OneGtuPixels = it->second;
	 Int_t CurrentGtuCounts(0);
	 Float_t CurrentGtuTheta(0);
	 Float_t CurrentGtuPhi(0);
	 for (UInt_t i(0); i < OneGtuPixels.size(); i++) {
//  	      RecoPixelData* pix = fEvent->GetRecoPixelData(fPixelsCluster[OneGtuPixels[i]]);
 	      RecoPixelData* pix = fAllActivePixels[OneGtuPixels[i]];
	      CurrentGtuCounts += pix->GetCounts();
	      CurrentGtuTheta  += pix->GetTheta()*pix->GetCounts();
	      CurrentGtuPhi    += pix->GetPhi()*pix->GetCounts();
         }
	 if (MaxCounts < CurrentGtuCounts) {
	     MaxCounts = CurrentGtuCounts;
	     fThetaFoVAtMax = CurrentGtuTheta/CurrentGtuCounts;
	     fPhiFoVAtMax   = CurrentGtuPhi/CurrentGtuCounts;
	 }
    }
    Msg(EsafMsg::Debug) << Form("FindThetaPhiAtMaximum(): Theta %3.2f deg, Phi %3.2f deg",
                                 fThetaFoVAtMax*RadToDeg(),fPhiFoVAtMax*RadToDeg()) << MsgDispatch;
   
}

//_____________________________________________________________________________
 Float_t EnergyModule::ThetaFoV(Float_t time) {
    //
    if (!fGraphTheta)                     return 0;
    Float_t fThetaFoVP0 = fGraphTheta->GetFunction("pol1")->GetParameter(0);
    Float_t fThetaFoVP1 = fGraphTheta->GetFunction("pol1")->GetParameter(1);
    return DegToRad()*(fThetaFoVP0 + time*fThetaFoVP1);
}

//_____________________________________________________________________________
 Float_t EnergyModule::PhiFoV(Float_t time) {
    //
    if (!fGraphPhi)                       return 0;
     // Get fit parameters to set up the functions
    Float_t fPhiFoVP0   = fGraphPhi->GetFunction("pol1")->GetParameter(0);
    Float_t fPhiFoVP1   = fGraphPhi->GetFunction("pol1")->GetParameter(1);
    return DegToRad()*(fPhiFoVP0 + time*fPhiFoVP1);
}

//_____________________________________________________________________________
 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.
    Int_t nSteps = 4*Int_t(fGtuMax  - fGtuMin); 
    Int_t CurrentStep = 0;
    Int_t GtuOffSet = 2; // extend the track time to (GtuMin - DeltaGtu*GtuOffSet, GtuMax + DeltaGtu*GtuOffSet)
    Float_t TimeStep = (fGtuMax  - fGtuMin + 2*GtuOffSet*fEvent->GetHeader().GetGtuLength()/microsecond)/nSteps;
    Float_t time = fGtuMin - GtuOffSet;
    map <Int_t, TPolyLine* >::const_iterator itPixels;
    Msg(EsafMsg::Info) << "BuildRecoveryFunction(): 0%" << MsgFlush;
    Bool_t next(0),subnext(0);
    Float_t left(0),right(10),subleft(0),subright(2.5);
    
    while (CurrentStep<=nSteps) {
         Float_t fThetaFoV = ThetaFoV(time);
	 Float_t fPhiFoV   = PhiFoV(time);
// 	 TH2F* Psf = fOpticsResponse->MakePsfHist(fThetaFoV,fPhiFoV,fEPspectrum);
         TH2F* Psf = fOpticsResponse->MakePsfHist(fThetaFoV,fPhiFoV,fMeanLambdaOnEntrancePupil);
	 Float_t IntegratedEfficacy(0);
         for (itPixels = fPixelLines.begin(); itPixels!=fPixelLines.end(); itPixels++) {
/*	      Int_t entry=-1;
	      for (UInt_t j(0); j<fPixelsCluster.size(); j++){
	           RecoPixelData* pix = fEvent->GetRecoPixelData(fPixelsCluster[j]);
		   if (itPixels->first == pix->GetPixelId()) {
		       entry = j;
		       break;
	           }
              }
	      if (entry == -1)                                       continue;
	      RecoPixelData* pix = fEvent->GetRecoPixelData(fPixelsCluster[entry]);
              if (pix->GetCounts() < fThreshold)                     continue;*/
              TPolyLine *p = itPixels->second;
	      Double_t *x = p->GetX();
              Double_t *y = p->GetY();	      
	      TCutG *PixelCut = new TCutG("PixelCut",p->GetN(),p->GetX(),p->GetY());
	      // find min and max of x and y to save time needed for integration
	      Float_t xmin(1.e20), xmax(-1.e20), ymin(1.e20), ymax(-1.e20);
	      for (Int_t i(0); i < p->GetN(); i++) {
	           if (xmin > x[i]) xmin = x[i];
		   if (xmax < x[i]) xmax = x[i];
	           if (ymin > y[i]) ymin = y[i];
		   if (ymax < y[i]) ymax = y[i];
	      }
	      Int_t xbinFirst = Psf->GetXaxis()->FindBin(xmin);
	      Int_t xbinLast  = Psf->GetXaxis()->FindBin(xmax);
	      Int_t ybinFirst = Psf->GetYaxis()->FindBin(ymin);
	      Int_t ybinLast  = Psf->GetYaxis()->FindBin(ymax);
	      
	      for (Int_t i = xbinFirst; i <= xbinLast; i++) {
	           Float_t x = Psf->GetXaxis()->GetBinCenter(i);
	           for (Int_t j = ybinFirst; j <=ybinLast; j++ ) {
	                Float_t y = Psf->GetYaxis()->GetBinCenter(j);
			if (PixelCut->IsInside(x,y)) {  
			    IntegratedEfficacy += Psf->GetBinContent(i,j);
			}  
		   }
	      }
	 }
	 // end of integration
	 fPsfOverlap->SetPoint(CurrentStep,time,IntegratedEfficacy);
	 fPsfIdealOverlap->SetPoint(CurrentStep,time,Psf->Integral());
	 SafeDelete(Psf);
	 time += TimeStep;
         CurrentStep++;
         Float_t fraction = 1.e2*CurrentStep/nSteps;
         if (left<fraction && fraction <=right) 
             next = kFALSE;
         else if (fraction>right) {
             next = kTRUE;
             left = right;
             right += 10;
         }
         if (subleft<fraction  && fraction <=subright) {
             subnext = kFALSE;
         }
         if (fraction > subright) {
             subnext = kTRUE;
             subleft = subright;
             subright +=2;
         }
         if (next) Msg(EsafMsg::Info) << left << "%" << MsgFlush;
         if (subnext) Msg(EsafMsg::Info)<< "." << MsgFlush;	 
    }
    Msg(EsafMsg::Info) << " -done. " << MsgDispatch;
    if(fCanvas) fCanvas->Clear();
    fPsfOverlap->Draw("A*X");
    fPsfOverlap->GetYaxis()->SetTitle("Psf Overlap with pixels, [mm#^2]");
    fPsfOverlap->GetXaxis()->SetTitle("Gtu Hit");
    SafeDelete(fRecovery);
    SafeDelete(fEffIdeal);
    fRecovery = new TSpline3("Recovery",fPsfOverlap->GetX(),fPsfOverlap->GetY(),fPsfOverlap->GetN());
    fRecovery->Draw("same");
    fEffIdeal = new TSpline3("Recovery",fPsfIdealOverlap->GetX(),fPsfIdealOverlap->GetY(),fPsfIdealOverlap->GetN());
    fEffIdeal->SetLineColor(kRed);
    fEffIdeal->Draw("same");
    if (fSaveEps == "yes")  
        if(fCanvas) fCanvas->Print(Form("output/check_recovery.run%d.event%d.eps",
	                 fEvent->GetHeader().GetRun(),fEvent->GetHeader().GetNum()));
    return kTRUE;
}
//_____________________________________________________________________________
 void EnergyModule::FillDataStructure() {
    // Fill data container object fData of type EnergyModuleData
    fData.fQuantumEff     = fQuantumEff;
    fData.fAttenuation    = fMeanAttenuation;
    fData.fYield          = fTotalYield;
    fData.fRmax           = fRmax;
    fData.fShowerLength   = fDeltaL;
    fData.fRecovery       = (TSpline3*)fRecovery->Clone();
    fData.fPhotoElectrons = (TGraphErrors*)fPhotoElectrons->Clone();
    fData.fMyPhotoElectrons = (TGraphErrors*)fMyPhotoElectrons->Clone();
    fData.fPsfIdeal       = (TSpline3*)fEffIdeal->Clone();
    fData.fFlux           = (TH1F*)fEPphotons->Clone();
}

//_____________________________________________________________________________
 Bool_t EnergyModule::ShowerTrack3D()  {
    // 
    // make the first point and get the shower track for energy of 1.e21 and reconstructed angles
    Float_t Length = 30*km; 
    EarthVector FirstPoint = fMaxPos + fShowerDir*Length;
    fShowerGenerator->SetShower(1.e21,fTheta,fPhi,FirstPoint);
    fShowerGenerator->SetQuiet(kTRUE);
    fShowerGenerator->SetDepthStep(0.5);
    fShowerTrack = (ShowerTrack*)fShowerGenerator->Get();
    
    // make estimated time (in GTU) distribution of photons on Entrance Pupil
    EarthVector ISS(0,0,fEusoAlt);
    Float_t tMin(1.e20), tMax(-1.e20);
    Double_t tempData[fShowerTrack->Size()][2];
    Double_t FluxMax(-1), tFluxMax(-1);
    for (UInt_t i(0); i < fShowerTrack->Size(); i++) {
         const ShowerStep& s = fShowerTrack->GetStep(i);
	 // get position and compute for it: fluo yield + attenuation + time at EUSO
	 EarthVector pos  = 0.5*(s.GetXYZi() + s.GetXYZf());
	 Float_t     step_length = (s.GetXYZi() - s.GetXYZf()).Mag();
	 Float_t     time = 0.5*(s.GetTimei() + s.GetTimef());
         EarthVector dir  = (ISS - pos);
	 Double_t    SolidAngle = 1/dir.Mag2()/(4*Pi());
	 Double_t    flux  = SolidAngle*s.GetNelectrons()*fMeanAttenuation*fTotalYield*step_length;
	 Float_t     EusoTime = time + dir.Mag()/EConst::C();
	 tempData[i][0] = EusoTime;
	 tempData[i][1] = flux;
	 if (tMin > EusoTime) tMin = EusoTime;
	 if (tMax < EusoTime) tMax = EusoTime;   
	 if (FluxMax < flux)  {
	     FluxMax  = flux;
	     tFluxMax = EusoTime;
	 }
    }
    SafeDelete(fEPphotons);
    // shift time distribution to have maximum close to that by p.e. distribution
    Float_t fGtuLength = fEvent->GetHeader().GetDetector()->GetElectronics().GetGtuLength();
    tMax -= (tFluxMax - fGtuLength*fGtuTimeAtMax);
    tMin -= (tFluxMax - fGtuLength*fGtuTimeAtMax);   
    Float_t bins = (tMax - tMin)/fGtuLength;
    fEPphotons = new TH1F("EPphotons","",Int_t(bins),tMin/fGtuLength,tMax/fGtuLength);
    for (UInt_t i(0); i < fShowerTrack->Size(); i++) {
         Float_t GtuTime = (tempData[i][0]-(tFluxMax - fGtuLength*fGtuTimeAtMax))/fGtuLength;
	 Float_t photons = tempData[i][1];
         fEPphotons->Fill(GtuTime,photons);
    }
    
    if (fSaveEps == "yes")  {
        if(fCanvas) {
	   fCanvas->Clear();
           fEPphotons->Draw();
	   fCanvas->Print(Form("output/check_shower_photons.run%d.event%d.eps",
	                 fEvent->GetHeader().GetRun(),fEvent->GetHeader().GetNum()));
	}
    }
    fShowerTrack->Reset();
    fShowerGenerator->Reset();
    return kTRUE;
}

//_____________________________________________________________________________
 void EnergyModule::FitTimeDistribution() {
    // Fit temporal distribution
    // define the data point and plot them for the check
    GtuPixels_t::iterator it;
    Pixels_t OneGtuPixels;
    Int_t point(0);
    for (it = fGtuPixels.begin(); it != fGtuPixels.end(); it++) {    
         OneGtuPixels.clear();
	 OneGtuPixels = it->second;
	 Int_t   CurrentGtuPixels(0);
	 Float_t CurrentGtuCounts(0);
	 for (UInt_t i(0); i < OneGtuPixels.size(); i++) {
 	      RecoPixelData* pix = fEvent->GetRecoPixelData(fPixelsCluster[OneGtuPixels[i]]);
	      if (pix->GetCounts() >= fThreshold) {
	          CurrentGtuCounts += pix->GetCounts();
		  CurrentGtuPixels++;
	      }
         }
	 // substruct the background
// 	 CurrentGtuCounts -= fBackground*CurrentGtuPixels;
	 fPhotoElectrons->SetPoint(point,it->first,CurrentGtuCounts);
	 fPhotoElectrons->SetPointError(point,0,Sqrt(1.e0/CurrentGtuCounts));
	 point++;
    }
    // make my cluster
    point = 0;
    for (it = fGtuMyPixels.begin(); it != fGtuMyPixels.end(); it++) {    
         OneGtuPixels.clear();
	 OneGtuPixels = it->second;
	 Int_t   CurrentGtuPixels(0);
	 Float_t CurrentGtuCounts(0);
	 for (UInt_t i(0); i < OneGtuPixels.size(); i++) {
	      RecoPixelData* pix = fAllActivePixels[OneGtuPixels[i]];
	      if (pix->GetCounts() >= fThreshold) {
	          CurrentGtuCounts += pix->GetCounts();
		  CurrentGtuPixels++;
	      }
         }
	 // substruct the background
// 	 CurrentGtuCounts -= fBackground*CurrentGtuPixels;
	 fMyPhotoElectrons->SetPoint(point,it->first,CurrentGtuCounts);
	 fMyPhotoElectrons->SetPointError(point,0,Sqrt(1.e0/CurrentGtuCounts));
	 point++;
    }
    
    if(fCanvas) fCanvas->Clear();
    fMyPhotoElectrons->SetMarkerColor(kBlue);
    fMyPhotoElectrons->SetMarkerStyle(22);
    fPhotoElectrons->SetMarkerStyle(21);
    fPhotoElectrons->SetMarkerColor(kBlack);
    fMyPhotoElectrons->Draw("APE");
    fPhotoElectrons->Draw("PE");
    fPhotoElectrons->GetYaxis()->SetTitle("p.e.");
    fPhotoElectrons->GetXaxis()->SetTitle("Gtu Hit");
    TLegend *legend=new TLegend(0.1,0.85,0.3,1);
    legend->SetTextFont(72);
    legend->SetTextSize(0.02);
    legend->AddEntry(fMyPhotoElectrons,"all pixels","lp");
    legend->AddEntry(fPhotoElectrons,"TDM cluster","lp");
    legend->Draw();
    
    // make the fit
    FillDataStructure();
    if (fFitMode == 1) FitMode1();
    if (fFitMode == 2) FitMode2();
    if (fFitMode == 3) {FitMode1(); FitMode2();}
    if (fFitMode == 4) FitMode2fumili();
    
}

//_____________________________________________________________________________
 void EnergyModule::FitMode1() {
    TMinuit *minuit = new TMinuit(3);
    minuit->SetPrintLevel(fMinuitOutputLevel);
    minuit->SetObjectFit(&fData);
    minuit->SetFCN(PhotoElectronsTimeFit);
    minuit->SetErrorDef(1.);
    Double_t par[3] = {1.e0,10,10};
    Double_t step[3] = {1.e-1,1.e-1,1.e-1};
    Double_t min[3] = {0.00, 0.0, 0.0};
    Double_t max[3] = {0.00, 1.e3, 1.e3};
    string cpar[3] = {"Energy","tMax","Sigma"};
    
    for( Int_t i=0; i<3; i++ ) {
        minuit->DefineParameter(i,cpar[i].c_str(),par[i],step[i],min[i],max[i]);
    }
    minuit->Migrad();
    
    // retrieve results
    Double_t outpar[3],err[3];
    for( Int_t i=0; i<3; i++ ) {
        minuit->GetParameter(i,outpar[i],err[i]);
    }
    fEnergy = outpar[0]*1e14; 
    fEnergyError = err[0]*1e14;
    Double_t amin,edm,errdef;
    Int_t nvpar,nparx,icstat;
    minuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
    fEnergyChi2  = amin/nvpar;
        
    TF1 *f1 = new TF1("f1",PhotoElectronsTotal,fEvent->GetStartGtu(),fEvent->GetEndGtu(),3);
    TF1 *f2 = new TF1("f2",PhotoElectronsRecovery,fEvent->GetStartGtu(),fEvent->GetEndGtu(),3);

    f1->SetParameters(outpar);
    f2->SetParameters(outpar);
    f1->Draw("same");
    f2->Draw("same");
    f1->SetLineStyle(1);   f1->SetLineColor(1);
    f2->SetLineStyle(2);   f2->SetLineColor(2);
    SaveFitInfo();
    SafeDelete(f1);
    SafeDelete(f2);
    SafeDelete(minuit);
}
//_____________________________________________________________________________
 void EnergyModule::FitMode2() {
    TMinuit *minuit = new TMinuit(2);
    minuit->SetPrintLevel(fMinuitOutputLevel);
    minuit->SetObjectFit(&fData);
    minuit->SetFCN(PhotoElectronsTimeFit2);
    minuit->SetErrorDef(1.);
    Double_t par[2] = {1.e-1,0.0};
    Double_t step[3] = {1.e-1,1.e0};
    Double_t min[3] = {0.00, 0.00};
    Double_t max[3] = {0.00, 0.00};
    string cpar[3] = {"Energy","tMax"};
    
    for( Int_t i=0; i<2; i++ ) {
        minuit->DefineParameter(i,cpar[i].c_str(),par[i],step[i],min[i],max[i]);
    }
//     minuit->FixParameter(1);
    minuit->Migrad();
    // retrieve results
    Double_t outpar[2],err[2];
    for( Int_t i=0; i<2; i++ ) {
        minuit->GetParameter(i,outpar[i],err[i]);
    }
    fEnergy = outpar[0]*1e15; 
    fEnergyError = err[0]*1e15;
    Double_t amin,edm,errdef;
    Int_t nvpar,nparx,icstat;
    minuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
    fEnergyChi2  = amin/nvpar;
    
    TF1 *f1 = new TF1("f1",PhotoElectrons2Total,fEvent->GetStartGtu(),fEvent->GetEndGtu(),2);
    TF1 *f2 = new TF1("f2",PhotoElectrons2Recovery,fEvent->GetStartGtu(),fEvent->GetEndGtu(),2);

    f1->SetParameters(outpar);
    f2->SetParameters(outpar);
    f1->Draw("same");
    f2->Draw("same");
    f1->SetLineStyle(1);   f1->SetLineColor(1);
    f2->SetLineStyle(2);   f2->SetLineColor(2);
    SaveFitInfo();
    SafeDelete(f1);
    SafeDelete(f2);
    SafeDelete(minuit);
}


//_____________________________________________________________________________
 void EnergyModule::FitMode2fumili() {
    TFumili *fumili = new TFumili(2);
    fumili->SetObjectFit(&fData);
    fumili->SetFCN(PhotoElectronsTimeFit2);
    fumili->SetErrorDef(1.);
    Double_t par[2] = {1.e-1,0.0};
    Double_t step[3] = {1.e-1,1.e0};
    Double_t min[3] = {0.00, 0.00};
    Double_t max[3] = {0.00, 0.00};
    string cpar[3] = {"Energy","tMax"};
    
    for( Int_t i=0; i<2; i++ ) {
        fumili->SetParameter(i,cpar[i].c_str(),par[i],step[i],min[i],max[i]);
    }
//     fumili->FixParameter(1);
//     fumili->Migrad();
    fumili->Minimize();
    // retrieve results
    Double_t outpar[2],err[2];
    char cname[100];
    Double_t vlow[2], vhigh[2];
    for( Int_t i=0; i<2; i++ ) {
	fumili->GetParameter(i,cname,outpar[i],err[i],vlow[i],vhigh[i]); 
    }
    fEnergy = outpar[0]*1e15; 
    fEnergyError = err[0]*1e15;
    Double_t amin,edm,errdef;
    Int_t nvpar,nparx;
    fumili->GetStats(amin,edm,errdef,nvpar,nparx);
    fEnergyChi2  = amin/nvpar;
    
    TF1 *f1 = new TF1("f1",PhotoElectrons2Total,fEvent->GetStartGtu(),fEvent->GetEndGtu(),2);
    TF1 *f2 = new TF1("f2",PhotoElectrons2Recovery,fEvent->GetStartGtu(),fEvent->GetEndGtu(),2);

    f1->SetParameters(outpar);
    f2->SetParameters(outpar);
    f1->Draw("same");
    f2->Draw("same");
    f1->SetLineStyle(1);   f1->SetLineColor(1);
    f2->SetLineStyle(2);   f2->SetLineColor(2);
    SaveFitInfo();
    SafeDelete(f1);
    SafeDelete(f2);
    SafeDelete(fumili);
}

//_____________________________________________________________________________
 void EnergyModule::SaveFitInfo() {

    TPaveText *cl = new TPaveText(0.78,0.6,1,1,"brNDC");
    cl->SetTextFont(72); cl->SetTextSize(0.02); cl->SetFillColor(18); cl->SetTextAlign(12);
    // Event Information
    cl->AddText("Truth Information");
    cl->AddText(Form("Energy = %.2E MeV",fEvent->GetHeader().GetTrueEnergy()));
    cl->AddText(Form("#Theta = %.2f deg",fEvent->GetHeader().GetTrueTheta()*RadToDeg()));
    cl->AddText(Form("#phi   = %.2f deg",fEvent->GetHeader().GetTruePhi()*RadToDeg()));
    cl->AddText("Reco Information");
    cl->AddText(Form("Energy = %.2E #pm %2.E MeV",fEnergy,fEnergyError));
    cl->AddText(Form("#Theta = %.2f deg",fTheta*RadToDeg()));
    cl->AddText(Form("#phi   = %.2f deg",fPhi*RadToDeg()));
    cl->Draw();
    if (fSaveEps == "yes")  
        if(fCanvas) fCanvas->Print(Form("output/check_timedistribution.run%d.event%d.eps",
	                 fEvent->GetHeader().GetRun(),fEvent->GetHeader().GetNum()));
    SafeDelete(cl);
}
//_____________________________________________________________________________
 void EnergyModule::EstimateBackground() {
    map <Int_t,Int_t> mapGtuCounts;
    map <Int_t,Int_t> mapGtuPixels;
    vector<RecoPixelData*> Bg;
    
    for (UInt_t i(0); i<fAllActivePixels.size(); i++) {
         // look if this pixel at the given time was taken as signal
         RecoPixelData* pix = fAllActivePixels[i];
	 Bool_t IsPixelFound = kFALSE;
	 for (UInt_t j(0); j<fPixelsCluster.size(); j++) {
 	      RecoPixelData* pix2 = fEvent->GetRecoPixelData(fPixelsCluster[j]);
	      if (pix->GetPixelId() == pix2->GetPixelId() && pix->GetGtu() == pix2->GetGtu()) {
	          if (pix->GetCounts() != pix2->GetCounts()) 
		  Msg(EsafMsg::Warning) << "EstimateBackground() two pixels with same Gtu and same Uid have different number of Counts. Is it possible?!" << MsgDispatch;
		  IsPixelFound = kTRUE;
                  break;  
	      }
	 }
	 if (!IsPixelFound) Bg.push_back(pix);
    }
    for (UInt_t i(0); i<Bg.size(); i++) {
         RecoPixelData* pix = Bg[i];
	 mapGtuCounts[pix->GetGtu()] += pix->GetCounts();
	 mapGtuPixels[pix->GetGtu()] += 1;
    } 
    map <Int_t,Int_t>::iterator it; 
    Int_t point(0);
    for (it = mapGtuCounts.begin(); it != mapGtuCounts.end(); it++) {
        Float_t noise = 1.e0*it->second/mapGtuPixels[it->first]; 
        fBgGraph->SetPoint(point++,it->first,noise);
    }
    if(fCanvas) fCanvas->Clear();
    fCanvas->Divide(1,2);
    fCanvas->cd(1);
    fBgGraph->Fit("pol0","Q");
    fBackground = fBgGraph->GetFunction("pol0")->GetParameter(0);
    fBgGraph->Draw("A*");
    fBgGraph->GetXaxis()->SetTitle("Gtu Hit");
    fBgGraph->GetYaxis()->SetTitle("Noise/pixel/GTU");
    fCanvas->cd(2);
    TH1F* hNoise = new TH1F("hNoise","Noise/pixel/GTU",100,1,2);
    Double_t *noise = fBgGraph->GetY();
    for (Int_t i(0); i < fBgGraph->GetN(); i++) {
         hNoise->Fill(noise[i]); 
    }
    hNoise->Draw();
    hNoise->GetXaxis()->SetTitle("Noise/pixel/GTU");
    fCanvas->cd(0);
    if (fSaveEps == "yes")  
        if(fCanvas) fCanvas->Print(Form("output/check_noisedistribution.run%d.event%d.eps",
	                 fEvent->GetHeader().GetRun(),fEvent->GetHeader().GetNum()));
    SafeDelete(hNoise);
}
//_____________________________________________________________________________
 Bool_t EnergyModule::Process(RecoEvent *ev) {
    // manager of the energy reconstruction chain
    fEvent      = ev;
    fRunpars    = fEvent->GetHeader().GetRunPars();
    fQuantumEff = fRunpars->GetPmtData().GetQuantum();
    CreateAuxilaryObjects();
    // Load information from previous modules
    LoadModules();
    // Make map <GtuHit_t,Pixels_t> 
    SortGtu();
    // Make maps <GtuHit_t,MeanTheta>,  <GtuHit_t,MeanPhi>
    BuildFoVAngleFunctions();
    // Make map of pmts and pixels to display 
    DisplayCluster();
    // Find Theta, Phi at Shower Maximum
    FindThetaPhiAtMaximum();
    // Estimate Fluorescence yield for the Shower Maximum
    MakeFluorYield();
    // Estimate Background
    EstimateBackground();
    // Build the recovery function along the track
    BuildRecoveryFunction();
    // Generate Shower Track, fill expected histo for flux of photons on EP
    ShowerTrack3D();
    // Fit the time distribution
    FitTimeDistribution();
    // Save data                                    
    MyData()->Add("Energy",fEnergy);
    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
    // 
    Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
    return kTRUE;
}

//_____________________________________________________________________________
 void EnergyModule::UserMemoryClean()  {
    // 
    // Here we clear any auxilary objects created during one event energy reconstruction
    //
    fTheta                     = -1; 
    fPhi                       = -1;
    fNumPoints                 = -1;
    fAltitude                  = -1; 
    fGtuMin                    =  1.e20;
    fGtuMax                    = -1.e20;
    fQuantumEff                = -1;
    fTotalYield                = -1; 
    fEnergy                    = -1;
    fEnergyError               = -1; 
    fEnergyChi2                = -1; 
    fFluorMaxPe                = -1; 
    fDeltaL                    = -1;
    fGtuAtMax                  = -1; 
    fGtuTimeAtMax              = -1;
    fThetaFoVAtMax             = -1; 
    fPhiFoVAtMax               = -1; 
    fMeanAttenuation           = -1;
    fMeanLambdaOnEntrancePupil = -1;
    fRmax                      = -1;
    fBackground                = -1;
    fXmax                      = -1.e20;
    fXmin                      =  1.e20;
    fYmax                      = -1.e20;
    fYmin                      =  1.e20;
    fData.Clear();
    fMaxPos.Clear(); 
    SafeDelete(fEPspectrum);
    SafeDelete(fEPphotons);
    fPixelsCluster.clear();
    fAllActivePixels.clear();
    fShowerDir.Clear();
    fPmtLines.clear();
    fPixelLines.clear();
    fGtuPixels.clear();
    fGtuMyPixels.clear();
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update 14-May-2005 21:31