--- calo/flight/CaloBragg/src/CaloBragg.cpp 2008/06/13 08:11:04 1.1.1.1 +++ 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 ){ @@ -109,7 +170,6 @@ // // - // Always calculate stdedx1 // Int_t view = 0; Int_t plane = 0; @@ -120,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; @@ -134,7 +198,7 @@ ptrack = L2->GetTrack(ntr); if ( ptrack ) track = ptrack->GetCaloTrack(); } else { - track = L2->GetCaloStoredTrack(ntr); + track = L2->GetCaloStoredTrack(ntr); //al momento e' vera solo questa riga }; // if ( !track && ntr >= 0 ){ @@ -155,168 +219,267 @@ // for(Int_t p=0; p<22; p++){ for(Int_t v=0; v<2; v++){ - cout<< L2->GetCaloLevel2()->cibar[p][v]; /*per usare traccia non del calo camboare cibar*/ - calorimetro[(2*p)+1-v][0] = L2->GetCaloLevel2()->cibar[p][v];//strip attraversata - calorimetro[(2*p)+1-v][1] = (epiano[p][v]); //energia del piano //(epiano[p][v])/0.89 + 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 - /*per ogni evento clacolo la conversione mip e w attraversato in equivalente Si*/ + /*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)); - conversione(); // out: 1) g/cm2 Si , 2) spessoreW equivalente in Si, 3)Mip corretta per inclinazione + for(Int_t ipla=0; ipla< 2*NPLA; ipla++) ordplane[ipla]=calorimetro[ipla][1]; //energia del piano - 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*/ + /*trova primo e ultimo piano attraversati*/ Int_t p = 0;//contatore piani - //primo parte da 0 e va in giù - while( estremi[0][1] == 0 && p<(2*NPLA) ){ - cout<<"piano "<

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.)){ - if( (calorimetro[p][0] != -1) && (calorimetro[p][1] >0.7)){ - estremi[1][0]=p; + 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; }; // - cout<<"estremi: in "< 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 rilasciata nel calo (aggiungendo quella 'teorica' nel W )*/ - - for(Int_t pl=0; pl<(2*NPLA); pl++){ - + /*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 - // cout<Print(); - if ( debug ) printf(" esci \n"); + if ( debug ) printf(" fine evento \n"); // }; +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], dEpianimean);//calcola rilascio energetico sui piani da media troncata - Enetrack(&tz1, &lpetot, &estremi[0][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 sp= spessore[0]*spessore[1]; - for(Int_t i=0;i<44;i++)Depth[i]=i*sp; + 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"); - //gStyle->SetLogy(); // + 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<44;i++)th->Fill(Depth[i],dEpianimean[i]); - for(Int_t i=0;i<44;i++)th2->Fill(Depth[i],calorimetro[i][1]*MIP); - th->Draw(); - th2->Draw("same"); - - tc->cd(2); - // - for(Int_t i=0;i<44;i++)th3->Fill(Depth[i],dEpianiloop[i]); - th3->Draw(); - th2->Draw("same"); - - tc->Modified(); - tc->Update(); + tc->Modified(); + tc->Update(); // gStyle->SetLabelSize(0); @@ -325,26 +488,19 @@ }; -void CaloBragg::SWAP( Float_t *A, Float_t *B ){ - Float_t Tmp = *A; - *A = *B; - *B = Tmp; -}; - - void CaloBragg::LoadParam(){ // elem[0] = 1.00794; //H 1 - elem[1] = 4.0026; //He 2 + elem[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 @@ -370,8 +526,7 @@ elem[31] = 72.61; //Ge 32 - -//parametri calorimetro + //parametri calorimetro NPLA = 22; NCHA = 96; nView = 2; @@ -392,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 */ @@ -420,7 +578,6 @@ Float_t ytgy = 0; Float_t a = 0.; - /*silicio*/ ytgx = ySi * L2->GetCaloLevel2()->tanx[0]; ytgy = ySi * L2->GetCaloLevel2()->tany[0]; @@ -428,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]; @@ -439,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 @@ -472,24 +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 @@ -499,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 @@ -510,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.; @@ -522,26 +672,28 @@ 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; -};//end ELOSS - +};//end ELOSS -void CaloBragg::Enetrack(Int_t* Z, Float_t* E0, Float_t* primo, Float_t out[]){ +void CaloBragg::Enetrack(Int_t* Z, Float_t* E0, Float_t* primo,Float_t* ultimo, Float_t out[]){ //calcola energia rilasciata sulla traccia (usa ELOSS) // in : Z =>carica @@ -560,16 +712,11 @@ memset(out, 0, 2*NPLA*sizeof(Float_t)); Float_t Massa = (elem[(*Z)-1] * MassP); - - - //loop piani (dal primo in cui entra) - for( Int_t ipla=((int)(*primo)); ipla< (2*NPLA); ipla++){ + 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; @@ -577,27 +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); - + 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{ Ezero = Ezero -dE;//energia residua - // cout<<"w calc "< 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.; + //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; - //do peso maggiore alle ultime 6 misure - // if((ipla >= estr2[1][0]-12)) wi=0.; - - Float_t arg = w*wi*(dE[ipla] - (calorimetro[ipla][1] * MIP)); - - sum += SQ(arg); // w*wi*(dEpiani[p][v]-(eplane[p][v]*MIP))));//( dEpiani[p][v] - (eplane[p][v]*MIP)); + 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.; - };//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 -//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) // Zstart => minimo z da cui patire // Zlimite => z a cui fermarsi (z al minimo di ionizz sul 1o piano) - // spessore => array conversione spessore Si, mip, W - // estr1 => array con primo[0][0] e ultimo[1][0] piano attraversati ed energie([x][1]) - // calo2[44][2]=> [0]strip attraversata [1]energia misurata su ognuno dei 44 piani // //out: array[4]=> chi2,Zbest,Ebest,piani saltati nel chi2 // - Float_t dEplan[2*NPLA];//energia rilasciata calcolata + memset(dEplan,0,2*NPLA*sizeof(Float_t)); Int_t Z = 0;// z iniziale Float_t Massa = 0.; - Float_t Stepint =(step)/1000.;//passo per il calcolo di energia + Float_t 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; @@ -713,121 +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; + 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 } } } -};//endloopze - +};//endloopze -void CaloBragg::mediatroncata(){ - //calcolo Z con media troncata e utilizzo questo Z per trovare l'energia migliore - //in: ordplane[44] => array con energia dei piani - // spess[3] => conversioni spessore di silicio, w, mip - // estr[2][2] => primo[0][0] e ultimo[1][0] piano attraversati ed energie[][1] - // calo[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano - // integrale => energia totale nel calorimetro considerando il W - // - // out[6] chi2,z,Etot,Pskip,flagmaxene,step - Float_t ordplane[44];//mi serve per la media troncat - memset(ordplane,0,44*sizeof(Float_t)); - for(Int_t ipla=0; ipla< 2*NPLA; ipla++) ordplane[ipla]=calorimetro[ipla][1]; //energia del piano - - //ordino tutte le energie dei piani in ordine crescente - Int_t ii=0; -// while( ii < NPLA-1 ){ - -// if(ordplane[ii+1] >= ordplane[ii]){ -// ii++; -// }else{ -// SWAP( &(ordplane[ii]), &(ordplane[ii+1])); -// ii=0; +// 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.; - //scelgo 4 piani minimo 'sensati' - Float_t sum4=0.; - Int_t pi=0; - cout<<"sum4="<0.7)&& sum4==0.){// (ordplane[pi] >50) - sum4=(ordplane[pi] + ordplane[pi+1] + ordplane[pi+2] + ordplane[pi+3]); - pi=2*NPLA; - }; - pi++; - }; - cout<<"sum4="<=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 @@ -839,42 +1047,53 @@ // calo1[44][2]=> energia[][1] e strip[][0] passaggio su ogni piano // integrale=> energia totale nel calorimetro considerando il W // - // out[6] chi2,z,Etot,Pskip + // out[4] chi2,z,Etot,Pskip /*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*/ - //Int_t step = 0; + Float_t zmin=1.; + Float_t zmax=32.; Float_t bestchitemp[4] = {0,0,0,0}; - // Float_t bestchi2[4] = {10000,0,0,0};//chi2,z,Etot,Pskip + bestchi2[0]=10000.; bestchi2[1]=0.; bestchi2[2]=0.; bestchi2[3]=0.; Float_t zero=0.; - //primo loop - // energia ezero, zstart zstop Si attrav 1 piano energie piani out - loopze(Integrale,zero,zmin,zmax); + //------------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(" Integrale %f , outene %f \n",Integrale,bestchi2[2]); - //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 = Integrale/1000.; - zero=bestchitemp[2]-step;//energia da 1 giro - 1step - zmin=bestchitemp[1]-2;// zda primo giro - if(zmin<1)zmin=1; + 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); // - loopze(step,zero,zmin,zmax); + if ( debug ) printf(" Integrale2 %f , outene %f step %f \n",Integrale,bestchi2[2],step); - // cout<<"z loop ok"<