CaloElectron/src/CaloElectron.cpp

#include <CaloElectron.h>

////////////////////////////////////////////////////////////
/// INITIALIZATIONS & GLOBAL PARAMETERS
////////////////////////////////////////////////////////////

CaloElectron_parameters * CaloElectron_parameters::_parameters = 0; 

/**
 * Set external parameters to default values
 */
void CaloElectron_parameters::SetDefault(){
    cout << "CaloElectron --> SET DEFAULT PARAMETERS" << endl;

    isimu = 0;
    dolongfit = 1;
    dolatcorr = 1;
    calibmode = 0;

    debug = false;

//     dodrawlong = false;
//     dodrawq0q1 = false;
//     dodrawlat = false;

    TAUMAX = 5.;
    RMAX   = 40.;//4.; //temporaneo???
    RMAXsp = 40.;
    
//     cout << "Shower cuts:"<<endl;
//     cout << "TAUMAX = "<<TAUMAX<<endl;
//     cout << "RMAX   = "<<RMAX<<" cm"<<endl;

    for(int ip=0; ip<22; ip++)
        for(int iv=0; iv<2; iv++)
            maskpl[ip][iv] = false; 
    //
    h_qtot[0] = NULL;
    h_qtot[1] = NULL;
    //
    h_qfit[0] = NULL;
    h_qfit[1] = NULL;
    //
    cl1 = NULL;

    //
    for(int i=0; i<NTAUBINMAX; i++){
        for(int j=0; j<3; j++){
            h_qtot_tau[i][j]     = NULL;
            h_qtot_tau_wav[i][j] = NULL;
            for(int ii=0; ii<CALIBNBIN; ii++){
                summW[i][j][ii] = 0.;
                summWF[i][j][ii] = 0.;
            }
        }
    }

    SetCalo(0,0);
    SetPar();
    
    Dump();

};
/**
 * Set external parameters to default values
 */
void CaloElectron_parameters::Dump(){
    cout << "CaloElectron --> DUMP PARAMETERS" << endl;
    cout << endl;
    cout << "Shower cuts (for energy determination):"<<endl;
    cout << " TAUMAX = "<<TAUMAX<<" TAU=t/T"<<endl;
    cout << " RMAX   = "<<RMAX<<" cm"<<endl;
//    cout << endl;
    cout << "MASKED PLANES:"<<endl;
    cout << " X -> ";
    for(int ip=0; ip<22; ip++)cout << maskpl[ip][0];  
    cout << endl;
    cout << " Y -> ";
    for(int ip=0; ip<22; ip++)cout << maskpl[ip][1];  
    cout << endl;
    cout << " ipmin = "<<ipmin<<endl;
    cout << " ipmax = "<<ipmax<<endl;
//    cout << endl;
    cout << "Flags:"<<endl;
    cout << " lateral correction ?      "<<dolatcorr<<endl;
    cout << " longitudinal correction ? "<<dolongfit<<endl;
    cout << " simulation ?              "<<isimu<<endl;
    cout << " calibration mode ?        "<<calibmode<<endl;
    cout << endl;
    
};
/**
 * Set calorimeter masks and book histograms.
 * @icalo Calorimeter configuration (icalo=0 --> all planes included)
 * @notused Num.of planes exluded (icalo=1 excludes bottom planes, icalo=2 exclude top planes)
 * NB: plane 19x is always excluded   
 * 
 */
void CaloElectron_parameters::SetCalo(int icalo,int notused){
    
    
    if(icalo==0){
        maskpl[18][0]=true;
        ipmin=0;
        ipmax=22;
    }
    if(icalo==1){
        maskpl[18][0]=true;
        ipmin=0;
        ipmax=22-notused;
        if(18-notused>=0)maskpl[18-notused][0]=true;
    }
    if(icalo==2){
        maskpl[18][0]=true;
        ipmin=notused;
        ipmax=22;
        if(18+notused<22)maskpl[18+notused][0]=true; 
    }

//     cout << "MASKED PLANES:"<<endl;
//     cout << "X -> ";
//     for(int ip=0; ip<22; ip++)cout << maskpl[ip][0];  
//     cout << endl;
//     cout << "Y -> ";
//     for(int ip=0; ip<22; ip++)cout << maskpl[ip][1];  
//     cout << endl;
//     cout << "ipmin = "<<ipmin<<endl;
//     cout << "ipmax = "<<ipmax<<endl;
    
    // histos for kolmogorov test 
//     if(h_qtot)h_qtot->Delete();
//     if(h_qfit)h_qfit->Delete();
//     h_qtot[0] = new TH2F("htot0","Matrix of energy releases (measured)",192,-0.5,191.5,22,-0.5,21.5);
//     h_qfit[0] = new TH2F("hfit0","Matrix of energy releases (fit)",     192,-0.5,191.5,22,-0.5,21.5);
//     h_qtot[1] = new TH2F("htot1","Matrix of energy releases (measured)",192,-0.5,191.5,22,-0.5,21.5);
//     h_qfit[1] = new TH2F("hfit1","Matrix of energy releases (fit)",     192,-0.5,191.5,22,-0.5,21.5);

    // binnaggio in funzione della coordinata, relativo al margine del piano
    int nbin = 100;
    double xbin[101];
    // retrieve the coordinate of the bottom-left corner of the plane
    xbin[0]=0.;
    int ibin = 0;
    for(int ise=0; ise<3; ise++){//sensors
        ibin++;
        if(ise==0)         xbin[ibin] = xbin[ibin-1] + DEAD;
        if(ise==1||ise==2) xbin[ibin] = xbin[ibin-1] + DEAD + GLUE + DEAD;
        for(int is=0; is<32; is++){//strip
            ibin++;
            xbin[ibin] = xbin[ibin-1] + PITCH;
        }       
    }
    ibin++;
    xbin[ibin] = xbin[ibin-1] + DEAD;

    if(h_qtot[0])h_qtot[0]->Delete();
    if(h_qtot[1])h_qtot[1]->Delete();
    if(h_qfit[0])h_qfit[0]->Delete();
    if(h_qfit[1])h_qfit[1]->Delete();

    h_qtot[0] = new TH2F("htot0","Matrix of energy releases (measured)",nbin,xbin,22,-0.5,21.5);
    h_qfit[0] = new TH2F("hfit0","Matrix of energy releases (fit)",     nbin,xbin,22,-0.5,21.5);
    h_qtot[1] = new TH2F("htot1","Matrix of energy releases (measured)",nbin,xbin,22,-0.5,21.5);
    h_qfit[1] = new TH2F("hfit1","Matrix of energy releases (fit)",     nbin,xbin,22,-0.5,21.5);

}
/**
 * Set external parameters to default values
 */
void CaloElectron_parameters::SetPar(){

    // ---------------
    // lateral profile (hardcoded parameters )
    // ---------------
    TString pathx,pathy,pathyd,pathyp;
    TString dir = gSystem->Getenv("PAM_CALIB");
    if(isimu==0){
        pathx=dir+"/cal-param/parxdati.txt";
        pathyd=dir+"/cal-param/paryddati.txt";
        pathyp=dir+"/cal-param/parypdati.txt";
    }else{
        pathx=dir+"/cal-param/parxsimu.txt";
        pathyd=dir+"/cal-param/parydsimu.txt";    
        pathyp=dir+"/cal-param/parypsimu.txt"; 
    }
    
    cout << "Load LATERAL PARAMETERS "<<endl;
    cout << " X        from "<<pathx << endl;
    cout << " Y-ODD    from "<<pathyd << endl;
    cout << " Y-EVEN   from "<<pathyp << endl;

    SetPar(pathx,pathyd,pathyp);
   
}
/**
 * Set external parameters from input files
 */
void CaloElectron_parameters::SetPar(TString pathx,TString pathyd,TString pathyp){

//these parameter are directly the values of rp,rt,p in each taubin
//they are obtained from the simulation in different taubin up to tau=3
    ifstream file;
    for(int iv=0;iv<2;iv++){
        if(iv==0){
            ifstream file(pathx);
            for(int i=0;i<12;i++)file>>file_tau[i]>>file_p[i][iv]>>file_rt[i][iv]>>file_rc[i][iv];
        } 
        if(iv==1){
            ifstream filep(pathyp);
            for(int i=0;i<12;i++)filep>>file_tau[i]>>file_p[i][iv]>>file_rt[i][iv]>>file_rc[i][iv];
            ifstream filed(pathyd);
            for(int i=0;i<12;i++)filed>>file_tau[i]>>file_p[i][iv+1]>>file_rt[i][iv+1]>>file_rc[i][iv+1];
        }
    }   
    
};

/**
 * Set calibration tools.
 * Create histos for lateral profile parameterization
 */
void CaloElectron_parameters::Calibration_SetTools(int ntau, float* taubin){

    if(!calibmode)return;

    cout << "CaloElectron --> Set calibration tools "<<endl;
    if( ntau > NTAUBINMAX ){
        ntau = NTAUBINMAX;
        cout << "(NB: max allowed n.tau-bins = "<<NTAUBINMAX<<")"<<endl;
    };
    par_ntau = ntau;    
    for(int i=0; i<=par_ntau; i++)par_taubin[i]=taubin[i];
    //
    cout << "Set TAU binning:"<<endl;
    cout << par_ntau<<" bins, from "<<par_taubin[0]<<" to "<<par_taubin[par_ntau]<<endl;
    cout << endl;
    //
    // create histos
    //
    TString htitle,hid;
    for(int i=0; i<par_ntau; i++){
        for(int j=0; j<3; j++){
            //
            hid.Form("hqtottau_%i_%i",i,j);
            htitle.Form("Lateral profile - view %i - tau-bin %i",j,i);
            cout <<hid<<" --> "<< htitle<<" "<<h_qtot_tau[i][j]<<endl;      
            if( h_qtot_tau[i][j] )h_qtot_tau[i][j]->Delete();
            h_qtot_tau[i][j] = new TH2F(hid.Data(),htitle.Data(),CALIBNBIN,-1*CALIBRANGE,1*CALIBRANGE,5000,0.,10.);
            //
            hid.Form("hqtottauwav_%i_%i",i,j);
            htitle.Form("Lateral profile - view %i - tau-bin %i (weighted average)",j,i);
            cout <<hid<<" --> "<< htitle<<" "<<h_qtot_tau_wav[i][j]<<endl;          
            if( h_qtot_tau_wav[i][j] )h_qtot_tau_wav[i][j]->Delete();
            h_qtot_tau_wav[i][j] = new TH1F(hid.Data(),htitle.Data(),CALIBNBIN,-1*CALIBRANGE,1*CALIBRANGE);
        }
    }
    
    return;
};
/**
 * Set calibration tools
 * Create histos for lateral profile parameterization
 */
void CaloElectron_parameters::Calibration_SetTools(){

    if(!calibmode)return;

    int ntau = 12;
    float fromtau = 0.;
    float totau   = 3.;
    float dtau = (totau-fromtau)/(float)ntau;
    vector<float> taubin;
    float *pt_taubin = taubin.get_allocator().allocate(ntau+1);
    pt_taubin[0]=fromtau;
    for(int i=1; i<=ntau; i++)pt_taubin[i]=(pt_taubin[0]+i*dtau);
    Calibration_SetTools(ntau,pt_taubin);
    return;
};
/**
 * Fit calibration histos
 */
void CaloElectron_parameters::Calibration_Fit(){

    if(!calibmode)return;

    TF1* frad = new TF1("frad",fradpro,-5.,5.,4);
    //
//    int nt_tau  =  CaloElectron_parameters::Get()->par_ntau;
    int ntau  =  par_ntau;
    //
    for(int itau = 0; itau<ntau; itau++){ //loop over tau-bins
        for(int j=0; j<3; j++){
            //
//          TProfile* hh = 0;
//          if( !h_qtot_tau[itau][j] )continue;
//          cout << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ "<<endl;
//          cout << " fitting profile histo "<<h_qtot_tau[itau][j]->GetName()<<endl;
//          cout << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ "<<endl;

//          hh = h_qtot_tau[itau][j]->ProfileX("hh");
//          frad->SetParameter(0,1.5);//rt
//          frad->SetParameter(1,0.5);//p
//          frad->SetParameter(2,0.4);//rc
//          frad->SetParameter(3,.01);//nrom
//          hh->Fit(frad,"R","",-2,2);
//          hh->Fit(frad,"R","",-5,5);
//          hh->Delete();
            //
            if( !h_qtot_tau_wav[itau][j] )continue;
            cout << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ "<<endl;
            cout << " fitting profile histo "<<h_qtot_tau_wav[itau][j]->GetName()<<endl;
            cout << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ "<<endl;

            TH1F *hh = h_qtot_tau_wav[itau][j];
            frad->SetParameter(0,1.5);//rt
            frad->SetParameter(1,0.5);//p
            frad->SetParameter(2,0.4);//rc
            frad->SetParameter(3,.01);//nrom
            hh->Fit(frad,"R","",-2,2);
            hh->Fit(frad,"R","",-5,5);
            hh->Delete();

        };
    };
    frad->Delete();

    return;
};
/**
 * Write calibration output to file
 */
