#include <CaloBragg.h>


ClassImp(CaloBragg);
//--------------------------------------
/**
 * Default constructor 
 */
CaloBragg::CaloBragg(){
  Clear();
};

CaloBragg::CaloBragg(PamLevel2 *l2p){  
  //
  Clear();
  LoadParam();
  //
  L2 = l2p;
  //
  if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
  //
  OBT = 0;
  PKT = 0;
  atime = 0;
  //
  debug = false;
  usetrack = false;
  //
};

void CaloBragg::Clear(){
  //
  tr = 0;
  sntr = 0;
  qtchi2 = 0.;
  qtz = 0.;
  qtetot = 0.;
  qtpskip = 0.;
  lpchi2 = 0.;
  lpz = 0.;
  lpetot = 0.;
  lppskip = 0.;
  memset(calorimetro,0,44*2*sizeof(Float_t));
  memset(spessore,0,3*sizeof(Float_t));
  memset(estremi,0,2*2*sizeof(Float_t));
  Integrale=0.;

  //
};

void CaloBragg::Print(){
  //

  if(!debug)  Process();
  //
  printf("========================================================================\n");
  printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack);
  printf(" chi 2 from truncated mean: %f \n", qtchi2);
  printf(" Z from truncated mean %f: \n", qtz);
  printf(" energy from truncated mean %f: \n", qtetot);
  printf(" plane not used for truncated mean %f: \n", qtpskip);
  printf(" chi 2 from loop %f:  \n", lpchi2);
  printf(" Z from loop %f: \n", lpz);
  printf(" energy from loop %f: \n", lpetot);
  printf(" plane not used for loop %f: \n", lppskip);
  printf("========================================================================\n");
  //
};

void CaloBragg::Delete(){
  Clear();
  //delete this;
};


void CaloBragg::Process(){
  Process(-1);
};

