#include #include #include #include //-------------------------------------- /** * Default constructor */ CaloNuclei::CaloNuclei(){ Clear(); }; CaloNuclei::CaloNuclei(PamLevel2 *l2p){ // Clear(); // L2 = l2p; // if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n"); // OBT = 0; PKT = 0; atime = 0; N = 5; R = 3; // debug = false; // debug = true; usetrack = true; usepl18x = false; // }; 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; // }; void CaloNuclei::Print(){ // 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"); // }; void CaloNuclei::Delete(){ Clear(); //delete this; }; 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"); return; }; // Bool_t newentry = false; // if ( L2->IsORB() ){ 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; }; // if ( !newentry ) return; // tr = ntr; // if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime); // Clear(); // 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); 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.; // // Calculate dedx1 and dedx3 // for ( Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ // mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); // 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; }; }; // }; // dedx3 /= 3.; // Float_t mesethr = dedx1 * 0.90; Float_t mesethr = 0.; if ( dedx1 > 0. ) mesethr = (sqrt(dedx1) - 0.50)*(sqrt(dedx1) - 0.50); Bool_t aldone = false; // retry: // if ( debug ) printf("retry\n"); // interplane = 0; // ethr = TMath::Max(defethr,mesethr); // // Find the interaction plane "interplane" // Int_t gapth = 3; Int_t nhit[2] = {0,0}; Int_t splane[2] = {-1,-1}; Int_t sview[2] = {-1,-1}; Int_t interpl[2] = {-1,-1}; Int_t interv[2] = {-1,-1}; Bool_t wmulthit[2] = {false,false}; Bool_t wgap[2] = {false,false}; Int_t ii = 0; while ( iiGetCaloLevel1()->istrip ){ // mip = L2->GetCaloLevel1()->DecodeEstrip(ii,view,plane,strip); // 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]; }; 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]++; }; }; }; // ii++; // }; // 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++){ // if ( nhit[view] > 1 && !wmulthit[view] && !wgap[view] ){ wmulthit[view] = true; interpl[view] = splane[view]; interv[view] = sview[view]; }; // if ( wmulthit[view] ) multhit = true; if ( wgap[view] ) gap = true; // // convert view and plane number of interaction plane into number of available dE/dx measurements before the interaction plane // if ( interpl[view] >= 0 ) { if ( interv[view] == 0 ){ winterplane[view] = (1 + interpl[view]) * 2; } else { winterplane[view] = (1 + interpl[view]) + (1 + interpl[view] - 1); }; if ( wmulthit[view] ) winterplane[view]--; }; }; if ( winterplane[0] > 0 && winterplane[1] > 0 ){ if ( multhit ){ interplane = TMath::Min(winterplane[0],winterplane[1]); } else { interplane = TMath::Max(winterplane[0],winterplane[1]); }; } else { if ( !winterplane[0] || !winterplane[1] ){ interplane = 0; } else { interplane = TMath::Max(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 ){ // 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; } else { ipl = (1 + plane) + (1 + plane - 1 ); }; if ( ipl > interplane ){ postq += mip; } else { preq += mip; 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; }; 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; }; }; }; }; // }; // ii++; // }; // // here we must calculate qpremean, order vector qme, select 3 lowest measurements and caculate the mean... // 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((Long64_t)interplane,qme,(Long64_t)ind,work); if ( qm >= qmt ){ if ( l < 3 ){ qpremean += qm; RN++; }; l++; if ( debug ) printf(" value no %i qm %f qmt %f \n",l,qm,qmt); }; ind++; }; // 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; } //////////////////////////////////// ////////////////////////////////// // 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); if ( mesethr2 < mesethr*0.90 ){ mesethr = (sqrt(dedx1) - 0.25)*(sqrt(dedx1) - 0.25); } else { mesethr = mesethr2; }; aldone = true; 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[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->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 =========================== //======================================================================= 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->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[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->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"); // };