void CaloElectron_parameters::Calibration_Save(TFile* file){

    if(!calibmode)return;

    if( !file || (file&&file->IsZombie()) )return;
    file->cd();
    //
//    int ntau        = CaloElectron_parameters::Get()->par_ntau;
    int ntau        = par_ntau;
    for(int itau=0; itau<ntau; itau++){
        for(int j=0; j<3; j++){

            //
            // fill weighted-average histogram
            //
            for(int ibin=0; ibin<CALIBNBIN; ibin++){
                if(summW[itau][j][ibin]==0.)continue;
                h_qtot_tau_wav[itau][j]->SetBinContent(ibin+1,summWF[itau][j][ibin]/summW[itau][j][ibin]);
                h_qtot_tau_wav[itau][j]->SetBinError(ibin+1,1./sqrt(summW[itau][j][ibin]));
            }
            
            if(h_qtot_tau[itau][j])     h_qtot_tau[itau][j]->Write();
            if(h_qtot_tau_wav[itau][j]) h_qtot_tau_wav[itau][j]->Write();
        }
    }    

    return;
};
/**
 * Load calibration output histos from existing file
 */
void CaloElectron_parameters::Calibration_Load(TFile* file){


    if( !file || (file&&file->IsZombie()) )return;
    file->cd();
    cout <<" File: "<<file->GetName();

    //
    TString hid;
    par_ntau=0;
    for(int i=0; i<100; i++){
        for(int j=0; j<3; j++){
            //
            hid.Form("hqtottau_%i_%i",i,j);
            TH2F* hh = (TH2F*)file->Get(hid.Data());
            if(!hh)continue;
            cout<<"histo --> "<<hid<<endl;
            h_qtot_tau[i][j]=hh;
            //
            hid.Form("hqtottauwav_%i_%i",i,j);
            TH1F* hhh = (TH1F*)file->Get(hid.Data());
            if(!hhh)continue;
            cout<<"histo --> "<<hid<<endl;
            h_qtot_tau_wav[i][j]=hhh;
            //
            par_ntau=i+1;
        }
    }

    calibmode = 1;

    return;
};


////////////////////////////////////////////////////////////
/// CLASS IMPLEMENTATION
////////////////////////////////////////////////////////////
/**
 * Reset event
 * it has to be called for each event
 */
void CaloElectron::Reset(){
        

    for(int ip=0; ip<22; ip++){
        for(int iv=0; iv<2; iv++){
            trkcoordx[ip][iv] = 0.;
            trkcoordy[ip][iv] = 0.;
//          trkcoordz[ip][iv] = 0.;
            trkstrip[ip][iv] = 0.;
            qplane0[ip][iv]  = 0.;
            qplane1[ip][iv]  = 0.;
            qplanefit[ip][iv]  = 0.;
            tplane[ip][iv]   = -100.;

            RMS[ip][iv]  = 0.;
            SKEW[ip][iv] = 0.;
            KURT[ip][iv] = 0.;
            TAU[ip][iv]  = 0.;
        }
    }    

    qtot = 0.;

    tg = 0.;
    tgx = 0.;
    tgy = 0.;

    corr__l = 1.;
    err__l = 0;
    chi2__l = 0.;
    chi2    = 0.;
    for(int i=0; i<4; i++)par__l[i]=0.; 

}
/**
 * \brief GetQ
 * Method to retrieve the total collected charge (within RMAX and TAUMAX).
 * @param icorr Correction level 
 * 0 = no correction
 * 1 = lateral correction 
 * 2 = longitudinal correction
 * 3 = ??
 */
float CaloElectron::GetQ(int icorr){
        
    float qx,qy,qtotc;
        
        
    qx = GetQView(icorr,0);
    qy = GetQView(icorr,1);
    qtotc = qx+qy;
    
//     if(qy<0.7*(qtotc/2.)){
//      cout<<"qx "<<qx<<" qy "<<qy<<" qtotc "<<qtotc<<endl;
//      cout<<"maybe this is image track???"<<endl;
//      qtotc=-100.;
//     }
    
    return qtotc;
        
};
/**
 * \brief GetQview
 * Method to retrieve the total collected charge in each view
 * @param icorr Correction level (0= no correction, 1=lateral correction, 2=longitudinal correction)
 * @param iv (0=xview, 1=yview) 
 */
float CaloElectron::GetQView(int icorr,int iv){
        
        
    float TAUMAX = CaloElectron_parameters::Get()->TAUMAX;
    int ipmin = CaloElectron_parameters::Get()->ipmin;
    int ipmax = CaloElectron_parameters::Get()->ipmax;

    float qqqq=0; 

//    if(onlyy && iv==0) return qqqq;
        

    float idmax,ipmaxl;   
    //se faccio così idmax deve essere rispetto al piano da cui inizio davvero
    idmax  = GetShowerMaximum(icorr);
    // ***BUG??***
    ipmaxl = idmax*TAUMAX+iv;
    // ***BUG??***
//    ipmaxl = idmax*TAUMAX;
    if( ipmaxl>(ipmax-ipmin) ) ipmaxl=ipmax-ipmin;
    
    //if(ipmaxl<22) cout<<"GetQView icorr "<<icorr<<" ipmaxl="<<ipmaxl<<endl;
    
    // ----------------------
    // total collected charge
    // ----------------------
    if(icorr==0){
        for(int ip=ipmin; ip<ipmin+ipmaxl; ip++){           
            qqqq += qplane0[ip][iv];
        }
    }
    // -------------------------------------
    // total charge after lateral correction
    // -------------------------------------
    else if(icorr==1){    
        for(int ip=ipmin; ip<ipmin+ipmaxl; ip++){
            qqqq += qplane1[ip][iv];
        }       
    }
    // ------------------------------------------
    // total charge after longitudinal correction
    // ------------------------------------------
    else if(icorr==2){
        for(int ip=ipmin; ip<ipmin+ipmaxl; ip++){
            qqqq += qplane1[ip][iv];
        }       
        qqqq=qqqq*corr__l;
    }
    // ------------------------------------------
    // total charge from longitudinal fit
    // ------------------------------------------
    else if(icorr==3){
        for(int ip=ipmin; ip<ipmin+ipmaxl; ip++){
            qqqq += qplanefit[ip][iv];      
        }       
    }
//    cout << "**** "<<ipmin << " "<<ipmin+ipmaxl<<" "<<qqqq<<" -- "<<qtot<<endl;
    
 
    return qqqq;
};
/**
 * \brief Evaluate shower maximum
 * It gives the number of W-layers traversed at shower maximum. 
 * The charge deposited in each layer is evaluated as the averege
 * between the adjacent x and y planes.
 * NB: the shower depth is evaluated relative to ipmin.
 * It goes from 0 to 22 (more generally, from 0 to ipmax-ipmin ).
 */
float CaloElectron::GetShowerMaximum(float qplane[][2]){
         
    int ipmin = CaloElectron_parameters::Get()->ipmin;
    int ipmax = CaloElectron_parameters::Get()->ipmax;
     
    int   maxlayer = 0;
    float maxq     = qplane[0][1];  //0Y  - first active layer
    float q;
    for(int il=1; il<22; il++){
        q=0.;
        int nxy=0;
        float qx = qplane[il-1][0];
        float qy = qplane[il][1];
        bool ISHIT = true;
        bool maskpl_x = CaloElectron_parameters::Get()->maskpl[il-1][0];
        bool maskpl_y = CaloElectron_parameters::Get()->maskpl[il][1];
        if( tplane[il-1][0] <0 || tplane[il][1]<0 )ISHIT=false;
        if(ISHIT){
//          if( !maskpl[il-1][0] )nxy++;
//          if( !maskpl[il][1]   )nxy++;
            if( !maskpl_x )nxy++;
            if( !maskpl_y   )nxy++;
            if(nxy>0)q=(qx+qy)/nxy;     
            if(q>maxq){
                maxq=q;
                maxlayer=il;
            }
        }
    }    
    if( tplane[21][0] > 0 ){
        q = qplane[21][0]; //21X - last active layer
        if(q>maxq){
            maxq=q;
            maxlayer=22;
        }
    }
    
    maxlayer=maxlayer-ipmin;
    if(maxlayer<0) maxlayer=0;
    if(maxlayer>ipmax-ipmin) maxlayer=ipmax-ipmin;
    return (float)maxlayer;
}

/**
 * \brief Evaluate shower maximum (id of shower maximum plane)
 */
float CaloElectron::GetShowerMaximum(int icorr){
    // ------------------------------------
    // shower maximum before any correction
    // ------------------------------------
    if(icorr==0) return GetShowerMaximum(qplane0);
    // ---------------------------------------
    // shower maximum after lateral correction
    // ---------------------------------------
    if(icorr==1) return GetShowerMaximum(qplane1);
    // ---------------------------------------

    // ---------------------------------------
    // shower maximum from longitudinal fit
    // ---------------------------------------
    if(icorr==2 || icorr==3){
        if( err__l==0 ){
            return par__l[2]*cos(atan(tg))/(WTICK/WX0);
        }else{
            return GetShowerMaximum(qplane1);
        }     
    }
    
    return 0.;
        
}
/**
 * \brief Evaluate shower maximum (in X0 unit, considering track inclination angle)
 */
float CaloElectron::GetShowerMaximumX0(int icorr){
        
        float maxX0=0;
        if(icorr==0 || icorr==1) maxX0=(GetShowerMaximum(icorr)/cos(atan(tg)))*(WTICK/WX0);
        if(icorr==2 || icorr==3) maxX0=par__l[2];
        
        return maxX0;
        
}

float CaloElectron::GetShowerMaximumViewX0(int view,int icorr){

        float maxX0=0;
        if(icorr==0 || icorr==1) maxX0=(GetShowerMaximumView(view,icorr)/cos(atan(tg)))*(WTICK/WX0);
        if(icorr==2 || icorr==3) maxX0=par__l[2];
        
        return maxX0;
        
}
float CaloElectron::GetShowerMaximumView(int view,float qplane[][2]){
    
    int ipmin    = CaloElectron_parameters::Get()->ipmin;
    int ipmax    = CaloElectron_parameters::Get()->ipmax;
    int   maxlayer = 0;
    float maxq     = qplane[0][view];
    for(int il=1; il<22; il++){
        if(qplane[il][view]>maxq){
            maxq     = qplane[il][view];
            maxlayer = il+1-view;
        }
    }
    maxlayer=maxlayer-ipmin;
    if(maxlayer<0) maxlayer=0;
    if(maxlayer>ipmax-ipmin) maxlayer=ipmax-ipmin;
    return (float)maxlayer;
    
}
float CaloElectron::GetShowerMaximumView(int view,int icorr){
    // ------------------------------------
    // shower maximum before any correction
    // ------------------------------------
    if(icorr==0) return GetShowerMaximumView(view,qplane0);
    // ---------------------------------------
    // shower maximum after lateral correction
    // ---------------------------------------
    if(icorr==1) return GetShowerMaximumView(view,qplane1);
    // ---------------------------------------

    // ---------------------------------------
    // shower maximum from longitudinal fit
    // ---------------------------------------
    if(icorr==2 || icorr==3){
        if(err__l==0){
            return par__l[2]*cos(atan(tg))/(WTICK/WX0);
        }else{
            return GetShowerMaximumView(view,qplane1);
        }     
    }
    
    return 0.;
        
}


/**
 * \brief Set event 
 */
