--- calo/flight/CaloBragg/src/CaloBragg.cpp 2008/11/28 15:59:25 1.5 +++ calo/flight/CaloBragg/src/CaloBragg.cpp 2011/06/17 11:56:42 1.12 @@ -25,23 +25,28 @@ // 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.; @@ -60,14 +65,15 @@ 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 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"); // }; @@ -82,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 ){ @@ -124,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; @@ -160,8 +220,8 @@ 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 }; }; @@ -177,27 +237,27 @@ //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((Long_64_t)44,ordplane,(Long_64_t)ind,work); - if (qm >= 0.7 ){ - if ( RN < 4 ){ - sum4 += qm; - RN++; - }; + 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.; + 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*/ @@ -212,7 +272,8 @@ }; 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]) ){ if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){ @@ -222,9 +283,37 @@ 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 @@ -232,13 +321,15 @@ Integrale+= 0.5*((calorimetro[pl-1][1] * MIP)+(calorimetro[pl][1] * MIP))*(spessore[1]); }; }; - Integrale=24000;//Integrale*1000; + //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*/ - Zdaloop(); // out: 1)chi2, 2)z, 3)Etot, 4)Pskip + if ( debug ) printf(" call Zdaloop with integrale %f \n",Integrale); + Zdaloop(); // out: 1)chi2, 2)z, 3)Etot, 4)Pskip if ( debug ) this->Print(); @@ -247,73 +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.); + +} + +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; + 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); -// }; - // - 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"); + 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"); - tc->Modified(); - tc->Update(); + tc->Modified(); + tc->Update(); // gStyle->SetLabelSize(0); @@ -327,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 @@ -360,7 +526,7 @@ elem[31] = 72.61; //Ge 32 -//parametri calorimetro + //parametri calorimetro NPLA = 22; NCHA = 96; nView = 2; @@ -381,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 */ @@ -416,7 +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]; @@ -426,21 +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 @@ -458,20 +625,23 @@ 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); - 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 @@ -480,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 @@ -502,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=(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); + // 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; @@ -543,7 +716,7 @@ 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; @@ -551,13 +724,22 @@ }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); + 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, &Ezero, &dE,II); if((Ezero-dE) <= Massa){//se l'energia depositata e' maggiore dell'energia della perticella stop return; }else{ @@ -587,66 +769,83 @@ 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 + // + 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 (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)){ // 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; + //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)); + 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 + 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){ - out[1] =79.; - break; - }; - - };//fine loop piani - //chi2,frammentato,pskip - out[0]=sum; - out[2]=badplanetot; + //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){ -// +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) @@ -657,21 +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)/1000.;//passo per il calcolo di energia + 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; @@ -679,43 +883,49 @@ 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 @@ -841,12 +1051,13 @@ /*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; + // 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.; @@ -856,8 +1067,13 @@ Float_t zero=0.; //------------primo loop ---------------------- // energia ezero, zstart zstop - loopze(Integrale,zero,zmin,zmax); + // 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(" Integrale %f , outene %f \n",Integrale,bestchi2[2]); //------------secondo loop ---------------------- for(Int_t i=0;i<4;i++) bestchitemp[i]=bestchi2[i]; @@ -867,9 +1083,16 @@ bestchi2[3] = 0.;//riazzero Float_t step = bestchitemp[2];// - zero=0; // qualsiasi altro valore peggiora le cose - zmin=zmax=bestchitemp[1]; - loopze(step,zero,zmin,zmax); // + 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); // + + if ( debug ) printf(" Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step); //chi2,z,Etot,Pskip