void CaloBragg::Process(Int_t ntr){
  //  
  if ( !L2 ){
    printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
    printf(" ERROR: CaloBragg variables not filled \n");
    return;
  };
  //
  Bool_t newentry = false;
  //
  if ( L2->IsORB() ){
    if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || ntr != sntr ){
      newentry = true;
      OBT = L2->GetOrbitalInfo()->OBT;
      PKT = L2->GetOrbitalInfo()->pkt_num;
      atime = L2->GetOrbitalInfo()->absTime;
      sntr = ntr;
    };
  } else {
    newentry = true;
  };
  //
  if ( !newentry ) return;
  //
  tr = ntr;
  //
  if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
  //
  Clear();

  //
  //
  // Always calculate stdedx1
  //
  Int_t view = 0;
  Int_t plane = 0;
  Int_t strip = 0;
  Float_t mip = 0.;
  Float_t epiano[22][2];
  memset(epiano,0,22*2*sizeof(Float_t));
  for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
    //
    mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
    epiano[plane][view]+=mip;
    //
    //
  };
  //
  //
  PamTrack *ptrack = 0;
  CaloTrkVar *track = 0;
  //
  if ( usetrack ){
    if ( ntr >= 0 ){
      ptrack = L2->GetTrack(ntr);
      if ( ptrack ) track = ptrack->GetCaloTrack();
    } else {
      track = L2->GetCaloStoredTrack(ntr);
    };
    //
    if ( !track && ntr >= 0 ){
      printf(" ERROR: cannot find any track!\n");
      printf(" ERROR: CaloBragg variables not completely filled \n");
      return;   
    };
  } else {
    if ( ntr >= 0 ){
      if ( debug ) printf(" ERROR: you asked not to use a track but you are looking for track number %i !\n",ntr);
      if ( debug ) printf(" ERROR: CaloBragg variables not completely filled \n");
      return;       
    };
  };
  //
  if(L2->GetCaloLevel2()->npcfit[0]==0 && L2->GetCaloLevel2()->npcfit[1]==0 && L2->GetCaloLevel2()->npcfit[2]==0 && L2->GetCaloLevel2()->npcfit[3]==0) return;// controllo sulla traccia nel calorimetro

  //
  for(Int_t p=0; p<22; p++){
    for(Int_t v=0; v<2; v++){
      cout<< L2->GetCaloLevel2()->cibar[p][v];    
      /*per usare traccia non del calo camboare cibar*/
      calorimetro[(2*p)+1-v][0] = L2->GetCaloLevel2()->cibar[p][v];//strip attraversata
      calorimetro[(2*p)+1-v][1] = (epiano[p][v]); //energia del piano //(epiano[p][v])/0.89
    };
  };


  /*per ogni evento clacolo la conversione mip e w attraversato in equivalente Si*/

  conversione(); // out: 1) g/cm2 Si , 2) spessoreW equivalente in Si, 3)Mip corretta per inclinazione

  cout<<"spessore"<<spessore[0]<<" gcm2si; "<<spessore[1]<<"W in Si; "<<spessore[2]<<"mip "<<endl;

  /*trova primo e ultimo piano attraversati*/

  Int_t p = 0;//contatore piani
  //primo parte da 0 e va in giù
  while( estremi[0][1] == 0  &&  p<(2*NPLA) ){
    cout<<"piano "<<p <<"; strip "<<calorimetro[p][0]<<"; en "<<calorimetro[p][1]<<endl;
    //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){
    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){ //0.7 soglia minima
      estremi[0][0]=p;
      estremi[0][1]=calorimetro[p][1] *MIP; //energia in MeV
    };
    p++;
  };
 //ultimo parte da 44 e sale
  p=43;
  while( (estremi[1][1] == 0)  &&  (p>(int)estremi[0][0]) ){
  cout<<"piano "<<p <<"; strip "<<calorimetro[p][0]<<"; en "<<calorimetro[p][1]<<endl;
  //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){
  if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){
      estremi[1][0]=p;
      estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV
    };
    p = p-1;
  };
  //
  cout<<"estremi: in "<<estremi[0][0]<<"piano, "<<estremi[0][1]<<"energia; out "<<estremi[1][0]<<"piano, "<<estremi[1][1]<<"energia."<<endl; 


  
  /*integrale: energia rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/

  for(Int_t pl=0; pl<(2*NPLA); pl++){

    //calcolo intergale in unita di spessori di silicio
    Integrale += calorimetro[pl][1] * MIP;//piano di silicio
    // cout<<Integrale<<" piani "<<pl
    cout<<"nel calo ho strip "<<calorimetro[pl][1]*MIP<<endl;
    // se non e'il 1o dopo l'Y (tutti i pari) c'e' il W          
    if(pl%2!=0){                                                             //equival W in Si
      Integrale+= 0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1]);
      // cout<<Integrale<<" W "<<pl<<endl;
      cout<<" W "<<0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1])<<endl;
    };
  };

  //  cout<<"integrale ok"<<endl;
  
  
  /*z se particella fosse al minimo*/  //energia1piano/mip corretta
  //Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));
  

  
  /*z ed energia con media troncata*/
  //  Float_t  bestchi2mean[4] = {0.,0.,0.,0.};//chi2,z,Etot,Pskip
  mediatroncata();  
  
  /*z ed energia con loop*/
  //  Float_t  bestchi2loop[4] = {0.,0.,0.,0.};//chi2,z,Etot,Pskip//chi2,z,Etot,Pskip
  Zdaloop();
  

  /*energia rilasciata da z migliore*/
 //  Float_t dEpianimean[2*NPLA];
  //  Int_t zet=(int)bestchi2mean[1];
  //Enetrack(&zet, &bestchi2mean[2], &estremi[0][0], dEpianimean);//calcola rilascio energetico sui piani

  // Float_t dEpianiloop[2*NPLA];
  //zet=(int)bestchi2loop[1];
  //Enetrack(&zet, &bestchi2loop[2], &estremi[0][0], dEpianiloop);//calcola rilascio energetico sui piani

  
  if ( debug ) this->Print();
  if ( debug ) printf(" esci \n");
  //
};