bool CaloElectron::Set(CaloLevel1 *cl1,float trackcoordinate[][2]){
         
    if(!cl1 || !trackcoordinate  ){
        cout <<" void CaloElectron::Set(CaloLevel1 *cl1, float tracoo[22][2]) -- ERROR -- input == NULL"<<endl;
        return false;
    }
    //
    CaloElectron_parameters::Get()->cl1 = cl1; //set CaloLevel1 event
    CaloStrip cstrip;


    // -------------------
    // retrieve parameters
    // -------------------
    bool debug    = CaloElectron_parameters::Get()->debug;
    float RMAX    = CaloElectron_parameters::Get()->RMAX;
    int dolatcorr = CaloElectron_parameters::Get()->dolatcorr;
    int dolongfit = CaloElectron_parameters::Get()->dolongfit;
    int ipmin     = CaloElectron_parameters::Get()->ipmin;
    int ipmax     = CaloElectron_parameters::Get()->ipmax;

    int calibmode = CaloElectron_parameters::Get()->calibmode;

//    gBenchmark->Start("conta");

    // --------------------------
    // reset energy-deposit plots
    // --------------------------
    TH2F* h_qtot[2];
    for(int iv=0; iv<2; iv++){
        h_qtot[iv] = CaloElectron_parameters::Get()->h_qtot[iv];
        if(!h_qtot[iv])return false;
        h_qtot[iv]->Reset();
    }

//     cout <<"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% "<< CaloElectron_parameters::Get() << endl;
//     cout <<"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% "<< CaloElectron_parameters::Get()->cl1 << endl;
//     cout <<"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% "<< cl1 << endl;


    if( cl1->istrip < 5 )return false;


//    float trkstrip[22][2]; // strip traversed by the track
    //
    for(int ip=0; ip<22; ip++){
        for(int iv=0; iv<2; iv++){          
            // --------------------------------------
            // determine strip traversed by the track
            // --------------------------------------
//          Float_t coo(0),lad(0),stri(0);
//          Float_t xcorner,ycorner;
//          GetCornerCoord(ip,iv,xcorner,ycorner);
//          if(iv==0) coo = (trackcoordinate[ip][iv]-xcorner)*10.;
//          if(iv==1) coo = (trackcoordinate[ip][iv]-ycorner)*10.;
//          if(coo<80.25) lad=1;
//          if(coo>80.25 && coo<160.75) lad=2;
//          if(coo>160.75) lad=3;
//          stri=32*(lad-1)+(coo-(lad-1)*80.5-0.96)/2.44;
//          Int_t trkstripnew = (Int_t)stri;
//          if(stri-trkstripnew>0.5) trkstripnew=trkstripnew+1;
//          trkstrip[ip][iv]=trkstripnew; 
//          if(coo<0. || coo> 241.)trkstrip[ip][iv]=-1.;
            //
            // non so perche` ma mi sembra che sia sbagliato...
            // ci riprovo a modo mio:
            //
            trkstrip[ip][iv]=-1.;
            float xy[2];
            for(int is=0; is<96; is++){     
                cstrip.Set(iv,ip,is);
                xy[0]   = cstrip.GetX();
                xy[1]   = cstrip.GetY();
                if( 
                    trackcoordinate[ip][iv] >= xy[iv] - 0.5*PITCH && 
                    trackcoordinate[ip][iv] <  xy[iv] + 0.5*PITCH && 
                    true){
                    trkstrip[ip][iv] = is;
                    break;
                }
            };
        }
    }


//    cout << " 1 ------------------ init"<<endl;
//    gBenchmark->Print("conta");

    
    //<><><><><><><><><><><><><><><><><><><><><><><><>
    // LOOP OVER HIT STRIPS
    // - to set charge collected in each plane
    // - to fill energy-deposit plot (if required)
    //<><><><><><><><><><><><><><><><><><><><><><><><>

    int isimu = CaloElectron_parameters::Get()->isimu;
    if(isimu==0)cstrip = CaloStrip(cl1,false);//conferma da Emiliano
    if(isimu==1)cstrip = CaloStrip(cl1,true);

    float xy[2];
    for(int ih=0;ih<cl1->istrip;ih++){
        int iv=-1;
        int ip=-1;
        int is=-1;
        cl1->DecodeEstrip(ih,iv,ip,is);
        cstrip.Set(iv,ip,is);
//      enestr[is][ip][iv] = cstrip.GetE();
        xy[0]              = cstrip.GetX();
        xy[1]              = cstrip.GetY();
        if( TMath::Abs(trackcoordinate[ip][iv]-xy[iv]) < RMAX )
            qplane0[ip][iv] += cstrip.GetE(); //<< QPLANE0
        bool maskpl = CaloElectron_parameters::Get()->maskpl[ip][iv];
        if( 
            ip>=ipmin && 
            ip<ipmax  && 
            !maskpl   &&
            true){
            qtot += cstrip.GetE();//<< QTOT
            //
//          h_qtot[iv]->Fill( (float)(is+iv*96), (float)ip, cstrip.GetE());
            // retrieve the coordinate of the bottom-left corner of the plane
            float corner[] = {0.,0.};
            GetCornerCoord(ip,iv,corner[0],corner[1]);
            h_qtot[iv]->Fill( xy[iv]-corner[iv], (float)ip, cstrip.GetE());
        }
//      cout <<(float)(is+iv*96)<<(float)ip<<cstrip.GetE()<<endl; 
    }
    
//    cout << " 2 ------------------ read charge"<<endl;
//    gBenchmark->Print("conta");
    
    // ------------------------------------------------
    // set track intersection coordinate in each plane
    // and evaluate shower depth
    // ------------------------------------------------

    float trkcoordz[22][2]; ///< track coordinates (PAMELA r.s.)

    for(int ip=0; ip<22; ip++){
        for(int iv=0; iv<2; iv++){
            cstrip.Set(iv,ip,0);
            trkcoordz[ip][iv] = cstrip.GetZ();
        }
    }    
    tgx = (trackcoordinate[0][0]-trackcoordinate[21][0])/(trkcoordz[0][0]-trkcoordz[21][0]);
    tgy = (trackcoordinate[0][1]-trackcoordinate[21][1])/(trkcoordz[0][1]-trkcoordz[21][1]);
    tg  = sqrt(pow(trackcoordinate[0][1]-trackcoordinate[21][1],2)+pow(trackcoordinate[0][0]-trackcoordinate[21][0],2))/(trkcoordz[0][1]-trkcoordz[21][1]);

    bool TRACK_OUT = false;
    for(int ip=0; ip<22; ip++){
        for(int iv=1; iv>=0; iv--){
            if(iv==0){
                trkcoordx[ip][0] = trackcoordinate[ip][0]; //track intersection coordinate
                trkcoordy[ip][0] = trackcoordinate[ip][1] + (trkcoordz[ip][0]-trkcoordz[ip][1])*tgy; //track intersection coordinate
            }else{           
                trkcoordx[ip][1] = trackcoordinate[ip][0] + (trkcoordz[ip][1]-trkcoordz[ip][0])*tgx; //track intersection coordinate
                trkcoordy[ip][1] = trackcoordinate[ip][1]; //track intersection coordinate
            }
            // ---------------------------------------------
            // if the track go out of the calorimeter, break        
            // ---------------------------------------------
            if( trkstrip[ip][iv] == -1 ){
                //TRACK_OUT = true;
                if(debug)cout << "ip="<<ip <<" iv="<<iv<<" --> TRACK OUT";
            }
            if( !TRACK_OUT ){
                // shower depth in units of radiation length (TO BE CHECKED)
                tplane[ip][iv]  = GetNWLayers(ip,iv) * sqrt(tgx*tgx+tgy*tgy+1) * WTICK/WX0;//sqrt(tgx*tgx+tgy*tgy+1)=1/cos(atan(tg)
                qplane1[ip][iv] = qplane0[ip][iv]; // initialize qplane1
            }else return false; 
        }
    }    

//    cout << " 3 ------------------ get track coordinates"<<endl;
//    gBenchmark->Print("conta");

    if(GetQ(0)<1){     
        cout<<"Don't go on ..GetQ(0) "<<GetQ(0)<<endl;
        return false;
    }
    if(GetShowerMaximum(0)<1 || GetShowerMaximum(0)>21){     
        cout<<"Don't go on ..GetShowerMaximum(0) "<<GetShowerMaximum(0)<<endl;
        return false;
    }
 
    // -----------------------------------------
    // fill calibration histograms (if required)
    // -----------------------------------------

    if( calibmode )Calibration_Fill();

    // -----------------------------------------
    // apply lateral correction
    // -----------------------------------------
    
    if(dolatcorr==0) return true;
    
    float maxlayer;
    int niter;
    float corr;
    
    maxlayer = 0.;
    corr     = 0.;
    niter    = 0;
    do{
        niter++;
        maxlayer = GetShowerMaximum(1);
        corr = ApplyLateralCorrection(maxlayer);
    }while( maxlayer != GetShowerMaximum(1) && niter < 3);

    if(debug && niter>1)cout<<"WARNING! - lateral correction iterated "<<niter<<" times"<<endl;

    if(debug)cout << "LATERAL      CORRECTION FACTOR            = "<<corr<<endl;

//    cout << " 4 ------------------ lateral correction"<<endl;
//    gBenchmark->Print("conta");

    if(GetQ(1)<1){     
        cout<<"Don't go on ..GetQ(1) "<<GetQ(1)<<endl;
        return false;
    }
    
    // -----------------------------------------
    // apply longitudinal correction
    //possible dolongfit values are:
    //-1 don't apply longitudinal correction
    //0 apply to q0 (charge collected by each plane inside RMAX)
    //1 apply to q1 (charge collected by each plane with lateral collection inside RMAX )
    // -----------------------------------------
    
    if(dolongfit<0) return true;
    
    int ret = ApplyLongitudinalFit(dolongfit);

//    cout << " 5 ------------------ longitudinal fit"<<endl;
//    gBenchmark->Print("conta");

    if( ret==0 ){ 
        corr    = ApplyLateralCorrection( GetShowerMaximum(2) );
        ret     = ApplyLongitudinalFit(dolongfit);
        corr__l = GetLongitudinalCorrection();  
        chi2    = ProfileTest();
    }
    if(debug){
        cout << "LATERAL      CORRECTION FACTOR (iterated) = "<<corr<<endl;
        cout << "LONGITUDINAL CORRECTION FACTOR            = "<<corr__l<<endl;
    }
//    cout << " 6 ------------------ iteration"<<endl;
//    gBenchmark->Print("conta");
     
//     if(ret==0 || ret==10){
//      if(debug)cout << ">> Longitudinal correction from level "<<dolongfit<<" corr__l "<<corr__l<<endl;
//      if(GetLongitudinalFcnout()>600){
//          cout<< "too big GetLongitudinalFcnout() "<<GetLongitudinalFcnout()<<endl;
//          //return false;
//      }
//     }else{
//      cout<< ">> problems with Longitudinal correction!"<<endl;
//      //return false;
//     }
     
    
    return true;

}
/**
 * \brief Set event 
 */
bool CaloElectron::Set(CaloLevel1 *cl1, CaloTrkVar* ctrk){
    
    if(!cl1 ||!ctrk){
        cout <<" void CaloElectron::Set(CaloLevel1 *cl1, CaloTrkVar* ctrk) -- ERROR -- input == NULL"<<endl;
        return false;
    }

    float trackcoordinate[22][2];

    for(int ip=0; ip<22; ip++){
      for(int iv=0; iv<2; iv++){
        trackcoordinate[ip][iv]=ctrk->tbar[ip][iv];
//      trkstrip[ip][iv]=ctrk->tibar[ip][iv];
      }
    }
    return Set(cl1,trackcoordinate);

};
 /**
 * \brief Set event 
 */
bool CaloElectron::Set(PamLevel2 *l2, int ntr){

    CaloLevel1* cl1 = l2->GetCaloLevel1();
    CaloTrkVar* ctrk = 0;
    if( ntr < l2->GetTrkLevel2()->GetNTracks() )ctrk = l2->GetTrack(ntr)->GetCaloTrack();    
    return Set(cl1,ctrk);

};


/**
 * \brief Method to apply lateral correction
 * @param maxl shower maximum, in units of W-layers
 */
