--- calo/flight/CaloNuclei/src/CaloNuclei.cpp 2007/11/26 08:48:24 1.11 +++ calo/flight/CaloNuclei/src/CaloNuclei.cpp 2008/11/28 16:00:19 1.15 @@ -1,4 +1,7 @@ #include +#include +#include +#include //-------------------------------------- /** @@ -23,6 +26,7 @@ R = 3; // debug = false; + // debug = true; usetrack = true; // }; @@ -41,6 +45,14 @@ 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; @@ -53,19 +65,25 @@ // 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(" 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"); // }; @@ -112,8 +130,9 @@ // if ( debug ) printf(" Always calculate stdedx1 \n"); // - // Always calculate stdedx1 + // Always calculate stdedx1 and maxrel // + Int_t cont=0; Int_t view = 0; Int_t plane = 0; Int_t strip = 0; @@ -142,15 +161,20 @@ // 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,stfpl); + 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]; + // 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<= qmt ){ if ( l < 3 ){ qpremean += qm; @@ -486,12 +514,11 @@ }; // 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((Long_64_t)interplane,qme2,(Long_64_t)ind,work); if ( qm2 >= qmt ){ if ( l < N ){ qpremeanN += qm2; @@ -502,9 +529,49 @@ }; 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; + } + //////////////////////////////////// + ////////////////////////////////// // if ( debug ) printf(" charge is %f \n",sqrt(qpremean)); // @@ -522,12 +589,284 @@ 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 =========================== + + + + + // if ( debug ) this->Print(); if ( debug ) printf(" esci \n");