--- calo/flight/CaloNuclei/src/CaloNuclei.cpp 2007/04/03 15:55:16 1.1 +++ calo/flight/CaloNuclei/src/CaloNuclei.cpp 2008/11/28 16:01:26 1.16 @@ -1,4 +1,7 @@ #include +#include +#include +#include //-------------------------------------- /** @@ -20,18 +23,36 @@ PKT = 0; atime = 0; N = 5; + R = 3; + // + debug = false; + // debug = true; + usetrack = true; // }; void CaloNuclei::Clear(){ // + UN = 0; + tr = 0; + sntr = 0; interplane = 0; preq = 0.; postq = 0.; + stdedx1 = 0.; + ethr = 0.; dedx1 = 0.; dedx3 = 0.; qpremean = 0.; qpremeanN = 0.; + maxrel = 0; + qNmin1 = 0; + qNmin1_w = 0; + charge_siegen1 = 0; + ZCalo_maxrel_b = 0; + ZCalo_dedx_b = 0; + ZCalo_dedx_defl= 0; + ZCalo_Nmin1_defl= 0; // multhit = false; gap = false; @@ -40,20 +61,30 @@ void CaloNuclei::Print(){ // - printf("==============================\n"); - printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime); - printf(" interplane: %i\n",interplane); - printf(" ethr: %f \n",ethr); - printf(" dedx1: %f \n",dedx1); - printf(" dedx3: %f \n",dedx3); - printf(" multhit: %i \n",multhit); - printf(" gap: %i \n",gap); - printf(" preq: %f \n",preq); - printf(" postq: %f \n",postq); - printf(" qpremean: %f \n",qpremean); - printf(" N: %i \n",N); - printf(" qpremeanN: %f \n",qpremeanN); - printf("==============================\n"); + Process(); + // + printf("========================================================================\n"); + printf(" OBT: %u PKT: %u ATIME: %u Track %i Use track %i \n",OBT,PKT,atime,tr,usetrack); + printf(" interplane [number of available dE/dx before interaction]:....... %i\n",interplane); + printf(" ethr [threshold used to determine interplane]:................... %f \n",ethr); + printf(" dedx1 [dE/dx from the first calorimeter plane]:.................. %f \n",dedx1); + printf(" stdedx1 [dE/dx from the first calorimeter plane standalone]:..... %f \n",stdedx1); + printf(" dedx3 [dE/dx (average) if the first 3 Si planes]:................ %f \n",dedx3); + printf(" multhit [true if interplane determined by multiple hits]:........ %i \n",multhit); + printf(" gap [true if interplane determined by a gap]:.................... %i \n",gap); + printf(" preq [total energy in MIP before the interaction plane]:......... %f \n",preq); + printf(" postq [total energy in MIP after the interaction plane]:......... %f \n",postq); + printf(" qpremean [truncated mean using 3 planes and 3 strips]:........... %f \n",qpremean); + printf(" N [no of used plane]:............................................ %i \n",N); + printf(" R [no strip used per plane ]:.................................... %i \n",R); + printf(" qpremeanN [truncated mean using N planes and R strips]:.......... %f \n",qpremeanN); + printf(" qNmin1 [truncated mean using N-1 planes and R strips]: .......... %f \n",qNmin1); + printf(" maxrel [dE/dx of strip with maximum release (I plane)]:.......... %f \n",maxrel); + printf(" ZCalo_maxrel_b [Z from maximum release in I Calo plane vs beta].. %f \n",ZCalo_maxrel_b); + printf(" ZCalo_dedx_b [Z from dedx in I Calo plane vs beta].. ............ %f \n",ZCalo_dedx_b); + printf(" ZCalo_dedx_defl [Z from dedx in I Calo plane vs deflection....... %f \n",ZCalo_dedx_defl); + printf(" ZCalo_Nmin1_defl [Z from truncated mean (N-1 pl) vs deflection].. %f \n",ZCalo_Nmin1_defl); + printf("========================================================================\n"); // }; @@ -64,7 +95,11 @@ void CaloNuclei::Process(){ - // + Process(0); +}; + +void CaloNuclei::Process(Int_t ntr){ + // if ( !L2 ){ printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n"); printf(" ERROR: CaloNuclei variables not filled \n"); @@ -74,11 +109,12 @@ Bool_t newentry = false; // if ( L2->IsORB() ){ - if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){ + if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || ntr != sntr ){ newentry = true; OBT = L2->GetOrbitalInfo()->OBT; PKT = L2->GetOrbitalInfo()->pkt_num; atime = L2->GetOrbitalInfo()->absTime; + sntr = ntr; }; } else { newentry = true; @@ -86,17 +122,89 @@ // if ( !newentry ) return; // - // printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime); + tr = ntr; + // + if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime); // Clear(); // - PamTrack *track = 0; - track = L2->GetTrack(0); + if ( debug ) printf(" Always calculate stdedx1 \n"); // + // Always calculate stdedx1 and maxrel + // + Int_t cont=0; Int_t view = 0; Int_t plane = 0; Int_t strip = 0; + Int_t indx = 0; + Float_t vfpl[96]; + Int_t stfpl[96]; + memset(vfpl, 0, 96*sizeof(Float_t)); + memset(stfpl, 0, 96*sizeof(Int_t)); Float_t mip = 0.; + for ( Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ + // + mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); + // + // put in vfpl vector the energy release on the first plane + // + if ( strip != -1 && view == 1 && plane == 0 ) { + stfpl[indx] = strip; + vfpl[indx] = mip; + indx++; + }; + // + }; + // + if ( debug ) printf(" find energy released along the strip of maximum on the first plane and on the two neighbour strips \n"); + // + // find energy released along the strip of maximum on the first plane and on the two neighbour strips + // + if ( indx > 0 ){ + Int_t mindx = (Int_t)TMath::LocMax(indx,vfpl); + for (Int_t ii=0; ii=0 && stfpl[ii] == (stfpl[mindx]-1) ) stdedx1 += vfpl[ii]; + if ( (mindx+1)<96 && stfpl[ii] == (stfpl[mindx]+1) ) stdedx1 += vfpl[ii]; + // if ( (mindx-1)>=0 && stfpl[ii] == stfpl[mindx-1] ) stdedx1 += vfpl[ii]; + // if ( (mindx+1)<96 && stfpl[ii] == stfpl[mindx+1] ) stdedx1 += vfpl[ii]; + }; + maxrel = vfpl[mindx]; + } else { + stdedx1 = 0.; + maxrel = 0.; + }; + // cout<= 0 ){ + ptrack = L2->GetTrack(ntr); + if ( ptrack ) track = ptrack->GetCaloTrack(); + } else { + track = L2->GetCaloStoredTrack(ntr); + }; + // + if ( !track && ntr >= 0 ){ + printf(" ERROR: cannot find any track!\n"); + printf(" ERROR: CaloNuclei variables not completely filled \n"); + return; + }; + } else { + if ( ntr >= 0 ){ + if ( debug ) printf(" ERROR: you asked not to use a track but you are looking for track number %i !\n",ntr); + if ( debug ) printf(" ERROR: CaloNuclei variables not completely filled \n"); + return; + }; + }; + // // Float_t defethr = 6. * 0.90; Float_t defethr = 6.25; // = (sqrt(9) - 0.5) ** 2.; // @@ -106,21 +214,42 @@ // mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); // - if ( strip != -1 && - view == 1 && - plane == 0 && - ( strip == (track->GetCaloTrack()->tibar[0][1]-1) || strip == (track->GetCaloTrack()->tibar[0][1]-2) || strip == (track->GetCaloTrack()->tibar[0][1]) ) - && true ){ - dedx1 += mip; - }; - if ( strip != -1 && - (( view == 1 && ( plane == 0 || plane == 1 ) ) || - ( view == 0 && plane == 0 )) && - (( view == 0 && ( strip == track->GetCaloTrack()->tibar[0][0] || strip == (track->GetCaloTrack()->tibar[0][0]-1) || strip == (track->GetCaloTrack()->tibar[0][0]-2) )) || - ( view == 1 && ( strip == track->GetCaloTrack()->tibar[0][1] || strip == (track->GetCaloTrack()->tibar[0][1]-1) || strip == (track->GetCaloTrack()->tibar[0][1]-2) )) || - ( view == 1 && ( strip == track->GetCaloTrack()->tibar[1][1] || strip == (track->GetCaloTrack()->tibar[1][1]-1) || strip == (track->GetCaloTrack()->tibar[1][1]-2) ))) && - true ){ - dedx3 += mip; + if ( ntr >= 0 ){ + // + if ( strip != -1 && + view == 1 && + plane == 0 && + ( strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) || strip == (track->tibar[0][1]) ) + && true ){ + dedx1 += mip; + }; + if ( strip != -1 && + (( view == 1 && ( plane == 0 || plane == 1 ) ) || + ( view == 0 && plane == 0 )) && + (( view == 0 && ( strip == track->tibar[0][0] || strip == (track->tibar[0][0]-1) || strip == (track->tibar[0][0]-2) )) || + ( view == 1 && ( strip == track->tibar[0][1] || strip == (track->tibar[0][1]-1) || strip == (track->tibar[0][1]-2) )) || + ( view == 1 && ( strip == track->tibar[1][1] || strip == (track->tibar[1][1]-1) || strip == (track->tibar[1][1]-2) ))) && + true ){ + dedx3 += mip; + }; + } else { + // + if ( strip != -1 && + view == 1 && + plane == 0 && + ( strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) || strip == (L2->GetCaloLevel2()->cibar[0][1]) ) + && true ){ + dedx1 += mip; + }; + if ( strip != -1 && + (( view == 1 && ( plane == 0 || plane == 1 ) ) || + ( view == 0 && plane == 0 )) && + (( view == 0 && ( strip == L2->GetCaloLevel2()->cibar[0][0] || strip == (L2->GetCaloLevel2()->cibar[0][0]-1) || strip == (L2->GetCaloLevel2()->cibar[0][0]-2) )) || + ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[0][1] || strip == (L2->GetCaloLevel2()->cibar[0][1]-1) || strip == (L2->GetCaloLevel2()->cibar[0][1]-2) )) || + ( view == 1 && ( strip == L2->GetCaloLevel2()->cibar[1][1] || strip == (L2->GetCaloLevel2()->cibar[1][1]-1) || strip == (L2->GetCaloLevel2()->cibar[1][1]-2) ))) && + true ){ + dedx3 += mip; + }; }; // }; @@ -133,7 +262,7 @@ // retry: // - // printf("retry\n"); + if ( debug ) printf("retry\n"); // interplane = 0; // @@ -154,30 +283,63 @@ // mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip); // - if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] && - ( strip == (track->GetCaloTrack()->tibar[plane][view]-1) || strip == (track->GetCaloTrack()->tibar[plane][view]-2) || strip == (track->GetCaloTrack()->tibar[plane][view]) ) - && true ){ - // printf(" inside loop: ii %i mip %f view %i plane %i strip %i tibar %i nhit %i splane %i sview %i \n",ii,mip,view,plane,strip,track->GetCaloTrack()->tibar[plane][view]-1,nhit[view],splane[view],sview[view]); - interpl[view] = plane; - interv[view] = view; - if ( splane[view] != plane || sview[view] != view ){ - if ( nhit[view] > 1 ){ - wmulthit[view] = true; - interpl[view] = splane[view]; - interv[view] = sview[view]; - // ii = L2->GetCaloLevel1()->istrip; + if ( ntr >= 0 ){ + if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] && + ( strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) ) + && true ){ + if ( debug ) printf(" inside loop: ii %i mip %f view %i plane %i strip %i tibar %i nhit %i splane %i sview %i \n",ii,mip,view,plane,strip,track->tibar[plane][view]-1,nhit[view],splane[view],sview[view]); + interpl[view] = plane; + interv[view] = view; + if ( splane[view] != plane || sview[view] != view ){ + if ( nhit[view] > 1 ){ + wmulthit[view] = true; + // if ( splane[view] == -1 ) splane[view] = 0; // + // if ( sview[view] == -1 ) sview[view] = view; // + interpl[view] = splane[view]; + interv[view] = sview[view]; }; - if ( plane > splane[view]+gapth ){ - wgap[view] = true; - interpl[view] = splane[view]; - interv[view] = sview[view]; - // ii = L2->GetCaloLevel1()->istrip; + if ( plane > splane[view]+gapth ){ + wgap[view] = true; + // if ( splane[view] == -1 ) splane[view] = 0;// + // if ( sview[view] == -1 ) sview[view] = view; // + interpl[view] = splane[view]; + interv[view] = sview[view]; + }; + splane[view] = plane; + sview[view] = view; + nhit[view] = 1; + } else { + nhit[view]++; + }; + }; + } else { + if ( strip != -1 && mip > ethr && !wmulthit[view] && !wgap[view] && + ( strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) ) + && true ){ + if ( debug ) printf(" inside loop: ii %i mip %f view %i plane %i strip %i cibar %i nhit %i splane %i sview %i \n",ii,mip,view,plane,strip,L2->GetCaloLevel2()->cibar[plane][view]-1,nhit[view],splane[view],sview[view]); + interpl[view] = plane; + interv[view] = view; + if ( splane[view] != plane || sview[view] != view ){ + if ( nhit[view] > 1 ){ + wmulthit[view] = true; + // if ( splane[view] == -1 ) splane[view] = 0; // + // if ( sview[view] == -1 ) sview[view] = view; // + interpl[view] = splane[view]; + interv[view] = sview[view]; + }; + if ( plane > splane[view]+gapth ){ + wgap[view] = true; + // if ( splane[view] == -1 ) splane[view] = 0;// + // if ( sview[view] == -1 ) sview[view] = view; // + interpl[view] = splane[view]; + interv[view] = sview[view]; + }; + splane[view] = plane; + sview[view] = view; + nhit[view] = 1; + } else { + nhit[view]++; }; - splane[view] = plane; - sview[view] = view; - nhit[view] = 1; - } else { - nhit[view]++; }; }; // @@ -185,7 +347,7 @@ // }; // - // printf("conversion interpl %i interv %i multhit %i interplane %i \n",interpl[0],interv[0],multhit,interplane); + if (debug ) printf("conversion interpl %i interv %i multhit %i interplane %i \n",interpl[0],interv[0],multhit,interplane); Int_t winterplane[2] = {-1,-1}; // for ( Int_t view = 0; view < 2; view++){ @@ -224,20 +386,26 @@ }; }; // - // printf("2conversion interpl %i interv %i multhit %i interplane %i \n",interpl[1],interv[1],multhit,interplane); - // printf("3conversion winterpl0 %i winterpl1 %i \n",winterplane[0],winterplane[1]); + if ( debug ) printf("2conversion interpl %i interv %i multhit %i interplane %i \n",interpl[1],interv[1],multhit,interplane); + if ( debug ) printf("3conversion winterpl0 %i winterpl1 %i \n",winterplane[0],winterplane[1]); // Int_t ipl = 0; if ( interplane > 0 ){ // // Calculate preq, postq, qpremean // + cont++; ii = 0; Int_t ind = -1; Int_t qsplane = -1; Int_t qsview = -1; + Int_t ind2 = -1; + Int_t qsplane2 = -1; + Int_t qsview2 = -1; Float_t qme[200]; memset(qme,0,200*sizeof(Float_t)); + Float_t qme2[2112]; + memset(qme2,0,2112*sizeof(Float_t)); // while ( iiGetCaloLevel1()->istrip ){ // @@ -253,15 +421,54 @@ postq += mip; } else { preq += mip; - if ( strip == (track->GetCaloTrack()->tibar[plane][view]-1) || strip == (track->GetCaloTrack()->tibar[plane][view]-2) || strip == (track->GetCaloTrack()->tibar[plane][view]) ){ - if ( qsplane != plane || qsview != view ){ - qsplane = plane; - qsview = view; - ind++; - if ( ind > 199 ) printf(" AAAGH!! \n"); - qme[ind] = 0.; + if ( ntr >= 0 ){ + if ( strip == (track->tibar[plane][view]-1) || strip == (track->tibar[plane][view]-2) || strip == (track->tibar[plane][view]) ){ + if ( qsplane != plane || qsview != view ){ + qsplane = plane; + qsview = view; + ind++; + if ( debug && ind > 199 ) printf(" AAAGH!! \n"); + qme[ind] = 0.; + }; + qme[ind] += mip; }; - qme[ind] += mip; + for ( Int_t ns = 0; ns < R ; ns++){ + Int_t ms = track->tibar[plane][view] - 1 - ns + (R - 1)/2; + if ( strip == ms ){ + if ( qsplane2 != plane || qsview2 != view ){ + qsplane2 = plane; + qsview2 = view; + ind2++; + if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n"); + qme2[ind2] = 0.; + }; + qme2[ind2] += mip; + }; + }; + } else { + if ( strip == (L2->GetCaloLevel2()->cibar[plane][view]-1) || strip == (L2->GetCaloLevel2()->cibar[plane][view]-2) || strip == (L2->GetCaloLevel2()->cibar[plane][view]) ){ + if ( qsplane != plane || qsview != view ){ + qsplane = plane; + qsview = view; + ind++; + if ( debug && ind > 199 ) printf(" AAAGH!! \n"); + qme[ind] = 0.; + }; + qme[ind] += mip; + }; + for ( Int_t ns = 0; ns < R ; ns++){ + Int_t ms = L2->GetCaloLevel2()->cibar[plane][view] - 1 - ns + (R - 1)/2; + if ( strip == ms ){ + if ( qsplane2 != plane || qsview2 != view ){ + qsplane2 = plane; + qsview2 = view; + ind2++; + if ( debug && ind2 > 2112 ) printf(" AAAGH!! \n"); + qme2[ind2] = 0.; + }; + qme2[ind2] += mip; + }; + }; }; }; // @@ -270,35 +477,103 @@ ii++; // }; - // + + + // // here we must calculate qpremean, order vector qme, select 3 lowest measurements and caculate the mean... // - // for (Int_t l=0; l < interplane; l++){ - // printf(" qme[%i] = %f \n",l,qme[l]); - // // if ( qme[ind] < 6.25 ) qme[ind] = 10000.; - // }; + if ( debug ){ + for (Int_t l=0; l < interplane; l++){ + printf(" qme[%i] = %f \n",l,qme[l]); + }; + }; + // Long64_t work[200]; ind = 0; Int_t l = 0; + Int_t RN = 0; Float_t qm = 0.; + Float_t qm2 = 0.; + // + Float_t qmt = ethr*0.8; // *0.9 + // Int_t uplim = TMath::Max(3,N); + Int_t uplim2 = interplane-1; + // while ( l < uplim && ind < interplane ){ - qm = TMath::KOrdStat(interplane,qme,ind,work); - if ( qm >= mesethr*0.9 ){ - if ( l < 3 ) qpremean += qm; - if ( l < N ) qpremeanN += qm; + qm = TMath::KOrdStat((Long64_t)interplane,qme,(Long64_t)ind,work); + if ( qm >= qmt ){ + if ( l < 3 ){ + qpremean += qm; + RN++; + }; l++; - // printf(" value no %i qm %f \n",l,qm); + if ( debug ) printf(" value no %i qm %f qmt %f \n",l,qm,qmt); }; ind++; }; // - qpremean /= l; - qpremeanN /= N; - // - // printf(" charge is %f \n",sqrt(qpremean)); + qpremean /= (Float_t)RN; + ind = 0; + l = 0; + RN = 0; + while ( l < uplim && ind < interplane ){ + qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work); + if ( qm2 >= qmt ){ + if ( l < N ){ + qpremeanN += qm2; + RN++; + }; + l++; + if ( debug ) printf(" qm2 value no %i qm %f qmt %f RN %i \n",l,qm2,qmt,RN); + }; + ind++; + }; + //////////////////////////////////// + //to calculate qNmin1/////////////// + /////////////////////////////////// + //values under threshold + qm2=0; + ind = 0; + l = 0; + RN = 0; + S2=0; + while ( l < uplim2 && ind= qmt ){ + if ( l < (interplane - 1 - S2)){ + qNmin1_w += qm2; + RN++; + }; + l++; + if ( debug ) printf(" qm2 value no %i qm %f qmt %f RN %i \n",l,qm2,qmt,RN); + }; + ind++; + }; + qpremeanN /= (Float_t)RN; + qNmin1_w /= (Float_t)RN; + UN = RN; + ///////set qNmin1 definition/////////// + if (interplane==1 || interplane==2){ + if (dedx1>0) qNmin1=dedx1; + else if (stdedx1>0) qNmin1=stdedx1; + } + else if (interplane > 2){ + qNmin1 = qNmin1_w; + } + //////////////////////////////////// + ////////////////////////////////// // - // printf("preq postq\n"); + if ( debug ) printf(" charge is %f \n",sqrt(qpremean)); // if ( mesethr != ethr && interplane >= 3 && !aldone ){ Float_t mesethr2 = (sqrt(qpremean) - 0.50)*(sqrt(qpremean) - 0.50); @@ -308,10 +583,292 @@ mesethr = mesethr2; }; aldone = true; - if ( mesethr > defethr ) goto retry; + if ( mesethr > defethr ){ + interplane = 0; + preq = 0.; + postq = 0.; + qpremean = 0.; + qpremeanN = 0.; + qNmin1 = 0; + multhit = false; + gap = false; + goto retry; + }; }; }; + + + +//======================================================================= +//=========== charge determination stdedx1 vs. beta =============== +//====================== Siegen method =========================== +//======================================================================= + +// Data from file Calo_Bands_New_7.dat +Float_t C0[9] = {0 , 1 , 2 , 3 , 4 , 5 , 6 , 8 , 90 }; +Float_t B0[9] = {0 , -2.03769 , 7.61781 , 19.7098 , 60.5598 , 57.9226 , 14.8368 , -1358.83 , 8200 }; +Float_t B1[9] = {0 , 0.0211274 , 9.32057e-010 , 4.47241e-07 , 1.44826e-06 , 2.6189e-05 , 0.00278178 , 55.5445 , 0 }; +Float_t B2[9] = {0 , -3.91742 , -20.0359 , -16.3043 , -16.9471 , -14.4622 , -10.9594 , -2.38014 , 0 }; +Float_t B3[9] = {0 , 11.1469 , -6.63105 , -27.8834 , -132.044 , -55.341 , 173.25 , 4115 , 0 }; +Float_t B4[9] = {0 , -14.3465 , -0.485215 , 18.8122 , 117.533 , -14.0898 , -325.269 , -4388.89 , 0 }; +Float_t B5[9] = {0 , 6.24281 , 3.96018 , 0 , -26.1881 , 42.9731 , 182.697 , 1661.01 , 0 }; + +Float_t x1[9],y1[9]; +Int_t n1 = 9; + +Float_t charge = 1000.; +Float_t beta = 100.; + +//------- First try track dependent beta + if( L2->GetTrkLevel2()->GetNTracks()>=1 ){ + PamTrack *TRKtrack = L2->GetTrack(0); + if (fabs(TRKtrack->GetToFTrack()->beta[12]) < 100.) beta = fabs(TRKtrack->GetToFTrack()->beta[12]); + } +//------- If no beta found, try standalone beta + if (beta == 100.) { + ToFTrkVar *ttrack = L2->GetToFStoredTrack(-1); + beta = fabs(ttrack->beta[12]); + } + + if (beta<2.) { // it makes no sense to allow beta=5 or so... + + Float_t mip=0; + mip=stdedx1 ; + + if (mip>0) { + + Float_t betahelp = pow(beta, 1.8); + Float_t ym = mip*betahelp; + Float_t xb = beta; + + for ( Int_t jj=0; jj<9; jj++ ){ + x1[jj] = B0[jj]+B1[jj]*pow(xb,B2[jj])+B3[jj]*xb+B4[jj]*xb*xb+B5[jj]*xb*xb*xb; + y1[jj] = C0[jj]*C0[jj] ; + } + + TGraph *gr1 = new TGraph(n1,x1,y1); + TSpline3 *spl1 = new TSpline3("grs",gr1); // use a cubic spline + Float_t chelp = spl1->Eval(ym); + charge = TMath::Sqrt(chelp); + gr1->Delete(); + spl1->Delete(); + + } // if (mip1>0) + } // beta < 100 + + + charge_siegen1 = charge; + +//======================= end charge Siegen =========================== + + +// //======================================================================= +// //=========== charge determination Maximum release vs. beta =============== +// //====================== Rome method =========================== +// //======================================================================= + + Float_t D0[7] = {0, 3 , 4 , 5 , 6, 8, 90}; + Float_t E1[7] = {0 ,923.553 , 659.842, 1113.97, 3037.25, 3034.84, 0}; + Float_t E2[7] = {0 ,6.92574 , 5.08865, 5.29349, 6.41442, 5.52969, 0}; + Float_t E3[7] = {0 ,9.7227 , 13.18, 23.5444, 38.2057, 63.6784, 80000}; + + Float_t xx1[7],yy1[7]; + n1 = 7; + + charge = 1000.; + mip=0; + + + if (beta<2.) { // it makes no sense to allow beta=5 or so... + + + mip=maxrel; + + if (mip>0) { + Float_t ym = mip; + Float_t xb = beta; + + for ( Int_t jj=0; jjEval(ym); + charge = TMath::Sqrt(chelp); + gr1->Delete(); + spl1->Delete(); + + + } // if (mip1>0) + } // beta < 100 + + + ZCalo_maxrel_b = charge; + + //======================= end charge Rome: maxril vs beta =========================== + + + +// ======================================================================= +// =========== charge determination dedx vs. beta =============== +// ====================== Rome method =========================== +// ======================================================================= + + Float_t F0[7] = {0.,3. ,4., 5. , 6., 8, 90}; + Float_t G1[7] = {0 ,642.935 , 848.684, 1346.05, 3238.82, 3468.6, 0}; + Float_t G2[7] = {0 ,6.2038 , 5.51723, 5.65265, 6.54089, 5.72723, 0}; + Float_t G3[7] = {0 ,9.2421 , 13.9858, 25.3912, 39.6332, 64.5674, 80000}; + + + charge = 1000.; + mip=0; + + + if (beta<2.) { // it makes no sense to allow beta=5 or so... + + + if( L2->GetTrkLevel2()->GetNTracks()>=1 ){ + mip=dedx1; + } + if (mip==0) mip=stdedx1; + + + if (mip>0) { + + Float_t ym = mip; + Float_t xb = beta; + + for ( Int_t jj=0; jjEval(ym); + charge = TMath::Sqrt(chelp); + gr1->Delete(); + spl1->Delete(); + + } //if (mip1>0) + } //beta < 100 + + ZCalo_dedx_b = charge; + + //======================= end charge Rome: dedx vs beta =========================== + + +//======================================================================= +//=========== charge determination dedx vs. defl =============== +//====================== Rome method =========================== +//======================================================================= + + //new + Float_t H0[7] = {0, 3 , 4 , 5 , 6, 8, 90 }; + Float_t I1[7] = {0 , 56.1019, 101.673, 109.225, 150.599, 388.531, 0}; + Float_t I2[7] = {0 , -12.5581, -22.5543, -15.9823, -28.2207, -93.6871, 0}; + Float_t I3[7] = {0 , 11.6218, 19.664, 32.1817, 45.7527, 84.5992, 80000}; + + +// Float_t H0[7] = {0, 3 , 4 , 5 , 6, 8, 90 }; +// Float_t I1[7] = {0 , 56.1019, 101.673, 155, 150.599, 388.531, 0}; +// Float_t I2[7] = {0 , -12.5581, -22.5543, -35.6217, -28.2207, -93.6871, 0}; +// Float_t I3[7] = {0 , 11.6218, 19.664, 34.3311, 45.7527, 84.5992, 8000}; + + + + charge = 1000.; + mip=0; + Float_t defl=0; + + + if (beta<2.) { // it makes no sense to allow beta=5 or so... + + if( L2->GetTrkLevel2()->GetNTracks()>=1 ){ + PamTrack *TRKtrack = L2->GetTrack(0); + mip=dedx1; + if (mip==0) mip=stdedx1; + defl=TRKtrack->GetTrkTrack()->al[4]; + + + if (mip>0 && defl<0.7 && defl>0) { + + Float_t ym = mip; + Float_t xb = defl; + + for ( Int_t jj=0; jjEval(ym); + charge = TMath::Sqrt(chelp); + gr1->Delete(); + spl1->Delete(); + + } // if (mip1>0 && defl<0.5 && defl>0) + }//Ntrack>=1 + } //beta < 100 + + ZCalo_dedx_defl = charge; + +//======================= end charge Rome: dedx vs defl =========================== + + +//============================================================================================ +//=========== charge determination Truncated mean (N-1 planes) vs. defl =================== +//================================ Rome method ======================================== +//============================================================================================ + + Float_t L0[7] = {0, 3 , 4 , 5 , 6, 8, 90}; + Float_t M1[7] = {0 , 63.0145, 120.504, 173.663, 245.33, 236.517, 0}; + Float_t M2[7] = {0 , -15.005, -31.0635, -39.4988, -60.5011, -46.3992, 0}; + Float_t M3[7] = {0 , 12.5037, 22.8652, 35.2907, 51.4678, 86.4155, 8000}; + + charge = 1000.; + mip=0; + + + if (beta<2.) { // it makes no sense to allow beta=5 or so... + + if( L2->GetTrkLevel2()->GetNTracks()>=1 ){ + mip=qNmin1; + + if (mip>0 && defl<0.7 && defl>0) { + + Float_t ym = mip; + Float_t xb = defl; + + for ( Int_t jj=0; jjEval(ym); + charge = TMath::Sqrt(chelp); + gr1->Delete(); + spl1->Delete(); + + } // if (mip1>0 && defl<0.5 && defl>0) + }//Ntrack>=1 + } //beta < 100 + + ZCalo_Nmin1_defl = charge; + +//======================= end charge Rome: Nmin1 vs defl =========================== + + + + + // - // printf(" esci \n"); + if ( debug ) this->Print(); + if ( debug ) printf(" esci \n"); // };