Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

HmaxByShapeMethodModule - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: HmaxByShapeMethodModule.cc,v 1.15 2005/10/02 14:18:08 thea Exp $
// Dmitry V. Naumov created Jun,  6 2004

#include "HmaxByShapeMethodModule.hh"
#include "RecoEvent.hh"
#include "RecoPixelData.hh"
#include "RecoRootEvent.hh"
#include "Atmosphere.hh"
#include "Config.hh"
#include "utils.hh"

#include <TProfile.h>
#include <TCanvas.h>
#include <TH1F.h>
#include <TF1.h>

using namespace sou;

ClassImp(HmaxByShapeMethodModule)

HmaxByShapeMethodModule* tmpHMaxpointer = NULL; // Points to the HmaxByShapeMethodModule class

// fFluorescence function
Double_t HmaxByShapeMethod_fFluorescence(Double_t *x, Double_t *par) {
    return par[0]*exp(-0.5*pow(x[0]-par[1],2)/pow(par[2],2));
}

// Sum of background and peak functions
Double_t HmaxByShapeMethod_fitFunction(Double_t *x, Double_t *par) {
    return HmaxByShapeMethod_fFluorescence(x,&par[0]) + tmpHMaxpointer->GetBackground();
}
//_____________________________________________________________________________
float HmaxByShapeMethod_ShowerMaximumDerivative(float t){
    float meanlambda = 360; // nm
    float xr = 2974/37.15;        // g/cm2
    float tmax = 1.7 + 0.76*log(tmpHMaxpointer->GetEnergy()/8.1e7);  
    float s = 2./(1.+tmax/t);
    
    float CosTheta = cos(tmpHMaxpointer->GetTheta());
    RecoPixelData* pix =             
        tmpHMaxpointer->GetClusters()[tmpHMaxpointer->PixelId(tmpHMaxpointer->GetFluorMaxTime())];
    float ct = cos(pix->GetTheta());
    return 1. - 2*log(s) - tmax/t*s - pow(400/meanlambda,4)/xr*CosTheta/ct;
}
//_____________________________________________________________________________
float HmaxByShapeMethod_DensityToAltitude(float h) {
    return Atmosphere::Get()->Air_Density(h) - tmpHMaxpointer->GetDensity();
}
//_____________________________________________________________________________
 HmaxByShapeMethodModule::HmaxByShapeMethodModule() : RecoModule("HmaxByShapeMethod") {
    // ctor
    tmpHMaxpointer = this;
    fEnergy        = 1.e20;
    fTemp           = NULL;
    fFitFcn         = NULL;
    fFluor          = NULL;

}

//_____________________________________________________________________________
 HmaxByShapeMethodModule::~HmaxByShapeMethodModule() {
    // dtor
    if (!fc)
        delete fc;
}