void CaloBragg::Draw(){

  Process();

 Float_t dEpianimean[44];
 Float_t dEpianiloop[44];
 Float_t Depth[44];
 Int_t tz=(Int_t)qtz;
 Int_t tz1=(Int_t)lpz;
 Enetrack(&tz, &qtetot, &estremi[0][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata
 Enetrack(&tz1, &lpetot, &estremi[0][0], dEpianiloop);//calcola rilascio energetico sui piani da loop

 Float_t sp= spessore[0]*spessore[1];
 for(Int_t i=0;i<44;i++)Depth[i]=i*sp;
  //
  gStyle->SetLabelSize(0.04);
  gStyle->SetNdivisions(510,"XY");
  //gStyle->SetLogy();
  //

	TString hid = Form("cCaloBragg");	
        TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));

	if ( tc ){
//	 tc->Clear();
	} else {
	 tc = new TCanvas(hid,hid);
	 tc->Divide(1,2);
	};
	//
	TString thid = Form("hCaloBragg");	
        TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
	if ( th ) th->Delete();
//	 th->Clear();
//	 th->Reset();
//	} else {
	th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.);
	th->SetMarkerStyle(20);
//	};
	//
	TString thid2 = Form("hCaloBragg2");	
        TH2F *th2  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid2));
	if ( th2 ) th2->Delete();
	th2 = new TH2F(thid2,thid2,300,-0.5,300.,1000,0.,150.);
	th2->SetMarkerStyle(20);
	th2->SetMarkerColor(kRed);
	//
	TString thid3 = Form("hCaloBragg3");	
        TH2F *th3  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid3));
	if ( th3 ) th3->Delete();
	th3 = new TH2F(thid3,thid3,300,-0.5,300.,1000,0.,150.);
	th3->SetMarkerStyle(20);
	th3->SetMarkerColor(kBlue);


	tc->cd(1);
	//
	for(Int_t i=0;i<44;i++)th->Fill(Depth[i],dEpianimean[i]);
	for(Int_t i=0;i<44;i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP);
	th->Draw();
	th2->Draw("same");

	tc->cd(2);
	//
	for(Int_t i=0;i<44;i++)th3->Fill(Depth[i],dEpianiloop[i]);
	th3->Draw();
	th2->Draw("same");

	tc->Modified();
	tc->Update();

  //
  gStyle->SetLabelSize(0);
  gStyle->SetNdivisions(1,"XY");
  //
};


void CaloBragg::SWAP( Float_t *A, Float_t *B ){
   Float_t  Tmp = *A;
   *A = *B;
   *B = Tmp;
};



