--- calo/flight/CaloNuclei/src/CaloNuclei.cpp 2007/04/05 13:38:32 1.4 +++ calo/flight/CaloNuclei/src/CaloNuclei.cpp 2007/11/26 08:48:24 1.11 @@ -23,15 +23,20 @@ R = 3; // debug = false; + 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.; @@ -47,19 +52,20 @@ 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("========================================================================\n"); // }; @@ -85,11 +91,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 +110,77 @@ // Clear(); // - PamTrack *track = 0; - track = L2->GetTrack(ntr); + if ( debug ) printf(" Always calculate stdedx1 \n"); + // + // Always calculate stdedx1 // 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,stfpl); + 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]; + }; + } else { + stdedx1 = 0.; + }; + // + if ( debug ) printf(" if ( !usetrack ) return: usetrack %i ntr %i \n",usetrack,ntr); + // + // + // if ( !usetrack ) return; + // + PamTrack *ptrack = 0; + CaloTrkVar *track = 0; + // + if ( usetrack ){ + if ( ntr >= 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.; // @@ -119,21 +190,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; + }; }; // }; @@ -167,32 +259,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 ){ - 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 ( 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]++; }; }; // @@ -273,29 +396,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; + }; + }; + }; }; // }; @@ -315,6 +464,7 @@ Long64_t work[200]; ind = 0; Int_t l = 0; + Int_t RN = 0; Float_t qm = 0.; Float_t qm2 = 0.; // @@ -325,28 +475,36 @@ while ( l < uplim && ind < interplane ){ qm = TMath::KOrdStat(interplane,qme,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); 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; + qpremeanN /= (Float_t)RN; + UN = RN; // if ( debug ) printf(" charge is %f \n",sqrt(qpremean)); // @@ -358,10 +516,20 @@ mesethr = mesethr2; }; aldone = true; - if ( mesethr > defethr ) goto retry; + if ( mesethr > defethr ){ + interplane = 0; + preq = 0.; + postq = 0.; + qpremean = 0.; + qpremeanN = 0.; + multhit = false; + gap = false; + goto retry; + }; }; }; // + if ( debug ) this->Print(); if ( debug ) printf(" esci \n"); // };