float CaloElectron::ApplyLateralCorrection( float maxl ){

    bool debug   = CaloElectron_parameters::Get()->debug;
    float TAUMAX = CaloElectron_parameters::Get()->TAUMAX;
//    float RMAX   = CaloElectron_parameters::Get()->RMAX;
    int ipmin    = CaloElectron_parameters::Get()->ipmin;
    int ipmax    = CaloElectron_parameters::Get()->ipmax;

    if( debug ) cout << "CaloElectron::ApplyLateralCorrection(maxl="<<maxl<<")"<<endl;

    if( maxl    <=0 )return 0.;
    if( GetQ(0) <=0 )return 0.;

    for(int ip=ipmin; ip<ipmax; ip++){
        for(int iv=1; iv>=0; iv--){

            double corr=1.;
            // evaluate the shower-depth of the plane
            float nlayer = (float)GetNWLayers(ip,iv);
            double tau   = (double)(nlayer-ipmin)/(double)maxl;
            int ibintau  = (int)(tau/0.25);
            if( ibintau>=12 ) ibintau=11;
            // retrieve the coordinate of the bottom-left corner of the plane
            float xcorner = 0;
            float ycorner = 0;
            GetCornerCoord(ip,iv,xcorner,ycorner);
            // evaluate track coordinate relative to the plane corner
            double xpl = (double)trkcoordx[ip][iv] - (double)xcorner;
            double ypl = (double)trkcoordy[ip][iv] - (double)ycorner;
//          trkcoordxpl[ip][iv] = xpl;
//          trkcoordypl[ip][iv] = ypl;
            //cout<<"ip "<<ip<<" iv "<<iv<<" xpl "<<xpl<<" ypl "<<ypl<<" trkcoordx "<<trkcoordx[ip][iv]<<" trkcoordy "<<trkcoordy[ip][iv]<<endl;
            // -----------------------------------------------------
            // evaluate lateral correction if: 
            // - the plane is traversed by the track (tplane!=-100)
            // - the plane is not the first active layer (tplane =0)
            // - the plane correspond to tau lower then TAUMAX
            // -----------------------------------------------------
            corr=1.;//correction has no effect if tau>TAUMAX                    
            bool maskpl = CaloElectron_parameters::Get()->maskpl[ip][iv];
//          if( !maskpl[ip][iv] ){
            if( !maskpl ){
                if( tplane[ip][iv] > 0 && tau<TAUMAX ){
                    int iiv=0;
                    if(iv==0) iiv=0;
                    if(iv==1 && ip%2) iiv=1;
                    if(iv==1 && !ip%2) iiv=2;
                    corr = GetLateralCorrection(xpl,ypl,ibintau,iiv);               
                }
            }            
            qplane1[ip][iv] = qplane0[ip][iv] / corr;    
//          if(debug && corr<0.5){ 
//              cout <<" qp0 "<<qplane0[ip][iv];
//              cout <<" qp1 "<<qplane1[ip][iv];
//              cout <<" ip  "<<ip<<" iv "<<iv<<" tau "<<tau<<endl;
//          }
        }
    }

    return GetQ(1)/GetQ(0);

//     corr = GetQ(1)/GetQ(0);
        
//      //questo disegna due grafici con la carica rilasciata nelle viste x e y prima e dopo la correzione 
//     if(dodrawq0q1){
    
//      Double_t tgx[22],tgy[22],qx[22],qy[22],qx1[22],qy1[22];
        
//      for(int ip=0; ip<22; ip++){
//          tgx[ip]=tplane[ip][0];
//          tgy[ip]=tplane[ip][1];
//          qx[ip]=qplane0[ip][0];
//          qy[ip]=qplane0[ip][1];
//          qx1[ip]=qplane1[ip][0];
//          qy1[ip]=qplane1[ip][1];
//      }
        
//      TGraph *gqx= new TGraph(22,tgx,qx);
//      TGraph *gqy= new TGraph(22,tgy,qy);
//      TGraph *gqx1= new TGraph(22,tgx,qx1);
//      TGraph *gqy1= new TGraph(22,tgy,qy1);
//      TCanvas a;a.Divide(2,1);a.cd(1);
//      a.cd(1);
//      gqx1->SetMarkerColor(2);gqx1->SetMarkerStyle(7);gqx1->Draw("AP");
//      gqx->SetMarkerColor(1);gqx->SetMarkerStyle(7);gqx->Draw("Psame");
//      a.cd(2);
//      gqy1->SetMarkerColor(2);gqy1->SetMarkerStyle(7);gqy1->Draw("AP");
//      gqy->SetMarkerColor(1);gqy->SetMarkerStyle(7);gqy->Draw("Psame");
//      a.Update();
//      int salvalo;
//      cout<<"save dodrawq0q1 (1 yes, 0 no)??";
//      cin>>salvalo;
//      if(salvalo>1)       a.Print(Form("./gq0q1_%d.root",salvalo));
//      else if(salvalo==1) a.Print("electronq0q1.root");
        
//     }// end of   if(dodrawq0q1)
    
//     //questo disegna la distribuzione laterale sul un certo piano e vista
//     //con la posizione dei gap e la funzione usata per correggere
//     if(dodrawq0q1){
    
//      int iv,ip,iiv(0);
//      float iwmax,iw;
//      cout<<"which plane and view?"<<endl;
//      cin>>ip>>iv;
//      if(iv==0)iiv=0;
//      if(iv==1 && ip%2) iiv=1;
//      if(iv==1 && !ip%2) iiv=2;
//      iwmax=GetShowerMaximum(0);
//      iw=ip+(iv-1);
//      int ibintau=(int)((iw/iwmax)/0.25);
//      cout<<"plane "<<ip<<" iv "<<iv<<" iiv "<<iiv<<" iwmax "<<iwmax<<" iw "<<iw<<" ibintau "<<(iw/iwmax)/0.25<<" "<<ibintau<<" ..setting parameter "<<endl;
//      cout<<file_rt[ibintau][iiv]<<" "<<file_p[ibintau][iiv]<<" "<<file_rc[ibintau][iiv]<<endl;
//      if(ibintau>=12) ibintau=11;
//      int nparam=5;
//      TF1 ff = TF1("fradpro",fradpro,-RMAXsp,RMAXsp,nparam);
//      Double_t param[5];
//      param[0]=file_rt[ibintau][iiv];
//      param[1]=file_p[ibintau][iiv];
//      param[2]=file_rc[ibintau][iiv];
//      param[3]=qplane0[ip][iv];
//      //param[4]=trkcoordx[ip][iv];
//      param[4]=0.;
//      cout<<"plane "<<ip<<" iv "<<iv<<" ibintau "<<ibintau<<" ..setting parameter "<<endl;
//      for(int i=0;i<nparam;i++){
//          ff.SetParameter(i,param[i]);
//          cout<<" "<<i<<" -> "<<param[i]<<endl;
//      }
    
    
//      CaloStrip st = CaloStrip();
//      CaloStrip st1 = CaloStrip();
//      CaloStrip st2 = CaloStrip();
//      if(isimu==0)st.UseStandardAlig();
//      if(isimu==1)st.UseMechanicalAlig();     
//      TH1F *hgap=new TH1F("hgap","hgap",34,-RMAXsp,RMAXsp);
//      float posgap1(0),posgap2(0);
//      st1.Set(iv,ip,31);st2.Set(iv,ip,32);
//      if (iv==0)     posgap1=(st1.GetX()+st2.GetX())/2.-trkcoordx[ip][iv];
//      if (iv==1)     posgap1=(st1.GetY()+st2.GetY())/2.-trkcoordy[ip][iv];
//      st1.Set(iv,ip,63);st2.Set(iv,ip,64);
//      if (iv==0)     posgap2=(st1.GetX()+st2.GetX())/2.-trkcoordx[ip][iv];
//      if (iv==1)     posgap2=(st1.GetY()+st2.GetY())/2.-trkcoordy[ip][iv];
        
//      hgap->Fill(posgap1,100.);
//      hgap->Fill(posgap2,100.);
        
        
//      TH1F *h=new TH1F("h","h",34,-RMAXsp,RMAXsp);
//      float binsize=8./34.;
//      double e(0),d(0),etot(0);
//      for(int is=0;is<96;is++){
//          e=enestr[is][ip][iv];
//          if(e>0.){
//              st.Set(iv,ip,is);
//              if (iv==0) d=st.GetX()-trkcoordx[ip][iv];
//              if (iv==1) d=st.GetY()-trkcoordy[ip][iv];
//              if(TMath::Abs(d)<RMAXsp){
//                  etot=etot+e;
//                  cout<<"filling is "<<is<<" d "<<d<<" e "<<e<<endl;
//                  h->Fill(d,e/binsize);
//              }
//          }
//      }
//      TCanvas a;a.Divide(2,1);a.cd(1);
//      h->Draw();//ff.Draw("same");
//      a.cd(2);ff.Draw(); h->Draw("same");hgap->SetFillColor(2);hgap->Draw("same");
//      a.Update();
//      cout<<"plane "<<ip<<" iv "<<iv<<" corr "<<qplane1[ip][iv]/qplane0[ip][iv]<<endl;
//      cout<<" h->Integral(width) "<<h->Integral("width")<<" h->Integral() "<<h->Integral()<<" ff.Integral() "<<ff.Integral(-12.,12.,param,1.e-12)<<" etot "<<etot<<endl;
        
//      int salvalo;
//      cout<<"save dodrawlat (1 yes, 0 no)??";
//      cin>>salvalo;
//      //potrei provare a calcolare il fattore di correzione su ogni piano da questo confronto e
//      //confrontarlo con quello che viene dal mio metodo
        
//     }
          
//     return corr;
    
}
////////////////////////////////////////////////////////////////////////            
/**
 * 1-dimension function describing lateral distribution of the 
 * shower as viewed by calorimeter
 * (projection of 3d function in one direction)
 *  
 *  xi[0] = x or y coordinate relative to shower axis
 *  parmin[0] = rt
 *  parmin[1] = p
 *  parmin[2] = rc
 *  parmin[3] = norm
 *  parmin[4] = x0,y0
 *
 */
////////////////////////////////////////////////////////////////////////            
double fradpro(Double_t *xi, Double_t *parmin) {

        double fradpromin2,p,rt,rc,es,x,pig,norm,c;
        x   = *xi;
        pig = acos(-1.);

        rt  = parmin[0];
        p   = parmin[1];
        rc  = parmin[2];

        norm = parmin[3];
        c    = parmin[4];
        x    = x-c;
        es   = 1.5;
        fradpromin2 = p*pig*pow(rc,2)/pow((pow(x,2)+pow(rc,2)),es);
        fradpromin2 = fradpromin2+(1-p)*pig*pow(rt,2)/pow((pow(x,2)+pow(rt,2)),es);
        fradpromin2 = norm*fradpromin2/(2*pig);     
        //cout<<"x,fradpromin2 "<< x<<" "<<fradpromin2  <<endl;
  return fradpromin2;
}
////////////////////////////////////////////////////////////////////////            
/**
 * 2-dimension function describing lateral distribution of the shower 
 *  
 *  xyi[0] = x coordinate relative to shower axis
 *  xyi[1] = y coordinate relative to shower axis
 *  parmin[0] = rt
 *  parmin[1] = p
 *  parmin[2] = rc
 *
 */
////////////////////////////////////////////////////////////////////////            
double frad(double *xyi, double *parmin) {
    
    double fradout,p,rt,rc,pig,r2;
        
    r2 = pow(xyi[0],2)+pow(xyi[1],2);
    pig = acos(-1.);

    rt = parmin[0];
    p = parmin[1];
    rc = parmin[2];

    fradout = p*2*pow(rc,2)/pow((r2+pow(rc,2)),2);
    fradout = fradout+(1-p)*2*pow(rt,2)/pow((r2+pow(rt,2)),2);
    fradout = fradout/(2*pig);
        
    return fradout;
};
////////////////////////////////////////////////////////////////////////            
/**
 * 2-dimension function describing lateral distribution of the shower, with 
 * the second variable implemented as a parameter 
 *  
 *  xyi[0] = x or y coordinate relative to shower axis
 *  parmin[0] = rt
 *  parmin[1] = p
 *  parmin[2] = rc
 *  parmin[3] = other coordinate
 *
 */
////////////////////////////////////////////////////////////////////////            
double fradx(double *xyi, double *parmin) {
    
    double fradout,p,rt,rc,pig,r2;
        
    r2 = pow(xyi[0],2)+pow(parmin[3],2);
    pig = acos(-1.);

    rt = parmin[0];
    p = parmin[1];
    rc = parmin[2];

    fradout = p*2*pow(rc,2)/pow((r2+pow(rc,2)),2);
    fradout = fradout+(1-p)*2*pow(rt,2)/pow((r2+pow(rt,2)),2);
    fradout = fradout/(2*pig);
        
    return fradout;
};


/**
 * \brief Evaluate lateral correction factor
 */
double CaloElectron::GetLateralCorrection(double xpl, double ypl, int ibintau, int iv){

    TF2 ff = TF2("frad",frad,-25.,25.,-25.,25.,3);
    
    int ivqui = iv;
    
    float file_rt = CaloElectron_parameters::Get()->file_rt[ibintau][ivqui];
    float file_p  = CaloElectron_parameters::Get()->file_p[ibintau][ivqui];
    float file_rc = CaloElectron_parameters::Get()->file_rc[ibintau][ivqui];

    ff.SetParameter(0,file_rt);
    ff.SetParameter(1,file_p);
    ff.SetParameter(2,file_rc);
    
    double corr = GetLateralCorrection(xpl,ypl,&ff,iv);
    return corr;
};
/**
 * \brief Evaluate lateral correction factor
 */
