CaloBragg/src/CaloBragg.cpp

#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;
  usepl18x = false;
  newchi2 = false;
  usenewBB = false;
  fzeta = -1.;
  //
};

void CaloBragg::Clear(){
  //
  ndf = 0;
  tr = 0;
  sntr = 0;
  //   qtchi2 = 0.;
  //   qtz = 0.;
  //   qtetot = 0.;
  //   qtpskip = 0.;
  lpchi2 = 0.;
  lpz = 0.;
  lpisotope= 0.;
  lpetot = 0.;
  lppskip = 0.;

  memset(calorimetro,0,44*2*sizeof(Float_t));
  memset(spessore,0,4*sizeof(Float_t));
  memset(estremi,0,2*2*sizeof(Float_t));
  Integrale=0.;
 
  for(Int_t l=0;l<44;l++){
    calorimetro[l][0]=-1.;
  }

};

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(" first plane: %f \n", estremi[0][0]);
  printf(" last plane: %f \n", estremi[1][0]);
  //   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(" isotope from loop %f: \n", lpisotope);
  printf(" energy from loop %f: \n", lpetot);
  printf(" plane not used for loop %f: \n", lppskip);
  printf(" ndf: %i \n",ndf);
  printf("========================================================================\n");
  //
};

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


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


void CaloBragg::CleanPlanes(Float_t epiano[22][2]){
  //  return;
  Int_t hitplanes = 0;
  for (Int_t i = 0; i<22; i++){
    for (Int_t j = 1; j>=0; j--){
      if ( epiano[i][j] > 0.7 ) hitplanes++;
    };
  };
  Float_t lowlim = 0.85;
  Float_t dedxone = 0.;
  Float_t step1 = 0.8*L2->GetCaloLevel2()->qtot/(Float_t)hitplanes;
  while ( dedxone < step1 ){
    for (Int_t i = 0; i<22; i++){
      for (Int_t j = 1; j>=0; j--){
        if ( epiano[i][j] >= step1 && dedxone < 0.7 ) dedxone = epiano[i][j];
      };
    };    
  }
  if ( dedxone < 0.7 ){
    for (Int_t i = 0; i<22; i++){
      for (Int_t j = 1; j>=0; j--){
        if ( epiano[i][j] > 0. && dedxone < 0.7 ) dedxone = epiano[i][j];
      };
    };    
  }
  //
  //  printf(" dedxone = %f step1 %f  \n",dedxone,step1);
  Bool_t revulsera = false;
  Bool_t nullius = false;
  Int_t nulliferus = 0;
  for (Int_t i = 0; i<22; i++){
    for (Int_t j = 1; j>=0; j--){
      if ( epiano[i][j] < dedxone*lowlim ){
        //        printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);
        epiano[i][j] = 0.;
      } else {
        //x        printf(" %i %i epiano %f limit %f nulliferus %i nullius %i \n",i,j,epiano[i][j],dedxone*lowlim,nulliferus,nullius);
        nulliferus = 0;
        revulsera = true;
      };
      if ( epiano[i][j] < 0.7 && revulsera ) nulliferus++;
      if ( nulliferus > 10 ) nullius = true;
      if ( nullius ) epiano[i][j] = 0.;
    };
  };

}


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();

  //
  //
  //
  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);
    //
    if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
    //
    epiano[plane][view]+=mip;
    //
    //
  };
  //
  this->CleanPlanes(*&epiano);
  //
  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);  //al momento e' vera solo questa riga
    };
    //
    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++){
      /*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 calcolo 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

  /*settaggio della soglia per il loop sulla determinazione del piano di partenza */
  Float_t ordplane[44];//mi serve per la media troncata
  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

  Long64_t work[200];
  Int_t ind = 0;
  //Int_t l = 0;
  Int_t RN = 0;
  Float_t sum4 = 0.;
  Float_t qm = 0.;
  while ( RN < 4 && ind < 44 ){
    qm = TMath::KOrdStat((Long64_t)44,ordplane,(Long64_t)ind,work);
    if (qm >= 0.7 ){    
      if ( RN < 4 ){
        sum4 += qm;
        RN++;
      };
    };
    ind++;
  };
  //
  //sum4 /= (Float_t)RN;
  Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));
  if(Zmean ==0.) Zmean=1.;
  if ( Zmean < 1. ) Zmean = 1.;


  /*trova primo e ultimo piano attraversati*/
  Int_t p = 0;//contatore piani
  //per il primo parte da 0 e va in giu'
  while( estremi[0][1] <= 0.  &&  p<(2*NPLA) ){ // era ==0 ma ricorda i problemi con Float == !!!!!
    //    if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >50.)){
    //   if( (calorimetro[p][0] >0) && (calorimetro[p][1]*MIP >0.3)){ //0.7 mip = 70MeV soglia minima
    if( (calorimetro[p][0] >0) && (calorimetro[p][1]*MIP >Zmean*0.7)){ // 70% della MIP
      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]) ){
    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;
  };
  //

  Float_t lastok = 0.;
  //  Bool_t goback = false;
  for ( int o = 0; o < estremi[1][0]; o++ ){
    //
    if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1];
    if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok;    
    //    if ( calorimetro[o][1] < 0.7 ) goback = true;
    //
  };
  lastok = 0.;
  //  if ( goback ){
  for ( int o = estremi[1][0]; o >= 0;  o-- ){
    //
    //    printf(" goback1: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
    if ( o < estremi[1][0] && calorimetro[o][1] > calorimetro[o+1][1]*1.2 && lastok > 0. ) calorimetro[o][1] = lastok;  
    if ( calorimetro[o][1] > 0.7 ) lastok = calorimetro[o][1];
    if ( calorimetro[o][1] < 0.7 && lastok > 0. ) calorimetro[o][1] = lastok;    
    //    printf(" goback2: o %i calo %f lastok %f \n",o,calorimetro[o][1],lastok);
    //
  };
  //  };

  if ( startZero ) {
    estremi[0][0] = 0.;
    //    estremi[0][1] = 0.;
  }
  
  /*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/
  for(Int_t pl=0; pl<(2*NPLA); pl++){ 
    //    printf(" integrale: calorimetro %f  \n",calorimetro[pl][1]);
    //calcolo intergale in unita di spessori di silicio
    Integrale += calorimetro[pl][1] * MIP;//piano di silicio
    // 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]);
    };
  };
  //Integrale=24000;//Integrale*1000;
  Integrale *= 1000.;
  
  /*z ed energia con media troncata*/
  //  mediatroncata();  // out:  1)chi2, 2)z, 3)Etot, 4)Pskip
  
  /*z ed energia con loop*/
  if ( debug ) printf(" call Zdaloop with integrale %f \n",Integrale);
  Zdaloop(); // out:  1)chi2, 2)z, 3)Etot, 4)Pskip

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