//_____________________________________________________________________________
 Bool_t HmaxByShapeMethodModule::Init() {
    Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
    Atmosphere::Get();
    Msg(EsafMsg::Info) << "Enabled Atmosphere: " << Atmosphere::Get()->GetType() << MsgDispatch;
    
    // Initialize the mean background
    // FIXME:need to be read from the Event and not from config file
    ConfigFileParser *pConfig = Config::Get()->GetCF("Electronics","EusoElectronics");
    fMeanBackground = pConfig->GetNum("EusoElectronics.fNightGlowRateOnAxis");
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t HmaxByShapeMethodModule::PreProcess() {
    
    fc = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("hmaxbyshape");
    if (!fc) fc = new TCanvas("hmaxbyshape","hmaxbyshape");

    fc->Clear();
    fc->cd();
    
    fc->SetFillColor(33);
    fc->SetFrameFillColor(41);
    fc->SetGrid();
    fTemp   = NULL;
    fFitFcn = NULL;
    fFluor  = NULL;
    fEv     = NULL;
    fAltitude = -1; 
    fFluorMaxAt = -1; 
    fFluorMaxPe = -1;
    fFluorMaxPeObs = -1;
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t HmaxByShapeMethodModule::Process(RecoEvent *ev) {
    fEv = ev;
    const RecoModuleData *td = ev->GetModuleData("TrackDirection");
    if ( td==NULL ) {
        Msg(EsafMsg::Panic) << "TrackDirectionModule is not found." << MsgDispatch;	    
        throw runtime_error( "HmaxByShapeMethodModule failed.");
    }
    fTheta    = td->GetDouble("Theta");
    fPhi      = td->GetDouble("Phi");
    fOmega.SetXYZ(TMath::Sin(fTheta)*TMath::Cos(fPhi), TMath::Sin(fTheta)*TMath::Sin(fPhi), 
		  TMath::Cos(fTheta));
    const RecoModuleData *gcm = ev->GetModuleData("GTUClustering");
    if ( gcm==NULL ) {
        Msg(EsafMsg::Panic) << "GTUClusteringModule is not found." << MsgDispatch;
        throw runtime_error("HmaxByShapeMethodModule failed.");
    }
    else {
        if ( gcm->GetObj("CluPixels") == NULL ) {
            Msg(EsafMsg::Panic) << "No data in event" << MsgDispatch;	    
            return kTRUE;
        }
        fPointsId = *(vector<Int_t>*)gcm->GetObj("CluPixels");

        if (fPointsId.size() == 0) {
            fNumPoints = 0;
            Msg(EsafMsg::Panic) << "No pixels selected: terminated." << MsgDispatch;	                
            return kTRUE;
        } else {
            fNumPoints = fPointsId.size();
            Msg(EsafMsg::Debug) << "Processing " << fNumPoints << " points" << MsgDispatch;
        }
    }
    fGtuLength = ev->GetHeader().GetGtuLength()/1000;
    fFirstHit = ev->GetStartGtu();
    fLastHit  = ev->GetEndGtu();
    Float_t Nbins = (Int_t) (fLastHit - fFirstHit)/fGtuLength;
    fTemp = (TH1F*)gDirectory->FindObject("fTemp_HMaxShapeModule");
    if ( fTemp != NULL) delete fTemp;
    
    fTemp = new TH1F("fTemp_HMaxShapeModule","",(Int_t)Nbins,fFirstHit,fLastHit);
    fTemp->GetXaxis()->SetTitle("time, #mu sec");
    fTemp->GetYaxis()->SetTitle("PhotoElectrons");

    // Save event unique tag
    fEventTag = Form("hmaxshape.run.%d.event%d",ev->GetHeader().GetRun(),ev->GetHeader().GetNum());
    for(Int_t i=0; i<fNumPoints; i++) {
        RecoPixelData *pix = ev->GetRecoPixelData(fPointsId[i]);
        fTemp->Fill(pix->GetGtu(),pix->GetCounts());
	cluster.push_back(pix);
    }
    Double_t par[3]={1,1,1};
    fFitFcn = new TF1("fFitFcn",HmaxByShapeMethod_fitFunction,ev->GetStartGtu(),ev->GetEndGtu(),3);
    fFluor  = new TF1("fFluor",HmaxByShapeMethod_fFluorescence,ev->GetStartGtu(),ev->GetEndGtu(),3);
    fFitFcn->SetNpx(500);
    fFitFcn->SetLineWidth(2);
    fFitFcn->SetLineColor(kMagenta);
    // find center of mass of fTemp
    // Try to find the fFluorescent peak
    float maxintfl(0);
    int FlWidth = 10;
    for(int i=0;i<(Int_t)fTemp->GetNbinsX();i++){
        if(maxintfl < fTemp->Integral(i,i+FlWidth) ){
        maxintfl = fTemp->Integral(i,i+FlWidth);
        double ftime_fTemp(0), fstat(0);
        for (int j=i;j<i+FlWidth;j++) {
            ftime_fTemp += fTemp->GetBinContent(j)*fTemp->GetBinCenter(j);
            fstat += fTemp->GetBinContent(j);
        }
        fFluorMaxTime  = ftime_fTemp/fstat;
        fFluorMaxPe    = fTemp->GetBinContent(fTemp->FindBin(fFluorMaxTime));
        }
    }
    fFitFcn->SetParLimits(0,0,fTemp->GetMaximum());
    fFitFcn->SetParLimits(1,0,ev->GetEndGtu());
    fFitFcn->SetParLimits(2,0,0.5*ev->GetEndGtu()); 
    fFitFcn->SetParameters(fFluorMaxPe,fFluorMaxTime,fTemp->GetRMS());
   
    fTemp->Fit("fFitFcn","RQ");
    fChi2 = fTemp->GetFunction("fFitFcn")->GetChisquare()/fTemp->GetFunction("fFitFcn")->GetNDF();
    fTemp->GetFunction("fFitFcn")->GetParameters(&par[0]);
    fFitFcn->SetParameters(par);
    fTemp->Draw("E");
    fFitFcn->Draw("same");
    fFluor->SetParameters(par);
    fFluorMaxPe   = par[0];
    fFluorMaxTime = par[1];
    if ( PixelId(fFluorMaxTime) == -1) return kFALSE;
    GetAltitude();
    MyData()->Add("Altitude",fAltitude*km);
    MyData()->Add("FluorMaxAt",fFluorMaxAt);
    MyData()->Add("FluorMaxPe",fFluorMaxPe);
    MyData()->Add("FluorMaxPeObs",fTemp->GetMaximum());
    MyData()->Add("DeltaL",DeltaL(fFluorMaxTime)*km);

    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t HmaxByShapeMethodModule::PostProcess() {
    if (fTemp!=(TH1F*)NULL) {         
        fTemp->SetMarkerStyle(21);
        fTemp->SetMarkerSize(0.8);
        fTemp->SetStats(0);
    }
    ConfigFileParser *pConfig = Config::Get()->GetCF("Reco","HmaxByShapeMethodModule");
    string saveEps = pConfig->GetStr("HmaxByShapeMethodModule.SaveEps");
 
    if (saveEps == "yes")
        fc->Print("output/"+fEventTag+".eps");
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t HmaxByShapeMethodModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
    fRecoRootEvent->GetRecoHmaxByShapeMethod().SetQuality(fChi2);
    fRecoRootEvent->GetRecoHmaxByShapeMethod().SetHmax(fAltitude);
    fRecoRootEvent->GetRecoHmaxByShapeMethod().SetTempHisto(fTemp);
    fRecoRootEvent->GetRecoHmaxByShapeMethod().SetErrorHmax(fAltitude -
                                                            fEv->GetHeader().GetTrueShowerMaxPos().Z()/km);
    return kTRUE;
}

//_____________________________________________________________________________
 Bool_t HmaxByShapeMethodModule::Done() {
     Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
     return kTRUE;
}

//_____________________________________________________________________________
 void HmaxByShapeMethodModule::UserMemoryClean()  {
    //
    //
    //
    SafeDelete(fFitFcn);
    SafeDelete(fFluor);
    SafeDelete(fTemp);
    cluster.clear();
}

//_____________________________________________________________________________
 Int_t HmaxByShapeMethodModule::PixelId(Double_t time) {
    fFluorMaxAt = -1;
    // This methods returns the pixel id number corresponding to the given time moment
    // returns -1 if no pixel found (WARNING)
    for(UInt_t i = 0; i < cluster.size() - 1; i++) 
	if (cluster[i]->GetGtu() < time && time <= cluster[i+1]->GetGtu()) {
	    fFluorMaxAt = i;
	    return i;
	}
    Msg(EsafMsg::Warning) << "PixelId can not determine pixel id corresponding to time " << time  << MsgDispatch;
    
    return fFluorMaxAt;
}

//_____________________________________________________________________________
 Double_t HmaxByShapeMethodModule::nOmega(Double_t time) {
    // This method computes the scalar product of the shower unit vector
    // Omega (sin(Theta)*cos(Phi), sin(Theta)*sin(Phi), cos(Theta))
    // and the given pixel unit vector aling its field of view direction
    Int_t GtuHit = PixelId(time);
    if (GtuHit == -1) 
	return -100;
    TVector3 pixelUV;
    RecoPixelData *pix = cluster[GtuHit];
    Double_t st = TMath::Sin(pix->GetTheta());
    Double_t ct = TMath::Cos(pix->GetTheta());
    Double_t sp = TMath::Sin(pix->GetPhi());
    Double_t cp = TMath::Cos(pix->GetPhi());
    pixelUV.SetXYZ(st*cp,st*sp,ct);
    return pixelUV*fOmega;
}

//_____________________________________________________________________________
 Double_t HmaxByShapeMethodModule::DeltaL(Double_t time)
{
    // This method returns the track length inclined with Theta and Phi zenith and azimuth angles
    // respectively seen by the given pixel in one GTU time interval

    return 0.3*fGtuLength/(1.+nOmega(time));
}

//_____________________________________________________________________________
 void HmaxByShapeMethodModule::GetAltitude()
{            
    Float_t x0 = 37.15*g/cm2;
    Float_t t1 = 10, t2 = 100;
    Float_t (*f)(float)= &HmaxByShapeMethod_ShowerMaximumDerivative;
    zbrac(f,&t1,&t2);
    Float_t t0 =  rtbis(f,t1,t2,0.01);
    Double_t FSecondDerivative = -1/(2*t0);
    Double_t FluorEntries = fFluor->Integral(fFirstHit,fLastHit)/fGtuLength;
    Double_t erfLeft  = TMath::Erf(sqrt(log(fFluorMaxPe/fFluor->Eval(fFirstHit))));
    Double_t erfRight = TMath::Erf(sqrt(log(fFluorMaxPe/fFluor->Eval(fLastHit))));
    fDensity = x0*fFluorMaxPe/FluorEntries*(erfLeft+erfRight)*
        TMath::Sqrt(-TMath::Pi()/FSecondDerivative/2)/(DeltaL(fFluorMaxTime)*km);
    DensityToAltitude();
}

//_____________________________________________________________________________

 void HmaxByShapeMethodModule::DensityToAltitude() {
    Float_t h1 = 0*km, h2 = 100*km; 
    Float_t (*f)(float)= &HmaxByShapeMethod_DensityToAltitude;
    zbrac(f,&h1,&h2);
    fAltitude = rtbis(f,h1,h2,0.01*km)/km;
}
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