double CaloElectron::GetLateralCorrection(double xpl, double ypl, TF2* ff,int iv){

    float RMAXsp = CaloElectron_parameters::Get()->RMAXsp;
    float RMAX   = CaloElectron_parameters::Get()->RMAX;

    double x1,x2,x3,x4,x5,ps[3];
    //    double x1,x2,x3,x4,x5,x6,ps[3];
    x1 = DEAD;
    x2 = x1 + PITCH*32;
    x3 = x2 + DEAD + GLUE + DEAD;
    x4 = x3 + PITCH*32;
    x5 = x4 + DEAD + GLUE + DEAD;
    //    x6 = x5 + PITCH*32;
    ps[0] = x1;
    ps[1] = x3;
    ps[2] = x5;
    double integgap(0),integgaptmp(0),integnogap(0),correggif(0);
    double epsil=0.0001;
    for (Int_t ii=0;ii<3;ii++){
        for (Int_t jj=0;jj<3;jj++){  
            double xmi = ps[ii] - xpl;
            double xma = xmi + PITCH*32;
            double ymi = ps[jj] - ypl;
            double yma = ymi + PITCH*32;
            int dointeg=1;
            if(iv==0){
                if(xmi<-RMAXsp && xma>-RMAXsp) xmi=-RMAXsp;
                if(xma>RMAXsp && xmi<RMAXsp) xma=RMAXsp;
                if(xmi>RMAXsp || xma<-RMAXsp)dointeg=0;
            }
            if(iv==1 || iv==2){
                if(ymi<-RMAXsp && yma>-RMAXsp) ymi=-RMAXsp;
                if(yma>RMAXsp && ymi<RMAXsp) yma=RMAXsp;
                if(ymi>RMAXsp || yma<-RMAXsp)dointeg=0;
            }
            if(dointeg==1){
                integgaptmp=ff->Integral(xmi,xma,ymi,yma,epsil);
                integgap += integgaptmp;
            }
        }
    }
    
    double RPMAX=40.;
    double xmin(0),xmax(0),ymin(0),ymax(0);
    if(iv==0){
        ymin=-RPMAX;
        ymax=RPMAX;
        xmin=-RMAX;
        xmax=RMAX;
    }
    if(iv==1 || iv==2){
        xmin=-RPMAX;
        xmax=RPMAX;
        ymin=-RMAX;
        ymax=RMAX;
    }
    integnogap=ff->Integral(xmin,xmax,ymin,ymax,epsil);
    
    correggif=integgap/integnogap;
    
    return correggif;

}
// ---------------------------------------------------------------------
// TMinuit does not allow fcn to be a member function, and the function
// arguments are fixed, so the one of the only ways to bring the data  
// into fcn is to declare the data as global.
// ---------------------------------------------------------------------
//
Int_t iuse;
Double_t zfit[44],tfit[44],errorzfit[44],errortfit[44];
//Float_t corr,m,q,lp1m,lp2m,sig1,sig2;Float_t parglobali[4];
//______________________________________________________________________________
Double_t func(Double_t *tin,Double_t *par)
{

    Double_t gaa,a,tm,et,value,t,t0;
    t  = *tin;
    et = par[0];
    a  = par[1];
    tm = par[2];
    t0 = par[3];
    gaa = TMath::Gamma(a);
    
    value = et*((a-1)/(tm*gaa))*pow(((a-1)*(t-t0)/tm),(a-1))*exp(-(a-1)*(t-t0)/tm);

    return value;
}
//______________________________________________________________________________
// fcn passes back f = chisq the function to be minimized.
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{

   int iusequi=0;
   Double_t chisq = 0.;
   for (int i=0;i<iuse; i++) {
      if(tfit[i]>=par[3]){
       iusequi++;
       Double_t fu=func(&tfit[i],par);
       Double_t delta  = (zfit[i]-fu)/errorzfit[i];
       chisq += delta*delta;
      }            
   }
   
   //cout<<"iusequi,par,fcn "<<iusequi<<" "<<par[0]<<" "<<par[1]<<" "<<par[2]<<" "<<chisq<<endl;
   f = chisq;

}

/**
 * \brief Method to apply longitudinal correction
 * @param longto Fit options:
 * -1 don't apply longitudinal correction
 *  0 apply to q0 (charge collected by each plane inside RMAX)
 *  1 apply to q1 (charge collected by each plane with lateral collection inside RMAX )
 */
int CaloElectron::ApplyLongitudinalFit(int longto){

    cout << "CaloElectron::ApplyLongitudinalFit("<<longto<<") ---> temporarily commented (problems with TMinuit!!)"<<endl;

//     bool debug   = CaloElectron_parameters::Get()->debug;
//     float TAUMAX = CaloElectron_parameters::Get()->TAUMAX;
// //    float RMAX   = CaloElectron_parameters::Get()->RMAX;
//     int ipmin    = CaloElectron_parameters::Get()->ipmin;
//     int ipmax    = CaloElectron_parameters::Get()->ipmax;

//     if( debug ) cout << "CaloElectron::ApplyLongitudinalFit("<<longto<<")"<<endl;

//     float ipmaxl = ipmax-ipmin;
    
//     /// ?????
//     if(GetShowerMaximum(longto)*TAUMAX<ipmax-ipmin){
//      corr__l        = 1; 
//      err__l = 10;
//      if(debug){
//          cout<<"GetShowerMaximum(longto) "<<GetShowerMaximum(longto)<<" TAUMAX "<<TAUMAX<<" ipmax-ipmin "<<ipmax-ipmin<<endl;
//      }else{
//          return err__l;     
//      }
//     }

// //  initialize TMinuit with a maximum of 5 params
//     TMinuit *minuit = new TMinuit(4); 

// //    cout << minuit << endl;
//     if(!debug) minuit->SetPrintLevel(-1);
//     if( debug) minuit->SetPrintLevel(-1);
//     minuit->SetFCN(fcn);
      

// //    cout << "<><><><><><><><><><><><><"<<endl;
// //  --------------------------------------------    
// //  Set layers to be used for the fit
// //  --------------------------------------------    
//     Double_t arglist[10];
//     Int_t ierflg = 0;
//     iuse = 0;//n.layers included in the fit
//     for(int ip=ipmin; ip<ipmin+ipmaxl; ip++){
//      //for(int iv=0; iv<2; iv++){
//      for(int iv=1; iv>=0; iv--){
//          bool maskpl = CaloElectron_parameters::Get()->maskpl[ip][iv];
//          if(!maskpl){
//              if(longto==0) zfit[iuse] = qplane0[ip][iv];
//              if(longto==1) zfit[iuse] = qplane1[ip][iv];
//              tfit[iuse] = (double)tplane[ip][iv]-(double)tplane[ipmin][1];     
//              errortfit[iuse] = 0.01;
//              errorzfit[iuse] = sqrt( zfit[iuse] );
//              if(zfit[iuse]<0.0001) errorzfit[iuse] = 1.;
// //           if(debug)cout<<"iuse,ip,iv,tplane,t,z "<<iuse<<" "<<ip<<" "<<iv<<" "<<tplane[ip][iv]<<" "<<tfit[iuse]<<" "<<zfit[iuse]<<endl;
//              iuse++;
//          }
//      }
//     }
//     if(debug)cout << "N.layers inlcuded in the fit: "<<iuse<<endl;

// //  --------------------------------------------    
// //  Set starting values and step sizes for param
// //  --------------------------------------------    
//     Double_t size     = (double)tplane[2][1]-(double)tplane[1][1];
//     Double_t initpar0 = 0;
//     for(int ip=ipmin+1; ip<ipmin+ipmaxl; ip++) 
//      initpar0 = initpar0 + size*(zfit[ip*2]+zfit[ip*2-1])/2;
//     Double_t initpar2 = GetShowerMaximumX0(longto);
//     Double_t initpar1 = pow(initpar2,0.9);
//     Double_t initpar3 = 0.;
    
//     Double_t vstart[4] = {initpar0, initpar1, initpar2, initpar3};
//     Double_t step[4]   = {1.,0.1,1.,0.1};
//     Double_t limmin[4] = {0,0.,0.,0};
//     Double_t limmax[4] = {0,20.,20.,0};
//     string parname[4]  = {"E","alpha","tm","t0"};
    
// //  --------------------------------------------    
// //  Now ready for minimization step with limits
// //  --------------------------------------------    
//     if(debug)cout << "--------------------- First step (with limits): "<<endl;
//     for (int i=0; i<4; i++)    
//      minuit->mnparm(i, parname[i].c_str(),vstart[i],step[i],limmin[i],limmax[i],ierflg);
//     minuit->FixParameter(3);//fix starting point   
//     arglist[0] = 500;
//     minuit->mnexcm("MIGRAD", arglist ,1,ierflg);
//     err__l = ierflg;
//     if(ierflg>0){
//      if(debug){
//          cout<<"ierflg "<<ierflg<<" ma prova ad andare avanti"<<endl;
//      }else{
//          return err__l;
//      }    
//     } 
   
//     Double_t vstart2[4];
//     Double_t vend[4],errvend[4];
    
//     for(int i=0; i<4; i++){
//      minuit->GetParameter(i,vend[i],errvend[i]);
//      if(debug)cout<<setw(7)<<parname[i].c_str()<<" v-start "<<setw(10)<<vstart[i]<<" --> v-end "<<setw(10)<<vend[i]<<" +- "<<errvend[i]<<endl;
        
//     }
        
// //     // ===========================================
// //     // NB NB NB !!!!!!!!!!!!!
// //     // temporaneamente eseguo il secondo step
// //     // solo in modalita` debug
// //     // (chiedere spiegazioni a elena)
// //     // ===========================================
// //     if( debug && err__l ){
        
// //  --------------------------------------------    
// //  Now ready for minimization step without limits
// //  --------------------------------------------    
//     if(debug)cout << "--------------------- Second step (without limits): "<<endl;
//     minuit->mnfree(0); //restore all fixed parameters
//     for (int i=0; i<4; i++){
//      vstart[i] = vend[i];
//      minuit->mnparm(i, parname[i].c_str(),vstart[i],step[i],0.,0.,ierflg);
//     }
//     minuit->FixParameter(3);//fix starting point   
//     arglist[0] = 500;
//     minuit->mnexcm("MIGRAD", arglist ,1,ierflg);
//     err__l = ierflg;
//     if(ierflg>0){
//      if(debug){
//          cout<<"ierflg "<<ierflg<<" ma prova ad andare avanti"<<endl;
//      }else{
//          return err__l;
//      }    
//     } 
//     for(int i=0; i<4; i++){
//      minuit->GetParameter(i,vend[i],errvend[i]);
//      if(debug)cout<<setw(7)<<parname[i].c_str()<<" v-start "<<setw(10)<<vstart[i]<<" --> v-end "<<setw(10)<<vend[i]<<" +- "<<errvend[i]<<endl;
//     }
    
// //     }

// //  --------------------------------------------    
// //  retrive parameters 
// //  --------------------------------------------    
// //    Double_t tup = TAUMAX*initpar2;
// //    TF1 *ffunc = new TF1("ffunc",func,0.,tup,4);
// //    ffunc->SetLineWidth(1);
//     for(int i=0; i<4; i++){
//      par__l[i] = vend[i];
//      if(par__l[i]<0 || par__l[2]>22 || par__l[1]>50){
//          cout<<" converge a valori assurdi... escluso"<<endl;
//          err__l = 5;
//          return err__l;
//      }
//     }

//     if(debug)cout << "--------------------- DONE! "<<endl;
    
    
//     for(int ip=0; ip<22; ip++){
//      for(int iv=0; iv<2; iv++){
//          bool maskpl = CaloElectron_parameters::Get()->maskpl[ip][iv];
//          if(!maskpl){
//              double tpla = (double)tplane[ip][iv];
//              Double_t fu1=func(&tpla,vend);
//              qplanefit[ip][iv] = fu1;
//          }
//      }
//     }
    
//     Double_t graddummy[4];
//     double outfcn;
//     minuit->Eval(4,graddummy,outfcn,vend,1);   
    
//     if(iuse>3)chi2__l = (float)outfcn/(float)(iuse-3);
//     else      chi2__l = 0.;
    
//     delete minuit; 
 
    return err__l;
    
}
/*
 * \brief Evaluate longitudinal correction 
 */