Float_t CaloBragg::Integral(){
  Process();

  Float_t dEpianiloop[44];
  Int_t tz1=(Int_t)lpz;
  Int_t ti1=(Int_t)lpisotope;

  Enetrack(&tz1, &ti1 , &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop


  Float_t integ = 0.;
  for(Int_t i=0;i<=estremi[1][0];i++){
    //    integ += dEplan[i];
    //printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]);
    integ += dEpianiloop[i];
    //    printf(" step %i integ %f deplan %f \n",i,integ,dEpianiloop[i]);
  }
  return integ;
}

Float_t CaloBragg::LastIntegral(){
  Process();

  Float_t integ = 0.;
  for(Int_t i=0;i<=estremi[1][0];i++){
    integ += dEplan[i];
    //printf(" step %i integ %f deplan %f \n",i,integ,dEplan[i]);
  }
  return integ;
}


void CaloBragg::Draw(){

  Process();

  this->Draw(0.,0.,0.);

}

void CaloBragg::Draw(Int_t Z, Int_t isotope, Float_t enetot){

  //  Float_t dEpianimean[44];
  Float_t dEpianiloop[44];
  Float_t Depth[44];
  //  Int_t tz=(Int_t)qtz;
  Int_t tz1= Z; 
  Int_t ti1= isotope; 
  Float_t enet = enetot;
  //  Float_t enet = lpetot;
  
  if ( Z > 0. && enetot > 0. ){
    estremi[0][0] = 0;
    estremi[1][0] = 43;


    Float_t    ytgx = 0.;
    Float_t    ytgy = 0.;
    
    //lunghezza effettiva di silicio attraversata (mm)
    Float_t    SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy));
    
    spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2
    
    /*tungsteno*/
    
    //rapporto tra rilasci energetici nei due materiali
    Float_t WCross = sqrt((yW*yW) + (ytgx*ytgx) + (ytgy*ytgy));//mm* rapporto lunghezze rad
    //gcm2W = WCross/10. * rhoW;
    
    //       (g/cm2W)/(g/cm2Si)
    spessore[3] =  (WCross/10.) * rhoW;
    Float_t a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);
    spessore[1] =  a;
    //riscala mip allo spessore attraversato
    spessore[2] = MIP*(SiCross/ySi);   

  } else {
    tz1=(Int_t)lpz;
    ti1=(Int_t)lpisotope;
    enet = lpetot;
    //  Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata
  
  }
  Enetrack(&tz1, &ti1, &enet, &estremi[0][0],&estremi[1][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");
  //
  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);
  //    };
  //
  tc->cd();
  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<=estremi[1][0];i++)th->Fill(Depth[i],dEpianimean[i]);
  for(Int_t i=0;i<=estremi[1][0];i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP);
  //    th->Draw();
  th2->Draw("same");

  //  tc->cd(2);
  tc->cd();
  //
  for(Int_t i=0;i<=estremi[1][0];i++){
    th3->Fill(Depth[i],dEpianiloop[i]);
    //    printf(" i %i Depth %f depianiloop %f \n",i,Depth[i],dEpianiloop[i]);
  }
  th3->Draw();
  th2->Draw("same");

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

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


