#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,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];
    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;
		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);