--- calo/flight/CaloBragg/src/CaloBragg.cpp 2008/09/22 20:01:19 1.3
+++ calo/flight/CaloBragg/src/CaloBragg.cpp 2008/11/28 16:01:38 1.6
@@ -3,7 +3,7 @@
ClassImp(CaloBragg);
//--------------------------------------
-/**
+/*
* Default constructor
*/
CaloBragg::CaloBragg(){
@@ -32,10 +32,10 @@
//
tr = 0;
sntr = 0;
- qtchi2 = 0.;
- qtz = 0.;
- qtetot = 0.;
- qtpskip = 0.;
+// qtchi2 = 0.;
+// qtz = 0.;
+// qtetot = 0.;
+// qtpskip = 0.;
lpchi2 = 0.;
lpz = 0.;
lpetot = 0.;
@@ -44,8 +44,11 @@
memset(spessore,0,3*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(){
@@ -55,10 +58,12 @@
//
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(" 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(" energy from loop %f: \n", lpetot);
@@ -162,17 +167,46 @@
/*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
- //cout<<"spessore= "<<spessore[0]<<" gcm2si; "<<spessore[1]<<"W in Si; "<<spessore[2]<<"mip "<<endl;
+ /*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 giù
- while( estremi[0][1] == 0 && p<(2*NPLA) ){
- //cout<<"piano "<<p <<"; strip "<<calorimetro[p][0]<<"; en "<<calorimetro[p][1]<<endl;
+ //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] != -1) && (calorimetro[p][1] >0.7)){ //0.7 soglia minima
+ // 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
};
@@ -180,17 +214,14 @@
};
//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.)){
+ while( (estremi[1][1] <= 0.) && (p>(int)estremi[0][0]) ){
if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){
- estremi[1][0]=p;//era p
+ estremi[1][0]=p;//
estremi[1][1]=calorimetro[p][1] *MIP;//energia in MeV
};
p = p-1;
};
//
-
/*integrale: energia totale rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/
for(Int_t pl=0; pl<(2*NPLA); pl++){
@@ -199,28 +230,15 @@
// 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<<" W "<<0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1])<<endl;
};
};
-
+ Integrale=24000;//Integrale*1000;
/*z ed energia con media troncata*/
-// cout<<"Media troncata"<<endl;
- mediatroncata(); // out: 1)chi2, 2)z, 3)Etot, 4)Pskip
+ // mediatroncata(); // out: 1)chi2, 2)z, 3)Etot, 4)Pskip
/*z ed energia con loop*/
-// cout<<"Zdaloop"<<endl;
- // Zdaloop(); // out: 1)chi2, 2)z, 3)Etot, 4)Pskip
-
-
- /*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
+ Zdaloop(); // out: 1)chi2, 2)z, 3)Etot, 4)Pskip
if ( debug ) this->Print();
@@ -233,12 +251,12 @@
Process();
- Float_t dEpianimean[44];
+// Float_t dEpianimean[44];
Float_t dEpianiloop[44];
Float_t Depth[44];
- Int_t tz=(Int_t)qtz;
+// Int_t tz=(Int_t)qtz;
Int_t tz1=(Int_t)lpz;
- Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata
+// Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata
Enetrack(&tz1, &lpetot, &estremi[0][0],&estremi[1][0], dEpianiloop);//calcola rilascio energetico sui piani da loop
Float_t sp= spessore[0]*spessore[1];
@@ -247,25 +265,23 @@
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);
+// tc->Divide(1,2);
};
//
- TString thid = Form("hCaloBragg");
- TH2F *th = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
- if ( th ) th->Delete();
+// 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);
+// th = new TH2F(thid,thid,300,-0.5,300.,1000,0.,150.);
+// th->SetMarkerStyle(20);
// };
//
TString thid2 = Form("hCaloBragg2");
@@ -285,9 +301,9 @@
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++)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();
+// th->Draw();
th2->Draw("same");
tc->cd(2);
@@ -393,7 +409,6 @@
Float_t ytgy = 0;
Float_t a = 0.;
-
/*silicio*/
ytgx = ySi * L2->GetCaloLevel2()->tanx[0];
ytgy = ySi * L2->GetCaloLevel2()->tany[0];
@@ -403,7 +418,6 @@
spessore[0] = SiCross/10. * rhoSi; //spessore silicio in g/cm2
-
/*tungsteno*/
ytgx = yW * L2->GetCaloLevel2()->tanx[0];
ytgy = yW * L2->GetCaloLevel2()->tany[0];
@@ -417,7 +431,6 @@
// (g/cm2W)/(g/cm2Si)
spessore[1] = a;
-
//riscala mip allo spessore attraversato
spessore[2] = MIP*(SiCross/ySi);
@@ -457,12 +470,8 @@
(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.);
+ Energia = (*x) * 0.307/28.09 * 14. *SQ(*z)/SQ(*Bet)*(0.5*log(lg) - SQ(*Bet) - C/14.);
*out =Energia;//out
@@ -482,8 +491,7 @@
Float_t Q=0.;
- Float_t v=0.;
- // Float_t Mass=0.;
+ Float_t v=0.;
Float_t gam=0.;
Float_t Bet=0.;
Float_t dEP=0.;
@@ -501,7 +509,7 @@
//perdita energia per un protone
Float_t protone =1.;
- Float_t Mass=MassP;
+ Float_t Mass=(elem[*Z-1]*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;
@@ -532,12 +540,8 @@
Float_t Massa = (elem[(*Z)-1] * MassP);
-
- //loop piani (dal primo in cui entra)
- // for( Int_t ipla=((int)(*primo)); ipla< (2*NPLA); ipla++){
- for( Int_t ipla=((int)(*primo)); ipla<= ((int)(*ultimo)); ipla++){
+ 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, &Ezero, &dE);//spessore in g/cm2!!
if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop
@@ -548,28 +552,23 @@
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);
- // cout<<"perdita per piano di W ="<<dE<<endl;
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
-// for(Int_t i=0;i< 44; i++)cout<<"deposito energetico (teorico) per il piano i"<<i<<", "<<out[i]<<endl;
-
-};//end Enetrack
+ };//fine loop piani
+};//end Enetrack
@@ -590,9 +589,8 @@
Float_t w,wi;
for(Int_t ipla=0; ipla<2*NPLA; ipla++){
-
//tutti i piani attraversati dalla traiettoria
- if(calorimetro[ipla][0] != -1.){
+ if(calorimetro[ipla][0] != -1.){ //
w=0.; //normalizzazione;
wi=1.;//peso
@@ -604,13 +602,13 @@
if((ipla>estremi[1][0])) wi=0.;
//normalizzazione
- if (calorimetro[ipla][1] != 0.) w=1./(calorimetro[ipla][1]* MIP);
+ 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)){
+ 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 è fermta) li considero non buoni
- if( (ipla <= estremi[1][0]) && (calorimetro[ipla][1] !=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.;
};
@@ -622,31 +620,31 @@
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.){
+ if(wi != 0.){//
PianoPrecedente= calorimetro[ipla][1];//tengo piano precedente
badplane = 0.;//azzero contatore piani scartati consecutivi
};
};
-
+
+ //da Emi
if(badplane > 2){
out[1] =79.;
break;
- };
-
+ };
+
};//fine loop piani
//chi2,frammentato,pskip
out[0]=sum;
out[2]=badplanetot;
-
-// if(out[1] ==79.)cout<<"frammentato !!!!!"<<endl;
};//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)
@@ -659,7 +657,6 @@
//
-
Float_t dEplan[2*NPLA];//energia rilasciata calcolata
memset(dEplan,0,2*NPLA*sizeof(Float_t));
@@ -667,7 +664,7 @@
Float_t Massa = 0.;
- Float_t Stepint =(step)/100.;//passo per il calcolo di energia //era 1000
+ Float_t Stepint =(step)/1000.;//passo per il calcolo di energia
Float_t energia =0.;//energia del loop
@@ -691,7 +688,7 @@
Massa = elem[inucl-1]*MassP;
//loop energia
- for(Int_t iene= 0; iene<1000; iene++){
+ 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 <=??
energia= Massa + (E0)+ iene*Stepint;//gli do un'energia totale (momento) massa+energia cinetica, aumentando la cinetica..
@@ -700,7 +697,7 @@
//calcolo chi2
chiquadro(dEplan,chi2);
- if( (chi2[1] != 79) ){//salto quelli che frammentano
+ 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 nel chi2
@@ -715,114 +712,120 @@
};//fine loop z
- // for (Int_t nu=0; nu<colmax; nu++){
- // for (Int_t en=0; en<rowmax; en++){
+
+ //Emi
for (Int_t nu=(int)(Zstart); nu<max; nu++){
for (Int_t en=0; en<1000; en++){
- if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] !=0.)){
+ if((matrixchi2[nu][en][0]<bestchi2[0]) && (matrixchi2[nu][en][0] >0.)){
bestchi2[0]= matrixchi2[nu][en][0];// chi2
bestchi2[1]= (Float_t)nu; // z
bestchi2[2]= matrixchi2[nu][en][1];//energia;
- bestchi2[3]= matrixchi2[nu][en][2];// totale piani saltati
+ bestchi2[3]= matrixchi2[nu][en][2];// totale piani saltati
}
}
}
- //==========================//
- // Int_t testz = (Int_t)bestchi2[1];
- // Enetrack(&testz, &bestchi2[2], &estremi[0][0],&estremi[1][0], dEplan);//calcola rilascio energetico sui piani
- // for(Int_t i=0;i<=estremi[1][0];i++){
- // cout<<"dEplan "<<dEplan[i]*10<<endl;
- // cout<<"calorimetro "<<calorimetro[i][1]<<endl;
- // }
- //==========================//
-};//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
- //Int_t ii=0;
- 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)/(4.*spessore[2])));
- if(Zmean ==0.) Zmean=1.;
- if ( Zmean < 1. ) Zmean = 1.;
- //
- // Zmean = 1.;
+};//endloopze
-// cout<<"sum4="<<sum4<<endl;
-// cout<<"Zmean="<<Zmean<<endl;
-
- //calcolo energia migliore per Z trovato con media troncata
-
- Float_t zmin=Zmean;
-
- bestchi2[0]=10000.;
- bestchi2[1]=0.;
- bestchi2[2]=0.;
- bestchi2[3]=0.;
- Float_t zero=0.;
- // cout<<"inizio media troncata"<<endl;
- // cout<<"input: step "<<Integrale<<", zero "<<zero<<", zmin "<<zmin<<", zmax"<<zmin<<endl;
- // step energia zstart zstop
- loopze(Integrale,zero,zmin,zmin);
-
- // cout<<"fine media troncata"<<endl;
- // cout<<"output: chi2 "<<bestchi2[0]<<", z"<<bestchi2[1]<<", Etot "<<bestchi2[2]<<", Pskip "<<bestchi2[3]<<endl;
+// 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);
- qtchi2=bestchi2[0];
- qtz=bestchi2[1];
- qtetot=bestchi2[2];
- qtpskip=bestchi2[3];
-};//end mediatroncata
- // cout<<"z media troncata ok"<<endl;
+// //======== 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
@@ -838,7 +841,7 @@
/*z se particella fosse al minimo*/ //energia1piano/mip corretta
- Float_t zmax = (sqrt(estremi[0][1]/spessore[2]));
+ 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*/
@@ -851,43 +854,23 @@
bestchi2[2]=0.;
bestchi2[3]=0.;
Float_t zero=0.;
- //primo loop
- // energia ezero, zstart zstop Si attrav 1 piano energie piani out
-
- // cout<<"inizio primo loop"<<endl;
- // cout<<"input: step "<<Integrale<<", zero "<<zero<<", zmin "<<zmin<<", zmax"<<zmax<<endl;
-
+ //------------primo loop ----------------------
+ // energia ezero, zstart zstop
loopze(Integrale,zero,zmin,zmax);
-// cout<<"fine primo loop"<<endl;
-// cout<<"output: chi2 "<<bestchi2[0]<<", z"<<bestchi2[1]<<", Etot "<<bestchi2[2]<<", Pskip "<<bestchi2[3]<<endl;
- //secondo loop
+ //------------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 = bestchitemp[2]/100.;//Integrale/10.;//era 1000
- //zero=bestchitemp[2]-step/2.;//energia da 1 giro - 1step
- //Float_t step = bestchitemp[2]/1000.;
- //zero=bestchitemp[2]-step;
- zero=0;
- zmin=bestchitemp[1]-2;// zda primo giro
- if(zmin<1)zmin=1;
- zmax=bestchitemp[1]+1;//
-
-// cout<<"inizio secondo loop"<<endl;
-// cout<<"input: step "<<step<<", zero "<<zero<<", zmin "<<zmin<<", zmax"<<zmax<<endl;
-
+ Float_t step = bestchitemp[2];//
+ zero=0; // qualsiasi altro valore peggiora le cose
+ zmin=zmax=bestchitemp[1];
loopze(step,zero,zmin,zmax); //
-// cout<<"fine secondo loop"<<endl;
-// cout<<"output: chi2 "<<bestchi2[0]<<", z"<<bestchi2[1]<<", Etot "<<bestchi2[2]<<", Pskip "<<bestchi2[3]<<endl;
-
-
-// cout<<"z loop ok"<<endl;
//chi2,z,Etot,Pskip
lpchi2=bestchi2[0];