--- calo/flight/CaloNuclei/src/CaloNuclei.cpp 2007/04/05 13:38:32 1.4 +++ calo/flight/CaloNuclei/src/CaloNuclei.cpp 2015/04/20 09:12:05 1.20 @@ -1,14 +1,17 @@ #include +#include +#include +#include //-------------------------------------- /** * Default constructor */ -CaloNuclei::CaloNuclei(){ - Clear(); -}; +// CaloNuclei::CaloNuclei(){ +// Clear(); +// }; -CaloNuclei::CaloNuclei(PamLevel2 *l2p){ +CaloNuclei::CaloNuclei(PamLevel2 *l2p,const char* alg){ // Clear(); // @@ -23,19 +26,35 @@ R = 3; // debug = false; + // debug = true; + usetrack = true; + usepl18x = false; + trkAlg = alg; // }; 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; @@ -47,19 +66,26 @@ Process(); // printf("========================================================================\n"); - printf(" OBT: %u PKT: %u ATIME: %u Track %i \n",OBT,PKT,atime,tr); - 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(" 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(" 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"); // }; @@ -85,11 +111,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; @@ -103,13 +130,87 @@ // Clear(); // - PamTrack *track = 0; - track = L2->GetTrack(ntr); + 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); + // + if ( !usepl18x && view==0 && plane==18 ) mip = 0.; + // + // + // 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,trkAlg); + if ( ptrack ) track = ptrack->GetCaloTrack(); + } else { + track = L2->GetCaloStoredTrack(ntr); + }; + // + if ( !track && ntr >= 0 ){ + printf(" ERROR: cannot find any track! \n"); + cout << "ERROR: trk.algorythm --> "<= 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.; // @@ -119,21 +220,44 @@ // 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 ( !usepl18x && view==0 && plane==18 ) mip = 0.; + // + 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; + }; }; // }; @@ -167,32 +291,66 @@ // 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 ){ - 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->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; - // if ( splane[view] == -1 ) splane[view] = 0; // - // if ( sview[view] == -1 ) sview[view] = view; // - interpl[view] = splane[view]; - interv[view] = sview[view]; + if ( !usepl18x && view==0 && plane==18 ) mip = 0.; + // + // + 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; - // 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]++; + }; + }; + } 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]++; }; }; // @@ -247,6 +405,7 @@ // // Calculate preq, postq, qpremean // + cont++; ii = 0; Int_t ind = -1; Int_t qsplane = -1; @@ -263,6 +422,9 @@ // mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip); // + if ( !usepl18x && view==0 && plane==18 ) mip = 0.; + // + // if ( strip != -1 ){ if ( view == 0 ){ ipl = (1 + plane) * 2; @@ -273,29 +435,55 @@ 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 ( debug && 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->GetCaloTrack()->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.; + 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; }; - 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; + }; + }; + }; }; // }; @@ -303,7 +491,9 @@ ii++; // }; - // + + + // // here we must calculate qpremean, order vector qme, select 3 lowest measurements and caculate the mean... // if ( debug ){ @@ -315,38 +505,87 @@ 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); + qm = TMath::KOrdStat((Long64_t)interplane,qme,(Long64_t)ind,work); if ( qm >= qmt ){ - if ( l < 3 ) qpremean += qm; + if ( l < 3 ){ + qpremean += qm; + RN++; + }; l++; if ( debug ) printf(" value no %i qm %f qmt %f \n",l,qm,qmt); }; ind++; }; // - qpremean /= l; - // + qpremean /= (Float_t)RN; ind = 0; l = 0; + RN = 0; while ( l < uplim && ind < interplane ){ - qm2 = TMath::KOrdStat(interplane,qme2,ind,work); + qm2 = TMath::KOrdStat((Long64_t)interplane,qme2,(Long64_t)ind,work); if ( qm2 >= qmt ){ - if ( l < N ) qpremeanN += qm2; + if ( l < N ){ + qpremeanN += qm2; + RN++; + }; l++; - if ( debug ) printf(" qm2 value no %i qm %f qmt %f \n",l,qm2,qmt); + if ( debug ) printf(" qm2 value no %i qm %f qmt %f RN %i \n",l,qm2,qmt,RN); }; ind++; }; - // - qpremeanN /= l; + //////////////////////////////////// + //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; + } + //////////////////////////////////// + ////////////////////////////////// // if ( debug ) printf(" charge is %f \n",sqrt(qpremean)); // @@ -358,10 +597,284 @@ 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->GetNTracks(trkAlg)>=1 ){ + PamTrack *TRKtrack = L2->GetTrack(0,trkAlg); + 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[9] = {0, 1, 2, 3 , 4 , 5 , 6, 8, 90}; + Float_t E1[9] = {0, 500, 500, 923.553 , 659.842, 1113.97, 3037.25, 3034.84, 0}; + Float_t E2[9] = {0, 11.0, 7.5, 6.92574 , 5.08865, 5.29349, 6.41442, 5.52969, 0}; + Float_t E3[9] = {0, 1.2, 4, 9.7227 , 13.18, 23.5444, 38.2057, 63.6784, 80000}; + + Float_t xx1[9],yy1[9]; + n1 = 9; + + 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[9] = {0., 1., 2., 3. ,4., 5. , 6., 8, 90}; + Float_t G1[9] = {0, 500, 500, 642.935 , 848.684, 1346.05, 3238.82, 3468.6, 0}; + Float_t G2[9] = {0, 11, 7.5, 6.2038 , 5.51723, 5.65265, 6.54089, 5.72723, 0}; + Float_t G3[9] = {0, 1.2, 4, 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->GetNTracks(trkAlg)>=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 =========================== +//======================================================================= + + Float_t H0[9] = {0, 1, 2, 3 , 4 , 5 , 6, 8, 90 }; + Float_t I1[9] = {0, 3.5, 40, 56.1019, 101.673, 109.225, 150.599, 388.531, 0}; + Float_t I2[9] = {0, -1, -13.6, -12.5581, -22.5543, -15.9823, -28.2207, -93.6871, 0}; + Float_t I3[9] = {0, 1, 5.3, 11.6218, 19.664, 32.1817, 45.7527, 84.5992, 80000}; + + + charge = 1000.; + mip=0; + Float_t defl=0; + + + if (beta<2.) { // it makes no sense to allow beta=5 or so... + + if( L2->GetNTracks(trkAlg)>=1 ){ + PamTrack *TRKtrack = L2->GetTrack(0,trkAlg); + 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[9] = {0, 1, 2, 3 , 4 , 5 , 6, 8, 90}; + Float_t M1[9] = {0, 3.5, 27, 63.0145, 120.504, 173.663, 245.33, 236.517, 0}; + Float_t M2[9] = {0, -1, -10.6, -15.005, -31.0635, -39.4988, -60.5011, -46.3992, 0}; + Float_t M3[9] = {0, 1, 7, 12.5037, 22.8652, 35.2907, 51.4678, 86.4155, 80000}; + + charge = 1000.; + mip=0; + + + if (beta<2.) { // it makes no sense to allow beta=5 or so... + + if( L2->GetNTracks(trkAlg)>=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 =========================== + + + + + // + if ( debug ) this->Print(); if ( debug ) printf(" esci \n"); // };