void CaloBragg::LoadParam(){

  //
  elem[0] = 1.00794; //H  1
  elem[1] = 4.0026;  //He 2
  elem[2] = 6.941;   //Li 3
  elem[3] = 9.012182;//Be 4
  elem[4] = 10.811;  //B  5
  elem[5] = 12.0107; //C  6
  elem[6] = 14.00674;//N  7
  elem[7] = 15.9994; //O  8
  elem[8] = 18.9984; //F  9
  elem[9] = 20.1797; //Ne 10
  elem[10] = 22.98977;//Na 11
  elem[11] = 24.3050; //Mg 12
  elem[12] = 26.9815; //Al 13 
  elem[13] = 28.0855; //Si 14
  elem[14] = 30.974;  //P  15
  elem[15] = 32.066;  //S  16
  elem[16] = 35.4527; //Cl 17
  elem[17] = 39.948;  //Ar 18
  elem[18] = 39.0983; //K  19
  elem[19] = 40.078;  //Ca 20
  elem[20] = 44.95591;//Sc 21
  elem[21] = 47.867;  //Ti 22
  elem[22] = 50.9415; //V  23
  elem[23] = 51.9961; //Cr 24
  elem[24] = 54.938049;//Mn 25
  elem[25] = 55.845;  //Fe 26
  elem[26] = 58.9332; //Co 27
  elem[27] = 58.6934; //Ni 28
  elem[28] = 63.546; //Cu 29
  elem[29] = 65.39; //Zn 30
  elem[30] = 69.723; //Ga 31
  elem[31] = 72.61; //Ge 32



//parametri calorimetro
  NPLA = 22;
  NCHA = 96;
  nView = 2;

  AA = 0.96;//mm larghezza strip
  ADIST = 80.5;//mm distanza tra pad
  PIANO = 8.59;//mm distanza 
  
  ySi = 0.38;//mm spessore silicio
  yW = 2.66;//mm spessore tungsteno
  rhoSi = 2.33;//g/cm3 densita' silicio
  rhoW = 19.3;//g/cm3 densita' tugsteno
  MIP = 0.106;//Mev g/cm2 energia al minimo nel silicio per 0.38 mm
  
  emin = 0.;
  
  //parametri bethe-bloch
  pigr = 3.1415;
  Na = 6.02e-23;
  ZA = 0.49; /*Z/A per Si*/
  ISi =182e-06; /*MeV*/
  Me = 0.511; /* MeV*/
  MassP = 931.27;/*MeV*/
  r2 = 7.95e-26; /*ro*ro in cm */
  
};


  
//
void CaloBragg::conversione(){

  //  calcolo spessore Si attraverato in funzione dell'inclinazione
  //  e conversione dello spessore di W in Si e correzione del valore
  //  della Mip pe lo spessore effettivo
  //
  //  in : evento
  //
  //  out: out[0] = gcm2Si = spessore silicio attraversato nel piano
  //       out[1] = WinSi = spessore equivalente in Si del W attraversato
  //       out[2] =  Mip = fattore conversione energia riscalato allo spessore attrversatonel piano

  Float_t SiCross=0.;
  Float_t WCross = 0.;
  Float_t ytgx = 0;
  Float_t ytgy = 0;
  Float_t a = 0.;


  /*silicio*/
  ytgx = ySi * L2->GetCaloLevel2()->tanx[0];
  ytgy = ySi * L2->GetCaloLevel2()->tany[0];
  
  //lunghezza effettiva di silicio attraversata (mm)
  SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy));
  
  spessore[0] = SiCross/10. * rhoSi; //spessore silicio in g/cm2
  
  
  /*tungsteno*/
  ytgx = yW * L2->GetCaloLevel2()->tanx[0];
  ytgy = yW * L2->GetCaloLevel2()->tany[0];
  
  //rapporto tra rilasci energetici nei due materiali
  WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad
  //gcm2W = WCross/10. * rhoW;

 a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);

  //       (g/cm2W)/(g/cm2Si)
  spessore[1] =  a;

  
  //riscala mip allo spessore attraversato
  spessore[2] = MIP*(SiCross/ySi);   

};//end conversione