void CaloBragg::LoadParam(){

  // elem[Z-1][isotop] 0 is the most common one
  // 
 
  elem[0][0] = 1.00782;   //H  1
  elem[0][1] = 2.01410;   // 2H (Isotope)
  elem[0][2] = -1.; 
  elem[0][3] = -1.; 
  elem[0][4] = -1.; 
  elem[0][5] = -1.; 
  elem[0][6] = -1.; 
 
  elem[1][0] = 4.002602;  //He 2
  elem[1][1] = 3.016029;  // 3He (Isotope)
  elem[1][2] = -1.; 
  elem[1][3] = -1.; 
  elem[1][4] = -1.; 
  elem[1][5] = -1.; 
  elem[1][6] = -1.; 

  elem[2][0] = 7.016004;  //Li 3
  elem[2][1] = 6.015123;  //6Li (Isotope)
  elem[2][2] = -1.; 
  elem[2][3] = -1.; 
  elem[2][4] = -1.; 
  elem[2][5] = -1.; 
  elem[2][6] = -1.; 

  elem[3][0] = 9.012182;  //Be 4
  elem[3][1] = 10.01353;  //10Be (Isotope) (most stable)
  elem[3][2] = -1.; 
  elem[3][3] = -1.; 
  elem[3][4] = -1.; 
  elem[3][5] = -1.; 
  elem[3][6] = -1.; 

  elem[4][0] = 11.0093;   //B 5
  elem[4][1] = 10.01294;  //10B  (Isotope)
  elem[4][2] = -1.; 
  elem[4][3] = -1.; 
  elem[4][4] = -1.; 
  elem[4][5] = -1.; 
  elem[4][5] = -1.; 

  elem[5][0] = 12.0107;   //C  6
  elem[5][1] = 13.00335;  //13C (Isotope)
  elem[5][2] = -1.; 
  elem[5][3] = -1.; 
  elem[5][4] = -1.; 
  elem[5][5] = -1.; 
  elem[5][5] = -1.; 

  elem[6][0] = 14.00674;  //N  7
  elem[6][1] = 15.00011;  //15N  (Isotope)
  elem[6][2] = -1.; 
  elem[6][3] = -1.; 
  elem[6][4] = -1.; 
  elem[6][5] = -1.; 
  elem[6][5] = -1.; 

  elem[7][0] = 15.99491;  //O  8
  elem[7][1] = 17.99916;  //18O  (Isotope)
  elem[7][2] = 16.99916;  //17O  (Isotope)
  elem[7][3] = -1.; 
  elem[7][4] = -1.; 
  elem[7][5] = -1.; 
  elem[7][5] = -1.; 

  elem[8][0] = 18.99840;  //F  9
  elem[8][1] = -1.; 
  elem[8][2] = -1.; 
  elem[8][3] = -1.; 
  elem[8][4] = -1.; 
  elem[8][5] = -1.; 
  elem[8][5] = -1.; 

  elem[9][0] = 19.99244;  //Ne 10
  elem[9][1] = 21.99138;  //22Ne (Isotope)
  elem[9][2] = 20.99384;  //21Ne 10
  elem[9][3] = -1.; 
  elem[9][4] = -1.; 
  elem[9][5] = -1.; 
  elem[9][6] = -1.; 

  elem[10][0] = 22.98977; //Na 11
  elem[10][1] = 21.99444; //22Na (Isotope) (most stable)
  elem[10][2] = -1.; 
  elem[10][3] = -1.; 
  elem[10][4] = -1.; 
  elem[10][5] = -1.; 
  elem[10][6] = -1.; 

  elem[11][0] = 23.98504; //Mg 12
  elem[11][1] = 25.98259; //26Mg (Isotope)
  elem[11][2] = 24.98504; //25Mg (Isotope)
  elem[11][3] = -1.; 
  elem[11][4] = -1.; 
  elem[11][5] = -1.; 
  elem[11][6] = -1.; 

  elem[12][0] = 26.98154; //Al 13 
  elem[12][1] = 25.98489; //26Al (Isotope) (most stable) 
  elem[12][2] = -1.; 
  elem[12][3] = -1.; 
  elem[12][4] = -1.; 
  elem[12][5] = -1.; 
  elem[12][6] = -1.; 

  elem[13][0] = 27.97692; //Si 14
  elem[13][1] = 28.97649; //29Si (Isotope) 
  elem[13][2] = 29.97377; //30Si (Isotope)
  elem[13][3] = -1.; 
  elem[13][4] = -1.; 
  elem[13][5] = -1.; 
  elem[13][6] = -1.; 

  elem[14][0] = 30.97376;  //P  15
  elem[14][1] = -1.; 
  elem[14][2] = -1.; 
  elem[14][3] = -1.; 
  elem[14][4] = -1.; 
  elem[14][5] = -1.; 
  elem[14][6] = -1.; 

  elem[15][0] = 31.97207;  //S  16
  elem[15][1] = 33.96787;  //34S  (Isotope)
  elem[15][2] = 32.97146;  //33S  (Isotope)
  elem[15][3] = 35.96708;  //36S  (Isotope)
  elem[15][4] = -1.; 
  elem[15][5] = -1.; 
  elem[15][6] = -1.; 

  elem[16][0] = 34.96885; //Cl 17
  elem[16][1] = 36.96831; //37Cl 17
  elem[16][2] = 35.96890; //36Cl (Isotope)
  elem[16][3] = -1.; 
  elem[16][4] = -1.; 
  elem[16][5] = -1.; 
  elem[16][6] = -1.; 

  elem[17][0] = 39.962383;  //Ar 18 
  elem[17][1] = 35.967545;  //36Ar (Isotope) 
  elem[17][2] = 37.962732;  //38Ar (Isotope) 
  elem[17][3] = 38.964313;  //39Ar (Isotope)
  elem[17][4] = -1.; 
  elem[17][5] = -1.; 
  elem[17][6] = -1.; 

  elem[18][0] = 38.963707; //K  19 
  elem[18][1] = 40.961825; //41K  (Isotope)
  elem[18][2] = 39.963998; //40K  (Isotope)
  elem[18][3] = -1.; 
  elem[18][4] = -1.; 
  elem[18][5] = -1.; 
  elem[18][6] = -1.; 

  elem[19][0] = 39.962590;  //Ca 20
  elem[19][1] = 43.955482;  //44Ca (Isotope)
  elem[19][2] = 41.958618;  //42Ca (Isotope)
  elem[19][3] = 42.958767;  //43Ca (Isotope)
  elem[19][4] = 45.953693;  //46Ca (Isotope)
  elem[19][5] = 40.962278;  //41Ca (Isotope)
  elem[19][6] = -1.; 

  elem[20][0] = 44.955912;//Sc 21
  elem[20][1] = -1.; 
  elem[20][2] = -1.; 
  elem[20][3] = -1.; 
  elem[20][4] = -1.; 
  elem[20][5] = -1.; 
  elem[20][6] = -1.; 

  elem[21][0] = 47.947946;  //Ti 22
  elem[21][1] = 45.952632;  //46Ti (Isotope)
  elem[21][2] = 46.951763;  //47Ti (Isotope)
  elem[21][3] = 48.947870;  //49Ti (Isotope)
  elem[21][4] = 49.944791;  //50Ti (Isotope)
  elem[21][5] = 43.959690;  //44Ti (Isotope) (half life 60y)
  elem[21][6] = -1.; 

  elem[22][0] = 50.943960; //V  23
  elem[22][1] = 49.947158; //50V  (Isotope)
  elem[22][2] = -1.; 
  elem[22][3] = -1.; 
  elem[22][4] = -1.; 
  elem[22][5] = -1.; 
  elem[22][6] = -1.; 

  elem[23][0] = 51.940507; //Cr 24
  elem[23][1] = 52.940649; //53Cr (Isotope)
  elem[23][2] = 49.946044; //50Cr (Isotope)
  elem[23][3] = 53.938880; //54Cr (Isotope)
  elem[23][4] = -1.; 
  elem[23][5] = -1.; 
  elem[23][6] = -1.; 

  elem[24][0] = 54.938049;//Mn 25
  elem[24][1] = 52.941290;//53Mn (Isotope)
  elem[24][2] = -1.; 
  elem[24][3] = -1.; 
  elem[24][4] = -1.; 
  elem[24][5] = -1.; 
  elem[24][6] = -1.; 

  elem[25][0] = 55.934937;  //Fe 26
  elem[25][1] = 53.939610;  //54Fe (Isotope)
  elem[25][2] = 56.935394;  //57Fe (Isotope)
  elem[25][3] = 57.933276;  //58Fe (Isotope)
  elem[25][4] = 59.934072;  //58Fe (Isotope)

  elem[26][0] = 58.933195; //Co 27
  elem[26][1] = 59.933817; //60Co (Isotope)
  elem[26][2] = -1.; 
  elem[26][3] = -1.; 
  elem[26][4] = -1.; 
  elem[26][5] = -1.; 
  elem[26][6] = -1.; 


  elem[27][0] = 57.935343; //Ni 28
  elem[27][1] = 61.928345; //62Ni (Isotope)
  elem[27][2] = 59.930786; //60Ni (Isotope)
  elem[27][3] = 60.931056; //61Ni (Isotope)
  elem[27][4] = 63.927966; //64Ni (Isotope)
  elem[27][5] = 58.934346; //59Ni (Isotope)
  elem[27][6] = -1.; 

  elem[28][0] = 62.929597; //Cu 29
  elem[28][1] = 64.927789; //65Cu (Isotope)
  elem[28][2] = -1.; 
  elem[28][3] = -1.; 
  elem[28][4] = -1.; 
  elem[28][5] = -1.; 
  elem[28][6] = -1.; 

  elem[29][0] = 63.929142; //Zn 30
  elem[29][1] = 65.926033; //66Zn (Isotope)
  elem[29][2] = 67.924844; //68Zn (Isotope)
  elem[29][3] = 66.927127; //67Zn (Isotope)
  elem[29][4] = 69.925319; //70Zn (Isotope)
  elem[29][5] = -1.; 
  elem[29][6] = -1.; 

  elem[30][0] = 68.925573; //Ga 31
  elem[30][1] = 70.924701; //71Ga (Isotope)
  elem[30][2] = -1.; 
  elem[30][3] = -1.; 
  elem[30][4] = -1.; 
  elem[30][5] = -1.; 
  elem[30][6] = -1.; 

  elem[31][0] = 73.921177; //Ge 32
  elem[31][1] = 71.922075; //72Ge (Isotope)
  elem[31][2] = 69.924247; //70Ge (Isotope)
  elem[31][3] = 75.921403; //76Ge (Isotope)
  elem[31][4] = 73.923459; //73Ge (Isotope)
  elem[31][5] = -1.; 
  elem[31][6] = -1.; 


  //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*/
  ISi = 171e-06; /*MeV*/
  IW  = 735e-06; /*MeV*/
  //  ISi =0.0001059994; /*GeV!!*/ no era giusto!!
  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;

  //       (g/cm2W)/(g/cm2Si)
  spessore[3] =  (WCross/10.) * rhoW;
  a=(WCross/SiCross)*(rhoW/rhoSi)*(1.145/1.664);  //(gcm2W)/(SiCross/10. * rhoSi)* (1.145/1.664);
  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, Float_t II){
 
  //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.;
  Float_t INo = ISi;

  if ( usenewBB ) INo = II;

  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(INo);
  //  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(INo,2.) +
    (3.858*pow(eta,-2.) - 0.1668*pow(eta,-4.) + 0.00158*pow(eta,-6.))*pow(10.,-9.)*pow(INo,3.);
  
  if(eta <= 0.13) C= C * log(eta/0.0653)/log(0.13/0.0653);      
 
  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, Int_t *isotope, Float_t *Etot, Float_t *out, Float_t II){
  
  /*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 gam=0.;
  Float_t Bet=0.;
  Float_t dEP=0.;
  
  // gamma  //  Mass = A * MassP; /*in Mev/c2*/ 
  gam =  (*Etot)/(elem[*Z-1][*isotope]*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.))));
  v= 121.4139*(Bet*pow((*Z),(2./3.))) + 0.0378*sin(190.7165*(Bet*pow((*Z),(2./3.)))); // EMI AAAAGGH!!

  //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=(elem[*Z-1]*MassP); //EMI
  //  BetheBloch(dx, &protone, &Mass, &gam, &Bet, &dEP);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP);

  BetheBloch(dx, &protone, &MassP, &gam, &Bet, &dEP, II);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP); //EMI
  
  *out= (SQ(Q)*(dEP));//*dx;


};//end ELOSS


