--- calo/flight/CaloBragg/src/CaloBragg.cpp 2008/09/22 20:01:19 1.3 +++ calo/flight/CaloBragg/src/CaloBragg.cpp 2011/06/17 11:56:42 1.12 @@ -3,7 +3,7 @@ ClassImp(CaloBragg); //-------------------------------------- -/** +/* * Default constructor */ CaloBragg::CaloBragg(){ @@ -25,27 +25,35 @@ // 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.; + // 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(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(){ @@ -55,14 +63,17 @@ // 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); printf(" plane not used for loop %f: \n", lppskip); + printf(" ndf: %i \n",ndf); printf("========================================================================\n"); // }; @@ -77,6 +88,56 @@ 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 ){ @@ -119,11 +180,15 @@ for ( Int_t i=0; iGetCaloLevel1()->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; @@ -155,72 +220,116 @@ 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 + 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 - //cout<<"spessore= "<= 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 "<

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 }; p++; }; - //ultimo parte da 44 e sale + + //ultimo parte da 44 e sale p=43; - while( (estremi[1][1] == 0) && (p>(int)estremi[0][0]) ){ - //cout<<"piano "<

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; }; // + 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]); - // cout<<" W "<<0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1])<Print(); @@ -229,75 +338,148 @@ }; +Float_t CaloBragg::Integral(){ + Process(); + + Float_t dEpianiloop[44]; + Int_t tz1=(Int_t)lpz; + Enetrack(&tz1, &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(); - 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],&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 + this->Draw(0.,0.); + +} - Float_t sp= spessore[0]*spessore[1]; - for(Int_t i=0;i<44;i++)Depth[i]=i*sp; +void CaloBragg::Draw(Int_t Z, 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; + 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; + enet = lpetot; + // Enetrack(&tz, &qtetot, &estremi[0][0],&estremi[1][0], dEpianimean);//calcola rilascio energetico sui piani da media troncata + + } + Enetrack(&tz1, &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(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(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(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(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"); - TString hid = Form("cCaloBragg"); - TCanvas *tc = dynamic_cast(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(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(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(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); - // - for(Int_t i=0;i<=estremi[1][0];i++)th3->Fill(Depth[i],dEpianiloop[i]); - th3->Draw(); - th2->Draw("same"); - - tc->Modified(); - tc->Update(); + tc->Modified(); + tc->Update(); // gStyle->SetLabelSize(0); @@ -311,14 +493,14 @@ // elem[0] = 1.00794; //H 1 - elem[1] = 4.0026; //He 2 + elem[1] = 4.002602; //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[8] = 18.9984032; //F 9 elem[9] = 20.1797; //Ne 10 elem[10] = 22.98977;//Na 11 elem[11] = 24.3050; //Mg 12 @@ -344,7 +526,7 @@ elem[31] = 72.61; //Ge 32 -//parametri calorimetro + //parametri calorimetro NPLA = 22; NCHA = 96; nView = 2; @@ -365,7 +547,10 @@ pigr = 3.1415; Na = 6.02e-23; ZA = 0.49; /*Z/A per Si*/ - ISi =182e-06; /*MeV*/ + // 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 */ @@ -393,7 +578,6 @@ Float_t ytgy = 0; Float_t a = 0.; - /*silicio*/ ytgx = ySi * L2->GetCaloLevel2()->tanx[0]; ytgy = ySi * L2->GetCaloLevel2()->tany[0]; @@ -401,8 +585,7 @@ //lunghezza effettiva di silicio attraversata (mm) SiCross = sqrt(SQ(ySi) + SQ(ytgx) + SQ(ytgy)); - spessore[0] = SiCross/10. * rhoSi; //spessore silicio in g/cm2 - + spessore[0] = (SiCross/10.) * rhoSi; //spessore silicio in g/cm2 /*tungsteno*/ ytgx = yW * L2->GetCaloLevel2()->tanx[0]; @@ -412,22 +595,19 @@ 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[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){ +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 @@ -445,23 +625,22 @@ 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(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.); + 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); - // 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 @@ -471,7 +650,7 @@ -void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Float_t *Etot, Float_t *out){ +void CaloBragg::ELOSS(Float_t *dx, Int_t *Z, Float_t *Etot, Float_t *out, Float_t II){ /*perdita di energia per ioni pesanti (come da routine geant)*/ // in : dx => spessore g/cm2 @@ -482,8 +661,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.; @@ -494,15 +672,18 @@ 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.)))); + 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=MassP; - BetheBloch(dx, &protone, &Mass, &gam, &Bet, &dEP);//ene non serve..go gamma.. BetheBloch(dx, 1, MassP, Etot/A, gam, Bet, &dEP); + // 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; @@ -532,14 +713,10 @@ 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!! + ELOSS(&spessore[0], Z, &Ezero, &dE, ISi);//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; @@ -547,29 +724,33 @@ }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= spessore[0]*spessore[1]; //((gcm2Si)*(WinSi))//spessore attraversato in g/cm2 - ELOSS(&sp, Z, &Ezero, &dE); - // cout<<"perdita per piano di W ="< 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 (iplaestremi[1][0]) && (calorimetro[ipla][1] >0.) ) wi=0.; - if((ipla>estremi[1][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); + //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; + //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; - 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)); + Float_t arg = w*wi*(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 + 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 + }; }; - }; - - 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 !!!!!"< 2){ + // printf(" AAAAAAAAAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG\n"); + out[1] =79.; + break; + }; + };//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){ -// +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) @@ -659,22 +856,26 @@ // - - 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)/100.;//passo per il calcolo di energia //era 1000 + 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 ( fzeta > 0. ){ + zmin = fzeta; + max = fzeta+1; + } Int_t colmax=32; Int_t rowmax=3000; @@ -682,147 +883,159 @@ Float_t matrixchi2[colmax][rowmax][3]; memset(matrixchi2, 0, colmax*rowmax*3*sizeof(Float_t)); + Int_t imin = 1-nostep/2; + Int_t imax = nostep/2; //loop elementi - for(Int_t inucl=(int)(Zstart); inucl0.)){ 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 "< 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="< 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"<=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,12 +1051,13 @@ /*z se particella fosse al minimo*/ //energia1piano/mip corretta - Float_t zmax = (sqrt(estremi[0][1]/spessore[2])); - if(zmax<31)zmax=zmax+1; + // 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[4] = {0,0,0,0}; bestchi2[0]=10000.; @@ -851,43 +1065,35 @@ 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"< loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,50); + loopze(Integrale*1.2/500.,Integrale/1000.,zmin,zmax,200); -// cout<<"fine primo loop"< loopze(step,step/2.,zmin,zmax,200); // + loopze(step,step/2.,zmin,zmax,500); // + if ( debug ) printf(" Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step); -// cout<<"z loop ok"<