void CaloBragg::BetheBloch(Float_t *x, Float_t *z, Float_t *Mass, Float_t *gam, Float_t *Bet, Float_t *out){
 
  //rilascio energetico con bethe bloch con correzioni 
  //in:    x: g/cm2
  //       z: carica
  //    Mass: Massa uma
  //     Ene: energia particella MeV//tolta
  //     gam: (etot/massa)
  //     Bet: rad((g2-1)/g2)
  //
  //out: energia rilasciata MeV

 
  Float_t eta =0.;
  Float_t Wmax =0.;
  Float_t lg =0.;
  Float_t Energia=0.;
  Float_t C=0.;
 

  eta = (*gam)*(*Bet);

  //Bet=3/gam;  SQ(*gam) * SQ(*Bet)
  Wmax = 2.* Me * SQ(eta) / (1. + 2.*(*gam)*Me/(*Mass) + SQ(Me)/SQ(*Mass));
  
  lg = 2.* Me * SQ(eta) * Wmax / SQ(ISi);
                                                                                               //  Energia = x* 2 * pigr * Na * r2 * Me * rhoSi *ZA*  SQ(z)/SQ(Bet) * lg;
  C=(0.42237*pow(eta,-2) + 0.0304*pow(eta,-4) - 0.00038*pow(eta,-6))*pow(10,-6)* pow(ISi,2) +
    (3.858*pow(eta,-2) - 0.1668*pow(eta,-4) + 0.00158*pow(eta,-6))*pow(10,-9)*pow(ISi,3);
  
  if(eta <= 0.13) C= C * log(eta/0.0653)/log(0.13/0.0653);      
  //   spessorecm  x ??/massSi x Zsi
  // Energia = (*x)* rhoSi * 0.307/28.09 * 14. *SQ(*z)/SQ(*Bet)*(0.5*log(lg) - SQ(*Bet) - C/14.);
  
 
  /*se ho gia' lo spessore in g/cm2 non se mejo???*/
  Energia = (*x) * 0.307/28.09 * 14. *SQ(*z)/SQ(*Bet)*(0.5*log(lg) - SQ(*Bet) - C/14.);

  *out =Energia;//out

};//end Bethebloch




void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Float_t *Etot, Float_t *out){
  
  /*perdita di energia per ioni pesanti (come da routine geant)*/ 
  //  in : dx    => spessore g/cm2
  //       Z     => carica
  //       Etot  => energia perticella
  //
  //  out:  energia persa


  Float_t Q=0.;
  Float_t v=0.;
                                                        //  Float_t Mass=0.;
  Float_t gam=0.;
  Float_t Bet=0.;
  Float_t dEP=0.;
  
  // gamma  //  Mass = A * MassP; /*in Mev/c2*/ 
  gam =  (*Etot)/(elem[*Z-1]*MassP); // E = gamma M c2
  

  Bet = sqrt((SQ(gam) -1.)/SQ(gam));
  
  v= 121.4139*(Bet/pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet/pow((*Z),(2./3.))));

  //carica effettiva
  Q= (*Z)*(1- (1.034 - 0.1777*exp(-0.08114*(*Z)))*exp(-v));

  //perdita energia per un protone
  Float_t protone =1.;
  Float_t Mass=MassP;
  BetheBloch(dx, &protone, &Mass, &gam, &Bet, &dEP);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP);
  
  *out= (SQ(Q)*(dEP));//*dx;

};//end ELOSS






void CaloBragg::Enetrack(Int_t* Z, Float_t* E0, Float_t* primo, Float_t out[]){
 
  //calcola energia rilasciata sulla traccia (usa ELOSS)
  //  in : Z             =>carica
  //       E0            =>energia    
  //       spess2[3]   => conversione spessore Si, Si in W, mip
  //       primo         => posizione primo piano attraversato
  //
  //  out: array[44]     =>rilasci energetici calcolati per ogni piano[44] dopo il primo(estremi[0][0])
  
  
  
  Float_t dE=0.; //energia rilasciata
  Float_t Ezero= *E0;//energia iniziale

  //azzero energia rilasciata sui piani
  memset(out, 0, 2*NPLA*sizeof(Float_t));
  
  Float_t Massa = (elem[(*Z)-1] * MassP);
  


  //loop piani      (dal primo in cui entra) 
  for( Int_t ipla=((int)(*primo)); ipla< (2*NPLA); ipla++){
    dE=0.;

    //spessore silicio corretto x inclinazione, z, energia, out:rilascio
    ELOSS(&spessore[0], Z, &Ezero, &dE);//spessore in g/cm2!!
    
    if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop
      out[ipla] = Ezero - Massa; //MeV
      return;
      
    }else{
      out[ipla] = dE; //MeV
      Ezero = Ezero - dE;//energia residua
    };

    //se sono su un piano Y (tutti i pari) dopo c'e' il tungsteno
    if(ipla%2 == 0){
      /*tungsteno*/
      dE=0.;
      Float_t sp= spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2
      ELOSS(&sp, Z, &Ezero, &dE);
      
      if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop
	return;
      }else{
	Ezero = Ezero -dE;//energia residua
	//	cout<<"w calc "<<dE<<endl;
      };
    };    
  };//fine loop piani
  
};//end Enetrack