float CaloElectron::GetLongitudinalCorrection(){

    bool debug   = CaloElectron_parameters::Get()->debug;
    float TAUMAX = CaloElectron_parameters::Get()->TAUMAX;
    int ipmin    = CaloElectron_parameters::Get()->ipmin;
    int ipmax    = CaloElectron_parameters::Get()->ipmax;
//  --------------------------------------------    
//  evaluate longitudinal 
//  --------------------------------------------    
    Double_t corr = 0.;
    Double_t t3,intecalo3,intecalo,tmaxcalo;
    if( err__l == 0 ){
        t3 = par__l[3] + TAUMAX*par__l[2];
        TF1 *ffunc = GetFunc_Longitudinal();
        intecalo3 = ffunc->Integral(par__l[3],t3);
//      if(debug)cout<<"integro da "<<par__l[3]<<" fino a t3 "<<t3<<" trovo intecalo3 "<<intecalo3<<endl;
        tmaxcalo = (double)tplane[ipmax-1][0]-(double)tplane[ipmin][1];
        if(t3<tmaxcalo){//if(t3<tplane[21][0]){
            intecalo=intecalo3;
//          if(debug)cout<<"siccome tmaxcalo "<<tmaxcalo<<" >t3 "<<" intecalo=intecalo3 "<<endl;
        }else{
            intecalo=ffunc->Integral(par__l[3],tmaxcalo);
//          if(debug)cout<<"integro da "<<par__l[3]<<" fino a tmaxcalo "<<tmaxcalo<<" trovo intecalo3 "<<intecalo<<endl;
        }
        corr = intecalo3/intecalo;   
    }else{
        corr=0;
    }

    corr__l = (float)corr;

    if(debug)  cout<<"Longitudinal correction:  "<<corr__l<<endl;   

    return corr__l;

}
/////////////////////////////////////////////////////////////////////
//
//
//
//
/////////////////////////////////////////////////////////////////////
TGraphErrors *CaloElectron::GetGraph_Longitudinal(){
    TGraphErrors *graph = new TGraphErrors(iuse ,tfit,zfit,errortfit,errorzfit);
    graph->SetMarkerSize(0.8);
    graph->SetMarkerStyle(20);    
    return graph;
}
//
TGraphErrors *CaloElectron::GetGraph_Longitudinal_Fit(){
    double zfunc[44];
    double errorzfunc[44];
    TF1 *fu = new TF1("fu",func,0.,tfit[iuse-1],4);
    for(int i=0; i<4; i++)fu->SetParameter(i,par__l[i]);
    for(int i=0; i<iuse; i++){
        zfunc[i] = fu->Eval(tfit[i]);;
        errorzfunc[i] = 0.;
        if(zfunc[i]>0.)errorzfunc[i] = sqrt(zfunc[i]);
    }    
    delete fu;
    TGraphErrors *graph = new TGraphErrors(iuse,tfit,zfunc,errortfit,errorzfunc);
    graph->SetMarkerSize(0.8);
    graph->SetMarkerColor(2);
    graph->SetMarkerStyle(20);    
    return graph;
}
//
TGraphErrors *CaloElectron::GetGraph_Longitudinal_Integral(){
    double zfitint[44];
    double errorzfitint[44];
    zfitint[0] = zfit[0];
    for(int i=1; i<iuse; i++){
        zfitint[i] = zfitint[i-1]+zfit[i];
        errorzfitint[i] = sqrt(zfitint[i]);
    }
    TGraphErrors *graph = new TGraphErrors(iuse,tfit,zfitint,errortfit,errorzfitint);
    graph->SetMarkerSize(0.8);
    graph->SetMarkerStyle(24);    
    return graph;
}
//
TGraphErrors *CaloElectron::GetGraph_Longitudinal_Integral_Fit(){
    double zfunc[44];
    double errorzfunc[44];
    TF1 *fu = new TF1("fu",func,0.,tfit[iuse-1],4);
    for(int i=0; i<4; i++)fu->SetParameter(i,par__l[i]);
    zfunc[0] = fu->Eval(tfit[0]);;
    for(int i=1; i<iuse; i++){
        zfunc[i] = zfunc[i-1] + fu->Eval(tfit[i]);;
        errorzfunc[i] = 0.;
        if(zfunc[i]>0.)errorzfunc[i] = sqrt(zfunc[i]);
    }    
    delete fu;
    TGraphErrors *graph = new TGraphErrors(iuse,tfit,zfunc,errortfit,errorzfunc);
    graph->SetMarkerSize(0.8);
    graph->SetMarkerStyle(24);    
    graph->SetMarkerColor(2);
    graph->SetLineColor(2);
    return graph;
}
//
TF1 *CaloElectron::GetFunc_Longitudinal(){
    TF1 *fu = new TF1("longitudinalprofile",func,0.,tfit[iuse-1],4);
    for(int i=0; i<4; i++)fu->SetParameter(i,par__l[i]);
    fu->SetLineColor(2);
//    fu->Draw();
//    if(err__l!=0)fu->SetLineStyle(12);
    return fu;
}
//
TGraphErrors *CaloElectron::GetGraph_Lateral(int ipp, int ivv){

    cout << "TGraphErrors *CaloElectron::GetGraph_Lateral("<<ipp<<","<<ivv<<") "<<endl;
//     cout <<"ooooooooo%%%%%%%%%%%%%%%%%%%%% "<< CaloElectron_parameters::Get() << endl;
//     cout <<"ooooooooo%%%%%%%%%%%%%%%%%%%%% "<< CaloElectron_parameters::Get()->cl1 << endl;
    CaloLevel1* cl1 = CaloElectron_parameters::Get()->cl1;
    if(!cl1)return NULL;

    float xx[96];
    float qq[96];
    float errxx[96];
    float errqq[96];
    CaloStrip cstrip;
    for(int is=0; is<96; is++){
        cstrip.Set(ivv,ipp,is);
        float xy[2];
        xy[0] = cstrip.GetX() - trkcoordx[ipp][ivv];
        xy[1] = cstrip.GetY() - trkcoordy[ipp][ivv];
        qq[is] = 0.;
        xx[is] = xy[ivv];
        errqq[is] = 0.;
        errxx[is] = 0.5*PITCH;  
    }
    // -----------------------------------
    // set charge collected in each plane
    // -----------------------------------
    int isimu = CaloElectron_parameters::Get()->isimu;
    if(isimu==0)cstrip = CaloStrip(cl1,false);
    if(isimu==1)cstrip = CaloStrip(cl1,true);
    for(int ih=0;ih<cl1->istrip;ih++){
        int iv=-1;
        int ip=-1;
        int is=-1;
        cl1->DecodeEstrip(ih,iv,ip,is);
        cstrip.Set(iv,ip,is);
        float xy[2];
        xy[0]              = cstrip.GetX() - trkcoordx[ip][iv];
        xy[1]              = cstrip.GetY() - trkcoordy[ip][iv];
        if( ip==ipp && iv==ivv ){
            qq[is] = cstrip.GetE();
            xx[is] = xy[iv];
            if(qq[is]>0)errqq[is] = sqrt( qq[is] );
        }
    }    
    TGraphErrors *graph = new TGraphErrors(96,xx,qq,errxx,errqq);
    graph->SetMarkerSize(0.8);
    graph->SetMarkerStyle(20);    
    return graph;
}
//
// TF1 *CaloElectron::GetFunc_Lateral(int ip, int iv, int icorr){


//     int iiv(0);
//     if(iv==0)iiv=0;
//     if(iv==1 && ip%2) iiv=1;
//     if(iv==1 && !ip%2) iiv=2;

//     float iwmax,iw;
//     iwmax = GetShowerMaximum(icorr);
//     iw    = ip+(iv-1);
//     int ibintau = (int)((iw/iwmax)/0.25);
//     if(ibintau>=12) ibintau=11;


//     float file_rt = CaloElectron_parameters::Get()->file_rt[ibintau][iiv];
//     float file_p  = CaloElectron_parameters::Get()->file_p[ibintau][iiv];
//     float file_rc = CaloElectron_parameters::Get()->file_rc[ibintau][iiv];
//     float RMAXsp = CaloElectron_parameters::Get()->RMAXsp;

//     int nparam=5;
//     TF1 *ff = new TF1("fradpro",fradpro,-RMAX,RMAX,nparam);
//     Double_t param[5];
//     param[0] = file_rt;
//     param[1] = file_p;
//     param[2] = file_rc;
//     param[3] = qplanefit[ip][iv];//qplane0[ip][iv];
//     param[4] = 0.;
//     for(int i=0;i<nparam;i++){
//      ff->SetParameter(i,param[i]);
// //   cout<<" "<<i<<" -> "<<param[i]<<endl;
//     }
//     return ff;

// }
/*
 * Return a graph with the integral of the 2-dimentional profile (normalized) over 
 * each strip (taking into account the gaps).
 */
TGraphErrors *CaloElectron::GetFunc_Lateral(int ip, int iv){

    cout << "TGraphErrors *CaloElectron::GetFunc_Lateral("<<ip<<","<<iv<<",) "<<endl;

    float xx[96];
    float qq[96];
    float errxx[96];
    float errqq[96];
    GetProfile(ip,iv,xx,qq,errxx,errqq);


    TGraphErrors *graph = new TGraphErrors(96,xx,qq,errxx,errqq);
//    graph->SetMarkerSize(0.8);
//    graph->SetMarkerStyle(1);    
    graph->SetLineColor(2);    
    return graph;

}
/**
 * evaluate the expected average profile, by integrating the expected charge distribution
 * over each strip
 *
 */
void CaloElectron::GetProfile(int ip, int iv, float xx[], float qq[],float errxx[], float errqq[] ){

//    bool debug = CaloElectron_parameters::Get()->debug;
//     float xx[96];
//     float qq[96];
//     float errxx[96];
//     float errqq[96];
    CaloLevel1* cl1 = CaloElectron_parameters::Get()->cl1;
    if(!cl1)cout << " ERROR -- CaloLevel1* cl1 == NULL"<<endl;
    if(!cl1)return;

    // x y-even or y-odd
    int iiv(0);
    if(iv==0)          iiv = 0;
    if(iv==1 && ip%2)  iiv = 1;
    if(iv==1 && !ip%2) iiv = 2;

    // shower maximum ==> tau-bin
//     float iwmax,iw;
//     iwmax = GetShowerMaximum(icorr);
// //    iw    = ip+(iv-1); //!!!!!!!!!!!!! BUG!?!?!
//     iw    = ip-(iv-1);
//     int ibintau = (int)((iw/iwmax)/0.25);
//     if(ibintau>=12) ibintau=11;
//     cout << " maximum (W)"<<iwmax<<endl;
//     cout << " ip      (W)"<<iw<<endl;
//     cout << " tau "<<iw/iwmax<<" tau-bin "<<ibintau<<endl;
    int ipmin    = CaloElectron_parameters::Get()->ipmin;
    float maxl   = GetShowerMaximum(2);
    float nlayer = (float)GetNWLayers(ip,iv);
    double tau   = (double)(nlayer-ipmin)/(double)maxl;
    int ibintau  = (int)(tau/0.25);
    if(ibintau>=12) ibintau=11;
//     if(debug){
//      cout << " maximum (W)"<<maxl<<endl;
//      cout << " ip      (W)"<<nlayer<<endl;
//      cout << " tau "<<tau<<" tau-bin "<<ibintau<<endl;
//     }
    // retrieve lateral-profile parameters
    float file_rt = CaloElectron_parameters::Get()->file_rt[ibintau][iiv];
    float file_p  = CaloElectron_parameters::Get()->file_p[ibintau][iiv];
    float file_rc = CaloElectron_parameters::Get()->file_rc[ibintau][iiv];
    float RMAXsp  = CaloElectron_parameters::Get()->RMAXsp;
//     if(debug){
//      cout << " RT "<<file_rt<<endl;
//      cout << " RC "<<file_rc<<endl;
//      cout << " p  "<<file_p<<endl;
//     }
    // 2-dimentional function
    TF2 ff = TF2("frad",frad,-25.,25.,-25.,25.,3);
    ff.SetParameter(0,file_rt);
    ff.SetParameter(1,file_p);
    ff.SetParameter(2,file_rc);

    // 1-dimentional function
    TF1 fx = TF1("fradx",fradx,-25.,25.,4);
    fx.SetParameter(0,file_rt);
    fx.SetParameter(1,file_p);
    fx.SetParameter(2,file_rc);
    
    // sensitive area
    //    double x1,x2,x3,x4,x5,x6,ps[3];
    double x1,x2,x3,x4,x5,ps[3];
    x1 = DEAD;
    x2 = x1 + PITCH*32;
    x3 = x2 + DEAD + GLUE + DEAD;
    x4 = x3 + PITCH*32;
    x5 = x4 + DEAD + GLUE + DEAD;
    //    x6 = x5 + PITCH*32;
    ps[0] = x1;
    ps[1] = x3;
    ps[2] = x5;

    // retrieve the coordinate of the bottom-left corner of the plane
    float corner[] = {0.,0.};
    GetCornerCoord(ip,iv,corner[0],corner[1]);
    

    // integral over each strip
    double integgap;
    double epsil = 0.0001;
    CaloStrip cstrip;
    
    //
    for(int is=0; is<96; is++){ //loop over strips
        cstrip.Set(iv,ip,is);
        float xy[2];
        float trkcoord[2];
        trkcoord[0] =  trkcoordx[ip][iv]; //track coordinate
        trkcoord[1] =  trkcoordy[ip][iv];
        xy[0] = cstrip.GetX() - trkcoordx[ip][iv];//strip coordinate (relative to track)
        xy[1] = cstrip.GetY() - trkcoordy[ip][iv];
        //
        qq[is] = 0.;
        xx[is] = xy[iv];
        errqq[is] = 0.;
        errxx[is] = 0.5*PITCH;
        //
        integgap = 0.;
        for (Int_t jj=0;jj<3;jj++){   //loop over the 3 sensors 
            //
            double ymi = ps[jj] + (double)corner[iv] - (double)trkcoord[iv];
            double yma = ymi + PITCH*32;
            double xmi = (double)xx[is]-(double)errxx[is];//strip size
            double xma = (double)xx[is]+(double)errxx[is];
            //
            int dointeg=1;
            if(iv==0){
                if(xmi<-RMAXsp && xma>-RMAXsp) xmi = -RMAXsp;
                if(xma>RMAXsp && xmi<RMAXsp)   xma =  RMAXsp;
                if(xmi>RMAXsp || xma<-RMAXsp)  dointeg = 0;
            }
            if(iv==1 || iv==2){
                if(ymi<-RMAXsp && yma>-RMAXsp) ymi = -RMAXsp;
                if(yma>RMAXsp && ymi<RMAXsp)   yma =  RMAXsp;
                if(ymi>RMAXsp || yma<-RMAXsp)  dointeg = 0;
            }
            if(dointeg==0) continue;

            //  
            bool DOFAST = false;
            if(fabs(xx[is])>2.)DOFAST=true;
        
            double integral = 0.;
            
            if(DOFAST){
            //
            // integrale veloce
            //

              //                double fmi,fma,k,c;
                double fmi,fma,k;
                fx.SetParameter(3,xmi);
                fmi = fx.Integral(ymi,yma);
                fx.SetParameter(3,xma);
                fma = fx.Integral(ymi,yma);
                k = 0.;
                if(fmi>0.&&fma>0.&&(xma-xmi)>0.)k=log(fma/fmi)/(xma-xmi);
                //c = log(fmi)-k*xmi;
                if(fabs(k)>1.e-6)integral += (fma-fmi)/k;
                    
            }else{
            //
            // integrale lento ma piu` preciso
            //
                integral = ff.Integral(xmi,xma,ymi,yma,epsil);

            }
            
            integgap += integral;           

//              cout << setw(2)<<jj<<setw(15)<<xmi<<setw(15)<<ymi<< setw(15)<<integral << setw(15)<<integgaptmp<<endl;
        }

        qq[is]    = integgap * qplanefit[ip][iv];
//      if( qq[is]>0 )errqq[is] = sqrt(qq[is]);
    }
}
/**
 *
 */