void CaloBragg::Enetrack(Int_t* Z, Int_t* isotope, Float_t* E0, Float_t* primo,Float_t* ultimo, 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][*isotope] * MassP);

  for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){
    dE=0.;
    //spessore silicio corretto x inclinazione, z, energia, out:rilascio
    ELOSS(&spessore[0], Z , isotope , &Ezero, &dE, ISi);//spessore in g/cm2!!

        if(dE!=dE) return; //controlla che non sia un NaN
    
        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
      if ( debug ) printf(" zompa %i out %f dE %f ezero %f \n",ipla,out[ipla],dE,Ezero);
    };
    //se sono su un piano Y (tutti i pari) dopo c'e' il tungsteno
    if(ipla%2 == 0){
      /*tungsteno*/
      dE=0.;
      Float_t sp = 0.;
      Float_t II = ISi;
      if ( usenewBB ){
        sp = spessore[3];        
        II = IW;
      } else {
        sp = spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato  in g/cm2
      }
      //      printf(" sp %f II %f \n",sp,II);
      ELOSS(&sp, Z, isotope , &Ezero, &dE,II);
      if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop
        return;
      }else{
        Ezero = Ezero -dE;//energia residua
      };
    };    

  };//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 badplane=0.;
  Float_t badplanetot=0.;
  Float_t w,wi;
  //
  if ( newchi2 ){
    ndf = 0;
    sum = 0.;
    for( Int_t ipla=((int)(estremi[0][0])); ipla<= ((int)(estremi[1][0])); ipla++){
      sum += pow((dE[ipla] - (calorimetro[ipla][1] * spessore[2]))/(0.05*dE[ipla]),2.);
      //      printf(" quiqui: dE %f calor %f spessore[2] %f \n",dE[ipla],spessore[2]*calorimetro[ipla][1],spessore[2]);
      ndf++;
    }
    ndf -= 2;
    if ( ndf > 0 ) sum /= (float)ndf;
    out[0] = sum;
    out[1] = 0.;
    out[2] = (int)(estremi[1][0])-ndf;
    //    printf(" sum %f ndf %i \n ",sum,ndf);
  } else {
    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.;
        if((ipla>estremi[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)){ // cosi' i piani senza rilascio non vengono considerati nel calcolo del chi2
          wi=0.;
          //se sono piani intermedi (non si e' 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;
        if(calorimetro[ipla][1] > 1200.) wi=0.5;
           if(debug) printf("chiquadro start \n "); 
        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));
        if(debug){
          printf("dedx  calcolata  %f e reale  %f  \n",dE[ipla],(calorimetro[ipla][1] * 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
        };
      };
          
      //da Emi
      if(badplane > 2){      
        //      printf(" AAAAAAAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG\n");
        out[1] =79.;
        break;
      };
  
    };//fine loop piani
    //chi2,frammentato,pskip  
    out[0]=sum;
    out[2]=badplanetot;
  }  
};//end chiquadro