void CaloBragg::chiquadro(Float_t dE[], Float_t out[]){
 
  // calcola chi2 tra energia calcolata e misurata
  // in : dE[44]       =>energia calcolata
  //      calo3[44][2]=> [0]strip attraversata [1]energia misurata per ogni piano
  //      estr2       => array con primo[0][0] e ultimo[1][0] piano attraversati ed energie[][1]
  //
  // out:  array[3]=> (chi2; piani scartati consecutivi(79= >3 quindi frammentato); piani scartati totale)


  Float_t sum = 0.;
  Float_t PianoPrecedente=0.;
  //  Float_t differenza =0.;
  Float_t badplane=0.;
  Float_t badplanetot=0.;
  Float_t w,wi;

  for(Int_t ipla=0; ipla<2*NPLA; ipla++){
   
    //tutti i piani attraversati dalla traiettoria
    if(calorimetro[ipla][0] != -1.){
      w=0.; //normalizzazione;
      wi=1.;//peso
 
      //tolgo piani attraversati dalla traccia ma precedenti il piano individuato come ingresso
      if (ipla<estremi[0][0])  wi=0.;
      
      //tolgo piani attraversati da traccia ma successivi all'ultimo se sono diversi da 0
      //    if((ipla>estremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.;

      //normalizzazione	
      if (calorimetro[ipla][1] != 0.)  w=1./(calorimetro[ipla][1]* MIP);    	
      
      //tolgo piani con rilasci inferiori al 30% del precedente
      if(calorimetro[ipla][1] <= (0.7*PianoPrecedente)){
	
	wi=0.;
	//se sono piani intermedi (non si è fermta) li considero non buoni
	if( (ipla < estremi[1][0]) && (calorimetro[ipla][1] !=0.)){
	  badplane+=1.;
	  badplanetot+=1.;
	};
      };
      

  
      //meno peso ai piani con rilasci maggiori di 1000 MIP
        if(calorimetro[ipla][1] > 1000) wi=0.5;
	    
	//do peso maggiore alle ultime 6 misure
	//	if((ipla >= estr2[1][0]-12)) wi=0.;

      Float_t arg  = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP));

      sum += SQ(arg); // w*wi*(dEpiani[p][v]-(eplane[p][v]*MIP))));//( dEpiani[p][v] - (eplane[p][v]*MIP));

      //se trovo piano non buono (tolto quindi wi=0) non modifico il piano precedente
      if(wi != 0.){
	PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente
	badplane = 0.;//azzero contatore piani scartati consecutivi
      };
    };
	  
    if(badplane > 2) out[1] =79.;
	  
  };//fine loop piani
  //chi2,frammentato,pskip  
  out[0]=sum;

  out[2]=badplanetot;
  
};//end chiquadro