float CaloElectron::ProfileTest(){
    
    int ipmin    = CaloElectron_parameters::Get()->ipmin;
    int ipmax    = CaloElectron_parameters::Get()->ipmax;
    float RMAX   = CaloElectron_parameters::Get()->RMAX;
//    float TAUMAX = CaloElectron_parameters::Get()->TAUMAX;

    TH2F* h_qfit[2];
    TH2F* h_qtot[2];
    //
    // retrieve histrograms of shower lateral profile
    //
    for(int iv=0; iv<2; iv++){
        h_qfit[iv] = CaloElectron_parameters::Get()->h_qfit[iv]; //average
        if(!h_qfit[iv])return 0.;
        h_qfit[iv]->Reset();                                     //-->reset

        h_qtot[iv] = CaloElectron_parameters::Get()->h_qtot[iv]; //data
        if(!h_qtot[iv])return 0.;
    }

//    float maxl = GetShowerMaximum(2);
    //
    double chi2tmp = 0.;
    int    ndof = 0;
//     cout << "<><><><><><><><><><><><><><>"<<endl;
//     cout << setw(3)  << "IP" << setw(3) << "IV" ;
//     cout << setw(10) << "TAU";
//     cout << setw(10) << "RMS";
//     cout << setw(10) << "SKEW";
//     cout << setw(10) << "KURT";
//     cout << endl;
    //
    for(int ip=0; ip<22; ip++){           // loop over planes
        for(int iv=1; iv>=0; iv--){        // loop over views
            bool maskpl = CaloElectron_parameters::Get()->maskpl[ip][iv];
            float corner[] = {0.,0.};
            GetCornerCoord(ip,iv,corner[0],corner[1]);
            if( 
                ip >= ipmin && 
                ip <  ipmax  && 
                !maskpl   &&
                true){
                //
                float xx[96];
                float qq[96];
                float errxx[96];
                float errqq[96];
                //
                float trkcoord[2];
                trkcoord[0] = trkcoordx[ip][iv]; 
                trkcoord[1] = trkcoordy[ip][iv]; 
                //
                GetProfile(ip,iv,xx,qq,errxx,errqq); //<<< average profile
//              float nlayer = (float)GetNWLayers(ip,iv);
//              double tau   = (double)(nlayer-ipmin)/(double)maxl;
                for(int is=0; is<96; is++) { 
                    if( 
                        fabs(xx[is]) < RMAX &&
//                      tau<TAUMAX &&
                        true){
//                      h_qfit[iv]->Fill( (float)(is+iv*96), (float)ip, qq[is]); 
                        h_qfit[iv]->Fill( xx[is]+trkcoord[iv]-corner[iv], (float)ip, qq[is]); 
                        float qqq = h_qtot[iv]->GetBinContent((is+iv*96+1),ip+1);
//                  cout << setw(15)<< qqq<<setw(15)<<qq[is]<<setw(15)<<pow(qq[is]-qqq,2.)<< endl;
                        if(qq[is]>0) chi2tmp += pow((double)(qq[is]-qqq),2.)/sqrt((double)qq[is]);
                        ndof++;
                    }
                }
//              cout << "<><><><><><><><><><><><><><>"<<endl;
//              cout << ip << " - "<<iv << endl;
//              TH1D* hfitpro = h_qfit[iv]->ProjectionX("hfitpx",ip+1,ip+2);
                TH1D* htotpro = h_qtot[iv]->ProjectionX("htotpx",ip+1,ip+2);
//              cout <<"Mean "<< setw(10)<<htotpro->GetMean()     << setw(10)<<hfitpro->GetMean() << endl;
//              cout <<"RMS  "<< setw(10)<<htotpro->GetRMS()      << setw(10)<<hfitpro->GetRMS() << endl;
//              cout <<"Skew "<< setw(10)<<htotpro->GetSkewness() << setw(10)<<hfitpro->GetSkewness() << endl;
//              cout <<"Kurt "<< setw(10)<<htotpro->GetKurtosis() << setw(10)<<hfitpro->GetKurtosis() << endl;
                if(GetShowerMaximumX0(2)>0.)
                    TAU[ip][iv]  = tplane[ip][iv]/GetShowerMaximumX0(2);
                if(htotpro->GetEntries()>4){
                    RMS[ip][iv]  = htotpro->GetRMS();
                    SKEW[ip][iv] = htotpro->GetSkewness();
                    KURT[ip][iv] = htotpro->GetKurtosis();
                }
//              cout << setw(3)  << ip << setw(3) << iv ;
//              cout << setw(10) << TAU[ip][iv];
//              cout << setw(10) << RMS[ip][iv];
//              cout << setw(10) << SKEW[ip][iv];
//              cout << setw(10) << KURT[ip][iv];
//              cout << endl;
//              hfitpro->Delete();
                htotpro->Delete();
//              cout << "<><><><><><><><><><><><><><>"<<endl;
            }
        }
    }    
//     cout << "<><><><><><><><><><><><><><>"<<endl;

//     double kolp = h_qfit->KolmogorovTest(h_qtot,"ND");
//     double kold = h_qfit->KolmogorovTest(h_qtot,"NDM");
//     h_qfit->Sumw2();
//     h_qtot->Sumw2();
//     double kolp = h_qfit->KolmogorovTest(h_qtot,"D");
//     double kold = h_qfit->KolmogorovTest(h_qtot,"DM");
//     double chi2tmp = h_qfit->Chi2Test(h_qtot,"WWP");
//     cout << "KKKKK "<<kolp << endl;
//     cout << "KKKKK "<<kold << endl;
//    cout << "KKKKK "<<chi2 << " "<<ndof<<" "<<chi2/(double)(ndof)<< endl;
    if( ndof>3 )chi2 = (float)(chi2tmp/(double)(ndof-3));
    else        chi2 = 0.;

    return chi2;
    
}
/**
 *
 */
float CaloElectron::GetMaxResidual(){
    float max = 0;
    for(int iv=0; iv<2; iv++){
        for(int ip=0; ip<22; ip++){
            if(
//              qplane1[ip][iv] > 0. &&
//              fabs(qplane1[ip][iv]-qplanefit[ip][iv])/sqrt(qplane1[ip][iv]) > fabs(max) && 
                fabs(qplane1[ip][iv]-qplanefit[ip][iv]) > fabs(max) && 
                true){
                max = qplane1[ip][iv]-qplanefit[ip][iv];
            }
        }
    }    
    return max;
}
int CaloElectron::GetNRun(){
    bool POS = true;
    bool NEG = true;
    int nrun = 0;
    for(int ip=0; ip<22; ip++){
        for(int iv=1; iv>=0; iv--){
            if( qplane1[ip][iv] > qplanefit[ip][iv] && NEG ){
                nrun++;
                POS = true; NEG = !POS; 
            }
            if( qplane1[ip][iv] <= qplanefit[ip][iv] && POS ){
                nrun++;
                NEG = true; POS = !NEG; 
            }
        }    
    }
    return nrun;
}
float CaloElectron::GetMaxRun(){
    bool POS = true;
    bool NEG = true;
    int nrun = 0;
    float sum[44]; for(int ip=0; ip<44; ip++)sum[ip]=0.;
    //
    for(int ip=0; ip<22; ip++){
        for(int iv=1; iv>=0; iv--){
            if( qplane1[ip][iv] > qplanefit[ip][iv] && NEG ){
                nrun++;
                POS = true; NEG = !POS;                 
            }
            if( qplane1[ip][iv] <= qplanefit[ip][iv] && POS ){
                nrun++;
                NEG = true; POS = !NEG; 
            }
            sum[nrun-1] += qplane1[ip][iv]-qplanefit[ip][iv];
        }    
    }
    //
    float max=0;
    for(int ip=0; ip<44; ip++)
        if( fabs(sum[ip]) > fabs(max) )max=sum[ip];

    return max;
}
/**
 * \brief Method to retrieve the coordinates of the plane bottom-left corner
 * (l'angolo è dell'inizio del rivelatore compreso il dead)
 */
void CaloElectron::GetCornerCoord(int ip, int iv, float& xorig, float& yorig){

// coordinate meccaniche dei piani in cm

    CaloStrip st = CaloStrip();
    int isimu = CaloElectron_parameters::Get()->isimu;
    if(isimu==0)st.UseStandardAlig();
    if(isimu==1)st.UseMechanicalAlig();
        
    yorig = -0.1*st.GetYalig()+0.1;
    xorig = -0.1*st.GetXalig();
    
    if(ip%2!=1){
        //ip 0,2,4,6,...,20
        if(iv==0){
            yorig = yorig;
            xorig = xorig+0.05;
        }else{
            yorig = yorig-0.05;
            xorig = xorig+0.1;
        }
    }else{      
        //ip 1,3,5,..,21
        if(iv==0){
            yorig = yorig-0.2;
            xorig = xorig-0.05;
        }else{
            yorig = yorig-0.15;
            xorig = xorig-0.1;
        }
    }
    
    
}
/**
 * \brief Print the event
 */
