#include <AlgAtNuReco.h>
Inheritance diagram for AlgAtNuReco:
Public Member Functions | |
| AlgAtNuReco () | |
| ~AlgAtNuReco () | |
| void | RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx) |
| void | Trace (const char *c) const |
Private Attributes | |
| TObjArray | hitplnbnk [500] |
| TObjArray | xtalkplnbnk [500] |
| TObjArray | clrplnbnk [500] |
| TObjArray | segplnbnk [500] |
| TObjArray | hitbnk |
| TObjArray | clrbnk |
| TObjArray | trksegbnk |
| TObjArray | shwsegbnk |
| TObjArray | trkbnk [2] |
| TObjArray | shwbnk [2] |
| TObjArray | trkcxt |
| TObjArray | shwcxt |
|
|
Definition at line 49 of file AlgAtNuReco.cxx. 00050 {
00051
00052 }
|
|
|
Definition at line 54 of file AlgAtNuReco.cxx. 00055 {
00056
00057 }
|
|
||||||||||||||||
|
Implements AlgBase. Definition at line 59 of file AlgAtNuReco.cxx. References ShowerAtNu::AddAssocTrack(), TrackSegmentAtNu::AddAssocTriToBeg(), TrackSegmentAtNu::AddAssocTriToEnd(), ShowerSegmentAtNu::AddCluster(), CandHandle::AddDaughterLink(), ShowerAtNu::AddHit(), TrackAtNu::AddHit(), ClusterAtNu::AddHit(), ShowerSegmentAtNu::AddSegment(), TrackSegmentAtNu::AddSegment(), TrackSegmentAtNu::AddSegToBeg(), TrackSegmentAtNu::AddSegToEnd(), CandAtNuRecoHandle::AddShower(), CandAtNuRecoHandle::AddTrack(), TrackSegmentAtNu::AddTriToBeg(), TrackSegmentAtNu::AddTriToEnd(), clrbnk, clrplnbnk, ShowerSegmentAtNu::ContainsClr(), AlgFactory::GetAlgHandle(), TrackSegmentAtNu::GetAssocTriBegAt(), TrackSegmentAtNu::GetAssocTriBegLast(), TrackSegmentAtNu::GetAssocTriEndLast(), ObjShowerAtNu::GetBegPlane(), ObjTrackAtNu::GetBegPlane(), ObjShowerAtNu::GetBegStrip(), ObjTrackAtNu::GetBegStrip(), ClusterAtNu::GetBegStrip(), ObjShowerAtNu::GetBegTime(), ObjTrackAtNu::GetBegTime(), CandContext::GetCandRecord(), HitAtNu::GetCharge(), ClusterAtNu::GetCharge(), CandStripHandle::GetCharge(), Registry::GetCharString(), ShowerSegmentAtNu::GetClrAt(), TrackSegmentAtNu::GetClrAt(), ShowerSegmentAtNu::GetClrLast(), TrackSegmentAtNu::GetClrLast(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), CandDeMuxDigitHandle::GetDeMuxDigitFlagWord(), VldContext::GetDetector(), HitAtNu::GetDigits(), ClusterAtNu::GetDigits(), Registry::GetDouble(), ObjTrackAtNu::GetEndDir(), ObjShowerAtNu::GetEndPlane(), ObjTrackAtNu::GetEndPlane(), ObjShowerAtNu::GetEndStrip(), ObjTrackAtNu::GetEndStrip(), ClusterAtNu::GetEndStrip(), ObjShowerAtNu::GetEndTime(), ObjTrackAtNu::GetEndTime(), ObjAtNu::GetHitAt(), ObjAtNu::GetHitLast(), AlgFactory::GetInstance(), Registry::GetInt(), CandContext::GetMom(), CandHandle::GetName(), CandHandle::GetNDaughters(), TrackSegmentAtNu::GetNPlanes(), ObjShowerAtNu::GetPartner(), ObjTrackAtNu::GetPartner(), ClusterAtNu::GetPlane(), HitAtNu::GetPlane(), CandStripHandle::GetPlane(), HitAtNu::GetPlaneView(), ObjAtNu::GetPlaneView(), ObjShowerAtNu::GetReseedFlag(), ObjTrackAtNu::GetReseedFlag(), ShowerSegmentAtNu::GetScore(), TrackSegmentAtNu::GetScore(), TrackSegmentAtNu::GetSegBegAt(), TrackSegmentAtNu::GetSegBegLast(), TrackSegmentAtNu::GetSegEndAt(), TrackSegmentAtNu::GetSegEndLast(), TrackSegmentAtNu::GetSegment(), HitAtNu::GetShwFlag(), ClusterAtNu::GetShwFlag(), ClusterAtNu::GetShwPlnFlag(), HitAtNu::GetStrip(), HitAtNu::GetTime(), TrackSegmentAtNu::GetTmpTrkFlag(), TrackSegmentAtNu::GetTriBegAt(), TrackSegmentAtNu::GetTriBegLast(), TrackSegmentAtNu::GetTriEndAt(), TrackSegmentAtNu::GetTriEndLast(), HitAtNu::GetTrkFlag(), TrackSegmentAtNu::GetTrkFlag(), ClusterAtNu::GetTrkFlag(), ClusterAtNu::GetTrkPlnFlag(), ObjAtNu::GetUID(), HitAtNu::GetUID(), CandHandle::GetVldContext(), ClusterAtNu::GetZPos(), HitAtNu::GetZPos(), hitbnk, hitplnbnk, TrackSegmentAtNu::IsBegAssoc(), HitAtNu::IsDiffuseShwAssoc(), ShowerAtNu::IsDiffuseShwAssoc(), ShowerSegmentAtNu::IsDiffuseShwAssoc(), TrackSegmentAtNu::IsEndAssoc(), ClusterAtNu::IsHitAssoc(), HitAtNu::IsShwAssoc(), ShowerSegmentAtNu::IsShwAssoc(), ClusterAtNu::IsShwAssoc(), ClusterAtNu::IsTrkAssoc(), CandShowerAtNuList::MakeCandidate(), CandTrackAtNuList::MakeCandidate(), MSG, CandRecord::SecureCandHandle(), segplnbnk, CandAtNuRecoHandle::SetAtNuRecoScore(), ShowerAtNu::SetBegVtxFlag(), TrackAtNu::SetBegVtxShower(), ClusterAtNu::SetClrPlnFlag(), CandAtNuRecoHandle::SetCPUTime(), CandShowerAtNuListHandle::SetCPUTime(), CandTrackAtNuListHandle::SetCPUTime(), ShowerAtNu::SetEndVtxFlag(), TrackAtNu::SetEndVtxShower(), CandHandle::SetName(), CandAtNuRecoHandle::SetNBlocks(), TrackSegmentAtNu::SetNPlanes(), ObjShowerAtNu::SetPartner(), ObjTrackAtNu::SetPartner(), ObjShowerAtNu::SetReseedFlag(), ObjTrackAtNu::SetReseedFlag(), TrackSegmentAtNu::SetSegment(), HitAtNu::SetShwFlag(), ClusterAtNu::SetShwFlag(), ClusterAtNu::SetShwPlnFlag(), CandHandle::SetTitle(), TrackSegmentAtNu::SetTmpTrkFlag(), HitAtNu::SetTrkFlag(), TrackSegmentAtNu::SetTrkFlag(), ClusterAtNu::SetTrkFlag(), ClusterAtNu::SetTrkPlnFlag(), ObjAtNu::SetUID(), HitAtNu::SetUID(), HitAtNu::SetXtalkFlag(), shwbnk, shwcxt, shwsegbnk, trkbnk, trkcxt, trksegbnk, and xtalkplnbnk. 00060 {
00061
00062 MSG("AlgAtNuReco", Msg::kDebug) << "AlgAtNuReco::RunAlg(...)" << endl;
00063
00064 // Get CandAtNuRecoHandle
00065 CandAtNuRecoHandle& atnu_handle = dynamic_cast<CandAtNuRecoHandle&>(ch);
00066
00067 // Unpack AlgConfig
00068 Double_t fPeCut,fTimeWin;
00069 Int_t fVtxFlag,fObjMax,fObjMaxFlag;
00070 TString fListOutTrk,fListOutShw;
00071 Int_t NEWOBJCTR,DELOBJCTR;
00072
00073 fListOutTrk = ac.GetCharString("ListOutTrk");
00074 fListOutShw = ac.GetCharString("ListOutShw");
00075 fPeCut = ac.GetDouble("PeCut");
00076 fVtxFlag = ac.GetInt("VtxFlag");
00077 fTimeWin=100.0;
00078 fObjMax=100000;
00079 NEWOBJCTR=0; DELOBJCTR=0;
00080
00081 // Unpack CandContext
00082 CandSliceListHandle* slice_list = (CandSliceListHandle*)(cx.GetDataIn());
00083
00084 // Switch on detector type
00085 Int_t itr;
00086 Int_t modtype,modnum,modpln,modstr;
00087 Int_t fEndPln,fEndStr;
00088 switch(slice_list->GetVldContext()->GetDetector()){
00089 case DetectorType::kFar:
00090 modtype=1; modnum=2; modpln=248; modstr=192; break;
00091 case DetectorType::kNear:
00092 modtype=2; modnum=1; modpln=282; modstr=96; break;
00093 case DetectorType::kCalib:
00094 modtype=3; modnum=1; modpln=60; modstr=24; break;
00095 default:
00096 modtype=-1; modnum=0; modpln=0; break;
00097 }
00098 fEndPln = modnum*(modpln+1); fEndStr = modstr;
00099
00100 MSG("AlgAtNuReco", Msg::kDebug)
00101 << " *** PeCut=" << fPeCut << " VtxFlag=" << fVtxFlag << endl;
00102
00103 MSG("AlgAtNuReco", Msg::kDebug)
00104 << " *** ModType=" << modtype << " ModNum=" << modnum << " "
00105 << " ModPln=" << modpln << " ModStr=" << modstr << endl;
00106
00107 MSG("AlgAtNuReco", Msg::kDebug)
00108 << " *** EndPln=" << fEndPln << " fEndStr=" << fEndStr << endl;
00109
00110
00111 // Clear timebanks
00112 Double_t dtime;
00113 Double_t fTime_TrackAtNu,fTime_ShowerAtNu,fTime_AtNuReco;
00114 fTime_TrackAtNu=0.0; fTime_ShowerAtNu=0.0; fTime_AtNuReco=0.0;
00115
00116 struct timeval tpbefore;
00117 struct timeval tpafter;
00118
00119 struct timeval atnubefore;
00120 struct timeval atnuafter;
00121
00122 gettimeofday(&atnubefore,0);
00123
00124 // Look for tracks and showers
00125 TIter sliceitr(slice_list->GetDaughterIterator());
00126 while(CandSliceHandle* slice = dynamic_cast<CandSliceHandle*>(sliceitr())){
00127 atnu_handle.AddDaughterLink(*slice);
00128 fObjMaxFlag=0;
00129
00130 MSG("AlgAtNuReco", Msg::kDebug)
00131 << " *** SLICE (" << atnu_handle.GetNDaughters() << ") strips= " << slice->GetNDaughters() << endl;
00132
00133
00134 // *************************************
00135 // * F O R M A T I O N O F H I T S *
00136 // *************************************
00137
00138 gettimeofday(&tpbefore,0);
00139
00140 // Sort strips into hit banks
00141 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of hits *** " << endl;
00142 Int_t i,j,k,l;
00143 Bool_t cont;
00144 Double_t chg,scr;
00145 Int_t id,pln,npln,mod,vuw,ctr,xtalk,assoc;
00146 Int_t Nrealhits=0,Nxtalkhits=0;
00147 TObjArray tmp;
00148 TIter stritr(slice->GetDaughterIterator());
00149 while(CandStripHandle* strip = dynamic_cast<CandStripHandle*>(stritr())){
00150 pln=strip->GetPlane(); chg=strip->GetCharge(); xtalk=0;
00151 if(pln>0 && pln<fEndPln && pln<500){
00152 TIter digitr(strip->GetDaughterIterator());
00153 while(CandDeMuxDigitHandle* digit = (CandDeMuxDigitHandle*)(digitr())){
00154 if( ( digit->GetDeMuxDigitFlagWord()<8 )
00155 && (digit->GetDeMuxDigitFlagWord() & CandDeMuxDigit::kXTalk)==(CandDeMuxDigit::kXTalk) ){
00156 xtalk=1;
00157 }
00158 }
00159 HitAtNu* hit = new HitAtNu(strip); hitbnk.Add(hit); NEWOBJCTR++;
00160 if(!xtalk && chg>fPeCut){ hitplnbnk[pln].Add(hit); Nrealhits++; }
00161 else{ hit->SetXtalkFlag(1); xtalkplnbnk[pln].Add(hit); Nxtalkhits++; }
00162 }
00163 atnu_handle.AddDaughterLink(*strip);
00164 }
00165 MSG("AlgAtNuReco", Msg::kDebug) << " *** real hits = " << Nrealhits << " xtalk hits = " << Nxtalkhits << endl;
00166
00167
00168 // ***************************************************
00169 // * F O R M A T I O N O F 1 D C L U S T E R S *
00170 // ***************************************************
00171
00172 // Form clusters
00173 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of clusters *** " << endl;
00174 for(i=1;i<fEndPln;i++){
00175 for(j=0;j<1+hitplnbnk[i].GetLast();j++){
00176 HitAtNu* hitj = (HitAtNu*)(hitplnbnk[i].At(j));
00177 if(hitj->GetUID()==0){
00178 ClusterAtNu* clr = new ClusterAtNu(hitj); NEWOBJCTR++;
00179 clrplnbnk[i].Add(clr); clrbnk.Add(clr);
00180 hitj->SetUID(1); ctr=1;
00181 while(ctr>0){
00182 ctr=0;
00183 for(k=0;k<1+hitplnbnk[i].GetLast();k++){
00184 HitAtNu* hitk = (HitAtNu*)(hitplnbnk[i].At(k));
00185 if(hitk->GetUID()==0 && clr->IsHitAssoc(hitk)){
00186 clr->AddHit(hitk); hitk->SetUID(1); ctr++;
00187 }
00188 }
00189 }
00190 }
00191 }
00192 }
00193
00194 // Identify track+shower clusters
00195 Int_t k0,k1,k2;
00196 Int_t nhits0,nhits1,nclrs0,nclrs1;
00197 TObjArray tmp0,tmp1,tmp2;
00198
00199 for(mod=0;mod<modnum;mod++){
00200 for(pln=0;pln<modpln+1;pln++){
00201 i=pln+mod*(modpln+1);
00202 for(j=0;j<1+clrplnbnk[i].GetLast();j++){
00203 ClusterAtNu* clr0 = (ClusterAtNu*)(clrplnbnk[i].At(j));
00204 if( clr0->GetCharge()>2.0 || clr0->GetDigits()>1.0 ){
00205 clr0->SetTrkFlag(1);
00206 }
00207 }
00208 }
00209 }
00210
00211 for(mod=0;mod<modnum;mod++){
00212 for(pln=1;pln<1+modpln;pln++){
00213 i=pln+mod*(modpln+1);
00214 for(j=0;j<1+clrplnbnk[i].GetLast();j++){
00215 ClusterAtNu* clr0 = (ClusterAtNu*)(clrplnbnk[i].At(j));
00216 nhits0=0; nhits1=0; nclrs0=0; nclrs1=0;
00217 tmp0.Clear(); tmp1.Clear();
00218 for(k=-1;k<2;k++){
00219 k0=i+2*k;
00220 if(k0>mod*(modpln+1) && k0<(mod+1)*(modpln+1)){
00221 for(k1=0;k1<1+clrplnbnk[k0].GetLast();k1++){
00222 ClusterAtNu* clr1 = (ClusterAtNu*)(clrplnbnk[k0].At(k1));
00223 if(clr0->GetCharge()>1.0 && clr1->GetCharge()>1.0){
00224 assoc=clr0->IsShwAssoc(clr1);
00225 if(assoc>0){ nhits0+=1+clr1->GetHitLast(); nclrs0+=1; tmp0.Add(clr1); }
00226 if(assoc>1){ nhits1+=1+clr1->GetHitLast(); nclrs1+=1; tmp1.Add(clr1); }
00227 }
00228 }
00229 }
00230 }
00231 if(nclrs1>1 && nhits1>4){
00232 for(k=0;k<1+tmp1.GetLast();k++){
00233 ClusterAtNu* clrtmp = (ClusterAtNu*)(tmp1.At(k));
00234 clrtmp->SetShwFlag(1);
00235 }
00236 }
00237 if(nclrs0>4 && nhits0>5){
00238 for(k=0;k<1+tmp0.GetLast();k++){
00239 ClusterAtNu* clrtmp = (ClusterAtNu*)(tmp0.At(k));
00240 clrtmp->SetShwFlag(1);
00241 }
00242 }
00243 if(clr0->GetShwFlag()<1 && clr0->GetCharge()>50.0){
00244 clr0->SetShwFlag(1);
00245 }
00246 }
00247 }
00248 }
00249
00250 // Identify track+shower planes
00251 for(mod=0;mod<modnum;mod++){
00252 for(pln=0;pln<modpln+1;pln++){
00253 i=pln+mod*(modpln+1); chg=0.0;
00254 for(j=0;j<1+clrplnbnk[i].GetLast();j++){
00255 ClusterAtNu* clr0 = (ClusterAtNu*)(clrplnbnk[i].At(j));
00256 chg+=clr0->GetCharge();
00257 }
00258 if(chg>0.0){
00259 for(j=0;j<1+clrplnbnk[i].GetLast();j++){
00260 ClusterAtNu* clr1 = (ClusterAtNu*)(clrplnbnk[i].At(j));
00261 if(1+clr1->GetHitLast()<3 && clr1->GetCharge()/chg>0.5) clr1->SetTrkPlnFlag(1);
00262 if(1+clr1->GetHitLast()>1 && clr1->GetCharge()/chg>0.5) clr1->SetShwPlnFlag(1);
00263 if(clr1->GetCharge()/chg>0.5) clr1->SetClrPlnFlag(1);
00264 }
00265 }
00266 }
00267 }
00268
00269 // Print out list of hits + 1D clusters
00270 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF HITS/CLUSTERS *** " << endl;
00271 for(i=0;i<fEndPln;i++){
00272 for(j=0;j<1+hitplnbnk[i].GetLast();j++){
00273 HitAtNu* hit = (HitAtNu*)(hitplnbnk[i].At(j));
00274 MSG("AlgAtNuReco", Msg::kVerbose)
00275 << " HIT:"
00276 << " pln=" << hit->GetPlane()
00277 << " strp=" << hit->GetStrip()
00278 << " chg=" << hit->GetCharge()
00279 << " time=" << hit->GetTime()
00280 << " dig=" << hit->GetDigits() << endl;
00281 }
00282 for(j=0;j<1+clrplnbnk[i].GetLast();j++){
00283 ClusterAtNu* clr = (ClusterAtNu*)(clrplnbnk[i].At(j));
00284 MSG("AlgAtNuReco", Msg::kVerbose)
00285 << " CLUSTER:"
00286 << " plane=" << clr->GetPlane()
00287 << " begstrip=" << clr->GetBegStrip()
00288 << " endstrip=" << clr->GetEndStrip()
00289 << " trkpln=" << clr->GetTrkPlnFlag()
00290 << " shwpln=" << clr->GetShwPlnFlag() << endl;
00291 }
00292 }
00293
00294 gettimeofday(&tpafter,0);
00295 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
00296 fTime_TrackAtNu+=0.5*dtime;
00297 fTime_ShowerAtNu+=0.5*dtime;
00298
00299
00300 // *********************************************
00301 // * F O R M A T I O N O F T R I P L E T S *
00302 // *********************************************
00303
00304 gettimeofday(&tpbefore,0);
00305
00306 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of triplets *** " << endl;
00307
00308 // Form triplets
00309 Int_t jnflag,jnflagtmp;
00310 Int_t km,kp,ktmp,ktmp1;
00311 for(mod=0;mod<modnum;mod++){
00312 for(pln=2;pln<modpln-1;pln++){
00313 i=pln+mod*(modpln+1);
00314 for(k0=0;k0<1+clrplnbnk[i].GetLast();k0++){
00315 ClusterAtNu* clr0 = (ClusterAtNu*)(clrplnbnk[i].At(k0));
00316 if(clr0->GetTrkFlag()>0){
00317 jnflag=0;
00318
00319 // look for triplets of form m 0 p
00320 if( 1+clrplnbnk[i-2].GetLast()>0 && 1+clrplnbnk[i+2].GetLast()>0 ){
00321 for(km=0;km<1+clrplnbnk[i-2].GetLast();km++){
00322 ClusterAtNu* clrm = (ClusterAtNu*)(clrplnbnk[i-2].At(km));
00323 for(kp=0;kp<1+clrplnbnk[i+2].GetLast();kp++){
00324 ClusterAtNu* clrp = (ClusterAtNu*)(clrplnbnk[i+2].At(kp));
00325 if(clrm->GetTrkFlag()>0 && clrp->GetTrkFlag()>0){
00326 assoc=clr0->IsTrkAssoc(clrm,clrp);
00327 if(assoc>0){
00328 TrackSegmentAtNu* seg0 = new TrackSegmentAtNu(clrm,clr0,clrp); NEWOBJCTR++;
00329 segplnbnk[i].Add(seg0); trksegbnk.Add(seg0); jnflag=1;
00330 if(assoc>1){
00331 clr0->SetTrkFlag(2);
00332 }
00333 }
00334 }
00335 }
00336 }
00337 }
00338
00339 // look for triplets of form m <-> 0 p
00340 if( 1+clrplnbnk[i+2].GetLast()>0 ){
00341
00342 // look for triplets of form m X 0 p
00343 if( pln>3 && 1+clrplnbnk[i-4].GetLast()>0 ){
00344 for(kp=0;kp<1+clrplnbnk[i+2].GetLast();kp++){
00345 ClusterAtNu* clrp = (ClusterAtNu*)(clrplnbnk[i+2].At(kp));
00346 for(km=0;km<1+clrplnbnk[i-4].GetLast();km++){
00347 ClusterAtNu* clrm = (ClusterAtNu*)(clrplnbnk[i-4].At(km));
00348 if(clrm->GetTrkFlag()>0 && clrp->GetTrkFlag()>0){
00349 assoc = clr0->IsTrkAssoc(clrm,clrp);
00350 if(assoc>0){
00351 jnflagtmp=0;
00352 if(jnflagtmp==0 && 1+clrplnbnk[i-2].GetLast()>0){
00353 for(ktmp=0;ktmp<1+clrplnbnk[i-2].GetLast();ktmp++){
00354 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i-2].At(ktmp));
00355 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00356 && clrtmp->IsTrkAssoc(clrm,clr0)>0 && clr0->IsTrkAssoc(clrtmp,clrp)>0 ) jnflagtmp=1;
00357 }
00358 }
00359 if(jnflagtmp==0){
00360 for(j=0;j<1;j++){
00361 TrackSegmentAtNu* seg0 = new TrackSegmentAtNu(clrm,clr0,clrp); NEWOBJCTR++;
00362 segplnbnk[i-2*j].Add(seg0); trksegbnk.Add(seg0); jnflag=1;
00363 }
00364 if(assoc>1){
00365 clr0->SetTrkFlag(2);
00366 }
00367 }
00368 }
00369 }
00370 }
00371 }
00372 }
00373
00374 // look for triplets of form m X X 0 p
00375 if( pln>5 && 1+clrplnbnk[i-6].GetLast()>0 ){
00376 for(kp=0;kp<1+clrplnbnk[i+2].GetLast();kp++){
00377 ClusterAtNu* clrp = (ClusterAtNu*)(clrplnbnk[i+2].At(kp));
00378 for(km=0;km<1+clrplnbnk[i-6].GetLast();km++){
00379 ClusterAtNu* clrm = (ClusterAtNu*)(clrplnbnk[i-6].At(km));
00380 if(clrm->GetTrkFlag()>0 && clrp->GetTrkFlag()>0){
00381 assoc = clr0->IsTrkAssoc(clrm,clrp);
00382 if(assoc>0){
00383 jnflagtmp=0;
00384 if(jnflagtmp==0 && 1+clrplnbnk[i-2].GetLast()>0){
00385 for(ktmp=0;ktmp<1+clrplnbnk[i-2].GetLast();ktmp++){
00386 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i-2].At(ktmp));
00387 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00388 && clrtmp->IsTrkAssoc(clrm,clr0)>0 && clr0->IsTrkAssoc(clrtmp,clrp)>0 ) jnflagtmp=1;
00389 }
00390 }
00391 if(jnflagtmp==0 && 1+clrplnbnk[i-4].GetLast()>0){
00392 for(ktmp=0;ktmp<1+clrplnbnk[i-4].GetLast();ktmp++){
00393 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i-4].At(ktmp));
00394 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00395 && clrtmp->IsTrkAssoc(clrm,clr0)>0 && clr0->IsTrkAssoc(clrtmp,clrp)>0 ) jnflagtmp=1;
00396 }
00397 }
00398 if(jnflagtmp==0 && 1+clrplnbnk[i-4].GetLast()>0 && 1+clrplnbnk[i-2].GetLast()>0){
00399 for(ktmp=0;ktmp<1+clrplnbnk[i-4].GetLast();ktmp++){
00400 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i-4].At(ktmp));
00401 for(ktmp1=0;ktmp1<1+clrplnbnk[i-2].GetLast();ktmp1++){
00402 ClusterAtNu* clrtmp1 = (ClusterAtNu*)(clrplnbnk[i-2].At(ktmp1));
00403 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0 && clrtmp1->GetTrkFlag()>0
00404 && clrtmp->IsTrkAssoc(clrm,clrtmp1)>0 && clrtmp1->IsTrkAssoc(clrtmp,clr0)>0
00405 && clr0->IsTrkAssoc(clrtmp1,clrp)>0 ) jnflagtmp=1;
00406 }
00407 }
00408 }
00409 if(jnflagtmp==0){
00410 for(j=0;j<1;j++){
00411 TrackSegmentAtNu* seg0 = new TrackSegmentAtNu(clrm,clr0,clrp); NEWOBJCTR++;
00412 segplnbnk[i-2*j].Add(seg0); trksegbnk.Add(seg0); jnflag=1;
00413 }
00414 if(assoc>1){
00415 clr0->SetTrkFlag(2);
00416 }
00417 }
00418 }
00419 }
00420 }
00421 }
00422 }
00423 }
00424
00425 // look for triplets of form m 0 <-> p
00426 if( 1+clrplnbnk[i-2].GetLast()>0 ){
00427
00428 // look for triplets of form m 0 X p
00429 if( pln<modpln-3 && 1+clrplnbnk[i+4].GetLast()>0 ){
00430 for(km=0;km<1+clrplnbnk[i-2].GetLast();km++){
00431 ClusterAtNu* clrm = (ClusterAtNu*)(clrplnbnk[i-2].At(km));
00432 for(kp=0;kp<1+clrplnbnk[i+4].GetLast();kp++){
00433 ClusterAtNu* clrp = (ClusterAtNu*)(clrplnbnk[i+4].At(kp));
00434 if(clrm->GetTrkFlag()>0 && clrp->GetTrkFlag()>0){
00435 assoc = clr0->IsTrkAssoc(clrm,clrp);
00436 if(assoc>0){
00437 jnflagtmp=0;
00438 if(jnflagtmp==0 && 1+clrplnbnk[i+2].GetLast()>0){
00439 for(ktmp=0;ktmp<1+clrplnbnk[i+2].GetLast();ktmp++){
00440 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i+2].At(ktmp));
00441 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00442 && clrtmp->IsTrkAssoc(clr0,clrp)>0 && clr0->IsTrkAssoc(clrm,clrtmp)>0 ) jnflagtmp=1;
00443 }
00444 }
00445 if(jnflagtmp==0){
00446 for(j=0;j<2;j++){
00447 TrackSegmentAtNu* seg0 = new TrackSegmentAtNu(clrm,clr0,clrp); NEWOBJCTR++;
00448 segplnbnk[i+2*j].Add(seg0); trksegbnk.Add(seg0); jnflag=1;
00449 }
00450 if(assoc>1){
00451 clr0->SetTrkFlag(2);
00452 }
00453 }
00454 }
00455 }
00456 }
00457 }
00458 }
00459
00460 // look for triplets of form m 0 X X p
00461 if( pln<modpln-5 && 1+clrplnbnk[i+6].GetLast()>0 ){
00462 for(km=0;km<1+clrplnbnk[i-2].GetLast();km++){
00463 ClusterAtNu* clrm = (ClusterAtNu*)(clrplnbnk[i-2].At(km));
00464 for(kp=0;kp<1+clrplnbnk[i+6].GetLast();kp++){
00465 ClusterAtNu* clrp = (ClusterAtNu*)(clrplnbnk[i+6].At(kp));
00466 if(clrm->GetTrkFlag()>0 && clrp->GetTrkFlag()>0){
00467 assoc = clr0->IsTrkAssoc(clrm,clrp);
00468 if(assoc>0){
00469 jnflagtmp=0;
00470 if(jnflagtmp==0 && 1+clrplnbnk[i+2].GetLast()>0){
00471 for(ktmp=0;ktmp<1+clrplnbnk[i+2].GetLast();ktmp++){
00472 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i+2].At(ktmp));
00473 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00474 && clrtmp->IsTrkAssoc(clr0,clrp)>0 && clr0->IsTrkAssoc(clrm,clrtmp)>0 ) jnflagtmp=1;
00475 }
00476 }
00477 if(jnflagtmp==0 && 1+clrplnbnk[i+4].GetLast()>0){
00478 for(ktmp=0;ktmp<1+clrplnbnk[i+4].GetLast();ktmp++){
00479 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i+4].At(ktmp));
00480 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00481 && clrtmp->IsTrkAssoc(clr0,clrp)>0 && clr0->IsTrkAssoc(clrm,clrtmp)>0 ) jnflagtmp=1;
00482 }
00483 }
00484 if(jnflagtmp==0 && 1+clrplnbnk[i+2].GetLast()>0 && 1+clrplnbnk[i+4].GetLast()>0){
00485 for(ktmp=0;ktmp<1+clrplnbnk[i+2].GetLast();ktmp++){
00486 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i+2].At(ktmp));
00487 for(ktmp1=0;ktmp1<1+clrplnbnk[i+4].GetLast();ktmp1++){
00488 ClusterAtNu* clrtmp1 = (ClusterAtNu*)(clrplnbnk[i+4].At(ktmp1));
00489 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0 && clrtmp1->GetTrkFlag()>0
00490 && clr0->IsTrkAssoc(clrm,clrtmp)>0 && clrtmp->IsTrkAssoc(clr0,clrtmp1)>0
00491 && clrtmp1->IsTrkAssoc(clrtmp,clrp) ) jnflagtmp=1;
00492 }
00493 }
00494 }
00495 if(jnflagtmp==0){
00496 for(j=0;j<3;j++){
00497 TrackSegmentAtNu* seg0 = new TrackSegmentAtNu(clrm,clr0,clrp); NEWOBJCTR++;
00498 segplnbnk[i+2*j].Add(seg0); trksegbnk.Add(seg0); jnflag=1;
00499 }
00500 if(assoc>1){
00501 clr0->SetTrkFlag(2);
00502 }
00503 }
00504 }
00505 }
00506 }
00507 }
00508 }
00509 }
00510
00511 // look for triplets of form m X 0 X p
00512 if( pln>3 && 1+clrplnbnk[i-4].GetLast()>0
00513 && pln<modpln-3 && 1+clrplnbnk[i+4].GetLast()>0 ){
00514 for(km=0;km<1+clrplnbnk[i-4].GetLast();km++){
00515 ClusterAtNu* clrm = (ClusterAtNu*)(clrplnbnk[i-4].At(km));
00516 for(kp=0;kp<1+clrplnbnk[i+4].GetLast();kp++){
00517 ClusterAtNu* clrp = (ClusterAtNu*)(clrplnbnk[i+4].At(kp));
00518 if(clrm->GetTrkFlag()>0 && clrp->GetTrkFlag()>0){
00519 assoc=clr0->IsTrkAssoc(clrm,clrp);
00520 if(assoc>0){
00521 jnflagtmp=0;
00522 if(jnflagtmp==0 && 1+clrplnbnk[i-2].GetLast()>0){
00523 for(ktmp=0;ktmp<1+clrplnbnk[i-2].GetLast();ktmp++){
00524 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i-2].At(ktmp));
00525 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00526 && clrtmp->IsTrkAssoc(clrm,clr0)>0 && clr0->IsTrkAssoc(clrtmp,clrp)>0 ) jnflagtmp=1;
00527 }
00528 }
00529 if(jnflagtmp==0 && 1+clrplnbnk[i+2].GetLast()>0){
00530 for(ktmp=0;ktmp<1+clrplnbnk[i+2].GetLast();ktmp++){
00531 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i+2].At(ktmp));
00532 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0
00533 && clrtmp->IsTrkAssoc(clr0,clrp)>0 && clr0->IsTrkAssoc(clrm,clrtmp)>0 ) jnflagtmp=1;
00534 }
00535 }
00536 if(jnflagtmp==0 && 1+clrplnbnk[i-2].GetLast()>0 && 1+clrplnbnk[i+2].GetLast()>0){
00537 for(ktmp=0;ktmp<1+clrplnbnk[i-2].GetLast();ktmp++){
00538 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrplnbnk[i-2].At(ktmp));
00539 for(ktmp1=0;ktmp1<1+clrplnbnk[i+2].GetLast();ktmp1++){
00540 ClusterAtNu* clrtmp1 = (ClusterAtNu*)(clrplnbnk[i+2].At(ktmp1));
00541 if( jnflagtmp==0 && clrtmp->GetTrkFlag()>0 && clrtmp1->GetTrkFlag()>0
00542 && clrtmp->IsTrkAssoc(clrm,clr0)>0 && clr0->IsTrkAssoc(clrtmp,clrtmp1)>0
00543 && clrtmp1->IsTrkAssoc(clr0,clrp) ) jnflagtmp=1;
00544 }
00545 }
00546 }
00547 if(jnflagtmp==0){
00548 for(j=0;j<2;j++){
00549 TrackSegmentAtNu* seg0 = new TrackSegmentAtNu(clrm,clr0,clrp); NEWOBJCTR++;
00550 segplnbnk[i+2*j].Add(seg0); trksegbnk.Add(seg0); jnflag=1;
00551 }
00552 if(assoc>1){
00553 clr0->SetTrkFlag(2);
00554 }
00555 }
00556 }
00557 }
00558 }
00559 }
00560 }
00561
00562 if( jnflag==1
00563 && 1+clr0->GetEndStrip()-clr0->GetBegStrip()<2 ){
00564 clr0->SetTrkFlag(2);
00565 }
00566
00567 }
00568 }
00569 }
00570 }
00571
00572 // Print out list of triplets
00573 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF TRIPLETS *** " << endl;
00574 for(i=0;i<fEndPln;i++){
00575 for(j=0;j<1+segplnbnk[i].GetLast();j++){
00576 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(segplnbnk[i].At(j));
00577 ClusterAtNu* clr0 = (ClusterAtNu*)(seg->GetClrAt(0));
00578 ClusterAtNu* clr1 = (ClusterAtNu*)(seg->GetClrAt(1));
00579 ClusterAtNu* clr2 = (ClusterAtNu*)(seg->GetClrAt(2));
00580 MSG("AlgAtNuReco", Msg::kVerbose)
00581 << " (" << clr0->GetPlane() << "," << clr0->GetBegStrip() << "->" << clr0->GetEndStrip() << ") "
00582 << " (" << clr1->GetPlane() << "," << clr1->GetBegStrip() << "->" << clr1->GetEndStrip() << ") "
00583 << " (" << clr2->GetPlane() << "," << clr2->GetBegStrip() << "->" << clr2->GetEndStrip() << ") " << endl;
00584 }
00585 }
00586
00587
00588 // *********************************************************
00589 // * F O R M A T I O N O F T R A C K S E G M E N T S *
00590 // *********************************************************
00591
00592 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of track segments *** " << endl;
00593 for(mod=0;mod<modnum;mod++){
00594 for(pln=3;pln<modpln-1;pln++){
00595 i=pln+mod*(modpln+1);
00596 if(1+segplnbnk[i].GetLast()>0 && 1+segplnbnk[i+2].GetLast()>0){
00597 for(k0=0;k0<1+segplnbnk[i].GetLast();k0++){
00598 TrackSegmentAtNu* seg0 = (TrackSegmentAtNu*)(segplnbnk[i].At(k0));
00599 for(kp=0;kp<1+segplnbnk[i+2].GetLast();kp++){
00600 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(segplnbnk[i+2].At(kp));
00601 if( seg0->IsEndAssoc(segp)>0 || segp->IsBegAssoc(seg0)>0 ){
00602 seg0->AddTriToEnd(segp); segp->AddTriToBeg(seg0);
00603 }
00604 }
00605 }
00606 }
00607 }
00608 }
00609
00610
00611 // Select preferred associations between triplets
00612 Int_t jnflagm,jnflagp;
00613 TObjArray tmpm,tmpp;
00614 for(mod=0;mod<modnum;mod++){
00615 for(pln=3;pln<modpln-1;pln++){
00616 i=pln+mod*(modpln+1);
00617
00618 // previous adjacent triplets ( i-2 <-> ) i <-> i+2
00619 tmpm.Clear();
00620 for(j=0;j<1+segplnbnk[i].GetLast();j++){
00621 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(segplnbnk[i].At(j));
00622 if(1+segm->GetTriEndLast()>0){
00623 tmpm.Add(segm);
00624 if(1+segm->GetTriBegLast()>0){
00625 segm->SetTmpTrkFlag(1);
00626 for(k=0;k<1+segm->GetTriBegLast();k++){
00627 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segm->GetTriBegAt(k));
00628 if( segm->GetTmpTrkFlag()<2 && segtmp->GetBegPlane()<segm->GetBegPlane() ){
00629 segm->SetTmpTrkFlag(2);
00630 }
00631 }
00632 }
00633 }
00634 }
00635
00636 // subsequent adjacent triplets i <-> i+2 ( <-> i+4 )
00637 tmpp.Clear();
00638 for(j=0;j<1+segplnbnk[i+2].GetLast();j++){
00639 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(segplnbnk[i+2].At(j));
00640 if(1+segp->GetTriBegLast()>0){
00641 tmpp.Add(segp);
00642 if(1+segp->GetTriEndLast()>0){
00643 segp->SetTmpTrkFlag(1);
00644 for(k=0;k<1+segp->GetTriEndLast();k++){
00645 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segp->GetTriEndAt(k));
00646 if( segp->GetTmpTrkFlag()<2 && segtmp->GetEndPlane()>segp->GetEndPlane() ){
00647 segp->SetTmpTrkFlag(2);
00648 }
00649 }
00650 }
00651 }
00652 }
00653
00654 // look for preferred joins ( i-2 <-> ) i <-> i+2
00655 for(j=0;j<1+tmpp.GetLast();j++){
00656 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(tmpp.At(j));
00657 jnflag=0;
00658 for(k=0;k<1+segp->GetTriBegLast();k++){
00659 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segp->GetTriBegAt(k));
00660 if(segtmp->GetTmpTrkFlag()>1) jnflag=1;
00661 }
00662 if(jnflag>0){
00663 for(k=0;k<1+segp->GetTriBegLast();k++){
00664 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segp->GetTriBegAt(k));
00665 if(segtmp->GetTmpTrkFlag()>0 && segtmp->GetTmpTrkFlag()<2) segtmp->SetTmpTrkFlag(2);
00666 }
00667 }
00668 }
00669
00670 // look for preferred joins i <-> i+2 ( <-> i+4 )
00671 for(j=0;j<1+tmpm.GetLast();j++){
00672 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(tmpm.At(j));
00673 jnflag=0;
00674 for(k=0;k<1+segm->GetTriEndLast();k++){
00675 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segm->GetTriEndAt(k));
00676 if(segtmp->GetTmpTrkFlag()>1) jnflag=1;
00677 }
00678 if(jnflag>0){
00679 for(k=0;k<1+segm->GetTriEndLast();k++){
00680 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segm->GetTriEndAt(k));
00681 if(segtmp->GetTmpTrkFlag()>0 && segtmp->GetTmpTrkFlag()<2) segtmp->SetTmpTrkFlag(2);
00682 }
00683 }
00684 }
00685
00686 // make joins
00687 if(1+tmpm.GetLast()>0 && 1+tmpp.GetLast()>0){
00688
00689 // look for doubly preferred joins ( i-2 <-> ) i <-> i+2 ( <-> i+4 )
00690 for(j=0;j<1+tmpm.GetLast();j++){
00691 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(tmpm.At(j));
00692 for(k=0;k<1+segm->GetTriEndLast();k++){
00693 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(segm->GetTriEndAt(k));
00694 if(segm->GetTmpTrkFlag()>1 && segp->GetTmpTrkFlag()>1){
00695 jnflagm=0;
00696 for(l=0;l<1+segm->GetTriBegLast();l++){
00697 TrackSegmentAtNu* segm2 = (TrackSegmentAtNu*)(segm->GetTriBegAt(l));
00698 if( jnflagm<1 && (segm2->IsEndAssoc(segp)||segp->IsBegAssoc(segm2)) ) jnflagm=1;
00699 }
00700 jnflagp=0;
00701 for(l=0;l<1+segp->GetTriEndLast();l++){
00702 TrackSegmentAtNu* segp2 = (TrackSegmentAtNu*)(segp->GetTriEndAt(l));
00703 if( jnflagp<1 && (segm->IsEndAssoc(segp2)||segp2->IsBegAssoc(segm)) ) jnflagp=1;
00704 }
00705 if(jnflagm>0 && jnflagp>0){
00706 segm->AddAssocTriToEnd(segp); segp->AddAssocTriToBeg(segm);
00707 }
00708 }
00709 }
00710 }
00711
00712 // look for singly preferred joins (1) ( i-2 <-> ) i <-> i+2
00713 for(j=0;j<1+tmpm.GetLast();j++){
00714 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(tmpm.At(j));
00715 if(segm->GetTmpTrkFlag()>1){
00716 jnflag=0;
00717 for(k=0;k<1+segm->GetTriEndLast();k++){
00718 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(segm->GetTriEndAt(k));
00719 if(segp->GetTmpTrkFlag()>1) jnflag=1;
00720 }
00721 if(jnflag<1){
00722 for(k=0;k<1+segm->GetTriEndLast();k++){
00723 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(segm->GetTriEndAt(k));
00724 jnflagm=0;
00725 for(l=0;l<1+segm->GetTriBegLast();l++){
00726 TrackSegmentAtNu* segm2 = (TrackSegmentAtNu*)(segm->GetTriBegAt(l));
00727 if( jnflagm<1 && (segm2->IsEndAssoc(segp)||segp->IsBegAssoc(segm2)) ) jnflagm=1;
00728 }
00729 if(jnflagm>0){
00730 segm->AddAssocTriToEnd(segp); segp->AddAssocTriToBeg(segm);
00731 }
00732 }
00733 }
00734 }
00735 }
00736
00737 // look for singly preferred joins (2) i <-> i+2 ( <-> i+4 )
00738 for(j=0;j<1+tmpp.GetLast();j++){
00739 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(tmpp.At(j));
00740 if(segp->GetTmpTrkFlag()>1){
00741 jnflag=0;
00742 for(k=0;k<1+segp->GetTriBegLast();k++){
00743 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(segp->GetTriBegAt(k));
00744 if(segm->GetTmpTrkFlag()>1) jnflag=1;
00745 }
00746 if(jnflag<1){
00747 for(k=0;k<1+segp->GetTriBegLast();k++){
00748 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(segp->GetTriBegAt(k));
00749 jnflagp=0;
00750 for(l=0;l<1+segp->GetTriEndLast();l++){
00751 TrackSegmentAtNu* segp2 = (TrackSegmentAtNu*)(segp->GetTriEndAt(l));
00752 if( jnflagp<1 && (segm->IsEndAssoc(segp2)||segp2->IsBegAssoc(segm)) ) jnflagp=1;
00753 }
00754 if(jnflagp>0){
00755 segm->AddAssocTriToEnd(segp); segp->AddAssocTriToBeg(segm);
00756 }
00757 }
00758 }
00759 }
00760 }
00761
00762 // look for other joins i <-> i+2
00763 for(j=0;j<1+tmpm.GetLast();j++){
00764 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(tmpm.At(j));
00765 for(k=0;k<1+segm->GetTriEndLast();k++){
00766 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(segm->GetTriEndAt(k));
00767 if(segm->GetTmpTrkFlag()<2 && segp->GetTmpTrkFlag()<2){
00768 jnflagm=0;
00769 for(l=0;l<1+segm->GetTriEndLast();l++){
00770 TrackSegmentAtNu* segmp = (TrackSegmentAtNu*)(segm->GetTriEndAt(l));
00771 if(segmp->GetTmpTrkFlag()>1) jnflagm=1;
00772 }
00773 jnflagp=0;
00774 for(l=0;l<1+segp->GetTriBegLast();l++){
00775 TrackSegmentAtNu* segpm = (TrackSegmentAtNu*)(segp->GetTriBegAt(l));
00776 if(segpm->GetTmpTrkFlag()>1) jnflagp=1;
00777 }
00778 if(jnflagm<1 && jnflagp<1){
00779 segm->AddAssocTriToEnd(segp); segp->AddAssocTriToBeg(segm);
00780 }
00781 }
00782 }
00783 }
00784
00785 }
00786
00787 // clear flags
00788 for(j=0;j<1+tmpm.GetLast();j++){
00789 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(tmpm.At(j));
00790 segm->SetTmpTrkFlag(0);
00791 }
00792
00793 for(j=0;j<1+tmpp.GetLast();j++){
00794 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(tmpp.At(j));
00795 segp->SetTmpTrkFlag(0);
00796 }
00797
00798 }
00799 }
00800
00801
00802 // Form segments from triplets
00803 TObjArray segbnk[2];
00804 for(mod=0;mod<modnum;mod++){
00805 for(pln=3;pln<modpln-1;pln++){
00806 i=pln+mod*(modpln+1);
00807
00808 for(j=0;j<1+segplnbnk[i].GetLast();j++){
00809 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(segplnbnk[i].At(j));
00810
00811 // isolated segments
00812 if(1+seg->GetTriBegLast()==0 && 1+seg->GetTriEndLast()==0){
00813 vuw=seg->GetPlaneView(); segbnk[vuw].Add(seg); seg->SetSegment(seg);
00814 }
00815
00816 // associated segments
00817 if(1+seg->GetAssocTriBegLast()>0 || 1+seg->GetAssocTriEndLast()>0){
00818 jnflag=0;
00819 if(1+seg->GetAssocTriBegLast()==1){
00820 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(seg->GetAssocTriBegAt(0));
00821 if(1+segtmp->GetAssocTriEndLast()==1){
00822 TrackSegmentAtNu* seg0 = (TrackSegmentAtNu*)(segtmp->GetSegment());
00823 seg0->AddSegment(seg); seg->SetSegment(seg0); seg->SetUID(-1);
00824 jnflag=1;
00825 }
00826 }
00827 if(jnflag<1){
00828 vuw=seg->GetPlaneView(); segbnk[vuw].Add(seg); seg->SetSegment(seg);
00829 for(k=0;k<1+seg->GetAssocTriBegLast();k++){
00830 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(seg->GetAssocTriBegAt(k));
00831 TrackSegmentAtNu* segtmp0 = (TrackSegmentAtNu*)(segtmp->GetSegment());
00832 seg->AddSegToBeg(segtmp0); segtmp0->AddSegToEnd(seg);
00833 }
00834 }
00835 }
00836
00837 }
00838 }
00839 }
00840
00841
00842 // Form associations between non-adjacent segments
00843 Int_t npts,npts0=4;
00844 Double_t inv_npts;
00845 for(vuw=0;vuw<2;vuw++){
00846 tmp1.Clear(); tmp2.Clear();
00847 for(i=0;i<1+segbnk[vuw].GetLast();i++){
00848 TrackSegmentAtNu* seg1 = (TrackSegmentAtNu*)(segbnk[vuw].At(i));
00849 if(1+seg1->GetSegEndLast()==0){
00850 jnflag=0; mod=(Int_t)((seg1->GetEndPlane())/(modpln+1)); npts=npts0;
00851 if( (modtype==1 && seg1->GetEndStrip()>80 && seg1->GetEndStrip()<=110)
00852 || (1+seg1->GetEndPlane()-seg1->GetBegPlane()/2>5) ){
00853 inv_npts = seg1->GetEndDir()/2.0;
00854 if(inv_npts<0) inv_npts=-inv_npts; if(inv_npts<0.01) inv_npts=0.01;
00855 npts=(Int_t)(1.0/inv_npts); if(npts<8) npts=8; if(npts>18) npts=18;
00856 }
00857 for(j=1;j<npts;j++){
00858 pln=seg1->GetEndPlane()+2*j;
00859 if(jnflag==0 && pln<(modpln+1)*(mod+1)){
00860 for(k=0;k<1+segplnbnk[pln].GetLast();k++){
00861 TrackSegmentAtNu* seg2 = (TrackSegmentAtNu*)(segplnbnk[pln].At(k));
00862 if( ( seg2->GetSegment()==seg2 && 1+seg2->GetSegBegLast()==0 )
00863 && ( seg2->GetBegPlane()-seg1->GetEndPlane()>0
00864 || seg1->GetClrLast()>5 && seg2->GetClrLast()>5 )
00865 && ( j<npts0
00866 || ( modtype==1
00867 && ( ( seg2->GetBegStrip()>80 && seg2->GetBegStrip()<=100
00868 && seg1->GetEndStrip()>90 && seg1->GetEndStrip()<=110 )
00869 || ( seg2->GetBegStrip()>90 && seg2->GetBegStrip()<=110
00870 && seg1->GetEndStrip()>80 && seg1->GetEndStrip()<=100 ) ) )
00871 || ( 1+seg1->GetEndPlane()-seg1->GetBegPlane()/2>5 ) ) ){
00872 if( seg1->IsEndAssoc(seg2) ){
00873 tmp1.Add(seg1); tmp2.Add(seg2); jnflag=1;
00874 MSG("AlgAtNuReco", Msg::kVerbose) << " NON-ADJ ASSOC: "
00875 << seg1->GetBegPlane() << " " << seg1->GetEndPlane() << " "
00876 << seg2->GetBegPlane() << " " << seg2->GetEndPlane() << endl;
00877 }
00878 }
00879 }
00880 }
00881 }
00882 }
00883 }
00884 for(i=0;i<1+tmp1.GetLast();i++){
00885 TrackSegmentAtNu* seg1 = (TrackSegmentAtNu*)(tmp1.At(i));
00886 TrackSegmentAtNu* seg2 = (TrackSegmentAtNu*)(tmp2.At(i));
00887 seg1->AddSegToEnd(seg2); seg2->AddSegToBeg(seg1);
00888 }
00889 }
00890
00891 // Print out list of segments
00892 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF TRACK SEGMENTS *** " << endl;
00893 for(vuw=0;vuw<2;vuw++){
00894 MSG("AlgAtNuReco", Msg::kVerbose) << "view: " << vuw << endl;
00895 for(i=0;i<1+segbnk[vuw].GetLast();i++){
00896 TrackSegmentAtNu* segment = (TrackSegmentAtNu*)(segbnk[vuw].At(i));
00897 MSG("AlgAtNuReco", Msg::kVerbose) << " i=" << i
00898 << " begpln: " << segment->GetBegPlane()
00899 << " begstr: " << segment->GetBegStrip()
00900 << " endpln: " << segment->GetEndPlane()
00901 << " endstr: " << segment->GetEndStrip() << endl;
00902 MSG("AlgAtNuReco", Msg::kVerbose) << " beg: " << endl;
00903 for(j=0;j<1+segment->GetSegBegLast();j++){
00904 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segment->GetSegBegAt(j));
00905 MSG("AlgAtNuReco", Msg::kVerbose) << " j=" << j
00906 << " begpln: " << segtmp->GetBegPlane()
00907 << " begstr: " << segtmp->GetBegStrip()
00908 << " endpln: " << segtmp->GetEndPlane()
00909 << " endstr: " << segtmp->GetEndStrip() << endl;
00910 }
00911 MSG("AlgAtNuReco", Msg::kVerbose) << " end: " << endl;
00912 for(j=0;j<1+segment->GetSegEndLast();j++){
00913 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(segment->GetSegEndAt(j));
00914 MSG("AlgAtNuReco", Msg::kVerbose) << " j=" << j
00915 << " begpln: " << segtmp->GetBegPlane()
00916 << " begstr: " << segtmp->GetBegStrip()
00917 << " endpln: " << segtmp->GetEndPlane()
00918 << " endstr: " << segtmp->GetEndStrip() << endl;
00919 }
00920 }
00921 }
00922
00923 gettimeofday(&tpafter,0);
00924 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
00925 fTime_TrackAtNu+=dtime;
00926
00927
00928 // ***********************************************************
00929 // * F O R M A T I O N O F S H O W E R S E G M E N T S *
00930 // ***********************************************************
00931
00932 gettimeofday(&tpbefore,0);
00933
00934 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of shower segments *** " << endl;
00935
00936 // Join associated clusters
00937 Int_t pln0;
00938 TObjArray segbnk2[2];
00939 for(mod=0;mod<modnum;mod++){
00940 for(pln=1;pln<1+modpln;pln++){
00941 i=pln+mod*(modpln+1);
00942 for(j=0;j<1+clrplnbnk[i].GetLast();j++){
00943 ClusterAtNu* clr0 = (ClusterAtNu*)(clrplnbnk[i].At(j));
00944 if(clr0->GetShwFlag()>0 && clr0->GetShwFlag()<2){
00945 ShowerSegmentAtNu* seg0 = new ShowerSegmentAtNu(clr0); NEWOBJCTR++;
00946 vuw = seg0->GetPlaneView(); segbnk2[vuw].Add(seg0);
00947 clr0->SetShwFlag(2); shwsegbnk.Add(seg0); ctr=1;
00948 while(ctr>0){
00949 ctr=0;
00950 pln0 = -2 + seg0->GetBegPlane();
00951 npln = 3 + (seg0->GetEndPlane()-seg0->GetBegPlane())/2;
00952 for(k=0;k<npln;k++){
00953 k0=pln0+2*k;
00954 if(k0>mod*(modpln+1) && k0<(mod+1)*(modpln+1)){
00955 for(k1=0;k1<1+clrplnbnk[k0].GetLast();k1++){
00956 ClusterAtNu* clr1 = (ClusterAtNu*)(clrplnbnk[k0].At(k1));
00957 if(clr1->GetShwFlag()>0 && clr1->GetShwFlag()<2){
00958 assoc=seg0->IsShwAssoc(clr1);
00959 if(assoc>0){
00960 seg0->AddCluster(clr1);
00961 clr1->SetShwFlag(2); ctr++;
00962 }
00963 }
00964 }
00965 }
00966 }
00967 }
00968 }
00969 }
00970 }
00971 }
00972
00973 // Print out list of shower segments
00974 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF SHOWER SEGMENTS *** " << endl;
00975 for(vuw=0;vuw<2;vuw++){
00976 for(i=0;i<1+segbnk2[vuw].GetLast();i++){
00977 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(segbnk2[vuw].At(i));
00978 MSG("AlgAtNuReco", Msg::kVerbose)
00979 << "vuw=" << vuw << " i=" << i << " "
00980 << seg->GetBegPlane() << " " << seg->GetEndPlane() << " "
00981 << seg->GetBegStrip() << " " << seg->GetEndStrip() << " "
00982 << 1+seg->GetHitLast() << endl;
00983 }
00984 }
00985
00986 gettimeofday(&tpafter,0);
00987 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
00988 fTime_ShowerAtNu+=dtime;
00989
00990
00991 // ******************************************
00992 // * C O M P A R E U / V V I E W S (1) *
00993 // ******************************************
00994
00995 Int_t shwctr,trkctr,plnctr;
00996 TObjArray mysegbnk[2],mysegbnk2[2];
00997 MSG("AlgAtNuReco", Msg::kDebug) << " *** comparing views (1) *** " << endl;
00998
00999 // determine track-like track segments
01000 gettimeofday(&tpbefore,0);
01001
01002 for(vuw=0;vuw<2;vuw++){
01003 for(i=0;i<1+segbnk[vuw].GetLast();i++){
01004 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(segbnk[vuw].At(i));
01005 shwctr=0; trkctr=0; plnctr=0;
01006 for(j=0;j<1+seg->GetClrLast();j++){
01007 ClusterAtNu* clr = (ClusterAtNu*)(seg->GetClrAt(j));
01008 if(clr->GetTrkFlag()>1){ trkctr++; if(clr->GetShwFlag()<2) shwctr++; }
01009 if(clr->GetTrkPlnFlag()>0){ plnctr++; }
01010 }
01011 if( plnctr>1 ){
01012 if(shwctr>0) seg->SetUID(1); if(shwctr>2) seg->SetUID(2);
01013 }
01014 }
01015 }
01016
01017 for(i=0;i<1+segbnk[0].GetLast();i++){
01018 TrackSegmentAtNu* segu = (TrackSegmentAtNu*)(segbnk[0].At(i));
01019 for(j=0;j<1+segbnk[1].GetLast();j++){
01020 TrackSegmentAtNu* segv = (TrackSegmentAtNu*)(segbnk[1].At(j));
01021 if( segu->GetBegPlane()+1<segv->GetEndPlane() && segu->GetEndPlane()>segv->GetBegPlane()+1
01022 && segu->GetUID()>0 && segv->GetUID()>0 ){
01023 if(segv->GetUID()>1) mysegbnk[0].Add(segu); if(segu->GetUID()>1) mysegbnk[1].Add(segv);
01024 }
01025 }
01026 }
01027
01028 for(vuw=0;vuw<2;vuw++){
01029 for(i=0;i<1+mysegbnk[vuw].GetLast();i++){
01030 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(mysegbnk[vuw].At(i));
01031 seg->SetUID(2);
01032 }
01033 }
01034
01035 for(vuw=0;vuw<2;vuw++){
01036 for(i=0;i<1+segbnk[vuw].GetLast();i++){
01037 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(segbnk[vuw].At(i));
01038 if( seg->GetUID()>0 && seg->GetUID()<2 ) seg->SetReseedFlag(1);
01039 }
01040 }
01041
01042 gettimeofday(&tpafter,0);
01043 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
01044 fTime_TrackAtNu+=dtime;
01045
01046 // determine shower-like shower segments
01047 gettimeofday(&tpbefore,0);
01048
01049 for(vuw=0;vuw<2;vuw++){
01050 for(i=0;i<1+segbnk2[vuw].GetLast();i++){
01051 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(segbnk2[vuw].At(i));
01052 trkctr=0; shwctr=0; plnctr=0;
01053 for(j=0;j<1+seg->GetClrLast();j++){
01054 ClusterAtNu* clr = (ClusterAtNu*)(seg->GetClrAt(j));
01055 if(clr->GetShwFlag()>1){ shwctr++; if(clr->GetTrkFlag()<2) trkctr++; }
01056 if(clr->GetShwPlnFlag()>0){ plnctr++; }
01057 }
01058 if( plnctr>-1 ){
01059 if(trkctr>-1) seg->SetUID(1); if(trkctr>0) seg->SetUID(2);
01060 }
01061 }
01062 }
01063
01064 for(i=0;i<1+segbnk2[0].GetLast();i++){
01065 ShowerSegmentAtNu* segu = (ShowerSegmentAtNu*)(segbnk2[0].At(i));
01066 for(j=0;j<1+segbnk2[1].GetLast();j++){
01067 ShowerSegmentAtNu* segv = (ShowerSegmentAtNu*)(segbnk2[1].At(j));
01068 if( segu->GetBegPlane()-2<segv->GetEndPlane() && segu->GetEndPlane()>segv->GetBegPlane()-2
01069 && segu->GetUID()>0 && segv->GetUID()>0 ){
01070 if(segv->GetUID()>1) mysegbnk2[0].Add(segu); if(segu->GetUID()>1) mysegbnk2[1].Add(segv);
01071 }
01072 }
01073 }
01074
01075 for(vuw=0;vuw<2;vuw++){
01076 for(i=0;i<1+mysegbnk2[vuw].GetLast();i++){
01077 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(mysegbnk2[vuw].At(i));
01078 seg->SetUID(2);
01079 }
01080 }
01081
01082 for(vuw=0;vuw<2;vuw++){
01083 for(i=0;i<1+segbnk2[vuw].GetLast();i++){
01084 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(segbnk2[vuw].At(i));
01085 if( seg->GetUID()>0 && seg->GetUID()<2 ) seg->SetReseedFlag(1);
01086 }
01087 }
01088
01089 gettimeofday(&tpafter,0);
01090 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
01091 fTime_ShowerAtNu+=dtime;
01092
01093
01094 // ***********************************************
01095 // * F O R M A T I O N O F 2 D T R A C K S *
01096 // ***********************************************
01097
01098 gettimeofday(&tpbefore,0);
01099
01100 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of 2D tracks *** " << endl;
01101 Int_t bpln,epln;
01102 Int_t objctr,objarrctr;
01103 Double_t dscr,topscr;
01104 Double_t strw,strdw;
01105 Int_t idm,idp;
01106 Int_t ntrks;
01107 Int_t tmpflag,trkflag;
01108 TObjArray tmpb,tmpe;
01109 TObjArray begbnk,endbnk;
01110 TObjArray tmpbegbnk,tmpendbnk;
01111 TObjArray tmpbnk,tmpbnk2;
01112
01113 if(1+segbnk[0].GetLast()>0 && 1+segbnk[1].GetLast()>0){
01114 for(vuw=0;vuw<2;vuw++){
01115
01116 cont=1; ntrks=0;
01117 while(cont){
01118 cont=0;
01119 begbnk.Clear(); endbnk.Clear();
01120
01121 // select seed segments
01122 tmp1.Clear(); tmp2.Clear();
01123 for(i=0;i<1+segbnk[vuw].GetLast();i++){
01124 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(segbnk[vuw].At(i));
01125 if( seg->GetUID()>0 && seg->GetUID()<3 ){
01126 ctr=0;
01127 for(j=0;j<1+seg->GetClrLast();j++){
01128 ClusterAtNu* clr = (ClusterAtNu*)(seg->GetClrAt(j));
01129 if( clr->GetTrkPlnFlag()>0 && clr->GetTrkFlag()>1 && clr->GetTrkFlag()<3 ) ctr++;
01130 }
01131 if( ntrks>0 && ctr<3 ) seg->SetUID(0);
01132 }
01133 if(seg->GetUID()==2) tmp1.Add(seg); if(seg->GetUID()==1) tmp2.Add(seg);
01134 }
01135
01136 if(1+tmp1.GetLast()==0){
01137 for(i=0;i<1+tmp2.GetLast();i++){
01138 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmp2.At(i));
01139 tmp1.Add(seg);
01140 }
01141 }
01142
01143 tmp2.Clear();
01144 if(1+tmp1.GetLast()>0){
01145 npln=-1;
01146 for(i=0;i<1+tmp1.GetLast();i++){
01147 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmp1.At(i));
01148 if( 1+seg->GetEndPlane()-seg->GetBegPlane()>npln ){
01149 npln=1+seg->GetEndPlane()-seg->GetBegPlane();
01150 }
01151 }
01152 if(npln>0){
01153 for(i=0;i<1+tmp1.GetLast();i++){
01154 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmp1.At(i));
01155 if( 1+seg->GetEndPlane()-seg->GetBegPlane()>12
01156 || 1+seg->GetEndPlane()-seg->GetBegPlane()>npln-5 ){
01157 tmp2.Add(seg);
01158 }
01159 }
01160 }
01161 }
01162
01163 // select seed tracks
01164 tmp1.Clear();
01165 if(1+tmp2.GetLast()>0){
01166 npln=-1;
01167 for(i=0;i<1+tmp2.GetLast();i++){
01168 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmp2.At(i));
01169 tmpm.Clear(); tmpp.Clear();
01170
01171 tmpb.Clear(); tmpb.Add(seg); tmpbnk.Clear();
01172 while(1+tmpb.GetLast()>0){
01173 tmp.Clear();
01174 for(j=0;j<1+tmpb.GetLast();j++){
01175 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmpb.At(j));
01176 tmp.Add(segtmp);
01177 }
01178 tmpb.Clear();
01179 for(j=0;j<1+tmp.GetLast();j++){
01180 tmpflag=0;
01181 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmp.At(j));
01182 for(k=0;k<1+segtmp->GetSegBegLast();k++){
01183 TrackSegmentAtNu* segb = (TrackSegmentAtNu*)(segtmp->GetSegBegAt(k));
01184 if(segb->GetUID()>-1 && segb->GetUID()<3){
01185 tmpflag=1;
01186 if(segb->GetTmpTrkFlag()<1){
01187 tmpb.Add(segb); segb->SetTmpTrkFlag(1); tmpbnk.Add(segb);
01188 }
01189 }
01190 }
01191 if(tmpflag==0){
01192 tmpm.Add(segtmp);
01193 }
01194 }
01195 }
01196 for(j=0;j<1+tmpbnk.GetLast();j++){
01197 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmpbnk.At(j));
01198 segtmp->SetTmpTrkFlag(0);
01199 }
01200
01201 tmpe.Clear(); tmpe.Add(seg); tmpbnk.Clear();
01202 while(1+tmpe.GetLast()>0){
01203 tmp.Clear();
01204 for(j=0;j<1+tmpe.GetLast();j++){
01205 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmpe.At(j));
01206 tmp.Add(segtmp);
01207 }
01208 tmpe.Clear();
01209 for(j=0;j<1+tmp.GetLast();j++){
01210 tmpflag=0;
01211 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmp.At(j));
01212 for(k=0;k<1+segtmp->GetSegEndLast();k++){
01213 TrackSegmentAtNu* sege = (TrackSegmentAtNu*)(segtmp->GetSegEndAt(k));
01214 if(sege->GetUID()>-1 && sege->GetUID()<3){
01215 tmpflag=1;
01216 if(sege->GetTmpTrkFlag()<1){
01217 tmpe.Add(sege); sege->SetTmpTrkFlag(1); tmpbnk.Add(sege);
01218 }
01219 }
01220 }
01221 if(tmpflag==0){
01222 tmpp.Add(segtmp);
01223 }
01224 }
01225 }
01226 for(j=0;j<1+tmpbnk.GetLast();j++){
01227 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmpbnk.At(j));
01228 segtmp->SetTmpTrkFlag(0);
01229 }
01230
01231 npts=-1;
01232 for(j=0;j<1+tmpm.GetLast();j++){
01233 TrackSegmentAtNu* segm = (TrackSegmentAtNu*)(tmpm.At(j));
01234 for(k=0;k<1+tmpp.GetLast();k++){
01235 TrackSegmentAtNu* segp = (TrackSegmentAtNu*)(tmpp.At(k));
01236 if(1+segp->GetEndPlane()-segm->GetBegPlane()>npts){
01237 npts=1+segp->GetEndPlane()-segm->GetBegPlane();
01238 }
01239 }
01240 }
01241 seg->SetNPlanes(npts);
01242 if(npts>npln) npln=npts;
01243 }
01244 if(npln>0){
01245 for(i=0;i<1+tmp2.GetLast();i++){
01246 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmp2.At(i));
01247 if( seg->GetNPlanes()>npln-5 ){
01248 tmp1.Add(seg);
01249 }
01250 }
01251 }
01252 }
01253
01254 // propagate segments
01255 tmpbnk2.Clear();
01256 for(i=0;i<1+tmp1.GetLast();i++){
01257 TrackSegmentAtNu* seg0 = (TrackSegmentAtNu*)(tmp1.At(i));
01258
01259 MSG("AlgAtNuReco", Msg::kVerbose)
01260 << " PROPAGATE SEED SEGMENT " << seg0->GetBegPlane() << " " << seg0->GetEndPlane() << endl;
01261
01262 seg0->SetTmpTrkFlag(3); trkflag=0;
01263 tmpbnk.Clear(); tmpbnk.Add(seg0);
01264
01265 // track backwards (1)
01266 tmpm.Clear(); tmpb.Clear(); tmpb.Add(seg0);
01267 while(1+tmpb.GetLast()>0){
01268 tmp.Clear();
01269 for(j=0;j<1+tmpb.GetLast();j++){
01270 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmpb.At(j));
01271 tmp.Add(segtmp);
01272 }
01273 tmpb.Clear();
01274 for(j=0;j<1+tmp.GetLast();j++){
01275 tmpflag=0;
01276 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmp.At(j));
01277 for(k=0;k<1+segtmp->GetSegBegLast();k++){
01278 TrackSegmentAtNu* segb = (TrackSegmentAtNu*)(segtmp->GetSegBegAt(k));
01279 if(segb->GetUID()>-1 && segb->GetUID()<3){
01280 tmpflag=1;
01281 if(segb->GetTmpTrkFlag()<1){
01282 tmpb.Add(segb); segb->SetTmpTrkFlag(1); tmpbnk.Add(segb);
01283 }
01284 }
01285 }
01286 if(tmpflag==0){
01287 tmpm.Add(segtmp);
01288 }
01289 }
01290 }
01291
01292 npln=999;
01293 for(j=0;j<1+tmpm.GetLast();j++){
01294 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpm.At(j));
01295 if(seg->GetBegPlane()<npln) npln=seg->GetBegPlane();
01296 }
01297
01298 tmpb.Clear();
01299 for(j=0;j<1+tmpm.GetLast();j++){
01300 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpm.At(j));
01301 if(seg->GetBegPlane()<npln+5){
01302 tmpb.Add(seg);
01303 }
01304 }
01305
01306 while(1+tmpb.GetLast()>0){
01307 tmp.Clear();
01308 for(j=0;j<1+tmpb.GetLast();j++){
01309 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpb.At(j));
01310 if(seg->GetTmpTrkFlag()<2){
01311 seg->SetTmpTrkFlag(2); tmp.Add(seg);
01312 }
01313 if(seg->GetTrkFlag()<1){
01314 seg->SetTrkFlag(1); tmpbnk2.Add(seg); trkflag=1;
01315 }
01316 }
01317 tmpb.Clear();
01318 for(j=0;j<1+tmp.GetLast();j++){
01319 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmp.At(j));
01320 for(k=0;k<1+seg->GetSegEndLast();k++){
01321 TrackSegmentAtNu* segb = (TrackSegmentAtNu*)(seg->GetSegEndAt(k));
01322 if(segb->GetTmpTrkFlag()>0){
01323 tmpb.Add(segb);
01324 }
01325 }
01326 }
01327 }
01328
01329 // track forwards (1)
01330 tmpp.Clear(); tmpe.Clear(); tmpe.Add(seg0);
01331 while(1+tmpe.GetLast()>0){
01332 tmp.Clear();
01333 for(j=0;j<1+tmpe.GetLast();j++){
01334 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmpe.At(j));
01335 tmp.Add(segtmp);
01336 }
01337 tmpe.Clear();
01338 for(j=0;j<1+tmp.GetLast();j++){
01339 tmpflag=0;
01340 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(tmp.At(j));
01341 for(k=0;k<1+segtmp->GetSegEndLast();k++){
01342 TrackSegmentAtNu* sege = (TrackSegmentAtNu*)(segtmp->GetSegEndAt(k));
01343 if(sege->GetUID()>-1 && sege->GetUID()<3){
01344 tmpflag=1;
01345 if(sege->GetTmpTrkFlag()<1){
01346 tmpe.Add(sege); sege->SetTmpTrkFlag(1); tmpbnk.Add(sege);
01347 }
01348 }
01349 }
01350 if(tmpflag==0){
01351 tmpp.Add(segtmp);
01352 }
01353 }
01354 }
01355
01356 npln=-999;
01357 for(j=0;j<1+tmpp.GetLast();j++){
01358 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpp.At(j));
01359 if(seg->GetEndPlane()>npln) npln=seg->GetEndPlane();
01360 }
01361
01362 tmpe.Clear();
01363 for(j=0;j<1+tmpp.GetLast();j++){
01364 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpp.At(j));
01365 if(seg->GetEndPlane()>npln-5){
01366 tmpe.Add(seg);
01367 }
01368 }
01369
01370 while(1+tmpe.GetLast()>0){
01371 tmp.Clear();
01372 for(j=0;j<1+tmpe.GetLast();j++){
01373 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpe.At(j));
01374 if(seg->GetTmpTrkFlag()<2){
01375 seg->SetTmpTrkFlag(2); tmp.Add(seg);
01376 }
01377 if(seg->GetTrkFlag()<1){
01378 seg->SetTrkFlag(1); tmpbnk2.Add(seg); trkflag=1;
01379 }
01380 }
01381 tmpe.Clear();
01382 for(j=0;j<1+tmp.GetLast();j++){
01383 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmp.At(j));
01384 for(k=0;k<1+seg->GetSegBegLast();k++){
01385 TrackSegmentAtNu* sege = (TrackSegmentAtNu*)(seg->GetSegBegAt(k));
01386 if(sege->GetTmpTrkFlag()>0){
01387 tmpe.Add(sege);
01388 }
01389 }
01390 }
01391 }
01392
01393 if(trkflag){
01394
01395 // track backwards (2)
01396 tmpbegbnk.Clear(); objarrctr=0;
01397 TObjArray* tmpbegsegbnk = new TObjArray(); NEWOBJCTR++;
01398 tmpbegsegbnk->Add(seg0); tmpbegbnk.Add(tmpbegsegbnk); objctr=1;
01399 while(objctr>0){
01400 objctr=0; npts=1+tmpbegbnk.GetLast();
01401 for(j=0;j<npts;j++){
01402 TObjArray* objtmp = (TObjArray*)(tmpbegbnk.At(j));
01403 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(objtmp->Last());
01404 tmpm.Clear();
01405 for(k=0;k<1+segtmp->GetSegBegLast();k++){
01406 TrackSegmentAtNu* segbeg = (TrackSegmentAtNu*)(segtmp->GetSegBegAt(k));
01407 if(segbeg->GetTmpTrkFlag()>1){
01408 tmpm.Add(segbeg);
01409 }
01410 }
01411 if(1+tmpm.GetLast()>0){
01412 for(k=1;k<1+tmpm.GetLast();k++){
01413 if(objarrctr<fObjMax){
01414 TObjArray* objnew = new TObjArray(*objtmp); NEWOBJCTR++;
01415 TrackSegmentAtNu* segbeg = (TrackSegmentAtNu*)(tmpm.At(k));
01416 objnew->Add(segbeg); tmpbegbnk.Add(objnew);
01417 objarrctr++;
01418 }
01419 else{
01420 if(!fObjMaxFlag){
01421 fObjMaxFlag=1;
01422 MSG("AlgAtNuReco", Msg::kWarning) << " *** CRAZY MEMORY GUZZLING *** " << endl;
01423 }
01424 }
01425 }
01426 TrackSegmentAtNu* segbeg = (TrackSegmentAtNu*)(tmpm.First());
01427 objtmp->Add(segbeg); objctr++;
01428 }
01429 }
01430 }
01431
01432 id=-1; topscr=-1.0;
01433 for(j=0;j<1+tmpbegbnk.GetLast();j++){
01434 TObjArray* objtmp = (TObjArray*)(tmpbegbnk.At(j));
01435 scr = seg0->GetScore(objtmp,0);
01436 if(scr>topscr){ topscr=scr; id=j; }
01437 }
01438
01439 TObjArray* objbegtmp = (TObjArray*)(tmpbegbnk.At(id));
01440 TObjArray* objbegnew = new TObjArray(*objbegtmp); begbnk.Add(objbegnew); NEWOBJCTR++;
01441 DELOBJCTR+=1+tmpbegbnk.GetLast(); tmpbegbnk.Delete();
01442
01443 // track forwards (2)
01444 tmpendbnk.Clear(); objarrctr=0;
01445 TObjArray* tmpendsegbnk = new TObjArray(); NEWOBJCTR++;
01446 tmpendsegbnk->Add(seg0); tmpendbnk.Add(tmpendsegbnk); objctr=1;
01447 while(objctr>0){
01448 objctr=0; npts=1+tmpendbnk.GetLast();
01449 for(j=0;j<npts;j++){
01450 TObjArray* objtmp = (TObjArray*)(tmpendbnk.At(j));
01451 TrackSegmentAtNu* segtmp = (TrackSegmentAtNu*)(objtmp->Last());
01452 tmpp.Clear();
01453 for(k=0;k<1+segtmp->GetSegEndLast();k++){
01454 TrackSegmentAtNu* segend = (TrackSegmentAtNu*)(segtmp->GetSegEndAt(k));
01455 if(segend->GetTmpTrkFlag()>1){
01456 tmpp.Add(segend);
01457 }
01458 }
01459 if(1+tmpp.GetLast()>0){
01460 for(k=1;k<1+tmpp.GetLast();k++){
01461 if(objarrctr<fObjMax){
01462 TObjArray* objnew = new TObjArray(*objtmp); NEWOBJCTR++;
01463 TrackSegmentAtNu* segend = (TrackSegmentAtNu*)(tmpp.At(k));
01464 objnew->Add(segend); tmpendbnk.Add(objnew);
01465 objarrctr++;
01466 }
01467 else{
01468 if(!fObjMaxFlag){
01469 fObjMaxFlag=1;
01470 MSG("AlgAtNuReco", Msg::kWarning) << " *** CRAZY MEMORY GUZZLING *** " << endl;
01471 }
01472 }
01473 }
01474 TrackSegmentAtNu* segend = (TrackSegmentAtNu*)(tmpp.First());
01475 objtmp->Add(segend); objctr++;
01476 }
01477 }
01478 }
01479
01480 id=-1; topscr=-1.0;
01481 for(j=0;j<1+tmpendbnk.GetLast();j++){
01482 TObjArray* objtmp = (TObjArray*)(tmpendbnk.At(j));
01483 scr = seg0->GetScore(0,objtmp);
01484 if(scr>topscr){ topscr=scr; id=j; }
01485 }
01486
01487 TObjArray* objendtmp = (TObjArray*)(tmpendbnk.At(id));
01488 TObjArray* objendnew = new TObjArray(*objendtmp); endbnk.Add(objendnew); NEWOBJCTR++;
01489 DELOBJCTR+=1+tmpendbnk.GetLast(); tmpendbnk.Delete();
01490 }
01491
01492 for(j=0;j<1+tmpbnk.GetLast();j++){
01493 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpbnk.At(j));
01494 seg->SetTmpTrkFlag(0);
01495 }
01496 }
01497
01498 for(j=0;j<1+tmpbnk2.GetLast();j++){
01499 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmpbnk2.At(j));
01500 seg->SetTrkFlag(0);
01501 }
01502
01503 // select best track
01504 if(1+begbnk.GetLast()>0 && 1+endbnk.GetLast()>0){
01505
01506 id=-1; topscr=-1.0;
01507 for(i=0;i<1+begbnk.GetLast();i++){
01508 TObjArray* objbeg = (TObjArray*)(begbnk.At(i));
01509 TObjArray* objend = (TObjArray*)(endbnk.At(i));
01510 TrackSegmentAtNu* seg0 = (TrackSegmentAtNu*)(objbeg->First());
01511 scr = seg0->GetScore(objbeg,objend);
01512 if(scr>topscr){ topscr=scr; id=i; }
01513 }
01514
01515 if(1+id>0){
01516 TObjArray* objbeg = (TObjArray*)(begbnk.At(id));
01517 TObjArray* objend = (TObjArray*)(endbnk.At(id));
01518 TrackSegmentAtNu* seg0 = (TrackSegmentAtNu*)(objbeg->First());
01519 for(i=1;i<1+objbeg->GetLast();i++){
01520 TrackSegmentAtNu* segbeg = (TrackSegmentAtNu*)(objbeg->At(i));
01521 seg0->AddSegment(segbeg); segbeg->SetUID(-1);
01522 }
01523 for(i=1;i<1+objend->GetLast();i++){
01524 TrackSegmentAtNu* segend = (TrackSegmentAtNu*)(objend->At(i));
01525 seg0->AddSegment(segend); segend->SetUID(-1);
01526 }
01527 seg0->SetUID(3);
01528 for(i=0;i<1+seg0->GetClrLast();i++){
01529 ClusterAtNu* clr = (ClusterAtNu*)(seg0->GetClrAt(i));
01530 clr->SetTrkFlag(3);
01531 }
01532 cont=1; ntrks++;
01533 }
01534 }
01535
01536 MSG("AlgAtNuReco", Msg::kVerbose) << " ... about to delete beg/end... " << endl;
01537 DELOBJCTR+=1+begbnk.GetLast(); DELOBJCTR+=1+endbnk.GetLast();
01538 begbnk.Delete(); endbnk.Delete();
01539 MSG("AlgAtNuReco", Msg::kVerbose) << " beg/end deleted ... " << endl;
01540 }
01541 }
01542 }
01543
01544 // Print out list of 2D tracks
01545 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF 2D TRACKS *** " << endl;
01546 for(vuw=0;vuw<2;vuw++){
01547 MSG("AlgAtNuReco", Msg::kVerbose) << "view: " << vuw << endl;
01548 for(i=0;i<1+segbnk[vuw].GetLast();i++){
01549 TrackSegmentAtNu* segment = (TrackSegmentAtNu*)(segbnk[vuw].At(i));
01550 if(segment->GetUID()>2){
01551 MSG("AlgAtNuReco", Msg::kVerbose) << " i=" << i
01552 << " begpln: " << segment->GetBegPlane()
01553 << " begstr: " << segment->GetBegStrip()
01554 << " endpln: " << segment->GetEndPlane()
01555 << " endstr: " << segment->GetEndStrip() << endl;
01556 }
01557 }
01558 }
01559
01560 gettimeofday(&tpafter,0);
01561 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
01562 fTime_TrackAtNu+=dtime;
01563
01564
01565 // *************************************************
01566 // * F O R M A T I O N O F 2 D S H O W E R S *
01567 // *************************************************
01568
01569 gettimeofday(&tpbefore,0);
01570
01571 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of 2D showers *** " << endl;
01572 Int_t nshws,maxshws;
01573 if(1+segbnk2[0].GetLast()>0 && 1+segbnk2[1].GetLast()>0){
01574 for(vuw=0;vuw<2;vuw++){
01575
01576 cont=1; nshws=0; maxshws=4;
01577 while(cont){
01578 cont=0;
01579
01580 // select segments
01581 tmp1.Clear(); tmp2.Clear();
01582 for(i=0;i<1+segbnk2[vuw].GetLast();i++){
01583 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(segbnk2[vuw].At(i));
01584 if(seg->GetUID()>0 && seg->GetUID()<3){
01585 ctr=0;
01586 for(j=0;j<1+seg->GetClrLast();j++){
01587 ClusterAtNu* clr = (ClusterAtNu*)(seg->GetClrAt(j));
01588 if( clr->GetShwPlnFlag()>0 && clr->GetShwFlag()>1 && clr->GetShwFlag()<3 ) ctr++;
01589 }
01590 if( nshws>0 && ctr<1 ) seg->SetUID(0);
01591 }
01592 if(seg->GetUID()==2) tmp1.Add(seg); if(seg->GetUID()==1) tmp2.Add(seg);
01593 }
01594
01595 if(1+tmp1.GetLast()==0){
01596 for(i=0;i<1+tmp2.GetLast();i++){
01597 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(tmp2.At(i));
01598 tmp1.Add(seg);
01599 }
01600 }
01601
01602 // determine primary segment
01603 npln=0; id=-1; topscr=0.0;
01604 for(i=0;i<1+tmp1.GetLast();i++){
01605 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(tmp1.At(i));
01606 npts = 1+(seg->GetEndPlane()-seg->GetBegPlane())/2;
01607 if(npts>npln-3){
01608 scr=seg->GetScore();
01609 if(npts>npln || scr>topscr){
01610 npln=npts; topscr=scr;
01611 id=i;
01612 }
01613 }
01614 }
01615
01616 // propagate primary segment
01617 if(1+id>0){
01618 ShowerSegmentAtNu* seg0 = (ShowerSegmentAtNu*)(tmp1.At(id));
01619
01620 mod=(Int_t)((seg0->GetEndPlane())/(modpln+1));
01621 seg0->SetUID(3);
01622 cont=1; nshws++;
01623
01624 ctr=1;
01625 while(ctr>0){
01626 ctr=0;
01627
01628 // add associated clusters
01629 pln0 = -4 + seg0->GetBegPlane();
01630 npln = 5 + (seg0->GetEndPlane()-seg0->GetBegPlane())/2;
01631 for(k=0;k<npln;k++){
01632 k0 = pln0+2*k;
01633 if(k0>mod*(modpln+1) && k0<(mod+1)*(modpln+1)){
01634 for(k1=0;k1<1+clrplnbnk[k0].GetLast();k1++){
01635 ClusterAtNu* clr1 = (ClusterAtNu*)(clrplnbnk[k0].At(k1));
01636 if(clr1->GetShwFlag()<2){
01637 assoc=seg0->IsShwAssoc(clr1);
01638 if( ( assoc>0 && clr1->GetTrkFlag()<3 )
01639 || ( assoc>1 && clr1->GetPlane()-seg0->GetBegPlane()>0 && seg0->GetEndPlane()-clr1->GetPlane()>0 ) ){
01640 seg0->AddCluster(clr1);
01641 clr1->SetShwFlag(2); ctr++;
01642 }
01643 }
01644 }
01645 }
01646 }
01647
01648 // add associated segments
01649 for(i=0;i<1+segbnk2[vuw].GetLast();i++){
01650 ShowerSegmentAtNu* segtmp = (ShowerSegmentAtNu*)(segbnk2[vuw].At(i));
01651 if( segtmp->GetUID()>0 && segtmp->GetUID()<3
01652 && (Int_t)((segtmp->GetEndPlane())/(modpln+1))==mod ){
01653 assoc=seg0->IsShwAssoc(segtmp);
01654 if(assoc>0){
01655 seg0->AddSegment(segtmp);
01656 segtmp->SetUID(-1); ctr++;
01657 }
01658 }
01659 }
01660
01661 }
01662
01663 for(i=0;i<1+seg0->GetClrLast();i++){
01664 ClusterAtNu* clr = (ClusterAtNu*)(seg0->GetClrAt(i));
01665 clr->SetShwFlag(3);
01666 }
01667
01668 }
01669
01670 }
01671 }
01672 }
01673
01674 // Print out list of 2D showers
01675 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF 2D SHOWERS *** " << endl;
01676 for(vuw=0;vuw<2;vuw++){
01677 MSG("AlgAtNuReco", Msg::kVerbose) << "view: " << vuw << endl;
01678 for(i=0;i<1+segbnk2[vuw].GetLast();i++){
01679 ShowerSegmentAtNu* segment = (ShowerSegmentAtNu*)(segbnk2[vuw].At(i));
01680 if(segment->GetUID()>2){
01681 MSG("AlgAtNuReco", Msg::kVerbose) << " i=" << i
01682 << " begpln: " << segment->GetBegPlane()
01683 << " begstr: " << segment->GetBegStrip()
01684 << " endpln: " << segment->GetEndPlane()
01685 << " endstr: " << segment->GetEndStrip() << endl;
01686 }
01687 }
01688 }
01689
01690 gettimeofday(&tpafter,0);
01691 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
01692 fTime_ShowerAtNu+=dtime;
01693
01694
01695 // ******************************************
01696 // * C O M P A R E U / V V I E W S (2) *
01697 // ******************************************
01698
01699 MSG("AlgAtNuReco", Msg::kDebug) << " *** matching views (2) *** " << endl;
01700
01701 gettimeofday(&tpbefore,0);
01702
01703 // match tracks
01704 TObjArray tmptrk[2],tmptrkseg[2];
01705 Int_t dstr,dpln,dpln1,dpln2;
01706 Int_t bpln2,bpln1,epln1,epln2;
01707 Int_t jnflag2[2];
01708 for(i=0;i<1+segbnk[0].GetLast();i++){
01709 TrackSegmentAtNu* segu = (TrackSegmentAtNu*)(segbnk[0].At(i));
01710 if(segu->GetUID()>2){
01711 for(j=0;j<1+segbnk[1].GetLast();j++){
01712 TrackSegmentAtNu* segv = (TrackSegmentAtNu*)(segbnk[1].At(j));
01713 if(segv->GetUID()>2){
01714 if( ( (segv->GetBegPlane()-segu->GetBegPlane()>-10
01715 && segv->GetBegPlane()-segu->GetBegPlane()<10)
01716 || (segv->GetEndPlane()-segu->GetEndPlane()>-10
01717 && segv->GetEndPlane()-segu->GetEndPlane()<10) )
01718 && ( segv->GetEndPlane()-segu->GetBegPlane()>2
01719 && segu->GetEndPlane()-segv->GetBegPlane()>2 )
01720 && ( segv->GetEndTime()-segu->GetBegTime()>-fTimeWin
01721 && segu->GetEndTime()-segv->GetBegTime()>-fTimeWin ) ){
01722 tmptrkseg[0].Add(segu); tmptrkseg[1].Add(segv);
01723 }
01724 }
01725 }
01726 }
01727 }
01728
01729 if(modnum>1){
01730 TObjArray* tmpmod = new TObjArray[2*(modnum-1)]; NEWOBJCTR+=2*(modnum-1);
01731 for(i=0;i<2*(modnum-1);i++){
01732 for(j=0;j<2;j++){
01733 TObjArray* tmparr = new TObjArray(); NEWOBJCTR++;
01734 tmpmod[i].Add(tmparr);
01735 }
01736 }
01737
01738 for(i=0;i<1+tmptrkseg[0].GetLast();i++){
01739 TrackSegmentAtNu* segu = (TrackSegmentAtNu*)(tmptrkseg[0].At(i));
01740 TrackSegmentAtNu* segv = (TrackSegmentAtNu*)(tmptrkseg[1].At(i));
01741 mod=(Int_t)((segu->GetEndPlane())/(modpln+1));
01742 if( segu->GetBegPlane()<mod*(modpln+1)+20
01743 && segv->GetBegPlane()<mod*(modpln+1)+20
01744 && mod>0 ){
01745 TObjArray* tmparrum = (TObjArray*)(tmpmod[2*mod-1].At(0));
01746 TObjArray* tmparrvm = (TObjArray*)(tmpmod[2*mod-1].At(1));
01747 tmparrum->Add(segu); tmparrvm->Add(segv);
01748 }
01749 if( segu->GetEndPlane()>(mod+1)*(modpln+1)-20
01750 && segv->GetEndPlane()>(mod+1)*(modpln+1)-20
01751 && mod<modnum-1 ){
01752 TObjArray* tmparru = (TObjArray*)(tmpmod[2*mod].At(0));
01753 TObjArray* tmparrv = (TObjArray*)(tmpmod[2*mod].At(1));
01754 tmparru->Add(segu); tmparrv->Add(segv);
01755 }
01756 }
01757
01758 for(mod=0;mod<modnum-1;mod++){
01759 TObjArray* tmparru = (TObjArray*)(tmpmod[2*mod].At(0));
01760 TObjArray* tmparrv = (TObjArray*)(tmpmod[2*mod].At(1));
01761 TObjArray* tmparrup = (TObjArray*)(tmpmod[2*mod+1].At(0));
01762 TObjArray* tmparrvp = (TObjArray*)(tmpmod[2*mod+1].At(1));
01763
01764 if( 1+tmparru->GetLast()>0 && 1+tmparrup->GetLast()>0 ){
01765
01766 cont=1;
01767 while(cont){
01768 cont=0;
01769 idm=-1; idp=-1; topscr=0;
01770 for(i=0;i<1+tmparru->GetLast();i++){
01771 TrackSegmentAtNu* segum = (TrackSegmentAtNu*)(tmparru->At(i));
01772 TrackSegmentAtNu* segvm = (TrackSegmentAtNu*)(tmparrv->At(i));
01773 if( (segum->GetPartner()==segvm && segvm->GetPartner()==segum)
01774 || (segum->GetPartner()==0 && segvm->GetPartner()==0
01775 && segum->GetUID()>2 && segvm->GetUID()>2) ){
01776 for(j=0;j<1+tmparrup->GetLast();j++){
01777 TrackSegmentAtNu* segup = (TrackSegmentAtNu*)(tmparrup->At(j));
01778 TrackSegmentAtNu* segvp = (TrackSegmentAtNu*)(tmparrvp->At(j));
01779 if(segup->GetUID()>2 && segvp->GetUID()>2){
01780 if( segum->IsEndAssoc(segup) && segvm->IsEndAssoc(segvp) ){
01781 if(segum->GetBegPlane()<=segvm->GetBegPlane()){
01782 bpln2=segum->GetBegPlane(); bpln1=segvm->GetBegPlane(); }
01783 if(segvm->GetBegPlane()<=segum->GetBegPlane()){
01784 bpln2=segvm->GetBegPlane(); bpln1=segum->GetBegPlane(); }
01785 if(segup->GetEndPlane()>=segvp->GetEndPlane()){
01786 epln2=segup->GetEndPlane(); epln1=segvp->GetEndPlane(); }
01787 if(segvp->GetEndPlane()>=segup->GetEndPlane()){
01788 epln2=segvp->GetEndPlane(); epln1=segup->GetEndPlane(); }
01789 dpln1=epln2-bpln2; dpln2=epln1-bpln1;
01790 scr=dpln2*dpln2/dpln1;
01791 if(scr>topscr){
01792 idm=i; idp=j; topscr=scr;
01793 }
01794 }
01795 }
01796 }
01797 }
01798 }
01799 if(1+idm>0 && 1+idp>0){
01800 TrackSegmentAtNu* segum = (TrackSegmentAtNu*)(tmparru->At(idm));
01801 TrackSegmentAtNu* segvm = (TrackSegmentAtNu*)(tmparrv->At(idm));
01802 TrackSegmentAtNu* segup = (TrackSegmentAtNu*)(tmparrup->At(idp));
01803 TrackSegmentAtNu* segvp = (TrackSegmentAtNu*)(tmparrvp->At(idp));
01804 segup->AddSegment(segum); segup->SetPartner(segvp);
01805 segum->SetUID(-1); segum->SetPartner(0);
01806 segvp->AddSegment(segvm); segvp->SetPartner(segup);
01807 segvm->SetUID(-1); segvm->SetPartner(0);
01808 cont=1;
01809 }
01810 }
01811 }
01812 }
01813
01814 for(i=0;i<2*(modnum-1);i++){
01815 DELOBJCTR+=1+tmpmod[i].GetLast(); tmpmod[i].Delete();
01816 }
01817 delete [] tmpmod; DELOBJCTR+=2*(modnum-1);
01818 }
01819
01820 // determine other best matched tracks
01821 cont=1;
01822 while(cont){
01823 cont=0;
01824 id=-1; topscr=0.0;
01825 for(i=0;i<1+tmptrkseg[0].GetLast();i++){
01826 TrackSegmentAtNu* segu = (TrackSegmentAtNu*)(tmptrkseg[0].At(i));
01827 TrackSegmentAtNu* segv = (TrackSegmentAtNu*)(tmptrkseg[1].At(i));
01828 if( segu->GetPartner()==segv && segv->GetPartner()==segu
01829 && segu->GetUID()<4 && segv->GetUID()<4 ){
01830 tmptrk[0].Add(segu); tmptrk[1].Add(segv);
01831 segu->SetUID(4); segv->SetUID(4);
01832 }
01833 else{
01834 if( segu->GetPartner()==0 && segv->GetPartner()==0
01835 && segu->GetUID()>2 && segv->GetUID()>2 ){
01836 if(segu->GetBegPlane()<=segv->GetBegPlane()){
01837 bpln2=segu->GetBegPlane(); bpln1=segv->GetBegPlane(); }
01838 if(segv->GetBegPlane()<=segu->GetBegPlane()){
01839 bpln2=segv->GetBegPlane(); bpln1=segu->GetBegPlane(); }
01840 if(segu->GetEndPlane()>=segv->GetEndPlane()){
01841 epln2=segu->GetEndPlane(); epln1=segv->GetEndPlane(); }
01842 if(segv->GetEndPlane()>=segu->GetEndPlane()){
01843 epln2=segv->GetEndPlane(); epln1=segu->GetEndPlane(); }
01844 dpln1=epln2-bpln2; dpln2=epln1-bpln1;
01845 scr=dpln2*dpln2/dpln1;
01846 if(scr>topscr){
01847 topscr=scr; id=i;
01848 }
01849 }
01850 }
01851 }
01852 if(1+id>0){
01853 TrackSegmentAtNu* segu = (TrackSegmentAtNu*)(tmptrkseg[0].At(id));
01854 TrackSegmentAtNu* segv = (TrackSegmentAtNu*)(tmptrkseg[1].At(id));
01855 segu->SetPartner(segv); segv->SetPartner(segu);
01856 cont=1;
01857 }
01858 }
01859
01860 ntrks=1+tmptrk[0].GetLast();
01861 for(vuw=0;vuw<2;vuw++){
01862 for(i=0;i<1+tmptrk[vuw].GetLast();i++){
01863 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmptrk[vuw].At(i));
01864 for(j=0;j<1+seg->GetClrLast();j++){
01865 ClusterAtNu* clr = (ClusterAtNu*)(seg->GetClrAt(j));
01866 clr->SetTrkFlag(4);
01867 }
01868 }
01869 }
01870
01871 gettimeofday(&tpafter,0);
01872 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
01873 fTime_TrackAtNu+=dtime;
01874
01875 gettimeofday(&tpbefore,0);
01876
01877 // match showers
01878 TObjArray tmpshw[2],tmpshwseg[2];
01879 for(i=0;i<1+segbnk2[0].GetLast();i++){
01880 ShowerSegmentAtNu* segu = (ShowerSegmentAtNu*)(segbnk2[0].At(i));
01881 if(segu->GetUID()>2){
01882 for(j=0;j<1+segbnk2[1].GetLast();j++){
01883 ShowerSegmentAtNu* segv = (ShowerSegmentAtNu*)(segbnk2[1].At(j));
01884 if(segv->GetUID()>2){
01885 if( ( segu->GetEndPlane()-segv->GetBegPlane()>-2
01886 && segv->GetEndPlane()-segu->GetBegPlane()>-2)
01887 && ( segu->GetEndTime()-segv->GetBegTime()>-fTimeWin
01888 && segv->GetEndTime()-segu->GetBegTime()>-fTimeWin ) ){
01889 tmpshwseg[0].Add(segu); tmpshwseg[1].Add(segv);
01890 }
01891 }
01892 }
01893 }
01894 }
01895
01896 // match showers that intersect tracks
01897 TObjArray tmpclr[2],tmpseg[2];
01898 if(1+tmptrk[0].GetLast()>0){
01899 for(i=0;i<1+tmptrk[0].GetLast();i++){
01900
01901 for(vuw=0;vuw<2;vuw++){
01902 tmpclr[vuw].Clear(); tmpseg[vuw].Clear();
01903 }
01904
01905 for(vuw=0;vuw<2;vuw++){
01906 TrackSegmentAtNu* segtrk = (TrackSegmentAtNu*)(tmptrk[vuw].At(i));
01907 for(j=0;j<1+segtrk->GetClrLast();j++){
01908 ClusterAtNu* clr = (ClusterAtNu*)(segtrk->GetClrAt(j));
01909 if(clr->GetShwFlag()>2) tmpclr[vuw].Add(clr);
01910 }
01911 }
01912
01913 if(1+tmpclr[0].GetLast()>0 && 1+tmpclr[1].GetLast()>0){
01914 for(j=0;j<1+tmpshwseg[0].GetLast();j++){
01915 for(vuw=0;vuw<2;vuw++){
01916 jnflag2[vuw]=0;
01917 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(tmpshwseg[vuw].At(j));
01918 if(seg->GetPartner()==0){
01919 for(k=0;k<1+tmpclr[vuw].GetLast();k++){
01920 ClusterAtNu* clr = (ClusterAtNu*)(tmpclr[vuw].At(k));
01921 if(seg->ContainsClr(clr)) jnflag2[vuw]=1;
01922 }
01923 }
01924 }
01925 if(jnflag2[0]>0 && jnflag2[1]>0){
01926 for(vuw=0;vuw<2;vuw++){
01927 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(tmpshwseg[vuw].At(j));
01928 tmpseg[vuw].Add(seg);
01929 }
01930 }
01931 }
01932 }
01933
01934 if(1+tmpseg[0].GetLast()>0 && 1+tmpseg[1].GetLast()>0){
01935 cont=1;
01936 while(cont){
01937 cont=0;
01938 id=-1; topscr=0.0;
01939 for(j=0;j<1+tmpseg[0].GetLast();j++){
01940 ShowerSegmentAtNu* segshwu = (ShowerSegmentAtNu*)(tmpseg[0].At(j));
01941 ShowerSegmentAtNu* segshwv = (ShowerSegmentAtNu*)(tmpseg[1].At(j));
01942 if( segshwu->GetPartner()==0 && segshwv->GetPartner()==0
01943 && segshwu->GetUID()>2 && segshwv->GetUID()>2 ){
01944 dpln1=segshwu->GetEndPlane()-segshwv->GetBegPlane();
01945 dpln2=segshwv->GetEndPlane()-segshwu->GetBegPlane();
01946 if(dpln1<=dpln2) scr=dpln1; if(dpln2<=dpln1) scr=dpln2;
01947 if(scr>topscr){
01948 id=j; scr=topscr;
01949 }
01950 }
01951 }
01952 if(1+id>0){
01953 ShowerSegmentAtNu* segshwu = (ShowerSegmentAtNu*)(tmpseg[0].At(id));
01954 ShowerSegmentAtNu* segshwv = (ShowerSegmentAtNu*)(tmpseg[1].At(id));
01955 segshwu->SetPartner(segshwv); segshwv->SetPartner(segshwu);
01956 cont=1;
01957 }
01958 }
01959 }
01960 }
01961 }
01962
01963 // determine other best matched showers
01964 cont=1;
01965 while(cont){
01966 cont=0;
01967 id=-1; topscr=0.0;
01968 for(i=0;i<1+tmpshwseg[1].GetLast();i++){
01969 ShowerSegmentAtNu* segu = (ShowerSegmentAtNu*)(tmpshwseg[0].At(i));
01970 ShowerSegmentAtNu* segv = (ShowerSegmentAtNu*)(tmpshwseg[1].At(i));
01971 if( segu->GetPartner()==segv && segv->GetPartner()==segu
01972 && segu->GetUID()<4 && segv->GetUID()<4 ){
01973 tmpshw[0].Add(segu); tmpshw[1].Add(segv);
01974 segu->SetUID(4); segv->SetUID(4);
01975 }
01976 else{
01977 if( segu->GetPartner()==0 && segv->GetPartner()==0
01978 && segu->GetUID()>2 && segv->GetUID()>2 ){
01979 dpln1 = segu->GetEndPlane()-segv->GetBegPlane();
01980 dpln2 = segv->GetEndPlane()-segu->GetBegPlane();
01981 if(dpln1<=dpln2) scr=dpln1; if(dpln2<=dpln1) scr=dpln2;
01982 if(scr>topscr){
01983 topscr=scr; id=i;
01984 }
01985 }
01986 }
01987 }
01988 if(1+id>0){
01989 ShowerSegmentAtNu* segu = (ShowerSegmentAtNu*)(tmpshwseg[0].At(id));
01990 ShowerSegmentAtNu* segv = (ShowerSegmentAtNu*)(tmpshwseg[1].At(id));
01991 segu->SetPartner(segv); segv->SetPartner(segu);
01992 cont=1;
01993 }
01994 }
01995
01996 nshws=1+tmpshw[0].GetLast();
01997 for(vuw=0;vuw<2;vuw++){
01998 for(i=0;i<1+tmpshw[vuw].GetLast();i++){
01999 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(tmpshw[vuw].At(i));
02000 for(j=0;j<1+seg->GetClrLast();j++){
02001 ClusterAtNu* clr = (ClusterAtNu*)(seg->GetClrAt(j));
02002 clr->SetShwFlag(4);
02003 }
02004 }
02005 }
02006
02007 gettimeofday(&tpafter,0);
02008 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02009 fTime_ShowerAtNu+=dtime;
02010
02011
02012 // *************************************
02013 // * R E S E E D I N G S H O W E R S *
02014 // *************************************
02015
02016 if( 1+tmpshw[0].GetLast()==0 && 1+tmpshw[1].GetLast()==0
02017 && 1+tmptrk[0].GetLast()==0 && 1+tmptrk[1].GetLast()==0 ){
02018
02019 gettimeofday(&tpbefore,0);
02020
02021 // low-multiplicity showers
02022 MSG("AlgAtNuReco",Msg::kDebug) << " reseeding low muliplicity showers " << endl;
02023
02024 for(vuw=0;vuw<2;vuw++){
02025 tmpclr[vuw].Clear(); tmpseg[vuw].Clear(); tmpshwseg[vuw].Clear();
02026 }
02027
02028 // reseed clusters
02029 for(mod=0;mod<modnum;mod++){
02030 for(pln=1;pln<1+modpln;pln++){
02031 i=pln+mod*(modpln+1);
02032 for(j=0;j<1+clrplnbnk[i].GetLast();j++){
02033 ClusterAtNu* clr0 = (ClusterAtNu*)(clrplnbnk[i].At(j));
02034 clr0->SetShwFlag(0);
02035 npts=0;
02036 for(k=-4;k<5;k++){
02037 k0=i+k;
02038 if(k0>mod*(modpln+1) && k0<(mod+1)*(modpln+1)){
02039 for(k1=0;k1<1+clrplnbnk[k0].GetLast();k1++){
02040 ClusterAtNu* clr1 = (ClusterAtNu*)(clrplnbnk[k0].At(k1));
02041 if(clr0->GetCharge()>1.0 && clr1->GetCharge()>1.0){
02042 npts+=1+clr1->GetHitLast();
02043 }
02044 }
02045 }
02046 }
02047 if(npts>3){
02048 vuw=clr0->GetPlaneView(); tmpclr[vuw].Add(clr0); clr0->SetShwFlag(1);
02049 }
02050 }
02051 }
02052 }
02053
02054 // reseed 2D showers
02055 if(1+tmpclr[0].GetLast()>0 && 1+tmpclr[1].GetLast()>0){
02056 for(vuw=0;vuw<2;vuw++){
02057 for(i=0;i<1+tmpclr[vuw].GetLast();i++){
02058 ClusterAtNu* clr0 = (ClusterAtNu*)(tmpclr[vuw].At(i));
02059 if(clr0->GetShwFlag()>0 && clr0->GetShwFlag()<2){
02060 ShowerSegmentAtNu* seg0 = new ShowerSegmentAtNu(clr0); NEWOBJCTR++;
02061 tmpseg[vuw].Add(seg0); shwsegbnk.Add(seg0);
02062 seg0->SetUID(2); seg0->SetReseedFlag(1);
02063 mod=(Int_t)((clr0->GetPlane())/(modpln+1)); clr0->SetShwFlag(2);
02064 ctr=1;
02065 while(ctr>0){
02066 ctr=0;
02067 pln0=-6+seg0->GetBegPlane();
02068 npln=7+(seg0->GetEndPlane()-seg0->GetBegPlane())/2;
02069 for(k=0;k<npln;k++){
02070 k0=pln0+2*k;
02071 if(k0>mod*(modpln+1) && k0<(mod+1)*(modpln+1)){
02072 for(k1=0;k1<1+clrplnbnk[k0].GetLast();k1++){
02073 ClusterAtNu* clr1 = (ClusterAtNu*)(clrplnbnk[k0].At(k1));
02074 if( clr1->GetShwFlag()>0 && clr1->GetShwFlag()<2
02075 && seg0->IsDiffuseShwAssoc(clr1)>0 ){
02076 seg0->AddCluster(clr1);
02077 clr1->SetShwFlag(2); ctr++;
02078 }
02079 }
02080 }
02081 }
02082 }
02083 }
02084 }
02085 }
02086 }
02087
02088 MSG("AlgAtNuReco",Msg::kVerbose) << " low multiplicity 2D showers - U=" << 1+tmpseg[0].GetLast()
02089 << " V=" << 1+tmpseg[1].GetLast() << endl;
02090
02091 // reseed 3D showers
02092 if( 1+tmpseg[0].GetLast()>0 && 1+tmpseg[1].GetLast()>0 ){
02093
02094 for(i=0;i<1+tmpseg[0].GetLast();i++){
02095 ShowerSegmentAtNu* segu = (ShowerSegmentAtNu*)(tmpseg[0].At(i));
02096 for(j=0;j<1+tmpseg[1].GetLast();j++){
02097 ShowerSegmentAtNu* segv = (ShowerSegmentAtNu*)(tmpseg[1].At(j));
02098 if( 1+segu->GetHitLast()>2 || 1+segv->GetHitLast()>2
02099 && segu->GetEndPlane()-segv->GetBegPlane()>-4
02100 && segv->GetEndPlane()-segu->GetBegPlane()>-4 ){
02101 tmpshwseg[0].Add(segu); tmpshwseg[1].Add(segv);
02102 }
02103 }
02104 }
02105
02106 id=-1; topscr=0.0;
02107 for(i=0;i<1+tmpshwseg[0].GetLast();i++){
02108 ShowerSegmentAtNu* segu = (ShowerSegmentAtNu*)(tmpshwseg[0].At(i));
02109 ShowerSegmentAtNu* segv = (ShowerSegmentAtNu*)(tmpshwseg[1].At(i));
02110 scr=1+segu->GetHitLast()+segv->GetHitLast();
02111 if(scr>topscr){
02112 topscr=scr; id=i;
02113 }
02114 }
02115
02116 if(1+id>0){
02117 ShowerSegmentAtNu* segu = (ShowerSegmentAtNu*)(tmpshwseg[0].At(id));
02118 ShowerSegmentAtNu* segv = (ShowerSegmentAtNu*)(tmpshwseg[1].At(id));
02119 segu->SetPartner(segv); segv->SetPartner(segu);
02120 segu->SetUID(4); segv->SetUID(4);
02121 tmpshw[0].Add(segu); tmpshw[1].Add(segv);
02122 }
02123
02124 }
02125
02126 gettimeofday(&tpafter,0);
02127 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02128 fTime_ShowerAtNu+=dtime;
02129 }
02130
02131
02132 // ***********************************************
02133 // * F O R M A T I O N O F 3 D T R A C K S *
02134 // ***********************************************
02135
02136 gettimeofday(&tpbefore,0);
02137
02138 // Form tracks
02139 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of 3D tracks *** " << endl;
02140 Int_t km1,kp1;
02141 Double_t swx,swy,swx2,swxy,sw;
02142 Double_t w,x,y,m,c;
02143 Int_t ustrps,vstrps,trkstrps,trkplns;
02144 Int_t bplnvuw[2],eplnvuw[2];
02145 for(i=0;i<1+tmptrk[0].GetLast();i++){
02146
02147 MSG("AlgAtNuReco", Msg::kVerbose) << " making track " << i << endl;
02148
02149 // beg/end planes
02150 bpln2=-1; epln2=-1;
02151 for(vuw=0;vuw<2;vuw++){
02152 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmptrk[vuw].At(i));
02153 if(bpln2<0 || seg->GetBegPlane()<bpln2) bpln2=seg->GetBegPlane();
02154 if(epln2<0 || seg->GetEndPlane()>epln2) epln2=seg->GetEndPlane();
02155 }
02156
02157 bpln1=-1; epln1=-1;
02158 for(vuw=0;vuw<2;vuw++){
02159 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmptrk[vuw].At(i));
02160 if(bpln1<0 || seg->GetBegPlane()>bpln2) bpln1=seg->GetBegPlane();
02161 if(epln1<0 || seg->GetEndPlane()<epln2) epln1=seg->GetEndPlane();
02162 }
02163
02164 for(vuw=0;vuw<2;vuw++){
02165 bplnvuw[vuw]=-1; eplnvuw[vuw]=-1;
02166 }
02167
02168 MSG("AlgAtNuReco", Msg::kVerbose)
02169 << bpln2 << " (" << bpln1 << ") -> (" << epln1 << ") " << epln2 << endl;
02170
02171 // unpack segments
02172 npln=1+epln2-bpln2;
02173 Int_t* strlst = new Int_t[npln]; NEWOBJCTR+=npln;
02174 Int_t* vuwlst = new Int_t[npln]; NEWOBJCTR+=npln;
02175 TObjArray* clrlst = new TObjArray[npln]; NEWOBJCTR+=npln;
02176 TObjArray* hitlst = new TObjArray[npln]; NEWOBJCTR+=npln;
02177 for(pln=0;pln<npln;pln++){
02178 strlst[pln]=-1; vuwlst[pln]=-1;
02179 }
02180
02181 for(vuw=0;vuw<2;vuw++){
02182 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmptrk[vuw].At(i));
02183 for(j=0;j<1+seg->GetClrLast();j++){
02184 ClusterAtNu* clr = (ClusterAtNu*)(seg->GetClrAt(j));
02185 pln = clr->GetPlane()-bpln2;
02186 if( pln>=0 && pln<npln ){
02187 clrlst[pln].Add(clr);
02188 }
02189 }
02190 }
02191
02192 // trim segments
02193 bpln1=bpln1-4;
02194 km=-1; pln=bpln1-bpln2; if(pln<0) pln=0;
02195 while(km<0 && pln<npln){
02196 if(1+clrlst[pln].GetLast()>0) km=pln; else pln++;
02197 }
02198 if(km>-1 && km+4<npln && 1+clrlst[km+2].GetLast()==0){
02199 if(1+clrlst[km+4].GetLast()>0){
02200 ClusterAtNu* clr = (ClusterAtNu*)(clrlst[km+4].First());
02201 if(clr->GetTrkPlnFlag()>0 && clr->GetShwPlnFlag()==0) km=bpln1-bpln2;
02202 }
02203 bpln1=bpln2+km;
02204 }
02205
02206 epln1=epln1+4;
02207 kp=-1; pln=epln1-bpln2; if(pln>npln-1) pln=npln-1;
02208 while(kp<0 && pln>-1){
02209 if(1+clrlst[pln].GetLast()>0) kp=pln; else pln--;
02210 }
02211 if(kp>-1 && kp-4>-1 && 1+clrlst[kp-2].GetLast()==0){
02212 if(1+clrlst[kp-4].GetLast()>0){
02213 ClusterAtNu* clr = (ClusterAtNu*)(clrlst[kp-4].First());
02214 if(clr->GetTrkPlnFlag()>0 && clr->GetShwPlnFlag()==0) kp=epln1-bpln2;
02215 }
02216 epln1=bpln2+kp;
02217 }
02218
02219 ustrps=0; vstrps=0; trkstrps=0;
02220 for(pln=0;pln<npln;pln++){
02221 if( bpln2+pln>bpln1 && bpln2+pln<epln1
02222 && 1+clrlst[pln].GetLast()>0 ){
02223 ClusterAtNu* clr = (ClusterAtNu*)(clrlst[pln].First());
02224 vuw = clr->GetPlaneView();
02225 vuwlst[pln]=vuw; strlst[pln]=1;
02226 if(bplnvuw[vuw]<0||bpln2+pln<bplnvuw[vuw]) bplnvuw[vuw]=bpln2+pln;
02227 if(eplnvuw[vuw]<0||bpln2+pln>eplnvuw[vuw]) eplnvuw[vuw]=bpln2+pln;
02228 if(vuw==0) ustrps++; if(vuw==1) vstrps++;
02229 if(clr->GetTrkPlnFlag()>0) trkstrps++;
02230 }
02231 }
02232
02233 MSG("AlgAtNuReco", Msg::kVerbose)
02234 << " U: (" << bplnvuw[0] << "," << eplnvuw[0] << ") "
02235 << " V: (" << bplnvuw[1] << "," << eplnvuw[1] << ") " << endl;
02236
02237 if( ustrps>2 && vstrps>2 && trkstrps>2 ){
02238
02239 // sort track planes
02240 for(vuw=0;vuw<2;vuw++){
02241 for(pln=0;pln<npln;pln++){
02242 if(vuwlst[pln]==vuw && strlst[pln]>0){
02243
02244 ClusterAtNu* clr0 = (ClusterAtNu*)(clrlst[pln].First());
02245 k=0; km1=-1;
02246 while(pln-k>0 && km1<0){
02247 k++; if(vuwlst[pln-k]==vuw) km1=k;
02248 }
02249 k=0; kp1=-1;
02250 while(pln+k<npln-1 && kp1<0){
02251 k++; if(vuwlst[pln+k]==vuw) kp1=k;
02252 }
02253
02254 if(km1>0 && kp1>0){
02255 ClusterAtNu* clrp1 = (ClusterAtNu*)(clrlst[pln+kp1].First());
02256 ClusterAtNu* clrm1 = (ClusterAtNu*)(clrlst[pln-km1].First());
02257 if( ( clrm1->GetPlane()>bplnvuw[vuw]
02258 || ( clr0->GetPlane()-clrm1->GetPlane()<3
02259 && clr0->GetTrkPlnFlag()>0 ) )
02260 && ( clrp1->GetPlane()<eplnvuw[vuw]
02261 || ( clrp1->GetPlane()-clr0->GetPlane()<3
02262 && clr0->GetTrkPlnFlag()>0 ) ) ){
02263 if(clr0->IsTrkAssoc(clrm1,clrp1)>1){
02264 strlst[pln]=3;
02265 if( ( clrp1->GetBegStrip()-clr0->GetEndStrip()>-1
02266 && clr0->GetBegStrip()-clrm1->GetEndStrip()>-1 )
02267 || ( clrp1->GetEndStrip()-clr0->GetBegStrip()<1
02268 && clr0->GetEndStrip()-clrm1->GetBegStrip()<1 ) ){
02269 if( clr0->GetPlane()-clrm1->GetPlane()<3
02270 && clrm1->GetTrkPlnFlag()>-1 ){ strlst[pln-km1]=3; }
02271 if( clrp1->GetPlane()-clr0->GetPlane()<3
02272 && clrp1->GetTrkPlnFlag()>-1 ){ strlst[pln+kp1]=3; }
02273 }
02274 }
02275 else{
02276 if( 1+clr0->GetEndStrip()-clr0->GetBegStrip()<3
02277 && clr0->GetTrkPlnFlag()>0 ){ strlst[pln]=3; }
02278 }
02279 }
02280 }
02281
02282 }
02283 }
02284 }
02285
02286 MSG("AlgAtNuReco", Msg::kVerbose) << " *** SORT TRACK HITS *** " << endl;
02287 for(pln=0;pln<npln;pln++){
02288 if( vuwlst[pln]>-1 ){
02289 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrlst[pln].First());
02290 MSG("AlgAtNuReco", Msg::kVerbose) << " " << pln << " " << clrtmp->GetPlane()
02291 << " " << clrtmp->GetBegStrip() << " " << clrtmp->GetEndStrip()
02292 << " " << strlst[pln] << endl;
02293 }
02294 }
02295
02296 // perform global linear fits
02297 for(vuw=0;vuw<2;vuw++){
02298 for(pln=0;pln<npln;pln++){
02299 if(vuwlst[pln]==vuw && strlst[pln]<2){
02300
02301 MSG("AlgAtNuReco", Msg::kVerbose) << " *** GLOBAL FIT *** pln=" << pln << endl;
02302
02303 k=0; kp=-1;
02304 while(pln+k<epln1-bpln2 && pln+k<npln-1 && kp<0){
02305 k++; if(vuwlst[pln+k]==vuw && strlst[pln+k]>2) kp=k;
02306 }
02307
02308 k=0; km=-1;
02309 while(pln-k>bpln1-bpln2 && pln-k>0 && km<0){
02310 k++; if(vuwlst[pln-k]==vuw && strlst[pln-k]>2) km=k;
02311 }
02312
02313 if(km>0){ km1=0-km-7; k1=0-km; } else{ km1=(bpln1-bpln2)-pln+1; k1=(bpln1-bpln2)-pln+1; }
02314 if(pln+km1<0){ km1=0-pln; } if(pln+k1<0){ k1=0-pln; }
02315 if(kp>0){ kp1=1+kp+7; k2=1+kp; } else{ kp1=(epln1-bpln2)-pln; k2=(epln1-bpln2)-pln; }
02316 if(pln+kp1>npln){ kp1=npln-pln; } if(pln+k2>npln){ k2=npln-pln; }
02317
02318 ClusterAtNu* clrseed = (ClusterAtNu*)(clrlst[pln].First());
02319 m=0.0; c=0.5*(clrseed->GetBegStrip()+clrseed->GetEndStrip());
02320 swx=0.0; swy=0.0; swx2=0.0; swxy=0.0; sw=0.0;
02321 for(k=km1;k<kp1;k++){
02322 if(vuwlst[pln+k]==vuw){
02323 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrlst[pln+k].First());
02324 for(l=0;l<1+clrtmp->GetHitLast();l++){
02325 HitAtNu* hittmp = (HitAtNu*)(clrtmp->GetHitAt(l));
02326 x=hittmp->GetZPos(); y=hittmp->GetStrip();
02327 w=hittmp->GetCharge()/clrtmp->GetCharge();
02328 if(strlst[pln+k]>2) w=5.0*w; else w=0.5*w;
02329 swx+=w*x; swy+=w*y; swx2+=w*x*x; swxy+=w*x*y; sw+=w;
02330 }
02331 }
02332 }
02333 if(sw>0.0){
02334 m=(sw*swxy-swx*swy)/(sw*swx2-swx*swx);
02335 c=(swy*swx2-swx*swxy)/(sw*swx2-swx*swx);
02336 }
02337
02338 for(k=k1;k<k2;k++){
02339 if(vuwlst[pln+k]==vuw && strlst[pln+k]<2){
02340 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrlst[pln+k].First());
02341 strw = m*clrtmp->GetZPos()+c;
02342 if(m>=0.0) strdw = 0.6+0.02*m; if(m<=0) strdw = 0.6-0.02*m;
02343 MSG("AlgAtNuReco", Msg::kVerbose) << " pln=" << pln+k << " strw=" << strw << " strdw=" << strdw << " (" << clrtmp->GetBegStrip() << " " << clrtmp->GetEndStrip() << ") " << endl;
02344 if( (km>0 && kp>0) || (km<0 && kp<0)
02345 || (clrtmp->GetEndStrip()>strw-1.6 && clrtmp->GetBegStrip()<strw+1.6) ){
02346 id=-1; dscr=0.0; scr=999.9;
02347 for(l=0;l<1+clrtmp->GetHitLast();l++){
02348 HitAtNu* hit = (HitAtNu*)(clrtmp->GetHitAt(l));
02349 dscr=hit->GetStrip()-strw; if(dscr<0) dscr=-dscr;
02350 if(dscr<scr){
02351 id=l; scr=dscr;
02352 }
02353 }
02354 if(1+id>0){
02355 HitAtNu* hit = (HitAtNu*)(clrtmp->GetHitAt(id));
02356 ClusterAtNu* clrnew = new ClusterAtNu(hit); clrbnk.Add(clrnew); NEWOBJCTR++;
02357 clrlst[pln+k].Add(clrnew); strlst[pln+k]=2;
02358 for(l=0;l<1+clrtmp->GetHitLast();l++){
02359 HitAtNu* hittmp = (HitAtNu*)(clrtmp->GetHitAt(l));
02360 if( ( hittmp->GetStrip()-hit->GetStrip()>0
02361 && hittmp->GetStrip()-hit->GetStrip()<strdw)
02362 || ( hittmp->GetStrip()-hit->GetStrip()<0
02363 && hittmp->GetStrip()-hit->GetStrip()>-strdw) ){
02364 clrnew->AddHit(hittmp);
02365 }
02366 }
02367 for(l=0;l<1+clrnew->GetHitLast();l++){
02368 HitAtNu* hittmp = (HitAtNu*)(clrnew->GetHitAt(l));
02369 MSG("AlgAtNuReco", Msg::kVerbose) << " ... add hit=" << hittmp->GetStrip() << " (" << pln+k << ") " << endl;
02370 }
02371 }
02372 }
02373 }
02374 }
02375 }
02376
02377 }
02378 }
02379
02380 // perform local linear fits
02381 for(vuw=0;vuw<2;vuw++){
02382 for(pln=0;pln<npln;pln++){
02383 if(vuwlst[pln]==vuw && strlst[pln]>1){
02384
02385 MSG("AlgAtNuReco", Msg::kVerbose) << " *** FIT (2) *** pln=" << pln << endl;
02386
02387 ClusterAtNu* clr = (ClusterAtNu*)(clrlst[pln].Last());
02388 strw=0.5*(clr->GetBegStrip()+clr->GetEndStrip());
02389 strdw=0.6;
02390
02391 if(1+clr->GetHitLast()>1){
02392 for(k=0;k<1+clr->GetHitLast();k++){
02393 HitAtNu* hit = (HitAtNu*)(clr->GetHitAt(k));
02394 hit->SetTrkFlag(1);
02395 }
02396 swx=0.0; swy=0.0; swx2=0.0; swxy=0.0; sw=0.0;
02397 for(k=-4;k<5;k++){
02398 if( pln+k>=0 && pln+k<npln
02399 && vuwlst[pln+k]==vuw && strlst[pln+k]>1 ){
02400 ClusterAtNu* clrtmp = (ClusterAtNu*)(clrlst[pln+k].Last());
02401 for(l=0;l<1+clrtmp->GetHitLast();l++){
02402 HitAtNu* hittmp = (HitAtNu*)(clrtmp->GetHitAt(l));
02403 x=hittmp->GetZPos(); y=hittmp->GetStrip();
02404 w=hittmp->GetCharge()/clrtmp->GetCharge();
02405 swx+=w*x; swy+=w*y; swx2+=w*x*x; swxy+=w*x*y; sw+=w;
02406 }
02407 }
02408 }
02409 if(sw>1.0){
02410 m=(sw*swxy-swx*swy)/(sw*swx2-swx*swx);
02411 c=(swy*swx2-swx*swxy)/(sw*swx2-swx*swx);
02412 strw=m*clr->GetZPos()+c;
02413 if(m>=0.0) strdw = 0.6+0.02*m; if(m<=0) strdw = 0.6-0.02*m;
02414 }
02415 }
02416
02417 MSG("AlgAtNuReco", Msg::kVerbose) << " pln=" << pln << " strw=" << strw << " strdw=" << strdw << endl;
02418
02419 id=-1; dscr=0.0; scr=999.0;
02420 for(k=0;k<1+clr->GetHitLast();k++){
02421 HitAtNu* hit = (HitAtNu*)(clr->GetHitAt(k));
02422 dscr=hit->GetStrip()-strw; if(dscr<0) dscr=-dscr;
02423 if(dscr<scr){
02424 id=k; scr=dscr;
02425 }
02426 }
02427 if(1+id>0){
02428 HitAtNu* hit = (HitAtNu*)(clr->GetHitAt(id));
02429 hitlst[pln].Add(hit); strlst[pln]=4;
02430 for(k=0;k<1+clr->GetHitLast();k++){
02431 HitAtNu* hittmp = (HitAtNu*)(clr->GetHitAt(k));
02432 if( ( hittmp->GetStrip()-hit->GetStrip()>0
02433 && hittmp->GetStrip()-hit->GetStrip()<strdw)
02434 || ( hittmp->GetStrip()-hit->GetStrip()<0
02435 && hittmp->GetStrip()-hit->GetStrip()>-strdw) ){
02436 hitlst[pln].Add(hittmp);
02437 }
02438 }
02439 }
02440 for(l=0;l<1+hitlst[pln].GetLast();l++){
02441 HitAtNu* hittmp = (HitAtNu*)(hitlst[pln].At(l));
02442 MSG("AlgAtNuReco", Msg::kVerbose) << " ... add hit=" << hittmp->GetStrip() << " (" << pln << ") " << endl;
02443 }
02444 }
02445 }
02446 }
02447
02448 // form tracks
02449 trkplns=0;
02450 for(pln=0;pln<npln;pln++){
02451 if( strlst[pln]>3 ){
02452 trkplns++;
02453 }
02454 }
02455
02456 if(trkplns>5){
02457 for(vuw=0;vuw<2;vuw++){
02458 TrackSegmentAtNu* seg = (TrackSegmentAtNu*)(tmptrk[vuw].At(i));
02459 TrackAtNu* trk = new TrackAtNu(slice); trkbnk[vuw].Add(trk); NEWOBJCTR++;
02460 if(seg->GetReseedFlag()||fObjMaxFlag) jnflag2[vuw]=1; else jnflag2[vuw]=0;
02461 for(pln=0;pln<npln;pln++){
02462 if(vuwlst[pln]==vuw && strlst[pln]>3){
02463 for(k=0;k<1+hitlst[pln].GetLast();k++){
02464 HitAtNu* hit = (HitAtNu*)(hitlst[pln].At(k));
02465 hit->SetTrkFlag(2); trk->AddHit(hit);
02466 }
02467 }
02468 }
02469 }
02470 TrackAtNu* trku = (TrackAtNu*)(trkbnk[0].Last());
02471 TrackAtNu* trkv = (TrackAtNu*)(trkbnk[1].Last());
02472 trku->SetPartner(trkv); trkv->SetPartner(trku);
02473
02474 if(jnflag2[0]>0 && jnflag2[1]>0){
02475 trku->SetReseedFlag(1); trkv->SetReseedFlag(1);
02476 }
02477 }
02478
02479 }
02480
02481 delete [] strlst; DELOBJCTR+=npln;
02482 delete [] vuwlst; DELOBJCTR+=npln;
02483 delete [] clrlst; DELOBJCTR+=npln;
02484 delete [] hitlst; DELOBJCTR+=npln;
02485 }
02486
02487 // Print out list of tracks
02488 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF TRACKS *** " << endl;
02489 for(vuw=0;vuw<2;vuw++){
02490 for(i=0;i<1+trkbnk[vuw].GetLast();i++){
02491 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02492 MSG("AlgAtNuReco", Msg::kVerbose)
02493 << " " << vuw << " " << i
02494 << " " << trk->GetBegPlane() << " " << trk->GetBegStrip()
02495 << " " << trk->GetEndPlane() << " " << trk->GetEndStrip() << endl;
02496 }
02497 }
02498
02499 gettimeofday(&tpafter,0);
02500 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02501 fTime_TrackAtNu+=dtime;
02502
02503
02504 // *************************************************
02505 // * F O R M A T I O N O F 3 D S H O W E R S *
02506 // *************************************************
02507
02508 gettimeofday(&tpbefore,0);
02509
02510 // Form showers
02511 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of 3D showers *** " << endl;
02512 for(i=0;i<1+tmpshw[0].GetLast();i++){
02513
02514 MSG("AlgAtNuReco", Msg::kVerbose) << " making shower " << i << endl;
02515
02516 for(vuw=0;vuw<2;vuw++){
02517 ShowerSegmentAtNu* seg = (ShowerSegmentAtNu*)(tmpshw[vuw].At(i));
02518 ShowerAtNu* shw = new ShowerAtNu(slice); shwbnk[vuw].Add(shw); NEWOBJCTR++;
02519 if(seg->GetReseedFlag()||fObjMaxFlag) jnflag2[vuw]=1; else jnflag2[vuw]=0;
02520 for(j=0;j<1+seg->GetHitLast();j++){
02521 HitAtNu* hit = (HitAtNu*)(seg->GetHitAt(j));
02522 shw->AddHit(hit); hit->SetShwFlag(2);
02523 }
02524 }
02525 ShowerAtNu* shwu = (ShowerAtNu*)(shwbnk[0].Last());
02526 ShowerAtNu* shwv = (ShowerAtNu*)(shwbnk[1].Last());
02527 shwu->SetPartner(shwv); shwv->SetPartner(shwu);
02528
02529 if(jnflag2[0]>0 && jnflag2[1]>0){
02530 shwu->SetReseedFlag(1); shwv->SetReseedFlag(1);
02531 }
02532 }
02533
02534
02535 // *********************************************************
02536 // * F O R M A T I O N O F V E R T E X S H O W E R S *
02537 // *********************************************************
02538
02539 if(fVtxFlag){
02540 MSG("AlgAtNuReco", Msg::kDebug) << " *** formation of vertex showers *** " << endl;
02541
02542 Int_t shwflag=0;
02543 for(i=0;i<1+shwbnk[0].GetLast();i++){
02544 if(shwflag==0){
02545 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[0].At(i));
02546 if(shw->GetReseedFlag()||1) shwflag=1;
02547 }
02548 }
02549
02550 TObjArray tmphittrk[2],tmphitbeg[2],tmphitend[2];
02551 for(i=0;i<1+trkbnk[0].GetLast();i++){
02552
02553 // beginning vertex shower (existing)
02554 for(vuw=0;vuw<2;vuw++){
02555 tmpclr[vuw].Clear();
02556 tmphittrk[vuw].Clear(); tmphitbeg[vuw].Clear();
02557 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02558 for(j=0;j<1+trk->GetHitLast();j++){
02559 HitAtNu* hit = (HitAtNu*)(trk->GetHitAt(j));
02560 if( hit->GetPlane()<trk->GetBegPlane()+9
02561 && 2*hit->GetPlane()<trk->GetBegPlane()+trk->GetEndPlane()+1
02562 && (Int_t)((hit->GetPlane())/(modpln+1))==(Int_t)((trk->GetBegPlane())/(modpln+1)) ){
02563 tmphittrk[vuw].Add(hit);
02564 }
02565 }
02566 }
02567
02568 for(j=0;j<1+shwbnk[0].GetLast();j++){
02569 for(vuw=0;vuw<2;vuw++){
02570 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02571 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[vuw].At(j));
02572 jnflag2[vuw]=0;
02573 if( shw->GetBegPlane()<trk->GetBegPlane()+9
02574 && 2*shw->GetBegPlane()<trk->GetBegPlane()+trk->GetEndPlane()+1
02575 && (Int_t)((shw->GetBegPlane())/(modpln+1))==(Int_t)((trk->GetBegPlane())/(modpln+1)) ){
02576 for(k=0;k<1+tmphittrk[vuw].GetLast();k++){
02577 if(jnflag2[vuw]==0){
02578 HitAtNu* hit = (HitAtNu*)(tmphittrk[vuw].At(k));
02579 if(shw->IsDiffuseShwAssoc(hit)>0) jnflag2[vuw]=1;
02580 }
02581 }
02582 }
02583 }
02584 if(jnflag2[0]>0 || jnflag2[1]>0){
02585 for(vuw=0;vuw<2;vuw++){
02586 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[vuw].At(j));
02587 tmpclr[vuw].Add(shw);
02588 }
02589 }
02590 }
02591
02592 if( 1+tmpclr[0].GetLast()>0 && 1+tmpclr[1].GetLast()>0 ){
02593 id=-1; topscr=9999.9;
02594 for(j=0;j<1+tmpclr[0].GetLast();j++){
02595 ShowerAtNu* shwu = (ShowerAtNu*)(tmpclr[0].At(j));
02596 ShowerAtNu* shwv = (ShowerAtNu*)(tmpclr[1].At(j));
02597 scr = shwu->GetBegPlane()+shwv->GetBegPlane();
02598 if(scr<topscr){
02599 topscr=scr; id=j;
02600 }
02601 }
02602 if(1+id>0){
02603 MSG("AlgAtNuReco", Msg::kVerbose) << " found existing beginning vertex shower " << endl;
02604 for(vuw=0;vuw<2;vuw++){
02605 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02606 ShowerAtNu* shw = (ShowerAtNu*)(tmpclr[vuw].At(id));
02607 trk->SetBegVtxShower(shw); shw->AddAssocTrack(trk); shw->SetBegVtxFlag(1);
02608 }
02609 }
02610 }
02611
02612 // beginning vertex shower (new)
02613 else{
02614 if( shwflag==0 && i==0 ){
02615 for(vuw=0;vuw<2;vuw++){
02616 tmphitbeg[vuw].Clear();
02617 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02618
02619 mod=(Int_t)((trk->GetBegPlane())/(modpln+1));
02620 for(j=-3;j<4;j++){
02621 pln=trk->GetBegPlane()+2*j;
02622 if( 2*pln<trk->GetBegPlane()+trk->GetEndPlane()+1
02623 && pln>mod*(modpln+1) && pln<(mod+1)*(modpln+1) ){
02624 for(k=0;k<1+hitplnbnk[pln].GetLast();k++){
02625 HitAtNu* hit = (HitAtNu*)(hitplnbnk[pln].At(k));
02626 if( hit->GetShwFlag()<2
02627 && ( hit->GetTrkFlag()<2 || hit->GetCharge()>50 ) ){
02628 jnflag=0;
02629 for(l=0;l<1+tmphittrk[vuw].GetLast();l++){
02630 if(jnflag==0){
02631 HitAtNu* hittrk = (HitAtNu*)(tmphittrk[vuw].At(l));
02632 if( hittrk->GetPlane()<trk->GetBegPlane()+7
02633 && hittrk->IsDiffuseShwAssoc(hit)>0 ){
02634 tmphitbeg[vuw].Add(hit); hit->SetShwFlag(1); jnflag=1;
02635 }
02636 }
02637 }
02638 }
02639 }
02640 }
02641 }
02642
02643 if(1+tmphitbeg[vuw].GetLast()>0){
02644 bpln=trk->GetBegPlane()-6; epln=trk->GetBegPlane()+6;
02645 ctr=1;
02646 while(ctr>0){
02647 ctr=0;
02648 pln0=bpln-4; npln=5+(epln-bpln)/2;
02649 for(j=0;j<npln;j++){
02650 pln=pln0+2*j;
02651 if(pln>mod*(modpln+1) && pln<(mod+1)*(modpln+1)){
02652 for(k=0;k<1+hitplnbnk[pln].GetLast();k++){
02653 HitAtNu* hit = (HitAtNu*)(hitplnbnk[pln].At(k));
02654 if(hit->GetShwFlag()<1 && hit->GetTrkFlag()<1){
02655 jnflag=0;
02656 for(l=0;l<1+tmphitbeg[vuw].GetLast();l++){
02657 if(jnflag==0){
02658 HitAtNu* hitshw = (HitAtNu*)(tmphitbeg[vuw].At(l));
02659 if( hitshw->IsDiffuseShwAssoc(hit)>0 ){
02660 tmphitbeg[vuw].Add(hit); hit->SetShwFlag(1);
02661 if(hit->GetPlane()<bpln) bpln=hit->GetPlane();
02662 if(hit->GetPlane()>epln) epln=hit->GetPlane();
02663 jnflag=1; ctr++;
02664 }
02665 }
02666 }
02667 }
02668 }
02669 }
02670 }
02671 }
02672 for(j=0;j<1+tmphitbeg[vuw].GetLast();j++){
02673 HitAtNu* hit = (HitAtNu*)(tmphitbeg[vuw].At(j));
02674 hit->SetShwFlag(0);
02675 }
02676 }
02677
02678 }
02679 }
02680 }
02681
02682 // end vertex shower (existing)
02683 for(vuw=0;vuw<2;vuw++){
02684 tmpclr[vuw].Clear();
02685 tmphittrk[vuw].Clear(); tmphitend[vuw].Clear();
02686 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02687 for(j=0;j<1+trk->GetHitLast();j++){
02688 HitAtNu* hit = (HitAtNu*)(trk->GetHitAt(j));
02689 if( hit->GetPlane()>trk->GetEndPlane()-9
02690 && 2*hit->GetPlane()>trk->GetBegPlane()+trk->GetEndPlane()-1
02691 && (Int_t)((hit->GetPlane())/(modpln+1))==(Int_t)((trk->GetEndPlane())/(modpln+1)) ){
02692 tmphittrk[vuw].Add(hit);
02693 }
02694 }
02695 }
02696
02697 for(j=0;j<1+shwbnk[0].GetLast();j++){
02698 for(vuw=0;vuw<2;vuw++){
02699 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02700 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[vuw].At(j));
02701 jnflag2[vuw]=0;
02702 if( shw->GetEndPlane()>trk->GetEndPlane()-9
02703 && 2*shw->GetEndPlane()>trk->GetBegPlane()+trk->GetEndPlane()-1
02704 && (Int_t)((shw->GetEndPlane())/(modpln+1))==(Int_t)((trk->GetEndPlane())/(modpln+1)) ){
02705 for(k=0;k<1+tmphittrk[vuw].GetLast();k++){
02706 if(jnflag2[vuw]==0){
02707 HitAtNu* hit = (HitAtNu*)(tmphittrk[vuw].At(k));
02708 if(shw->IsDiffuseShwAssoc(hit)>0) jnflag2[vuw]=1;
02709 }
02710 }
02711 }
02712 }
02713 if(jnflag2[0]>0 || jnflag2[1]>0){
02714 for(vuw=0;vuw<2;vuw++){
02715 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[vuw].At(j));
02716 tmpclr[vuw].Add(shw);
02717 }
02718 }
02719 }
02720
02721 if( 1+tmpclr[0].GetLast()>0 && 1+tmpclr[1].GetLast()>0 ){
02722 id=-1; topscr=-9999.9;
02723 for(j=0;j<1+tmpclr[0].GetLast();j++){
02724 ShowerAtNu* shwu = (ShowerAtNu*)(tmpclr[0].At(j));
02725 ShowerAtNu* shwv = (ShowerAtNu*)(tmpclr[1].At(j));
02726 scr = shwu->GetEndPlane()+shwv->GetEndPlane();
02727 if(scr>topscr){
02728 topscr=scr; id=j;
02729 }
02730 }
02731 if(1+id>0){
02732 MSG("AlgAtNuReco", Msg::kVerbose) << " found existing end vertex shower " << endl;
02733 for(vuw=0;vuw<2;vuw++){
02734 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02735 ShowerAtNu* shw = (ShowerAtNu*)(tmpclr[vuw].At(id));
02736 trk->SetEndVtxShower(shw); shw->AddAssocTrack(trk); shw->SetEndVtxFlag(1);
02737 }
02738 }
02739 }
02740
02741 // end vertex shower (new)
02742 else{
02743 if( shwflag==0 && i==0 ){
02744 for(vuw=0;vuw<2;vuw++){
02745 tmphitend[vuw].Clear();
02746 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02747
02748 mod=(Int_t)((trk->GetEndPlane())/(modpln+1));
02749 for(j=-3;j<4;j++){
02750 pln=trk->GetEndPlane()+2*j;
02751 if( 2*pln>trk->GetBegPlane()+trk->GetEndPlane()-1
02752 && pln>mod*(modpln+1) && pln<(mod+1)*(modpln+1) ){
02753 for(k=0;k<1+hitplnbnk[pln].GetLast();k++){
02754 HitAtNu* hit = (HitAtNu*)(hitplnbnk[pln].At(k));
02755 if( hit->GetShwFlag()<2
02756 && ( hit->GetTrkFlag()<2 || hit->GetCharge()>50 ) ){
02757 jnflag=0;
02758 for(l=0;l<1+tmphittrk[vuw].GetLast();l++){
02759 if(jnflag==0){
02760 HitAtNu* hittrk = (HitAtNu*)(tmphittrk[vuw].At(l));
02761 if( hittrk->GetPlane()>trk->GetEndPlane()-7
02762 && hittrk->IsDiffuseShwAssoc(hit)>0 ){
02763 tmphitend[vuw].Add(hit); hit->SetShwFlag(1); jnflag=1;
02764 }
02765 }
02766 }
02767 }
02768 }
02769 }
02770 }
02771
02772 if(1+tmphitend[vuw].GetLast()>0){
02773 bpln=trk->GetEndPlane()-6; epln=trk->GetEndPlane()+6;
02774 ctr=1;
02775 while(ctr>0){
02776 ctr=0;
02777 pln0=bpln-4; npln=5+(epln-bpln)/2;
02778 for(j=0;j<npln;j++){
02779 pln=pln0+2*j;
02780 if(pln>mod*(modpln+1) && pln<(mod+1)*(modpln+1)){
02781 for(k=0;k<1+hitplnbnk[pln].GetLast();k++){
02782 HitAtNu* hit = (HitAtNu*)(hitplnbnk[pln].At(k));
02783 if(hit->GetShwFlag()<1 && hit->GetTrkFlag()<1){
02784 jnflag=0;
02785 for(l=0;l<1+tmphitend[vuw].GetLast();l++){
02786 if(jnflag==0){
02787 HitAtNu* hitshw = (HitAtNu*)(tmphitend[vuw].At(l));
02788 if( hitshw->IsDiffuseShwAssoc(hit)>0 ){
02789 tmphitend[vuw].Add(hit); hit->SetShwFlag(1);
02790 if(hit->GetPlane()<bpln) bpln=hit->GetPlane();
02791 if(hit->GetPlane()>epln) epln=hit->GetPlane();
02792 jnflag=1; ctr++;
02793 }
02794 }
02795 }
02796 }
02797 }
02798 }
02799 }
02800 }
02801 for(j=0;j<1+tmphitend[vuw].GetLast();j++){
02802 HitAtNu* hit = (HitAtNu*)(tmphitend[vuw].At(j));
02803 hit->SetShwFlag(0);
02804 }
02805 }
02806
02807 }
02808 }
02809 }
02810
02811 if( 1+tmphitbeg[0].GetLast()>0 && 1+tmphitbeg[1].GetLast()>0 ){
02812 MSG("AlgAtNuReco", Msg::kVerbose) << " found new vertex shower at beginning of track " << endl;
02813 for(vuw=0;vuw<2;vuw++){
02814 ShowerAtNu* shw = new ShowerAtNu(slice); NEWOBJCTR++;
02815 for(j=0;j<1+tmphitbeg[vuw].GetLast();j++){
02816 HitAtNu* hit = (HitAtNu*)(tmphitbeg[vuw].At(j));
02817 shw->AddHit(hit); hit->SetShwFlag(2);
02818 }
02819 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02820 trk->SetBegVtxShower(shw); shw->AddAssocTrack(trk);
02821 shw->SetBegVtxFlag(1); shw->SetReseedFlag(1);
02822 shwbnk[vuw].Add(shw);
02823 }
02824 }
02825
02826 if( 1+tmphitend[0].GetLast()>0 && 1+tmphitend[1].GetLast()>0 ){
02827 MSG("AlgAtNuReco", Msg::kVerbose) << " found new vertex shower at end of track " << endl;
02828 for(vuw=0;vuw<2;vuw++){
02829 ShowerAtNu* shw = new ShowerAtNu(slice); NEWOBJCTR++;
02830 for(j=0;j<1+tmphitend[vuw].GetLast();j++){
02831 HitAtNu* hit = (HitAtNu*)(tmphitend[vuw].At(j));
02832 shw->AddHit(hit); hit->SetShwFlag(2);
02833 }
02834 TrackAtNu* trk = (TrackAtNu*)(trkbnk[vuw].At(i));
02835 trk->SetEndVtxShower(shw); shw->AddAssocTrack(trk);
02836 shw->SetEndVtxFlag(1); shw->SetReseedFlag(1);
02837 shwbnk[vuw].Add(shw);
02838 }
02839 }
02840
02841 }
02842 }
02843
02844
02845 // *********************************************
02846 // * M O P U P R E M A I N I N G H I T S *
02847 // *********************************************
02848
02849 // Mopping up remaining hits
02850 MSG("AlgAtNuReco", Msg::kDebug) << " *** mop up remaining hits *** " << endl;
02851 Double_t shwwinplns,shwwinstrps;
02852 if( 1+shwbnk[0].GetLast()>0 ){
02853 shwwinplns=8; shwwinstrps=12;
02854 if(1+trkbnk[0].GetLast()==0){
02855 shwwinplns=16; shwwinstrps=25; }
02856 for(mod=0;mod<modnum;mod++){
02857 for(pln=1;pln<1+modpln;pln++){
02858 i=pln+mod*(modpln+1);
02859 for(j=0;j<1+hitplnbnk[i].GetLast();j++){
02860 HitAtNu* hit = (HitAtNu*)(hitplnbnk[i].At(j));
02861 if(hit->GetShwFlag()<1 && hit->GetTrkFlag()<1){
02862 vuw = hit->GetPlaneView();
02863 id=-1; scr=shwwinplns+shwwinstrps;
02864 for(k=0;k<1+shwbnk[vuw].GetLast();k++){
02865 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[vuw].At(k));
02866 if( (Int_t)((shw->GetEndPlane())/(modpln+1)) == mod ){
02867 for(l=0;l<1+shw->GetHitLast();l++){
02868 if(1+shwbnk[vuw].GetLast()>1 || id<0){
02869 HitAtNu* hittmp = (HitAtNu*)(shw->GetHitAt(l));
02870 if( (hittmp->GetShwFlag()>1)
02871 && (hit->GetTrkFlag()<1||hittmp->IsShwAssoc(hit)>0) ){
02872 dpln=hit->GetPlane()-hittmp->GetPlane(); if(dpln<0) dpln=-dpln;
02873 dstr=hit->GetStrip()-hittmp->GetStrip(); if(dstr<0) dstr=-dstr;
02874 if( dpln<shwwinplns && dstr<shwwinstrps && dpln+dstr<scr ){
02875 id=k; scr=dpln+dstr;
02876 }
02877 }
02878 }
02879 }
02880 }
02881 }
02882 if(1+id>0){
02883 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[vuw].At(id));
02884 shw->AddHit(hit); hit->SetShwFlag(1);
02885 }
02886 }
02887 }
02888 }
02889 }
02890 }
02891
02892 // Print list of showers
02893 MSG("AlgAtNuReco", Msg::kVerbose) << " *** LIST OF SHOWERS *** " << endl;
02894 for(vuw=0;vuw<2;vuw++){
02895 MSG("AlgAtNuReco", Msg::kVerbose) << "view: " << vuw << endl;
02896 for(i=0;i<1+shwbnk[vuw].GetLast();i++){
02897 ShowerAtNu* shw = (ShowerAtNu*)(shwbnk[vuw].At(i));
02898 MSG("AlgAtNuReco", Msg::kVerbose)
02899 << " i=" << i << " "
02900 << shw->GetBegPlane() << " " << shw->GetBegStrip() << " "
02901 << shw->GetEndPlane() << " " << shw->GetEndStrip() << endl;
02902 }
02903 }
02904
02905 gettimeofday(&tpafter,0);
02906 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02907 fTime_ShowerAtNu+=dtime;
02908
02909 // ***************************************************
02910 // * D E L E T E T E M P O R A R Y O B J E C T S *
02911 // ***************************************************
02912
02913 // clear/delete hit/cluster objects
02914 gettimeofday(&tpbefore,0);
02915 for(i=0;i<fEndPln;i++){
02916 hitplnbnk[i].Clear(); // hit plane bank
02917 xtalkplnbnk[i].Clear();// xtalk plane bank
02918 clrplnbnk[i].Clear(); // cluster plane bank
02919 }
02920 DELOBJCTR+=1+clrbnk.GetLast(); clrbnk.Delete(); // cluster bank
02921 gettimeofday(&tpafter,0);
02922 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02923 fTime_TrackAtNu+=0.5*dtime;
02924 fTime_ShowerAtNu+=0.5*dtime;
02925
02926 // clear/delete track segment objects
02927 gettimeofday(&tpbefore,0);
02928 for(i=0;i<fEndPln;i++){
02929 segplnbnk[i].Clear(); // track segment plane bank
02930 }
02931 DELOBJCTR+=1+trksegbnk.GetLast(); trksegbnk.Delete(); // track segment bank
02932 gettimeofday(&tpafter,0);
02933 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02934 fTime_TrackAtNu+=dtime;
02935
02936 // clear/delete shower segment objects
02937 gettimeofday(&tpbefore,0);
02938 DELOBJCTR+=1+shwsegbnk.GetLast(); shwsegbnk.Delete(); // shower segment bank
02939 gettimeofday(&tpafter,0);
02940 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02941 fTime_ShowerAtNu+=dtime;
02942
02943 } // iteration over slices
02944
02945
02946 // *********************************************************
02947 // * F O R M A T I O N O F C A N D X X X H A N D L E S *
02948 // *********************************************************
02949 MSG("AlgAtNuReco", Msg::kDebug) << " Formation of CandXXXHandles " << endl;
02950
02951 // create Singleton Instance of AlgFactory
02952 AlgFactory &af = AlgFactory::GetInstance();
02953
02954 // create CandTrackAtNuListHandle
02955 gettimeofday(&tpbefore,0);
02956
02957 TObjArray* objtrk = 0;
02958 for(itr=0;itr<1+trkbnk[0].GetLast();itr++){
02959 TrackAtNu* trku = (TrackAtNu*)(trkbnk[0].At(itr));
02960 TrackAtNu* trkv = (TrackAtNu*)(trkbnk[1].At(itr));
02961 objtrk = new TObjArray(); NEWOBJCTR++;
02962 objtrk->Add(trku);
02963 objtrk->Add(trkv);
02964 trkcxt.Add(objtrk);
02965 }
02966
02967 MSG("AlgAtNuReco", Msg::kDebug) << "*** CREATE CandTrackAtNuListHandle *** " << endl;
02968 AlgHandle ah_trk = af.GetAlgHandle("AlgTrackAtNuList", "default");
02969 CandContext cx_trk(this, cx.GetMom());
02970 cx_trk.SetCandRecord(cx.GetCandRecord());
02971 cx_trk.SetDataIn(&trkcxt);
02972 CandTrackAtNuListHandle trk_list = CandTrackAtNuList::MakeCandidate(ah_trk, cx_trk);
02973 trk_list.SetName(fListOutTrk.Data());
02974 trk_list.SetTitle(TString("Created by AtNuFindModule from ").Append(slice_list->GetName()));
02975
02976 TIter trkitr(trk_list.GetDaughterIterator());
02977 while(CandTrackAtNuHandle* trk = dynamic_cast<CandTrackAtNuHandle*>(trkitr())){
02978 atnu_handle.AddTrack(trk);
02979 }
02980
02981 gettimeofday(&tpafter,0);
02982 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
02983 fTime_TrackAtNu+=dtime;
02984
02985 // create CandShowerAtNuListHandle
02986 gettimeofday(&tpbefore,0);
02987
02988 TObjArray* objshw = 0;
02989 for(itr=0;itr<1+shwbnk[0].GetLast();itr++){
02990 ShowerAtNu* shwu = (ShowerAtNu*)(shwbnk[0].At(itr));
02991 ShowerAtNu* shwv = (ShowerAtNu*)(shwbnk[1].At(itr));
02992 objshw = new TObjArray(); NEWOBJCTR++;
02993 objshw->Add(shwu);
02994 objshw->Add(shwv);
02995 shwcxt.Add(objshw);
02996 }
02997
02998 MSG("AlgAtNuReco", Msg::kDebug) << "*** CREATE CandShowerAtNuListHandle *** " << endl;
02999 AlgHandle ah_shw = af.GetAlgHandle("AlgShowerAtNuList", "default");
03000 CandContext cx_shw(this, cx.GetMom());
03001 cx_shw.SetCandRecord(cx.GetCandRecord());
03002 cx_shw.SetDataIn(&shwcxt);
03003 CandShowerAtNuListHandle shw_list = CandShowerAtNuList::MakeCandidate(ah_shw, cx_shw);
03004 shw_list.SetName(fListOutShw.Data());
03005 shw_list.SetTitle(TString("Created by AtNuFindModule from ").Append(slice_list->GetName()));
03006
03007 TIter shwitr(shw_list.GetDaughterIterator());
03008 while(CandShowerAtNuHandle* shw = dynamic_cast<CandShowerAtNuHandle*>(shwitr())){
03009 atnu_handle.AddShower(shw);
03010 }
03011
03012 gettimeofday(&tpafter,0);
03013 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
03014 fTime_ShowerAtNu+=dtime;
03015
03016 // ***************************************************
03017 // * D E L E T E R E M A I N I N G O B J E C T S *
03018 // ***************************************************
03019
03020 // delete hit objects
03021 gettimeofday(&tpbefore,0);
03022 DELOBJCTR+=1+hitbnk.GetLast(); hitbnk.Delete(); // hit bank
03023 gettimeofday(&tpafter,0);
03024 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
03025 fTime_TrackAtNu+=0.5*dtime;
03026 fTime_ShowerAtNu+=0.5*dtime;
03027
03028 // delete track objects
03029 gettimeofday(&tpbefore,0);
03030 for(itr=0;itr<2;itr++){
03031 DELOBJCTR+=1+trkbnk[itr].GetLast(); trkbnk[itr].Delete(); // track bank
03032 }
03033 DELOBJCTR+=1+trkcxt.GetLast(); trkcxt.Delete(); // candcontext for tracks
03034 gettimeofday(&tpafter,0);
03035 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
03036 fTime_TrackAtNu+=dtime;
03037
03038 // delete shower objects
03039 gettimeofday(&tpbefore,0);
03040 for(itr=0;itr<2;itr++){
03041 DELOBJCTR+=1+shwbnk[itr].GetLast(); shwbnk[itr].Delete(); // shower bank
03042 }
03043 DELOBJCTR+=1+shwcxt.GetLast(); shwcxt.Delete(); // candcontext for showers
03044 gettimeofday(&tpafter,0);
03045 dtime=1000.0*(tpafter.tv_sec-tpbefore.tv_sec)+0.001*(tpafter.tv_usec-tpbefore.tv_usec);
03046 fTime_ShowerAtNu+=dtime;
03047
03048
03049 // *************************************************
03050 // * S E T A T N U R E C O P A R A M E T E R S *
03051 // *************************************************
03052
03053 MSG("AlgAtNuReco", Msg::kDebug) << " *** CANDATNURECO PARAMETERS *** " << endl;
03054
03055
03056 Double_t atnurecoscore=0.0;
03057
03058
03059 atnu_handle.SetAtNuRecoScore(atnurecoscore);
03060
03061 gettimeofday(&atnuafter,0);
03062 dtime=1000.0*(atnuafter.tv_sec-atnubefore.tv_sec)+0.001*(atnuafter.tv_usec-atnubefore.tv_usec);
03063 fTime_AtNuReco+=dtime;
03064
03065
03066 // ***************************************************************
03067 // * G I V E C A N D H A N D L E S T O C A N D R E C O R D *
03068 // ***************************************************************
03069
03070 MSG("AlgAtNuReco", Msg::kDebug) << " *** MEMORY BOOK-KEEPING : "
03071 << " new=" << NEWOBJCTR
03072 << " delete=" << DELOBJCTR << endl;
03073
03074 atnu_handle.SetNBlocks(NEWOBJCTR);
03075
03076 trk_list.SetCPUTime(fTime_TrackAtNu);
03077 MSG("AlgAtNuReco", Msg::kDebug) << " *** TRACK RECO TIME = " << fTime_TrackAtNu << endl;
03078
03079 shw_list.SetCPUTime(fTime_ShowerAtNu);
03080 MSG("AlgAtNuReco", Msg::kDebug) << " *** SHOWER RECO TIME = " << fTime_ShowerAtNu << endl;
03081
03082 atnu_handle.SetCPUTime(fTime_AtNuReco);
03083 MSG("AlgAtNuReco", Msg::kDebug) << " *** TOTAL ATNURECO TIME = " << fTime_AtNuReco << endl;
03084
03085
03086 CandRecord* candrec = (CandRecord*)(cx.GetCandRecord());
03087 candrec->SecureCandHandle(trk_list);
03088 candrec->SecureCandHandle(shw_list);
03089
03090 return;
03091
03092 }
|
|
|
Reimplemented from AlgBase. Definition at line 3094 of file AlgAtNuReco.cxx. 03095 {
03096
03097 }
|
|
|
Definition at line 25 of file AlgAtNuReco.h. Referenced by RunAlg(). |
|
|
Definition at line 21 of file AlgAtNuReco.h. Referenced by RunAlg(). |
|
|
Definition at line 24 of file AlgAtNuReco.h. Referenced by RunAlg(). |
|
|
Definition at line 19 of file AlgAtNuReco.h. Referenced by RunAlg(). |
|
|
Definition at line 22 of file AlgAtNuReco.h. Referenced by RunAlg(). |
|
|
Definition at line 30 of file AlgAtNuReco.h. Referenced by RunAlg(). |
|
|
Definition at line 32 of file AlgAtNuReco.h. Referenced by RunAlg(). |
|
|
Definition at line 27 of file AlgAtNuReco.h. Referenced by RunAlg(). |
|
|
Definition at line 29 of file AlgAtNuReco.h. Referenced by RunAlg(). |
|
|
Definition at line 31 of file AlgAtNuReco.h. Referenced by RunAlg(). |
|
|
Definition at line 26 of file AlgAtNuReco.h. Referenced by RunAlg(). |
|
|
Definition at line 20 of file AlgAtNuReco.h. Referenced by RunAlg(). |
1.3.9.1