//loopze(Integrale,&zero,&zmin, &zmax, spess1 , estr3, calo1, bestchi2);
void CaloBragg::loopze( Float_t step, Float_t E0,Float_t Zstart, Float_t Zlimite){
//
  //loop su z ed energie per trovare miglior z (ed energia)
  //in:  nloop     => energia massima da provare (nloop x E0)
  //     E0        => energia iniziale (intergale)         
  //     Zstart    => minimo z da cui patire
  //     Zlimite   => z a cui fermarsi (z al minimo di ionizz sul 1o piano)
  //     spessore  => array conversione spessore Si, mip, W
  //     estr1   => array con primo[0][0] e ultimo[1][0] piano attraversati ed energie([x][1])
  //     calo2[44][2]=> [0]strip attraversata [1]energia misurata su ognuno dei 44 piani
  //
  //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2
  //

  Float_t dEplan[2*NPLA];//energia rilasciata calcolata
  memset(dEplan,0,2*NPLA*sizeof(Float_t));

  Int_t Z = 0;// z iniziale

  Float_t Massa = 0.;

  Float_t Stepint =(step)/1000.;//passo per il calcolo di energia

  Float_t energia =0.;//energia del loop

  Float_t chi2[3] = {0,0,0};//out dal calcolo chi2: chi2, piani consecutivi saltati, piani totali saltati
  
  Int_t max=32;//max z di cui so la massa :P
  if((Zlimite)<=31) max=(int)(Zlimite) + 1;
  
  Int_t colmax=32;
  Int_t rowmax=3000;

  Float_t matrixchi2[colmax][rowmax][3];
  memset(matrixchi2, 0, colmax*rowmax*3*sizeof(Float_t));


  //loop elementi  
  for(Int_t inucl=(int)(Zstart); inucl<max; inucl++){
    
    Z= inucl;
   
    Massa = elem[inucl-1]*MassP;
    
    //loop energia 
    for(Int_t iene= 0; iene<3000; iene++){
    
      energia=  Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..
      //cout<<"folse"<<estremi[0][0];
      Enetrack(&Z, &energia, &estremi[0][0], dEplan);//calcola rilascio energetico sui piani
      
      //calcolo chi2
      chiquadro(dEplan,chi2);
      
      if( (chi2[1] != 79) ){//salto quelli che frammentano
	matrixchi2[inucl][iene][0]=chi2[0];//valore chi2 per questo z a questa energia
	matrixchi2[inucl][iene][1]=energia;//energia per questo chi2
	matrixchi2[inucl][iene][2]=chi2[2];//piani saltati ne chi2
      } else {	
	matrixchi2[inucl][iene][0]=1000;//valore chi2 per questo z a questa energia
	matrixchi2[inucl][iene][1]=1000;//energia per questo chi2
	matrixchi2[inucl][iene][2]=1000;//piani saltati ne chi2
      }
    }//fine loop energia

 
 };//fine loop z

  for (Int_t nu=0; nu<colmax; nu++){
    for (Int_t en=0; en<rowmax; en++){
      if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] !=0.)){
 	bestchi2[0]= matrixchi2[nu][en][0];// chi2
 	bestchi2[1]= (Float_t)nu;
 	bestchi2[2]= matrixchi2[nu][en][1];//energia;
 	bestchi2[3]=	matrixchi2[nu][en][2];// totale piani saltati
      }
    }
  }

};//endloopze