void CaloElectron::Print(){
    
    int ipmin = CaloElectron_parameters::Get()->ipmin;
    int ipmax = CaloElectron_parameters::Get()->ipmax;
//    cout<<"Print idiev "<<idiev<<endl;
    cout << " n. W-layers -        qx0       qy0 ------->      qx1       qy1 ------->      qx(fit)   qy(fit)"<<endl;
    float qx0,qy0,qx1,qy1,qx2,qy2;
    qx0=0.;
    qy0=0.;
    qx1=0.;
    qy1=0.;
    qx2=0.;
    qy2=0.;
    for(int il=ipmin; il<=ipmax; il++){
        if(il>ipmin){
            qx0 = qplane0[il-1][0];
            qy0 = qplane0[il][1];
            qx1 = qplane1[il-1][0];
            qy1 = qplane1[il][1];
            qx2 = qplanefit[il-1][0];
            qy2 = qplanefit[il][1];
        }
        if(il==ipmin){
//          qx0 = 0;
            qy0 = qplane0[ipmin][1];
//          qx1 = 0;
            qy1 = qplane1[ipmin][1];
//          qx2 = 0;
            qy2 = qplanefit[ipmin][1];
        }
        if(il==ipmax){
            qx0 = qplane0[ipmax-1][0];
//          qy0 = 0.;
            qx1 = qplane1[ipmax-1][0];
//          qy1 = 0.;
            qx2 = qplanefit[ipmax-1][0];
//          qy2 = 0.;
        }
        cout << "          "<<setw(2)<<il<<" - "<<setw(10)<<qx0<<setw(10)<<qy0<<"       "<<setw(10)<<qx1<<setw(10)<<qy1<<"       "<<setw(10)<<qx2<<setw(10)<<qy2<<endl;
    }
    cout<<" "<<endl;
    float TAUMAX = CaloElectron_parameters::Get()->TAUMAX;
    float RMAX   = CaloElectron_parameters::Get()->RMAX;
    cout << "Total collectced charge (mip) "<<qtot<<endl;
    cout << "Charge outside cylinder ( RMAX="<<RMAX<<" ,TAUMAX="<<TAUMAX<<" )";
    cout << qtot - GetQ(0);
    if(qtot>0.)cout <<" ( "<< 100.*(qtot-GetQ(0))/qtot<<" % )"<<endl;    
    for(int icorr=0;icorr<4;icorr++){
        cout << " ----------------------------------------------   "<<icorr;
        if(icorr==0)cout <<"(no corr.)"<<endl;
        if(icorr==1)cout <<"(lateral corr.)"<<endl;
        if(icorr==2)cout <<"(lateral+longitudinal corr.)"<<endl;
        if(icorr==3)cout <<"(fit)"<<endl;
        cout << " - Collected charge (mip)       "<<GetQ(icorr)<<endl;
        cout << " - Shower maximum (W-layer, X0) "<<GetShowerMaximum(icorr)<<" "<<GetShowerMaximumX0(icorr)<<endl;
        cout << " - Energy (GeV)                 "<<GetEnergy(icorr)<<endl;     
    }
    cout << "Longitudinal parameters:"<<endl;
    cout << "- E             :"<<par__l[0]<<endl;
    cout << "- alpha         :"<<par__l[1]<<endl;
    cout << "- Tm            :"<<par__l[2]<<endl;
    cout << "- T0            :"<<par__l[3]<<endl;
    cout << "- chi2          :"<<chi2__l<<endl;
    cout << "- minuit status :"<<err__l<<endl;
    cout << endl ;
    cout << "- chi2 (global) :"<<chi2<<endl;
    cout << endl ;

//    cout<<" rigidity "<<rig<<" qtotlevel2 "<<qtotl2<<endl;
    
}
//______________________________________________________________________________
Double_t paraloss(Double_t *tin,Double_t *par)
{
    
    Double_t t,p0,p1,p2,p3,t0,value;
    t=*tin;
    p0=par[0];
    p1=par[1];
    p2=par[2];
    p3=par[3];
    t0=1.7;//t0=1.5;
    Double_t paralossdit0=1.00659;
    if(t<t0) return paralossdit0;
    value=p0+p1*t+p2*pow(t,2)+p3*pow(t,3);
    return value;
}

/**
 * \brief Collect corrected charge to energy
 */
float CaloElectron::GetEnergy(int icorr){
    
    float energy = 0;
    float collectedq = 0;
    float p0 = 0;
    float p1 = 0;
    
    if(icorr == 0){
        // ---------------------------------------
        // before any correction
        // ---------------------------------------
        
        //the same as energy 1.. not to be used 
        p0=24.14;
        p1=253.3;
        
    }else if(icorr == 1 || icorr == 3){
        // ---------------------------------------
        // after lateral correction 4 giugno 2007
        // ---------------------------------------
        
        //linear (only with fit of points 5-15gev where longitudinal containment is satisfactory)
        //we expect a sistematic error for greater energies that has to be corrected by longitudinal fit
        
        //10-9-07
        p0=24.14;
        p1=253.3;
        
    }else if(icorr == 2){
        // ---------------------------------------
        // after longitudinal fit
        // ---------------------------------------
        
        //10-9-07
        p0=-13.86;
        p1=264.4;
        
        
    }else{

        return energy;

    }

    collectedq=GetQ(icorr);
//    if(onlyy)collectedq=collectedq*2;
    
    energy = (collectedq-p0)/p1;

    if(GetQ(icorr)<=0) energy=0.;
    
    if(icorr == 3){
        TF1 *f4=new TF1("f4",paraloss,0.,3.,4);//f1->FixParameter(2,1.3);
        f4->SetParameter(0,-0.8013231);
        f4->SetParameter(1,2.276176);
        f4->SetParameter(2,-0.8295547);
        f4->SetParameter(3,0.06835571);
        float lnmax=TMath::Log(par__l[2]);
        float lnmax1=TMath::Log(GetShowerMaximumX0(1));
        if(lnmax>2.4 || GetIerflg()>0){
            cout<<"max "<<GetShowerMaximumX0(1)<<" "<<GetShowerMaximumX0(2)<<" lnmax "<<lnmax<<" ier "<<GetIerflg()<<endl;
            lnmax=lnmax1;
        }
        float corrparaloss=f4->Eval(lnmax);
        //cout<<"icorr 3...lnmax "<<lnmax<<" corrparaloss "<<corrparaloss<<endl;
        energy=energy/corrparaloss;
    }
    
    return energy;
    
};
    
/*
 defined cuts to have gaussian error from calorimeter reconstruction for alignement purpose
*/
bool CaloElectron::IsInsideFiducial(int level){
 
 
//     float ma=GetShowerMaximum(1);
//     if(ma>=20 || ma<=4)return false;
//     int ima=(int)ma+ipmin;
//     float strima=0;
//     for(int iv=0;iv<=1;iv++){
//      strima=GetStripW(ima,iv);
//      if(strima<11 || strima>86 || strima==33 ||  strima==32 || strima==65 || strima==64) return false;
//     }
    
//     if(level==1) return true;
    
//     if(outfcn>600) return false;
//     if(minuitierflg>0) return false;
//     float lma1=log(par__l[1]);
//     float lma2=log(par__l[2]);
//     float lma1expected=0.92*lma2;
//     if(lma2<1.3 || lma2>2.3)  return false;
//     if(TMath::Abs(lma1-lma1expected)>0.5)  return false;
    
    return true;
 
}

//////////////////////////////////////////////////////////////////////////

/**
 * Fill calibration histos
 */
void CaloElectron::Calibration_Fill(){

    bool debug      = CaloElectron_parameters::Get()->debug;
    float TAUMAX    = CaloElectron_parameters::Get()->TAUMAX;
    float RMAX      = CaloElectron_parameters::Get()->RMAX;
    int ntau        = CaloElectron_parameters::Get()->par_ntau;
    float *taubin   = CaloElectron_parameters::Get()->par_taubin;
    CaloLevel1* cl1 = CaloElectron_parameters::Get()->cl1;
    //
    if(!cl1)cout << "NB --> missing CaloLevel1, needed for calibration !!"<<endl;
    if(!cl1)return;
    //
    int isimu       = CaloElectron_parameters::Get()->isimu;

    // /////////////////////////////////////
    // some parameters to select good shower
    // /////////////////////////////////////
    float stoll_g = 5;//maximum n.strips from gaps
    float stoll_b = 10;//maximum n.strips from borders
    //    float toll = stoll_b*PITCH;//3.;//(cm) distance of shower axis from border on the last silicon plane
    float qmin = 10;//(mip)minimum energy deposit per plane
    //

    // ---------------
    // lateral profile
    // ---------------
    CaloStrip cstrip;
    if(isimu==0)cstrip = CaloStrip(cl1,false);
    if(isimu==1)cstrip = CaloStrip(cl1,true);
    if(debug)cout << "======  n.strip "<<cl1->istrip<<endl;
    if(debug)cout << "======  s.max (X0) "<<GetShowerMaximumX0(0);
    float showermax = GetShowerMaximumX0(0);

    // -------------------------------
    // containment conditions (global)
    // -------------------------------
    // shower axis must be contained within 3 cm (on the last plane)
    // (otherwise the shower maximum could be not properly evaluated)
    //    bool CONTAINED = false;
    //
    // last silicon plane
    //
    int ip = 21;
    int iv = 0;//x
    //    cstrip.Set(iv,ip,0);
    //    float xcorner_l = cstrip.GetX();
    //    float ycorner_l = cstrip.GetY();
    cstrip.Set(iv,ip,95);
    //    float xcorner_r = cstrip.GetX();
    //    float ycorner_r = cstrip.GetY();         
    //    
    //    if(
    //  trkcoordx[ip][iv] > xcorner_l + toll &&
    //  trkcoordx[ip][iv] < xcorner_r - toll &&
    //  trkcoordy[ip][iv] > ycorner_l + toll &&
    //  trkcoordy[ip][iv] < ycorner_r - toll &&
    //  true)CONTAINED = true;
    // 
//    if(!CONTAINED)return; // discard event if shower is not contained
    
    // -----------------------------------------------
    // first loop over planes to determine the tau-bin
    // and verify some general conditions:
    // - shower axis in each plane far from gaps
    // - minimal energy deposit
    // ----------------------------------------------- 
    int ibintau[22][2];
    //
    for(int ip=0; ip<22; ip++){
        for(int iv=0; iv<2; iv++){
            //
            ibintau[ip][iv]=-1;
            // ------------
            // shower depth
            // ------------         
            float tau = tplane[ip][iv]/showermax;
            if(debug)cout <<" --> plane "<<ip<<" view "<<iv<<" tau "<<tau<<endl;
            //
            if(tau>TAUMAX)continue; // discard plane if tau exceed TAUMAX
            //
            // ---------------------------------------
            // containment conditions (plane by plane)
            // ---------------------------------------
            // shower axis should be far from gaps
            //      
            if(qplane0[ip][iv]<qmin)continue;//discard plane if energy deposit is low
            //
            bool CONTAINED = false;
            //
            if(
                ( trkstrip[ip][iv] >= stoll_b      && trkstrip[ip][iv] <  32 - stoll_g )||
                ( trkstrip[ip][iv] >= 32 + stoll_g && trkstrip[ip][iv] <  64 - stoll_g )||
                ( trkstrip[ip][iv] >= 64 + stoll_g && trkstrip[ip][iv] <  96 - stoll_b )||
                false)CONTAINED = true;
            // 
            if(!CONTAINED)continue; // discard plane if shower traverses gaps
            //
            for(int itau = 0; itau<ntau; itau++){
                if(taubin[0]>tau)break;
                if(taubin[itau+1]>=tau){
                    ibintau[ip][iv]=itau;
                    break;
                }
            }
            if(debug)cout << " --> tau-bin "<<ibintau<<endl;
            //
        }
    }
    // ----------------------------------------------------
    // then loop over hit strips to fill profile histograms
    // ----------------------------------------------------
//     for(int ih=0;ih<cl1->istrip;ih++){//loop over hit strips
//      int iv=-1;
//      int ip=-1;
//      int is=-1;
//      cl1->DecodeEstrip(ih,iv,ip,is);
//
//  NBNB --> credo di dover riempire con TUTTE le strisce
//
    for(int ip=0; ip<22; ip++){ 
        for(int iv=0; iv<2; iv++){

            int ivv=iv;    //x (0) and yeven (1)
            if(ip%2)ivv=2; //yodd (2)
            //
            if(ibintau[ip][iv]==-1)continue; //outside tau range or plane discarded
            //
            for(int is=0; is<96; is++){

                cl1->GetEstrip(iv,ip,is);
                cstrip.Set(iv,ip,is);
                //
                float xy[2];
                xy[0]              = cstrip.GetX() - trkcoordx[ip][iv];
                xy[1]              = cstrip.GetY() - trkcoordy[ip][iv];
                //
                if( TMath::Abs(xy[iv]) >= RMAX )continue;// discard strip if  exceed TAUMAX
                //
                // FRACTION OF CHARGE COLLECTED BY THE STRIP
                //
                float pitch = 1.;//PITCH;
                float F = cstrip.GetE();
                F = F / qplane0[ip][iv]; // normalized to total charge
                F = F / pitch;   // per unit length
                //
                // weigths
                //
                double calo_noise = 0.2; //MIP
//              if(F==0.)calo_noise = 0.7; //MIP (threshold)            
                double DF2 = 0.;
                DF2 += pow(calo_noise/cstrip.GetE(),2.)*pow((double)F,2.);
//              DF2 += pow((double)F,2.)/(double)qplane0[ip][iv];//sbagliato
                double W = 1./DF2;
                // -----------
                // fill histos
                // -----------
                TH2F* h_qtot_tau = CaloElectron_parameters::Get()->h_qtot_tau[ibintau[ip][iv]][ivv];
                if( !h_qtot_tau )
                    CaloElectron_parameters::Get()->Calibration_SetTools();
                //
                if( !h_qtot_tau )return;;
                h_qtot_tau->Fill(xy[iv],F);
                // -------------------------
                // evaluate weighted average
                // -------------------------
                int ibin;               
                for(ibin=-1; ibin<CALIBNBIN; ibin++){               
                    if( -1*CALIBRANGE + (ibin+1)*2*CALIBRANGE/CALIBNBIN >= xy[iv] )break;
                }
                if(ibin==-1)continue;//strip out of range
                //
                // problema: quando il segnale e` 0. in realta` e` <0.7 
                // quindi si introduce un bias... che spero pero` sia piccolo
                //
                //
                CaloElectron_parameters::Get()->summW[ibintau[ip][iv]][ivv][ibin] 
                    += W;
                CaloElectron_parameters::Get()->summWF[ibintau[ip][iv]][ivv][ibin] 
                    += W*(double)F;

            }               
        }
    }


    return;
};

ClassImp(CaloElectron);
ClassImp(CaloElectron_parameters);