void CaloBragg::loopze( Float_t step, Float_t E0,Float_t Zstart, Float_t Zlimite, Int_t nostep = 1000){
  //
  //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)
  //
  //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2
  //

 //printf("entrato");
  memset(dEplan,0,2*NPLA*sizeof(Float_t));

  Int_t Z = 0;// z iniziale
  
  Int_t isotope=0;

  Float_t Massa = 0.;

  Float_t Stepint =(step)/(Float_t)nostep;//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 zmin = (int)Zstart;  
  Int_t max=32;//max z di cui so la massa :P
  if((Zlimite)<=31) max=(int)(Zlimite) + 1;

if(debug)  printf("loopze inizio max  %d \n",max); 
  if ( fzeta > 0. ){
    zmin = fzeta;
    max = fzeta+1;
  }
  
  Int_t colmax=32;
  Int_t rowmax=3000;
  Int_t isomax=7;

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

  Int_t imin = 1-nostep/2;
  Int_t imax = nostep/2;

  //loop elementi  
  for(Int_t inucl=zmin; inucl<max; inucl++){
    
    Z= inucl;

  //loop isotopi
  while ( elem[inucl-1][isotope] > 0. ){  
    Massa = elem[inucl-1][isotope]*MassP; 

    //loop energia 
    Int_t iene2 = -1;

    //    for(Int_t iene= 0; iene<1000; iene++){// da non cambiare in base a Stepint altrimenti cambia la matrice bestchi2!!!cosi' non raggiungo mai integrale!!!!! mettere <=??

    for(Int_t iene= imin; iene<imax; iene++){// da non cambiare in base a Stepint altrimenti cambia la matrice bestchi2!!!cosi' non raggiungo mai integrale!!!!! mettere <=??
         iene2++;
                 energia=  Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..
         
         
                 if( fene > 0. ) energia=fene;   //forza l'energia
                 if (debug) printf("loopze energia %f, z %d, isotopo %d ,iene %d\n",energia,Z,isotope,iene);
//              printf(" energia %f , forzata %f \n",energia,fene);
                 Enetrack(&Z, &isotope, &energia, &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani
      
         chiquadro(dEplan,chi2); //calcolo chi2
                 if (debug) printf("loopze chi %f \n",chi2[0]);
     if(debug && TMath::Finite(chi2[0])==1 && (TMath::IsNaN(chi2[0])!=1) )    printf("loopze    fin mat %f \n",chi2[0]); 
      //      printf(" last deplan from: Z = %i iene %i energia %f chi2 %f \n",inucl,iene,energia,chi2[0]);
         if( (chi2[1] != 79.) ){//salto quelli che frammentano
                matrixchi2[inucl][isotope][iene2][0]=chi2[0];//valore chi2 per questo z a questa energia
                matrixchi2[inucl][isotope][iene2][1]=energia;//energia per questo chi2
                        matrixchi2[inucl][isotope][iene2][2]=chi2[2];//piani saltati nel chi2
            if( fene > 0. ) break;
         } else {       
                        matrixchi2[inucl][isotope][iene2][0]=1000.;//valore chi2 per questo z a questa energia
                        matrixchi2[inucl][isotope][iene2][1]=1000.;//energia per questo chi2
                        matrixchi2[inucl][isotope][iene2][2]=1000.;//piani saltati nel chi2
                break;
      }

    }//fine loop energia 
 
  isotope++; //incremento il contatore isotopi
  }//fine loop isotopi
  isotope=0; //riazzero il contatore isotopi
 
  }//fine loop z

  isotope=0;//non dovrebbe servire

  //Emi
  for (Int_t nu=zmin; nu<max; nu++){
    while(elem[nu-1][isotope]> 0.){
        for (Int_t en=0; en<nostep; en++){
        if((matrixchi2[nu][isotope][en][0]<bestchi2[0]) && (matrixchi2[nu][isotope][en][0] >0.)){
          bestchi2[0]= matrixchi2[nu][isotope][en][0];// chi2
          bestchi2[1]= (Float_t)nu; // z
          bestchi2[2]= matrixchi2[nu][isotope][en][1];//energia;
          bestchi2[3]= matrixchi2[nu][isotope][en][2];// totale piani saltati
      bestchi2[4]= (Float_t)isotope; //isotopo
        }
      }
    isotope++;
    }
   isotope=0;
  }

};//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[4] chi2,z,Etot,Pskip

//   Float_t ordplane[44];//mi serve per la media troncata
//   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

//     Long64_t work[200];
//     Int_t ind = 0;
//     //Int_t l = 0;
//     Int_t RN = 0;
//     Float_t sum4 = 0.;
//     Float_t qm = 0.;
//     //
//     //Float_t qmt = ethr*0.8; // *0.9
//     //
//     //Int_t uplim = TMath::Max(3,N);
//     //
//     while ( RN < 4 && ind < 44 ){
//       qm = TMath::KOrdStat(44,ordplane,ind,work);
//       if (qm >= 0.7 ){       
//      if ( RN < 4 ){
//        sum4 += qm;
//        RN++;
//      };
// //   l++;
// //   if ( debug ) printf(" value no %i qm %f sum4 %f \n",l,qm,sum4); 
//       };
//       ind++;
//     };
//     //
//     sum4 /= (Float_t)RN;
//     Float_t Zmean = (sqrt((sum4*MIP)/(((Float_t)RN)*spessore[2])));//ma non e'/1??
//     if(Zmean ==0.) Zmean=1.;
//     if ( Zmean < 1. ) Zmean = 1.;

    
// //     Zmean =round(Zmean);
// //     if(Zmean <1.) Zmean=1.;

// //     if(Zmean >0.)Zmean =round(Zmean);

//     //======== per i nuclei=======
//     if (Zmean >=2.){
//     ind = 0;
//     RN = 0;
//     sum4 = 0.;
//     qm = 0.;
//     while ( RN < 4 && ind < 44 ){
//       qm = TMath::KOrdStat(44,ordplane,ind,work);
//       if (qm >= (Zmean*Zmean)-Zmean*Zmean*0.2 ){     
//      if ( RN < 4 ){
//        sum4 += qm;
//        RN++;
//      };
//       };
//       ind++;
//     };
//     //
//     sum4 /= (Float_t)RN;
//     Zmean = (sqrt((sum4*MIP)/(4.*spessore[2])));//ma non e' /1??
//     }

  
//   //calcolo energia migliore per Z trovato con media troncata
//   //  Float_t zmin=Zmean;
//     Float_t zmin=round(Zmean);
  
//   bestchi2[0]=10000.;
//   bestchi2[1]=0.;
//   bestchi2[2]=0.;
//   bestchi2[3]=0.;
//   Float_t zero=0.;

//   //    step   energia zstart zstop  
//   loopze(Integrale,zero,zmin,zmin);
  

//   qtchi2=bestchi2[0];
//   qtz=bestchi2[1];
//   qtetot=bestchi2[2];
//   qtpskip=bestchi2[3];
// };//end mediatroncata


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[4] chi2,z,Etot,Pskip


  /*z se particella fosse al minimo*/  //energia1piano/mip corretta
  //  Float_t zmax = round(sqrt(estremi[0][1]/spessore[2]));
  //  if(zmax<31)zmax=zmax+1;
  
  /*calcolo Z ed E con loop sui vari elementi ed energie*/

  Float_t zmin=1.;
  Float_t zmax=32.;
  Float_t bestchitemp[5] = {0,0,0,0,0};

  bestchi2[0]=numeric_limits<Float_t>::max();
  bestchi2[1]=0.;
  bestchi2[2]=0.;
  bestchi2[3]=0.;
  bestchi2[4]=0.;
  Float_t zero=0.;
  //------------primo loop   ----------------------
  //     energia   ezero, zstart  zstop 
  //  loopze(Integrale,zero,zmin,zmax);

  //->  loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,50);
  loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,200);

  //  loopze(Integrale*2.,Integrale/100.,zmin,zmax);
  if ( debug) printf("Zdaloop start Integrale %f , outene %f \n",Integrale,bestchi2[2]);
  
  //------------secondo loop  ----------------------
  for(Int_t i=0;i<5;i++) bestchitemp[i]=bestchi2[i];
  bestchi2[0]=numeric_limits<Float_t>::max();
  bestchi2[1] = 0.;
  bestchi2[2] = 0.;
  bestchi2[3] = 0.;
  bestchi2[4] = 0.;//riazzero

  Float_t step = bestchitemp[2];//
  zero=0.;  // qualsiasi altro valore peggiora le cose
  //  zmin=zmax=bestchitemp[1];
  zmin=bestchitemp[1]-1;
  zmax=bestchitemp[1]+1;
  //loopze(step,zero,zmin,zmax); //

  //->  loopze(step,step/2.,zmin,zmax,200); //
  loopze(step,step/2.,zmin,zmax,500); //

  //step = bestchitemp[2];//

  //loopze(step/2,3*step/4.,zmin,zmax,500); //

  if ( debug ) printf("Zdaloop Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step);

  //chi2,z,Etot,Pskip
  lpchi2=bestchi2[0];
  lpz=bestchi2[1];
  lpetot=bestchi2[2];
  lppskip=bestchi2[3];
  lpisotope=bestchi2[4];
};//endZdaloop