void CaloBragg::mediatroncata(){
  //calcolo Z con media troncata e utilizzo questo Z per trovare l'energia migliore
  //in: ordplane[44]   => array con energia dei piani
  //    spess[3]    => conversioni spessore di silicio, w,  mip
  //    estr[2][2]  => primo[0][0] e ultimo[1][0] piano attraversati ed energie[][1]
  //    calo[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano
  //    integrale      => energia totale nel calorimetro considerando il W
  //
  // out[6] chi2,z,Etot,Pskip,flagmaxene,step
  Float_t ordplane[44];//mi serve per la media troncat
  memset(ordplane,0,44*sizeof(Float_t));

  for(Int_t ipla=0; ipla< 2*NPLA; ipla++)  ordplane[ipla]=calorimetro[ipla][1]; //energia del piano

  
  //ordino tutte le energie dei piani in ordine crescente
  Int_t ii=0;
//   while( ii < NPLA-1 ){  
    
//     if(ordplane[ii+1] >= ordplane[ii]){
//       ii++;
//     }else{
//       SWAP( &(ordplane[ii]), &(ordplane[ii+1]));
//       ii=0;
//     };
//   };


  //scelgo 4 piani minimo 'sensati'
  Float_t sum4=0.;
  Int_t pi=0;
  cout<<"sum4="<<sum4<<endl;
  while(pi<2*NPLA){
    if ((ordplane[pi] >0.7)&& sum4==0.){// (ordplane[pi] >50)
      sum4=(ordplane[pi] + ordplane[pi+1] + ordplane[pi+2] + ordplane[pi+3]);
      pi=2*NPLA;
    };
    pi++;
  };
  cout<<"sum4="<<sum4<<endl;
  Float_t Zmean = (sqrt((sum4*MIP)/(4.*spessore[2])));
  cout<<"Zmean="<<Zmean<<endl;
  if(Zmean ==0.) Zmean=1.;
  
  /*z se particella fosse al minimo*/  //energia1piano/mip corretta
  //  Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));

  
  //calcolo energia migliore per Z trovato con media troncata
  //Int_t step = 0;
  Float_t zmin=Zmean;
  //  Float_t bestchi2[4] = {10000,0,0,0};//chi2,z,Etot,Pskip //era bestchi2[5]
  bestchi2[0]=10000.;
  bestchi2[1]=0.;
  bestchi2[2]=0.;
  bestchi2[3]=0.;
  Float_t zero=0.;
    //    step  energia   zstart  zstop Si attrav  1 piano  energie piani  out 
  loopze(Integrale,zero,zmin,zmin);
  
  qtchi2=bestchi2[0];
  qtz=bestchi2[1];
  qtetot=bestchi2[2];
  qtpskip=bestchi2[3];
};//end mediatroncata
  //  cout<<"z media troncata ok"<<endl;



void CaloBragg::Zdaloop(){
  //calcolo Z con un loop su tutti i possibli Z ed energie
  //in: ordplane[44]=> array con energia dei piani
  //    spess1[3]=> conversioni spessore di silicio, w e mip
  //    estr3[2][2]=>  primo[0][0] e ultimo[1][0] piano ed energie
  //    calo1[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano
  //    integrale=> energia totale nel calorimetro considerando il W
  //
  // out[6] chi2,z,Etot,Pskip


  /*z se particella fosse al minimo*/  //energia1piano/mip corretta
  Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));
  
  /*calcolo Z ed E con loop sui vari elementi ed energie*/
  //Int_t step = 0;
  Float_t zmin=1.;
  Float_t bestchitemp[4] = {0,0,0,0};
  //  Float_t bestchi2[4] = {10000,0,0,0};//chi2,z,Etot,Pskip
  bestchi2[0]=10000.;
  bestchi2[1]=0.;
  bestchi2[2]=0.;
  bestchi2[3]=0.;
  Float_t zero=0.;
  //primo loop   
  //     energia   ezero, zstart  zstop Si attrav  1 piano  energie piani  out 
  loopze(Integrale,zero,zmin,zmax);
  
  //secondo loop
  for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i];
  bestchi2[0] = 10000.;
  bestchi2[1] = 0.;
  bestchi2[2] = 0.;
  bestchi2[3] = 0.;//riazzero

  Float_t step = Integrale/1000.;
  zero=bestchitemp[2]-step;//energia da 1 giro - 1step 
  zmin=bestchitemp[1]-2;// zda primo giro
  if(zmin<1)zmin=1;
  zmax=bestchitemp[1]+1;

  loopze(step,zero,zmin,zmax);

  //  cout<<"z loop ok"<<endl;
  
  //chi2,z,Etot,Pskip
  lpchi2=bestchi2[0];
  lpz=bestchi2[1];
  lpetot=bestchi2[2];
  lppskip=bestchi2[3];

};//endZdaloop












