00001 #include "AlgTrackAtNu.h"
00002 #include "CandTrackAtNuHandle.h"
00003
00004 #include "Algorithm/AlgConfig.h"
00005 #include "Algorithm/AlgFactory.h"
00006 #include "Algorithm/AlgHandle.h"
00007
00008 #include "MessageService/MsgService.h"
00009 #include "JobControl/JobCModuleRegistry.h"
00010 #include "Calibrator/Calibrator.h"
00011 #include "RecoBase/CandStripHandle.h"
00012 #include "CandDigit/CandDeMuxDigitHandle.h"
00013 #include "CandDigit/CandDeMuxDigitListHandle.h"
00014 #include "Candidate/CandContext.h"
00015
00016 #include "UgliGeometry/UgliGeomHandle.h"
00017 #include "Validity/VldContext.h"
00018 #include "Plex/PlexPlaneId.h"
00019
00020 #include "HitAtNu.h"
00021 #include "ShowerAtNu.h"
00022 #include "TrackAtNu.h"
00023
00024 #include "TObjArray.h"
00025
00026
00027
00028
00029
00030 ClassImp(AlgTrackAtNu)
00031
00032 CVSID("$Id: AlgTrackAtNu.cxx,v 1.12 2006/05/22 16:44:42 rhatcher Exp $");
00033
00034 AlgTrackAtNu::AlgTrackAtNu()
00035 {
00036
00037 }
00038
00039 AlgTrackAtNu::~AlgTrackAtNu()
00040 {
00041
00042 }
00043
00044 void AlgTrackAtNu::RunAlg(AlgConfig& ac, CandHandle &ch, CandContext &cx)
00045 {
00046
00047 MSG("AlgTrackAtNu", Msg::kDebug) << "AlgTrackAtNu::RunAlg(...)" << endl;
00048
00049
00050 CandTrackAtNuHandle& track = dynamic_cast<CandTrackAtNuHandle&>(ch);
00051
00052
00053 Double_t fFibreIndex;
00054 Int_t fAtNuAnaOnOff;
00055 fFibreIndex = ac.GetDouble("FibreIndex");
00056 fAtNuAnaOnOff = ac.GetInt("AtNuAnaOnOff");
00057 MSG("AlgTrackAtNu", Msg::kDebug) << " FibreIndex=" << fFibreIndex << endl;
00058 MSG("AlgTrackAtNu", Msg::kDebug) << " AtNuAnaOnOff=" << fAtNuAnaOnOff << endl;
00059
00060
00061 const TObjArray* arr = dynamic_cast<const TObjArray*>(cx.GetDataIn());
00062 TrackAtNu* trku = (TrackAtNu*)(arr->At(0));
00063 TrackAtNu* trkv = (TrackAtNu*)(arr->At(1));
00064 track.SetCandSlice(trku->GetCandSliceHandle());
00065
00066 CandRecord* candrec = (CandRecord*)(cx.GetCandRecord());
00067 VldContext *vldc = (VldContext*)(candrec->GetVldContext());
00068 UgliGeomHandle ugh(*vldc);
00069
00070
00071 Calibrator& cal = Calibrator::Instance();
00072 cal.ReInitialise(*vldc);
00073
00074
00075 Int_t atmosflag,modtype,modnum,modpln,modstr;
00076 switch(candrec->GetVldContext()->GetDetector()){
00077 case Detector::kFar:
00078 atmosflag=1; modtype=1; modnum=2; modpln=248; modstr=192; break;
00079 case Detector::kNear:
00080 atmosflag=0; modtype=2; modnum=1; modpln=282; modstr=96; break;
00081 case Detector::kCalib:
00082 atmosflag=0; modtype=3; modnum=1; modpln=60; modstr=24; break;
00083 default:
00084 atmosflag=0; modtype=-1; modnum=0; modpln=0; modstr=0; break;
00085 }
00086
00087 Int_t i,j,k,l;
00088 Int_t pln,bpln,epln,npln;
00089 Double_t begZsm1=-9999.9,endZsm1=-9999.9;
00090 Double_t begZsm2=-9999.9,endZsm2=-9999.9;
00091 Int_t begsm1=-1,endsm1=-1,begsm2=-1,endsm2=-1;
00092
00093 for(j=0;j<486;j++){
00094 PlexPlaneId myplaneid(vldc->GetDetector(),j,1);
00095 UgliSteelPlnHandle myplanehandle = ugh.GetSteelPlnHandle(myplaneid);
00096 if(myplanehandle.IsValid()){
00097 if(j<=modpln && begZsm1<0.0){ begsm1=j; begZsm1=myplanehandle.GetZ0(); }
00098 if(j<=modpln && myplanehandle.GetZ0()>endZsm1){ endsm1=j; endZsm1=myplanehandle.GetZ0(); }
00099 if(j>modpln && begZsm2<0.0){ begsm2=j; begZsm2=myplanehandle.GetZ0(); }
00100 if(j>modpln && myplanehandle.GetZ0()>endZsm2){ endsm2=j; endZsm2=myplanehandle.GetZ0(); }
00101 }
00102 }
00103
00104 MSG("AlgTrackAtNu", Msg::kDebug) << " *** detector geometry *** " << endl
00105 << " SM1 : pln " << begsm1 << " -> " << endsm1
00106 << " Z " << begZsm1 << " -> " << endZsm1 << endl
00107 << " SM2 : pln " << begsm2 << " -> " << endsm2
00108 << " Z " << begZsm2 << " -> " << endZsm2 << endl;
00109
00110 if(trku->GetBegPlane()<trkv->GetBegPlane()) bpln=trku->GetBegPlane(); else bpln=trkv->GetBegPlane();
00111 if(trku->GetEndPlane()>trkv->GetEndPlane()) epln=trku->GetEndPlane(); else epln=trkv->GetEndPlane();
00112
00113
00114
00115
00116
00117 MSG("AlgTrackAtNu", Msg::kDebug) << " *** calculate strip co-ordinates *** " << endl;
00118
00119 npln = 1+epln-bpln;
00120 Double_t* U = new Double_t[npln];
00121 Double_t* V = new Double_t[npln];
00122 Double_t* Z = new Double_t[npln];
00123 Double_t* Q = new Double_t[npln];
00124 Double_t* dS = new Double_t[npln];
00125 Double_t* dSsteel = new Double_t[npln];
00126 Double_t* Qm = new Double_t[npln];
00127 Double_t* Qp = new Double_t[npln];
00128 Double_t* CTm = new Double_t[npln];
00129 Double_t* CTp = new Double_t[npln];
00130 Double_t* Wm = new Double_t[npln];
00131 Double_t* Wp = new Double_t[npln];
00132 Int_t* plnvuw = new Int_t[npln];
00133 Int_t* plnnum = new Int_t[npln];
00134 for(i=0;i<npln;i++){
00135 plnvuw[i]=-1; plnnum[i]=0;
00136 U[i]=0.0; V[i]=0.0; Z[i]=0.0; Q[i]=0.0;
00137 Qm[i]=0.0; Qp[i]=0.0; CTm[i]=0.0; CTp[i]=0.0;
00138 Wm[i]=0.0; Wp[i]=0.0;
00139 }
00140
00141
00142 for(i=0;i<1+trku->GetHitLast();i++){
00143 HitAtNu* hit = (HitAtNu*)(trku->GetHitAt(i));
00144 CandStripHandle* strip = (CandStripHandle*)(hit->GetCandStripHandle());
00145 pln = hit->GetPlane()-bpln;
00146 fHitArr[pln].Add(hit);
00147 if(plnvuw[pln]<0){
00148 plnvuw[pln]=hit->GetPlaneView();
00149 Z[pln]=hit->GetZPos();
00150 }
00151 U[pln]+=hit->GetCharge()*hit->GetTPos();
00152 Q[pln]+=hit->GetCharge();
00153 track.AddDaughterLink(*strip);
00154 track.SetInShower(strip,hit->GetShwFlag());
00155 plnnum[pln]++;
00156 }
00157
00158 for(i=0;i<1+trkv->GetHitLast();i++){
00159 HitAtNu* hit = (HitAtNu*)(trkv->GetHitAt(i));
00160 CandStripHandle* strip = (CandStripHandle*)(hit->GetCandStripHandle());
00161 pln = hit->GetPlane()-bpln;
00162 fHitArr[pln].Add(hit);
00163 if(plnvuw[pln]<0){
00164 plnvuw[pln]=hit->GetPlaneView();
00165 Z[pln]=hit->GetZPos();
00166 }
00167 V[pln]+=hit->GetCharge()*hit->GetTPos();
00168 Q[pln]+=hit->GetCharge();
00169 track.AddDaughterLink(*strip);
00170 track.SetInShower(strip,hit->GetShwFlag());
00171 plnnum[pln]++;
00172 }
00173
00174 Int_t nplanes=0,nstrips=0;
00175 for(i=0;i<npln;i++){
00176 if(plnvuw[i]>-1){
00177 if(Q[i]>0.0){
00178 U[i]=U[i]/Q[i];
00179 V[i]=V[i]/Q[i];
00180 nstrips+=plnnum[i];
00181 nplanes++;
00182 }
00183 }
00184 }
00185
00186
00187 Int_t flag,vuw;
00188 Int_t km1,kp1,km2,kp2;
00189 for(pln=0;pln<npln;pln++){
00190 if(plnvuw[pln]>-1){
00191 flag=0;
00192 vuw=plnvuw[pln];
00193 k=0; km1=-1; km2=-1;
00194 while(pln-k>0 && km2<0){
00195 k++;
00196 if(plnvuw[pln-k]>-1 && plnvuw[pln-k]!=vuw){
00197 if(km1<0) km1=k; else km2=k;
00198 }
00199 }
00200 k=0; kp1=-1; kp2=-1;
00201 while(pln+k<npln-1 && kp2<0){
00202 k++;
00203 if(plnvuw[pln+k]>-1 && plnvuw[pln+k]!=vuw){
00204 if(kp1<0) kp1=k; else kp2=k;
00205 }
00206 }
00207 if(km1>0 && kp1>0){
00208 if(vuw==0) V[pln]=(kp1*V[pln-km1]+km1*V[pln+kp1])/(km1+kp1);
00209 if(vuw==1) U[pln]=(kp1*U[pln-km1]+km1*U[pln+kp1])/(km1+kp1);
00210 flag=1;
00211 }
00212 else{
00213 if(km1>0){
00214 if(vuw==0) if(km2>0) V[pln]=V[pln-km1]+(Z[pln]-Z[pln-km1])*(V[pln-km1]-V[pln-km2])/(Z[pln-km1]-Z[pln-km2]); else V[pln]=V[pln-km1];
00215 if(vuw==1) if(km2>0) U[pln]=U[pln-km1]+(Z[pln]-Z[pln-km1])*(U[pln-km1]-U[pln-km2])/(Z[pln-km1]-Z[pln-km2]); else U[pln]=U[pln-km1];
00216 flag=1;
00217 }
00218 if(kp1>0){
00219 if(vuw==0) if(kp2>0) V[pln]=V[pln+kp1]+(Z[pln]-Z[pln+kp1])*(V[pln+kp1]-V[pln+kp2])/(Z[pln+kp1]-Z[pln+kp2]); else V[pln]=V[pln+kp1];
00220 if(vuw==1) if(kp2>0) U[pln]=U[pln+kp1]+(Z[pln]-Z[pln+kp1])*(U[pln+kp1]-U[pln+kp2])/(Z[pln+kp1]-Z[pln+kp2]); else U[pln]=U[pln+kp1];
00221 flag=1;
00222 }
00223 }
00224 if(flag==0){
00225 plnvuw[pln]=-1;
00226 }
00227 }
00228 }
00229
00230
00231 Double_t utmp =0., vtmp = 0., ztmp = 0.;
00232 Double_t dq,du,dv,dz,ds;
00233 Double_t range_thru_detector,range_thru_steel;
00234 Double_t sumz_thru_steel,sumz_thru_detector;
00235 flag=0;
00236 range_thru_detector=0.0; range_thru_steel=0.0;
00237 sumz_thru_steel=0.0; sumz_thru_detector=0.0;
00238 for(i=0;i<npln;i++){
00239 if(plnvuw[i]>-1){
00240 if(flag){
00241 du=U[i]-utmp; dv=V[i]-vtmp; dz=Z[i]-ztmp;
00242 ds=sqrt(du*du+dv*dv+dz*dz);
00243 if(dz>1.0) dq=2.0; else dq=1.0;
00244 range_thru_steel+=dq*ds*(0.0254/dz);
00245 sumz_thru_steel+=dq*0.0254;
00246 range_thru_detector+=ds;
00247 sumz_thru_detector+=dz;
00248 }
00249 dS[i]=range_thru_detector;
00250 dSsteel[i]=range_thru_steel;
00251 utmp=U[i]; vtmp=V[i]; ztmp=Z[i];
00252 flag=1;
00253 }
00254 }
00255
00256
00257 Double_t tpos = 0., fibre = 0.;
00258 Double_t ctime = 0., digchg = 0.;
00259 Double_t indx=fFibreIndex;
00260 for(pln=0;pln<npln;pln++){
00261 for(i=0;i<1+fHitArr[pln].GetLast();i++){
00262 HitAtNu* hit = (HitAtNu*)(fHitArr[pln].At(i));
00263
00264 CandStripHandle* strip = (CandStripHandle*)(hit->GetCandStripHandle());
00265 if( strip->GetPlaneView()==PlaneView::kU
00266 || strip->GetPlaneView()==PlaneView::kX ) tpos = -V[pln];
00267 if( strip->GetPlaneView()==PlaneView::kV
00268 || strip->GetPlaneView()==PlaneView::kY ) tpos = U[pln];
00269 PlexStripEndId stripid = strip->GetStripEndId();
00270 UgliStripHandle striphandle = ugh.GetStripHandle(stripid);
00271
00272 if(strip->GetNDigit(StripEnd::kPositive)>0){
00273 fibre = striphandle.ClearFiber(StripEnd::kPositive)+striphandle.WlsPigtail(StripEnd::kPositive)+striphandle.GetHalfLength()-tpos;
00274 ctime = 3.0e8*strip->GetTime(StripEnd::kPositive)-indx*fibre;
00275 digchg=strip->GetCharge(StripEnd::kPositive);
00276 Qp[pln]+=digchg; CTp[pln]+=digchg*ctime;
00277 }
00278
00279 if(strip->GetNDigit(StripEnd::kNegative)>0){
00280 fibre = striphandle.ClearFiber(StripEnd::kNegative)+striphandle.WlsPigtail(StripEnd::kNegative)+striphandle.GetHalfLength()+tpos;
00281 ctime = 3.0e8*strip->GetTime(StripEnd::kNegative)-indx*fibre;
00282 digchg=strip->GetCharge(StripEnd::kNegative);
00283 Qm[pln]+=digchg; CTm[pln]+=digchg*ctime;
00284 }
00285
00286 }
00287 }
00288
00289
00290 Int_t plnm,plnp;
00291 plnm=-1; plnp=-1;
00292 for(i=0;i<npln;i++){
00293 pln=bpln+i;
00294 if(plnvuw[i]>-1){
00295 track.SetU(pln,U[i]);
00296 track.SetV(pln,V[i]);
00297 if(Qm[i]>0.0){
00298 CTm[i]=CTm[i]/Qm[i];
00299 track.SetT(pln,StripEnd::kNegative,CTm[i]/3.0e8);
00300 }
00301 if(Qp[i]>0.0){
00302 CTp[i]=CTp[i]/Qp[i];
00303 track.SetT(pln,StripEnd::kPositive,CTp[i]/3.0e8);
00304 }
00305 if(plnm<0||i<plnm) plnm=i; if(plnp<0||i>plnp) plnp=i;
00306 }
00307 }
00308
00309 MSG("AlgTrackAtNu", Msg::kDebug) << " *** track hits *** " << endl;
00310 for(pln=0;pln<npln;pln++){
00311 if( plnvuw[pln]>-1 ){
00312 MSG("AlgTrackAtNu",Msg::kDebug)
00313 << " plane=" << pln << " view=" << plnvuw[pln]
00314 << " z=" << Z[pln] << " u=" << U[pln] << " v=" << V[pln]
00315 << " s=" << dS[pln] << endl;
00316 }
00317 }
00318
00319
00320
00321
00322
00323 MSG("AlgTrackAtNu", Msg::kDebug) << " *** vertex + direction *** " << endl;
00324
00325
00326 Double_t mu,cu,mv,cv,err;
00327 Double_t u,v,z,w;
00328 Double_t Uw,Uwz,Uwzz,Uwu,Uwuu,Uwzu;
00329 Double_t Vw,Vwz,Vwzz,Vwv,Vwvv,Vwzv;
00330
00331 Int_t Upts,Vpts,MAXpts=4;
00332 Double_t precou,precov,precos,precoz;
00333 Int_t ndf;
00334
00335 Double_t bdiru,bdirv,bdirx,bdiry,bdirz;
00336 Double_t bvtxu,bvtxv,bvtxx,bvtxy,bvtxz;
00337 Double_t bdirerr=-1.0;
00338 Int_t bndf=-1;
00339
00340 Double_t ediru,edirv,edirx,ediry,edirz;
00341 Double_t evtxu,evtxv,evtxx,evtxy,evtxz;
00342 Double_t edirerr=-1.0;
00343 Int_t endf=-1;
00344
00345 Double_t lindiru,lindirv,lindirx,lindiry,lindirz;
00346 Double_t lindirerr=-1.0;
00347 Int_t linndf=-1;
00348
00349
00350 Uw=0.0; Uwz=0.0; Uwzz=0.0; Uwu=0.0; Uwuu=0.0; Uwzu=0.0;
00351 Vw=0.0; Vwz=0.0; Vwzz=0.0; Vwv=0.0; Vwvv=0.0; Vwzv=0.0;
00352 Upts=0; Vpts=0;
00353
00354 for(i=0;i<npln;i++){
00355 if(i>=plnm && plnvuw[i]>-1){
00356 u=U[i]; v=V[i]; z=Z[i]; w=1.0;
00357
00358 if(plnvuw[i]==0 && Upts<MAXpts){
00359 Uw+=w;
00360 Uwu+=w*u; Uwz+=w*z; Uwuu+=w*u*u;
00361 Uwzu+=w*u*z; Uwzz+=w*z*z;
00362 Upts++;
00363 }
00364
00365 if(plnvuw[i]==1 && Vpts<MAXpts){
00366 Vw+=w;
00367 Vwv+=w*v; Vwz+=w*z; Vwvv+=w*v*v;
00368 Vwzv+=w*v*z; Vwzz+=w*z*z;
00369 Vpts++;
00370 }
00371 }
00372 }
00373
00374 mu=0.0; mv=0.0; err=-1.0; ndf=-1;
00375 if(Upts>1 && Vpts>1){
00376 mu=(Uw*Uwzu-Uwz*Uwu)/(Uw*Uwzz-Uwz*Uwz); cu=(Uwu*Uwzz-Uwz*Uwzu)/(Uw*Uwzz-Uwz*Uwz);
00377 mv=(Vw*Vwzv-Vwz*Vwv)/(Vw*Vwzz-Vwz*Vwz); cv=(Vwv*Vwzz-Vwz*Vwzv)/(Vw*Vwzz-Vwz*Vwz);
00378 err=(Uwuu-2.0*mu*Uwzu-2.0*cu*Uwu+mu*mu*Uwzz+2.0*mu*cu*Uwz+Uw*cu*cu)+(Vwvv-2.0*mv*Vwzv-2.0*cv*Vwv+mv*mv*Vwzz+2.0*mv*cv*Vwz+Vw*cv*cv);
00379 ndf=Upts+Vpts-4;
00380 }
00381
00382 precou=mu; precov=mv;
00383 precos=sqrt(precou*precou+precov*precov+1.0);
00384 precou=precou/precos; precov=precov/precos; precoz=1.0/precos;
00385 bdiru=precou; bdirv=precov; bdirz=precoz;
00386 bdirx=0.7071*(bdiru-bdirv); bdiry=0.7071*(bdiru+bdirv);
00387 bdirerr=err;
00388 bndf=ndf;
00389
00390 bvtxu=U[plnm]-0.02*(bdiru/bdirz); bvtxv=V[plnm]-0.02*(bdirv/bdirz); bvtxz=Z[plnm]-0.02;
00391 bvtxx=0.7071*(bvtxu-bvtxv); bvtxy=0.7071*(bvtxu+bvtxv);
00392
00393
00394 Uw=0.0; Uwz=0.0; Uwzz=0.0; Uwu=0.0; Uwuu=0.0; Uwzu=0.0;
00395 Vw=0.0; Vwz=0.0; Vwzz=0.0; Vwv=0.0; Vwvv=0.0; Vwzv=0.0;
00396 Upts=0; Vpts=0;
00397
00398 for(i=npln-1;i>-1;i--){
00399 if(i<=plnp && plnvuw[i]>-1){
00400 u=U[i]; v=V[i]; z=Z[i]; w=1.0;
00401
00402 if(plnvuw[i]==0 && Upts<MAXpts){
00403 Uw+=w;
00404 Uwu+=w*u; Uwz+=w*z; Uwuu+=w*u*u;
00405 Uwzu+=w*u*z; Uwzz+=w*z*z;
00406 Upts++;
00407 }
00408
00409 if(plnvuw[i]==1 && Vpts<MAXpts){
00410 Vw+=w;
00411 Vwv+=w*v; Vwz+=w*z; Vwvv+=w*v*v;
00412 Vwzv+=w*v*z; Vwzz+=w*z*z;
00413 Vpts++;
00414 }
00415 }
00416 }
00417
00418 mu=0.0; mv=0.0; err=-1.0; ndf=-1;
00419 if(Upts>1 && Vpts>1){
00420 mu=(Uw*Uwzu-Uwz*Uwu)/(Uw*Uwzz-Uwz*Uwz); cu=(Uwu*Uwzz-Uwz*Uwzu)/(Uw*Uwzz-Uwz*Uwz);
00421 mv=(Vw*Vwzv-Vwz*Vwv)/(Vw*Vwzz-Vwz*Vwz); cv=(Vwv*Vwzz-Vwz*Vwzv)/(Vw*Vwzz-Vwz*Vwz);
00422 err=(Uwuu-2.0*mu*Uwzu-2.0*cu*Uwu+mu*mu*Uwzz+2.0*mu*cu*Uwz+Uw*cu*cu)+(Vwvv-2.0*mv*Vwzv-2.0*cv*Vwv+mv*mv*Vwzz+2.0*mv*cv*Vwz+Vw*cv*cv);
00423 ndf=Upts+Vpts-4;
00424 }
00425
00426 precou=mu; precov=mv;
00427 precos=sqrt(precou*precou+precov*precov+1.0);
00428 precou=precou/precos; precov=precov/precos; precoz=1.0/precos;
00429 ediru=precou; edirv=precov; edirz=precoz;
00430 edirx=0.7071*(ediru-edirv); ediry=0.7071*(ediru+edirv);
00431 edirerr=err;
00432 endf=ndf;
00433
00434 evtxu=U[plnp]+0.04*(ediru/edirz); evtxv=V[plnp]+0.04*(edirv/edirz); evtxz=Z[plnp]+0.04;
00435 evtxx=0.7071*(evtxu-evtxv); evtxy=0.7071*(evtxu+evtxv);
00436
00437 MSG("AlgTrackAtNu", Msg::kDebug) << " beg vertex : UVXYZ = (" << bvtxu << "," << bvtxv << "," << bvtxx << "," << bvtxy << "," << bvtxz << ") " << endl;
00438
00439 MSG("AlgTrackAtNu", Msg::kDebug) << " end vertex : UVXYZ = (" << evtxu << "," << evtxv << "," << evtxx << "," << evtxy << "," << evtxz << ") " << endl;
00440
00441
00442 Uw=0.0; Uwz=0.0; Uwzz=0.0; Uwu=0.0; Uwuu=0.0; Uwzu=0.0;
00443 Vw=0.0; Vwz=0.0; Vwzz=0.0; Vwv=0.0; Vwvv=0.0; Vwzv=0.0;
00444 Upts=0; Vpts=0;
00445
00446 for(i=0;i<npln;i++){
00447 if(i>=plnm && i<=plnp && plnvuw[i]>-1){
00448 u=U[i]; v=V[i]; z=Z[i]; w=Q[i];
00449
00450 if(plnvuw[i]==0){
00451 Uw+=w;
00452 Uwu+=w*u; Uwz+=w*z; Uwuu+=w*u*u;
00453 Uwzu+=w*u*z; Uwzz+=w*z*z;
00454 Upts++;
00455 }
00456
00457 if(plnvuw[i]==1){
00458 Vw+=w;
00459 Vwv+=w*v; Vwz+=w*z; Vwvv+=w*v*v;
00460 Vwzv+=w*v*z; Vwzz+=w*z*z;
00461 Vpts++;
00462 }
00463 }
00464 }
00465
00466 mu=0.0; mv=0.0; err=-1.0; ndf=-1;
00467 if(Upts>1 && Vpts>1){
00468 mu=(Uw*Uwzu-Uwz*Uwu)/(Uw*Uwzz-Uwz*Uwz); cu=(Uwu*Uwzz-Uwz*Uwzu)/(Uw*Uwzz-Uwz*Uwz);
00469 mv=(Vw*Vwzv-Vwz*Vwv)/(Vw*Vwzz-Vwz*Vwz); cv=(Vwv*Vwzz-Vwz*Vwzv)/(Vw*Vwzz-Vwz*Vwz);
00470 err=(Uwuu-2.0*mu*Uwzu-2.0*cu*Uwu+mu*mu*Uwzz+2.0*mu*cu*Uwz+Uw*cu*cu)+(Vwvv-2.0*mv*Vwzv-2.0*cv*Vwv+mv*mv*Vwzz+2.0*mv*cv*Vwz+Vw*cv*cv);
00471 ndf=Upts+Vpts-4;
00472 }
00473
00474 precou=mu; precov=mv;
00475 precos=sqrt(precou*precou+precov*precov+1.0);
00476 precou=precou/precos; precov=precov/precos; precoz=1.0/precos;
00477 lindiru=precou; lindirv=precov; lindirz=precoz;
00478 lindirx=0.7071*(lindiru-lindirv); lindiry=0.7071*(lindiru+lindirv);
00479 lindirerr=err;
00480 linndf=ndf;
00481
00482
00483
00484
00485
00486 MSG("AlgTrackAtNu", Msg::kDebug) << " *** containment stuff *** " << endl;
00487
00488 Double_t begr=0.0,endr=0.0;
00489 Double_t begtrace=0.0,endtrace=0.0;
00490 Double_t begtraceZ=0.0,endtraceZ=0.0;
00491
00492 if( atmosflag ){
00493
00494
00495 Double_t xm,xp,ym,yp,um,up,vm,vp;
00496 Double_t rmin;
00497
00498 rmin=4.0;
00499 up=4.0-bvtxu; if(up<rmin) rmin=up;
00500 um=4.0+bvtxu; if(um<rmin) rmin=um;
00501 vp=4.0-bvtxv; if(vp<rmin) rmin=vp;
00502 vm=4.0+bvtxv; if(vm<rmin) rmin=vm;
00503 xp=4.0-bvtxx; if(xp<rmin) rmin=xp;
00504 xm=4.0+bvtxx; if(xm<rmin) rmin=xm;
00505 yp=4.0-bvtxy; if(yp<rmin) rmin=yp;
00506 ym=4.0+bvtxy; if(ym<rmin) rmin=ym;
00507 begr=rmin;
00508
00509 rmin=4.0;
00510 up=4.0-evtxu; if(up<rmin) rmin=up;
00511 um=4.0+evtxu; if(um<rmin) rmin=um;
00512 vp=4.0-evtxv; if(vp<rmin) rmin=vp;
00513 vm=4.0+evtxv; if(vm<rmin) rmin=vm;
00514 xp=4.0-evtxx; if(xp<rmin) rmin=xp;
00515 xm=4.0+evtxx; if(xm<rmin) rmin=xm;
00516 yp=4.0-evtxy; if(yp<rmin) rmin=yp;
00517 ym=4.0+evtxy; if(ym<rmin) rmin=ym;
00518 endr=rmin;
00519
00520 Int_t trkflag;
00521 Double_t posX,posY,posU,posV,posZ;
00522 Double_t dirX,dirY,dirU,dirV,dirZ;
00523 Double_t dX,dY,dU,dV,dR;
00524 Double_t trace,traceZ;
00525 Double_t tmpU,tmpV,tmpZ;
00526 Double_t tmpdirU,tmpdirV,tmpdirZ;
00527 Double_t gradu,gradv,grads;
00528 Double_t tmpgrad;
00529 Int_t jm,jp;
00530
00531
00532 tmpU = bvtxu; tmpV = bvtxv; tmpZ = bvtxz;
00533 tmpdirU = bdiru; tmpdirV = bdirv; tmpdirZ = bdirz;
00534 gradu = (tmpdirU/tmpdirZ); gradv = (tmpdirV/tmpdirZ);
00535
00536 for(j=0;j<10;j++){
00537 jm=plnm+j; jp=-1;
00538 if( jm<npln && plnvuw[jm]>-1 && Z[jm]-bvtxz<1.0 ){
00539 if( jp<0 && jm+2<10 && jm+2<npln && plnvuw[jm+2]>-1 && Z[jm+2]-bvtxz<1.0 ) jp=jm+2;
00540 if( jp<0 && jm+4<10 && jm+4<npln && plnvuw[jm+4]>-1 && Z[jm+4]-bvtxz<1.0 ) jp=jm+4;
00541 if( jp>=0 && plnvuw[jm]==plnvuw[jp] ){
00542 if( plnvuw[jm]==0 ){
00543 tmpgrad=(U[jp]-U[jm])/(Z[jp]-Z[jm]);
00544 if( (gradu<0&&tmpgrad<gradu)||(gradu>0.0&&tmpgrad>gradu) ){
00545 gradu=tmpgrad; tmpU=U[jm]+gradu*(tmpZ-Z[jm]);
00546 }
00547 }
00548 if( plnvuw[jm]==1 ){
00549 tmpgrad=(V[jp]-V[jm])/(Z[jp]-Z[jm]);
00550 if( (gradv<0&&tmpgrad<gradv)||(gradv>0.0&&tmpgrad>gradv) ){
00551 gradv=tmpgrad; tmpV=V[jm]+gradv*(tmpZ-Z[jm]);
00552 }
00553 }
00554 }
00555 }
00556 }
00557
00558 posU=tmpU; posV=tmpV; posZ=tmpZ; posX=0.7071*(posU-posV); posY=0.7071*(posU+posV);
00559 grads=sqrt(gradu*gradu+gradv*gradv+1.0);
00560 dirU=-gradu/grads; dirV=-gradv/grads; dirZ=-1.0/grads;
00561 dirX=0.7071*(dirU-dirV); dirY=0.7071*(dirU+dirV);
00562 trkflag=0; dR=0.0;
00563
00564 if(trkflag==0){
00565 if(dirX>0.0){
00566 dX = 4.0-posX; dY = dX*(dirY/dirX);
00567 if( posY+dY<4.0*(sqrt(2.0)-1.0) && posY+dY>-4.0*(sqrt(2.0)-1.0) ){
00568 trkflag=1; if(dX>0.0) dR=sqrt(dX*dX+dY*dY); if(dX<0.0) dR=-sqrt(dX*dX+dY*dY);
00569 }
00570 }
00571 if(dirX<0.0){
00572 dX = 4.0+posX; dY = -dX*(dirY/dirX);
00573 if( posY+dY<4.0*(sqrt(2.0)-1.0) && posY+dY>-4.0*(sqrt(2.0)-1.0) ){
00574 trkflag=1; if(dX>0.0) dR=sqrt(dX*dX+dY*dY); if(dX<0.0) dR=-sqrt(dX*dX+dY*dY);
00575 }
00576 }
00577 }
00578
00579 if(trkflag==0){
00580 if(dirY>0.0){
00581 dY = 4.0-posY; dX = dY*(dirX/dirY);
00582 if( posX+dX<4.0*(sqrt(2.0)-1.0) && posX+dX>-4.0*(sqrt(2.0)-1.0) ){
00583 trkflag=1; if(dY>0.0) dR=sqrt(dX*dX+dY*dY); if(dY<0.0) dR=-sqrt(dX*dX+dY*dY);
00584 }
00585 }
00586 if(dirY<0.0){
00587 dY = 4.0+posY; dX = -dY*(dirX/dirY);
00588 if( posX+dX<4.0*(sqrt(2.0)-1.0) && posX+dX>-4.0*(sqrt(2.0)-1.0) ){
00589 trkflag=1; if(dY>0.0) dR=sqrt(dX*dX+dY*dY); if(dY<0.0) dR=-sqrt(dX*dX+dY*dY);
00590 }
00591 }
00592 }
00593
00594 if(trkflag==0){
00595 if(dirU>0.0){
00596 dU = 4.0-posU; dV = dU*(dirV/dirU);
00597 if( posV+dV<4.0*(sqrt(2.0)-1.0) && posV+dV>-4.0*(sqrt(2.0)-1.0) ){
00598 trkflag=1; if(dU>0.0) dR=sqrt(dV*dV+dU*dU); if(dU<0.0) dR=-sqrt(dV*dV+dU*dU);
00599 }
00600 }
00601 if(dirU<0.0){
00602 dU = 4.0+posU; dV = -dU*(dirV/dirU);
00603 if( posV+dV<4.0*(sqrt(2.0)-1.0) && posV+dV>-4.0*(sqrt(2.0)-1.0) ){
00604 trkflag=1; if(dU>0.0) dR=sqrt(dV*dV+dU*dU); if(dU<0.0) dR=-sqrt(dV*dV+dU*dU);
00605 }
00606 }
00607 }
00608
00609 if(trkflag==0){
00610 if(dirV>0.0){
00611 dV = 4.0-posV; dU = dV*(dirU/dirV);
00612 if( posU+dU<4.0*(sqrt(2.0)-1.0) && posU+dU>-4.0*(sqrt(2.0)-1.0) ){
00613 trkflag=1; if(dV>0.0) dR=sqrt(dV*dV+dU*dU); if(dV<0.0) dR=-sqrt(dV*dV+dU*dU);
00614 }
00615 }
00616 if(dirV<0.0){
00617 dV = 4.0+posV; dU = -dV*(dirU/dirV);
00618 if( posU+dU<4.0*(sqrt(2.0)-1.0) && posU+dU>-4.0*(sqrt(2.0)-1.0) ){
00619 trkflag=1; if(dV>0.0) dR=sqrt(dV*dV+dU*dU); if(dV<0.0) dR=-sqrt(dV*dV+dU*dU);
00620 }
00621 }
00622 }
00623
00624 traceZ=0.0; trace=0.0;
00625 if(sqrt(dirX*dirX+dirY*dirY)>0.0){
00626 traceZ = -(dR*dirZ)/sqrt(dirX*dirX+dirY*dirY);
00627 }
00628 if(traceZ>0.0){
00629 if(begZsm1>-100.0 && posZ-begZsm1>-0.06 && traceZ>posZ-begZsm1) traceZ = posZ-begZsm1;
00630 if(begZsm2>-100.0 && posZ-begZsm2>-0.06 && traceZ>posZ-begZsm2) traceZ = posZ-begZsm2;
00631 }
00632 if(traceZ<0.0){
00633 if(begZsm1>-100.0 && posZ-begZsm1>-0.06 && traceZ<-(posZ-begZsm1)) traceZ = -(posZ-begZsm1);
00634 if(begZsm2>-100.0 && posZ-begZsm2>-0.06 && traceZ<-(posZ-begZsm2)) traceZ = -(posZ-begZsm2);
00635 }
00636 if(dirZ<0.0){
00637 trace=-traceZ/dirZ;
00638 }
00639 begtrace=trace;
00640 begtraceZ=traceZ;
00641
00642
00643 tmpU = evtxu; tmpV = evtxv; tmpZ = evtxz;
00644 tmpdirU = ediru; tmpdirV = edirv; tmpdirZ = edirz;
00645 gradu = (tmpdirU/tmpdirZ); gradv = (tmpdirV/tmpdirZ);
00646
00647 for(j=0;j<10;j++){
00648 jm=-1; jp=plnp-1-j;
00649 if( jp>-1 && plnvuw[jp]>-1 && Z[jp]-evtxz>-1.0 ){
00650 if( jm<0 && jp-2>npln-1-10 && jp-2>-1 && plnvuw[jp-2]>-1 && Z[jp-2]-evtxz>-1.0 ) jm=jp-2;
00651 if( jm<0 && jp-4>npln-1-10 && jp-4>-1 && plnvuw[jp-4]>-1 && Z[jp-4]-evtxz>-1.0 ) jm=jp-4;
00652 if( jm>=0 && plnvuw[jm]==plnvuw[jp] ){
00653 if( plnvuw[jp]==0 ){
00654 tmpgrad=(U[jp]-U[jm])/(Z[jp]-Z[jm]);
00655 if( (gradu<0&&tmpgrad<gradu)||(gradu>0.0&&tmpgrad>gradu) ){
00656 gradu=tmpgrad; tmpU=U[jp]+gradu*(tmpZ-Z[jp]);
00657 }
00658 }
00659 if( plnvuw[jp]==1 ){
00660 tmpgrad=(V[jp]-V[jm])/(Z[jp]-Z[jm]);
00661 if( (gradv<0&&tmpgrad<gradv)||(gradv>0.0&&tmpgrad>gradv) ){
00662 gradv=tmpgrad; tmpV=V[jp]+gradv*(tmpZ-Z[jp]);
00663 }
00664 }
00665 }
00666 }
00667 }
00668
00669 posU=tmpU; posV=tmpV; posZ=tmpZ; posX=0.7071*(posU-posV); posY=0.7071*(posU+posV);
00670 grads=sqrt(gradu*gradu+gradv*gradv+1.0);
00671 dirU=gradu/grads; dirV=gradv/grads; dirZ=1.0/grads;
00672 dirX=0.7071*(dirU-dirV); dirY=0.7071*(dirU+dirV);
00673 trkflag=0; dR=0.0;
00674
00675 if(trkflag==0){
00676 if(dirX>0.0){
00677 dX = 4.0-posX; dY = dX*(dirY/dirX);
00678 if( posY+dY<4.0*(sqrt(2.0)-1.0) && posY+dY>-4.0*(sqrt(2.0)-1.0) ){
00679 trkflag=1; if(dX>0.0) dR=sqrt(dX*dX+dY*dY); if(dX<0.0) dR=-sqrt(dX*dX+dY*dY);
00680 }
00681 }
00682 if(dirX<0.0){
00683 dX = 4.0+posX; dY = -dX*(dirY/dirX);
00684 if( posY+dY<4.0*(sqrt(2.0)-1.0) && posY+dY>-4.0*(sqrt(2.0)-1.0) ){
00685 trkflag=1; if(dX>0.0) dR=sqrt(dX*dX+dY*dY); if(dX<0.0) dR=-sqrt(dX*dX+dY*dY);
00686 }
00687 }
00688 }
00689
00690 if(trkflag==0){
00691 if(dirY>0.0){
00692 dY = 4.0-posY; dX = dY*(dirX/dirY);
00693 if( posX+dX<4.0*(sqrt(2.0)-1.0) && posX+dX>-4.0*(sqrt(2.0)-1.0) ){
00694 trkflag=1; if(dY>0.0) dR=sqrt(dX*dX+dY*dY); if(dY<0.0) dR=-sqrt(dX*dX+dY*dY);
00695 }
00696 }
00697 if(dirY<0.0){
00698 dY = 4.0+posY; dX = -dY*(dirX/dirY);
00699 if( posX+dX<4.0*(sqrt(2.0)-1.0) && posX+dX>-4.0*(sqrt(2.0)-1.0) ){
00700 trkflag=1; if(dY>0.0) dR=sqrt(dX*dX+dY*dY); if(dY<0.0) dR=-sqrt(dX*dX+dY*dY);
00701 }
00702 }
00703 }
00704
00705 if(trkflag==0){
00706 if(dirU>0.0){
00707 dU = 4.0-posU; dV = dU*(dirV/dirU);
00708 if( posV+dV<4.0*(sqrt(2.0)-1.0) && posV+dV>-4.0*(sqrt(2.0)-1.0) ){
00709 trkflag=1; if(dU>0.0) dR=sqrt(dV*dV+dU*dU); if(dU<0.0) dR=-sqrt(dV*dV+dU*dU);
00710 }
00711 }
00712 if(dirU<0.0){
00713 dU = 4.0+posU; dV = -dU*(dirV/dirU);
00714 if( posV+dV<4.0*(sqrt(2.0)-1.0) && posV+dV>-4.0*(sqrt(2.0)-1.0) ){
00715 trkflag=1; if(dU>0.0) dR=sqrt(dV*dV+dU*dU); if(dU<0.0) dR=-sqrt(dV*dV+dU*dU);
00716 }
00717 }
00718 }
00719
00720 if(trkflag==0){
00721 if(dirV>0.0){
00722 dV = 4.0-posV; dU = dV*(dirU/dirV);
00723 if( posU+dU<4.0*(sqrt(2.0)-1.0) && posU+dU>-4.0*(sqrt(2.0)-1.0) ){
00724 trkflag=1; if(dV>0.0) dR=sqrt(dV*dV+dU*dU); if(dV<0.0) dR=-sqrt(dV*dV+dU*dU);
00725 }
00726 }
00727 if(dirV<0.0){
00728 dV = 4.0+posV; dU = -dV*(dirU/dirV);
00729 if( posU+dU<4.0*(sqrt(2.0)-1.0) && posU+dU>-4.0*(sqrt(2.0)-1.0) ){
00730 trkflag=1; if(dV>0.0) dR=sqrt(dV*dV+dU*dU); if(dV<0.0) dR=-sqrt(dV*dV+dU*dU);
00731 }
00732 }
00733 }
00734
00735 traceZ = 0.0; trace=0.0;
00736 if(sqrt(dirX*dirX+dirY*dirY)>0.0){
00737 traceZ = (dR*dirZ)/sqrt(dirX*dirX+dirY*dirY);
00738 }
00739 if(traceZ>0.0){
00740 if(endZsm1>-100.0 && endZsm1-posZ>-0.06 && traceZ>endZsm1-posZ) traceZ = endZsm1-posZ;
00741 if(endZsm2>-100.0 && endZsm2-posZ>-0.06 && traceZ>endZsm2-posZ) traceZ = endZsm2-posZ;
00742 }
00743 if(traceZ<0.0){
00744 if(endZsm1>-100.0 && endZsm1-posZ>-0.06 && traceZ<-(endZsm1-posZ)) traceZ = -(endZsm1-posZ);
00745 if(endZsm2>-100.0 && endZsm2-posZ>-0.06 && traceZ<-(endZsm2-posZ)) traceZ = -(endZsm2-posZ);
00746 }
00747 if(dirZ>0.0){
00748 trace = traceZ/dirZ;
00749 }
00750 endtrace=trace;
00751 endtraceZ=traceZ;
00752
00753 }
00754
00755
00756
00757
00758
00759 MSG("AlgTrackAtNu", Msg::kDebug) << " *** direction using timing *** " << endl;
00760
00761
00762 Int_t ctflag;
00763 Double_t ctmin,ctmax,ctave,ctrms;
00764 Double_t sq,sqct,sqctct,swgt;
00765 Double_t wm,wp;
00766 for(i=0;i<npln;i++){
00767 if( plnvuw[i]>-1 ){
00768 wm=0.0; wp=0.0;
00769 if(Qm[i]>0.0) wm=1.0; if(Qp[i]>0.0) wp=1.0;
00770 ctflag=0; ctmin=0.0; ctmax=0.0; ctave=0.0; ctrms=0.0;
00771 sq=0.0; sqct=0.0; sqctct=0.0; swgt=0.0;
00772 for(j=-4;j<5;j++){
00773 if( j!=0 && i+j>-1 && i+j<npln && plnvuw[i+j]>-1 ){
00774 if( Qm[i+j]>1.0
00775 && !( Qp[i+j]>0.0 && !(CTm[i+j]-CTp[i+j]>-4.0&&CTm[i+j]-CTp[i+j]<4.0) ) ){
00776 if(ctflag){
00777 if(CTm[i+j]<ctmin) ctmin=CTm[i+j]; if(CTm[i+j]>ctmax) ctmax=CTm[i+j];
00778 }
00779 else{
00780 ctmin=CTm[i+j]; ctmax=CTm[i+j]; ctflag=1;
00781 }
00782 sq+=Qm[i+j]; sqct+=Qm[i+j]*CTm[i+j]; sqctct+=Qm[i+j]*CTm[i+j]*CTm[i+j];
00783 swgt+=1.0;
00784 }
00785 if( Qp[i+j]>1.0
00786 && !( Qm[i+j]>0.0 && !(CTp[i+j]-CTm[i+j]>-4.0&&CTp[i+j]-CTm[i+j]<4.0) ) ){
00787 if(ctflag){
00788 if(CTp[i+j]<ctmin) ctmin=CTp[i+j]; if(CTp[i+j]>ctmax) ctmax=CTp[i+j];
00789 }
00790 else{
00791 ctmin=CTp[i+j]; ctmax=CTp[i+j]; ctflag=1;
00792 }
00793 sq+=Qp[i+j]; sqct+=Qp[i+j]*CTp[i+j]; sqctct+=Qp[i+j]*CTp[i+j]*CTp[i+j];
00794 swgt+=1.0;
00795 }
00796 }
00797 }
00798 if( swgt>1.0 ){
00799 ctave=sqct/sq; ctrms=0.0;
00800 if((sqctct/sq)-(sqct/sq)*(sqct/sq)>0.0){
00801 ctrms=sqrt((sqctct/sq)-(sqct/sq)*(sqct/sq));
00802 }
00803 if( (Qm[i]>1.0)
00804 && (CTm[i]-ctmin>-2.0 && CTm[i]-ctmax<2.0)
00805 && (CTm[i]-ctave<ctrms+3.0 && CTm[i]-ctave>-ctrms-3.0) ){
00806 wm=Qm[i];
00807 }
00808 if( (Qp[i]>1.0)
00809 && (CTp[i]-ctmin>-2.0 && CTp[i]-ctmax<2.0)
00810 && (CTp[i]-ctave<ctrms+3.0 && CTp[i]-ctave>-ctrms-3.0) ){
00811 wp=Qp[i];
00812 }
00813 }
00814 Wm[i]=wm; Wp[i]=wp;
00815 }
00816 }
00817
00818 for(i=0;i<npln;i++){
00819 if( plnvuw[i]>-1 ){
00820 MSG("AlgTrackAtNu", Msg::kDebug) << " S=" << dS[i]
00821 << " CTm=" << CTm[i] << " CTp=" << CTp[i]
00822 << " Wm=" << Wm[i] << " Wp=" << Wp[i] << endl;
00823 }
00824 }
00825
00826
00827 Double_t ctimeslope=0.0,ctimeoffset=0.0,ctimeaverage=0.0,ctimescatter=0.0;
00828 Double_t Sw,Sws,Swss,Swt,Swtt,Swst;
00829 Double_t s,t;
00830 Int_t npts=0;
00831 Sw=0.0; Sws=0.0; Swss=0.0; Swt=0.0; Swtt=0.0; Swst=0.0;
00832 for(i=0;i<npln;i++){
00833 if( plnvuw[i]>-1 ){
00834 if( Wm[i]>0.0 ){
00835 s=dS[i]; t=CTm[i]; w=Wm[i];
00836 Sws+=w*s; Swss+=w*s*s; Swt+=w*t; Swtt+=w*t*t; Swst+=w*s*t;
00837 Sw+=w; npts++;
00838 }
00839 if( Wp[i]>0.0 ){
00840 s=dS[i]; t=CTp[i]; w=Wp[i];
00841 Sws+=w*s; Swss+=w*s*s; Swt+=w*t; Swtt+=w*t*t; Swst+=w*s*t;
00842 Sw+=w; npts++;
00843 }
00844 }
00845 }
00846 if(npts>2){
00847 ctimeslope=(Sw*Swst-Sws*Swt)/(Sw*Swss-Sws*Sws);
00848 ctimeoffset=(Swt*Swss-Sws*Swst)/(Sw*Swss-Sws*Sws);
00849 ctimeaverage=(Swt/Sw); ctimescatter=0.0;
00850 if( (Swtt/Sw)-(Swt/Sw)*(Swt/Sw)>0.0 && (Swss/Sw)-(Sws/Sw)*(Sws/Sw)>0.0 ){
00851 ctimescatter=((Swst/Sw)-(Sws/Sw)*(Swt/Sw))/((sqrt((Swss/Sw)-(Sws/Sw)*(Sws/Sw)))*(sqrt((Swtt/Sw)-(Swt/Sw)*(Swt/Sw))));
00852 }
00853 }
00854
00855
00856 Double_t C,Csigma,Ctrunc;
00857 Double_t cp,cm;
00858 Double_t chisqp=0.0,chisqm=0.0;
00859 Int_t chisqndfp=0,chisqndfm=0;
00860 Double_t dir,time_score,mytimeoffset;
00861 Int_t itr;
00862
00863 MSG("AlgTrackAtNu", Msg::kDebug) << " POSITIVE FIT: " << endl;
00864 chisqp=-99999.9; cp=-99999.9;
00865 C=ctimeaverage; Csigma=-99999.9; Ctrunc=-99999.9;
00866 for(itr=0;itr<2;itr++){
00867 Swtt=0.0; Swt=0.0; Sw=0.0; npts=0;
00868 for(pln=0;pln<npln;pln++){
00869 if( plnvuw[pln]>-1 ){
00870 if( Wm[pln]>0.0 ){
00871 w=Wm[pln]; t=CTm[pln]-dS[pln]+C;
00872 if(Ctrunc<0.0||(t>-Ctrunc&&t<Ctrunc)){
00873 Swtt+=w*t*t; Swt+=w*t; Sw+=w; npts++;
00874 }
00875 }
00876 if( Wp[pln]>0.0 ){
00877 w=Wp[pln]; t=CTp[pln]-dS[pln]+C;
00878 if(Ctrunc<0.0||(t>-Ctrunc&&t<Ctrunc)){
00879 Swtt+=w*t*t; Swt+=w*t; Sw+=w; npts++;
00880 }
00881 }
00882 }
00883 }
00884 if(npts>1){
00885 C=C-Swt/Sw; Csigma=0.0;
00886 if((Swtt/Sw)-(Swt/Sw)*(Swt/Sw)>0.0){
00887 Csigma=sqrt((Swtt/Sw)-(Swt/Sw)*(Swt/Sw));
00888 }
00889 chisqp=Csigma; chisqndfp=npts-1; cp=-C;
00890 Ctrunc=Csigma+3.0;
00891 }
00892 }
00893 MSG("AlgTrackAtNu", Msg::kDebug) << " *** C=" << C << " Csigma=" << Csigma << " *** " << endl;
00894
00895
00896 MSG("AlgTrackAtNu", Msg::kDebug) << " NEGATIVE FIT: " << endl;
00897 chisqm=-99999.9; cm=-99999.9;
00898 C=ctimeaverage; Csigma=-99999.9; Ctrunc=-99999.9;
00899 for(itr=0;itr<2;itr++){
00900 Swtt=0.0; Swt=0.0; Sw=0.0; npts=0;
00901 for(pln=0;pln<npln;pln++){
00902 if( plnvuw[pln]>-1 ){
00903 if( Wm[pln]>0.0 ){
00904 w=Wm[pln]; t=CTm[pln]+dS[pln]+C;
00905 if(Ctrunc<0.0||(t>-Ctrunc&&t<Ctrunc)){
00906 Swtt+=w*t*t; Swt+=w*t; Sw+=w; npts++;
00907 }
00908 }
00909 if( Wp[pln]>0.0 ){
00910 w=Wp[pln]; t=CTp[pln]+dS[pln]+C;
00911 if(Ctrunc<0.0||(t>-Ctrunc&&t<Ctrunc)){
00912 Swtt+=w*t*t; Swt+=w*t; Sw+=w; npts++;
00913 }
00914 }
00915 }
00916 }
00917 if(npts>1){
00918 C=C-Swt/Sw; Csigma=0.0;
00919 if((Swtt/Sw)-(Swt/Sw)*(Swt/Sw)>0.0){
00920 Csigma=sqrt((Swtt/Sw)-(Swt/Sw)*(Swt/Sw));
00921 }
00922 chisqm=Csigma; chisqndfm=npts-1; cm=-C;
00923 Ctrunc=Csigma+3.0;
00924 }
00925 }
00926 MSG("AlgTrackAtNu", Msg::kDebug) << " *** C=" << C << " Csigma=" << Csigma << " *** " << endl;
00927
00928 dir=0.0;
00929 time_score=0.0; mytimeoffset=0.0;
00930 if(chisqm>0.0 && chisqp>0.0){
00931 if(chisqp>chisqm){
00932 dir=-1.0; time_score=1-chisqm/chisqp; mytimeoffset=cm;
00933 }
00934 else{
00935 dir=+1.0; time_score=1-chisqp/chisqm; mytimeoffset=cp;
00936 }
00937 time_score=dir*time_score;
00938 }
00939
00940 MSG("AlgTrackAtNu", Msg::kDebug) << " timing parameters " << endl
00941 << " chisqp=" << chisqp << " (" << chisqndfp << ") "
00942 << " chisqm=" << chisqm << " (" << chisqndfm << ") " << endl
00943 << " score=" << time_score << endl;
00944
00945
00946
00947
00948
00949
00950
00951 MSG("AlgTrackAtNu", Msg::kDebug) << " *** vertex showers *** " << endl;
00952 Bool_t begvtxshw=0,endvtxshw=0;
00953 Int_t begvtxshwstrips=0,endvtxshwstrips=0;
00954 Bool_t begvtxreseeded=0,endvtxreseeded=0;
00955
00956 ShowerAtNu* begshwU = (ShowerAtNu*)(trku->GetBegVtxShower());
00957 ShowerAtNu* begshwV = (ShowerAtNu*)(trkv->GetBegVtxShower());
00958 if( begshwU && begshwV ){
00959 begvtxshw=1;
00960 begvtxshwstrips=(1+begshwU->GetHitLast())+(1+begshwV->GetHitLast());
00961 if(begshwU->GetReseedFlag() && begshwV->GetReseedFlag()) begvtxreseeded=1;
00962 }
00963
00964 ShowerAtNu* endshwU = (ShowerAtNu*)(trku->GetEndVtxShower());
00965 ShowerAtNu* endshwV = (ShowerAtNu*)(trkv->GetEndVtxShower());
00966 if( endshwU && endshwV ){
00967 endvtxshw=1;
00968 endvtxshwstrips=(1+endshwU->GetHitLast())+(1+endshwV->GetHitLast());
00969 if(endshwU->GetReseedFlag() && endshwV->GetReseedFlag()) endvtxreseeded=1;
00970 }
00971
00972 MSG("AlgTrackAtNu", Msg::kVerbose) << " vertex showers " << endl
00973 << " beg vtxshw strps = " << begvtxshwstrips
00974 << " end vtxshw strps = " << endvtxshwstrips << endl;
00975
00976
00977
00978
00979
00980
00981 MSG("AlgTrackAtNu", Msg::kDebug) << " *** slice stuff *** " << endl;
00982
00983 Int_t ntrkplns=0;
00984 Double_t trkph=0.0,shwph=0.0;
00985 Double_t assocph=0.0,totph=0.0;
00986 Double_t assocphlocal=0.0,totphlocal=0.0;
00987 Double_t begQmax=0.0,endQmax=0.0;
00988 Double_t begRmax=0.0,endRmax=0.0;
00989 Double_t begRdigits=0.0,endRdigits=0.0;
00990 Double_t begUmean=0.0,endUmean=0.0;
00991 Double_t begUwidth=0.0,endUwidth=0.0;
00992 Double_t begVmean=0.0,endVmean=0.0;
00993 Double_t begVwidth=0.0,endVwidth=0.0;
00994 Int_t begplndigits=bpln,endplndigits=epln;
00995 Int_t bstr,estr,xtalk;
00996 Double_t chg;
00997
00998 const CandSliceHandle* slice = dynamic_cast<const CandSliceHandle*>(track.GetCandSlice());
00999 TIter sliceitr(slice->GetDaughterIterator());
01000 while(CandStripHandle* strip = (CandStripHandle*)(sliceitr())){
01001 pln=strip->GetPlane(); chg=strip->GetCharge();
01002 xtalk=0;
01003 TIter digitr(strip->GetDaughterIterator());
01004 while(CandDeMuxDigitHandle* digit = dynamic_cast<CandDeMuxDigitHandle*>(digitr())){
01005 if( (digit->GetDeMuxDigitFlagWord()<8)
01006 && ( (digit->GetDeMuxDigitFlagWord() & CandDeMuxDigit::kXTalk)==(CandDeMuxDigit::kXTalk) ) ){
01007 xtalk=1;
01008 }
01009 }
01010 if( !xtalk && chg>2.0
01011 && pln>-1 && pln<500 ){
01012 fSliceStrpArr[pln].Add(strip);
01013 }
01014 }
01015
01016 TIter stripitr(track.GetDaughterIterator());
01017 while(CandStripHandle* strip = (CandStripHandle*)(stripitr())){
01018 pln=strip->GetPlane(); chg=strip->GetCharge();
01019 xtalk=0;
01020 TIter digitr(strip->GetDaughterIterator());
01021 while(CandDeMuxDigitHandle* digit = dynamic_cast<CandDeMuxDigitHandle*>(digitr())){
01022 if( (digit->GetDeMuxDigitFlagWord()<8)
01023 && ( (digit->GetDeMuxDigitFlagWord() & CandDeMuxDigit::kXTalk)==(CandDeMuxDigit::kXTalk) ) ){
01024 xtalk=1;
01025 }
01026 }
01027 if( !xtalk && chg>2.0
01028 && pln>-1 && pln<500 ){
01029 fTrackStrpArr[pln].Add(strip);
01030 }
01031 }
01032
01033
01034 for(i=0;i<500;i++){
01035 if( 1+fSliceStrpArr[i].GetLast()>0 ){
01036 bstr=-1; estr=-1;
01037 for(j=0;j<1+fTrackStrpArr[i].GetLast();j++){
01038 CandStripHandle* strip = (CandStripHandle*)(fTrackStrpArr[i].At(j));
01039 if(bstr<0||strip->GetStrip()<bstr) bstr=strip->GetStrip();
01040 if(estr<0||strip->GetStrip()>estr) estr=strip->GetStrip();
01041 trkph+=strip->GetCharge();
01042 }
01043 assocphlocal=0.0; totphlocal=0.0;
01044 for(j=0;j<1+fSliceStrpArr[i].GetLast();j++){
01045 CandStripHandle* strip = (CandStripHandle*)(fSliceStrpArr[i].At(j));
01046 if( bstr>-1 && estr>-1
01047 && strip->GetStrip()>bstr-2 && strip->GetStrip()<estr+2 ){
01048 assocphlocal+=strip->GetCharge(); assocph+=strip->GetCharge();
01049 }
01050 totphlocal+=strip->GetCharge(); totph+=strip->GetCharge();
01051 }
01052 if( totphlocal>0.0 && totphlocal<80.0
01053 && assocphlocal/totphlocal>0.8 ){
01054 ntrkplns++;
01055 }
01056 }
01057 }
01058 shwph=totph-trkph;
01059
01060
01061
01062
01063
01064
01065 if( atmosflag && fAtNuAnaOnOff ){
01066
01067 TObjArray tmpu,tmpv;
01068 Double_t dT,dU,dV,dR,dQ;
01069 Double_t myUwt2,myUwt,myUw;
01070 Double_t myVwt2,myVwt,myVw;
01071 Double_t Qmax,Rmax;
01072
01073
01074 Rmax=0.0;
01075 myUwt2=0.0; myUwt=0.0; myUw=0.0;
01076 myVwt2=0.0; myVwt=0.0; myVw=0.0;
01077 for(j=-4;j<=5;j++){
01078 tmpu.Clear(); tmpv.Clear();
01079
01080 if( bpln+j>0 && bpln+j<500
01081 && !(bpln<249&&bpln+j>249) && !(bpln>249&&bpln+j<249) ){
01082 for(k=0;k<1+fSliceStrpArr[bpln+j].GetLast();k++){
01083 CandStripHandle* strip = (CandStripHandle*)(fSliceStrpArr[bpln+j].At(k));
01084 dT=0.0; dQ=0.0;
01085 if(strip->GetPlaneView()==PlaneView::kU){
01086 dT=(strip->GetTPos())-(bvtxu+bdiru*(strip->GetZPos()-bvtxz));
01087 dQ=strip->GetCharge();
01088 myUwt2+=dQ*dT*dT; myUwt+=dQ*dT; myUw+=dQ;
01089 }
01090 if(strip->GetPlaneView()==PlaneView::kV){
01091 dT=(strip->GetTPos())-(bvtxv+bdirv*(strip->GetZPos()-bvtxz));
01092 dQ=strip->GetCharge();
01093 myVwt2+=dQ*dT*dT; myVwt+=dQ*dT; myVw+=dQ;
01094 }
01095 }
01096 }
01097
01098 if( bpln+j>0 && bpln+j<500
01099 && !(bpln<249&&bpln+j>249) && !(bpln>249&&bpln+j<249) ){
01100 for(k=0;k<1+fSliceStrpArr[bpln+j].GetLast();k++){
01101 CandStripHandle* strip = (CandStripHandle*)(fSliceStrpArr[bpln+j].At(k));
01102 if(strip->GetPlaneView()==PlaneView::kU){ tmpu.Add(strip); }
01103 if(strip->GetPlaneView()==PlaneView::kV){ tmpv.Add(strip); }
01104 }
01105 }
01106
01107 if( bpln+j-1>0 && bpln+j-1<500
01108 && !(bpln<249&&bpln+j-1>249) && !(bpln>249&&bpln+j-1<249) ){
01109 for(k=0;k<1+fSliceStrpArr[bpln+j-1].GetLast();k++){
01110 CandStripHandle* strip = (CandStripHandle*)(fSliceStrpArr[bpln+j-1].At(k));
01111 if(strip->GetPlaneView()==PlaneView::kU){ tmpu.Add(strip); }
01112 if(strip->GetPlaneView()==PlaneView::kV){ tmpv.Add(strip); }
01113 }
01114 }
01115
01116 if(1+tmpu.GetLast()>0 && 1+tmpv.GetLast()>0){
01117 for(k=0;k<1+tmpu.GetLast();k++){
01118 for(l=0;l<1+tmpv.GetLast();l++){
01119 CandStripHandle* strip = (CandStripHandle*)(tmpu.At(k));
01120 CandStripHandle* strip1 = (CandStripHandle*)(tmpv.At(l));
01121 dU=strip->GetTPos()-bvtxu; dV=strip1->GetTPos()-bvtxv;
01122 dR=sqrt(dU*dU+dV*dV); if(dR>Rmax) Rmax=dR;
01123 }
01124 }
01125 }
01126 }
01127
01128 begRmax=Rmax;
01129
01130 if(myUw>0.0){
01131 begUmean=myUwt/myUw;
01132 if((myUwt2/myUw)>(myUwt/myUw)*(myUwt/myUw)){
01133 begUwidth=sqrt((myUwt2/myUw)-(myUwt/myUw)*(myUwt/myUw));
01134 }
01135 }
01136
01137 if(myVw>0.0){
01138 begVmean=myVwt/myVw;
01139 if((myVwt2/myVw)>(myVwt/myVw)*(myVwt/myVw)){
01140 begVwidth=sqrt((myVwt2/myVw)-(myVwt/myVw)*(myVwt/myVw));
01141 }
01142 }
01143
01144
01145 Rmax=0.0;
01146 myUwt2=0.0; myUwt=0.0; myUw=0.0;
01147 myVwt2=0.0; myVwt=0.0; myVw=0.0;
01148 for(j=-4;j<=5;j++){
01149 tmpu.Clear(); tmpv.Clear();
01150
01151 if( epln+j>0 && epln+j<500
01152 && !(epln<249&&epln+j>249) && !(epln>249&&epln+j<249) ){
01153 for(k=0;k<1+fSliceStrpArr[epln+j].GetLast();k++){
01154 CandStripHandle* strip = (CandStripHandle*)(fSliceStrpArr[epln+j].At(k));
01155 dT=0.0; dQ=0.0;
01156 if(strip->GetPlaneView()==PlaneView::kU){
01157 dT=(strip->GetTPos())-(evtxu+ediru*(strip->GetZPos()-evtxz));
01158 dQ=strip->GetCharge();
01159 myUwt2+=dQ*dT*dT; myUwt+=dQ*dT; myUw+=dQ;
01160 }
01161 if(strip->GetPlaneView()==PlaneView::kV){
01162 dT=(strip->GetTPos())-(evtxv+edirv*(strip->GetZPos()-evtxz));
01163 dQ=strip->GetCharge();
01164 myVwt2+=dQ*dT*dT; myVwt+=dQ*dT; myVw+=dQ;
01165 }
01166 }
01167 }
01168
01169 if( epln+j>0 && epln+j<500
01170 && !(epln<249&&epln+j>249) && !(epln>249&&epln+j<249) ){
01171 for(k=0;k<1+fSliceStrpArr[epln+j].GetLast();k++){
01172 CandStripHandle* strip = (CandStripHandle*)(fSliceStrpArr[epln+j].At(k));
01173 if(strip->GetPlaneView()==PlaneView::kU){ tmpu.Add(strip); }
01174 if(strip->GetPlaneView()==PlaneView::kV){ tmpv.Add(strip); }
01175 }
01176 }
01177
01178 if( epln+j-1>0 && epln+j-1<500
01179 && !(epln<249&&epln+j-1>249) && !(epln>249&&epln+j-1<249) ){
01180 for(k=0;k<1+fSliceStrpArr[epln+j-1].GetLast();k++){
01181 CandStripHandle* strip = (CandStripHandle*)(fSliceStrpArr[epln+j-1].At(k));
01182 if(strip->GetPlaneView()==PlaneView::kU){ tmpu.Add(strip); }
01183 if(strip->GetPlaneView()==PlaneView::kV){ tmpv.Add(strip); }
01184 }
01185 }
01186
01187 if(1+tmpu.GetLast()>0 && 1+tmpv.GetLast()>0){
01188 for(k=0;k<1+tmpu.GetLast();k++){
01189 for(l=0;l<1+tmpv.GetLast();l++){
01190 CandStripHandle* strip = (CandStripHandle*)(tmpu.At(k));
01191 CandStripHandle* strip1 = (CandStripHandle*)(tmpv.At(l));
01192 dU=strip->GetTPos()-evtxu; dV=strip1->GetTPos()-evtxv;
01193 dR=sqrt(dU*dU+dV*dV); if(dR>Rmax) Rmax=dR;
01194 }
01195 }
01196 }
01197 }
01198
01199 endRmax=Rmax;
01200
01201 if(myUw>0.0){
01202 endUmean=myUwt/myUw;
01203 if((myUwt2/myUw)>(myUwt/myUw)*(myUwt/myUw)){
01204 endUwidth=sqrt((myUwt2/myUw)-(myUwt/myUw)*(myUwt/myUw));
01205 }
01206 }
01207
01208 if(myVw>0.0){
01209 endVmean=myVwt/myVw;
01210 if((myVwt2/myVw)>(myVwt/myVw)*(myVwt/myVw)){
01211 endVwidth=sqrt((myVwt2/myVw)-(myVwt/myVw)*(myVwt/myVw));
01212 }
01213 }
01214
01215
01216
01217 Qmax=0.0;
01218 for(j=-4;j<=5;j++){
01219 dQ=0.0;
01220 if( bpln+j>0 && bpln+j<500
01221 && !(bpln<249&&bpln+j>249) && !(bpln>249&&bpln+j<249) ){
01222 for(k=0;k<1+fSliceStrpArr[bpln+j].GetLast();k++){
01223 CandStripHandle* strip = (CandStripHandle*)(fSliceStrpArr[bpln+j].At(k));
01224 dQ+=strip->GetCharge();
01225 }
01226 }
01227 if(dQ>Qmax) Qmax=dQ;
01228 }
01229 begQmax=Qmax;
01230
01231 Qmax=0.0;
01232 for(j=-4;j<=5;j++){
01233 dQ=0.0;
01234 if( epln+j>0 && epln+j<500
01235 && !(epln<249&&epln+j>249) && !(epln>249&&epln+j<249) ){
01236 for(k=0;k<1+fSliceStrpArr[epln+j].GetLast();k++){
01237 CandStripHandle* strip = (CandStripHandle*)(fSliceStrpArr[epln+j].At(k));
01238 dQ+=strip->GetCharge();
01239 }
01240 }
01241 if(dQ>Qmax) Qmax=dQ;
01242 }
01243 endQmax=Qmax;
01244
01245
01246 Int_t plnM = 0, plnP = 0;
01247 Double_t totq,totqt;
01248 Double_t opos,tpos,upos,vpos,xpos,ypos,rpos;
01249 Double_t xm,xp,ym,yp,um,up,vm,vp;
01250
01251 begplndigits=bpln;
01252 begRdigits=4.0;
01253 if(bpln<249){
01254 plnM=begsm1; plnP=bpln+3; if(plnP>endsm1) plnP=endsm1;
01255 }
01256 if(bpln>249){
01257 plnM=begsm2; plnP=bpln+3; if(plnP>endsm2) plnP=endsm2;
01258 }
01259 for(i=plnM;i<plnP;i++){
01260 if( i-1>=0
01261 && 1+fSliceStrpArr[i].GetLast()>0 && 1+fSliceStrpArr[i-1].GetLast()>0 ){
01262 totqt=0.0; totq=0.0;
01263 for(j=0;j<1+fSliceStrpArr[i-1].GetLast();j++){
01264 CandStripHandle* strip = (CandStripHandle*)(fSliceStrpArr[i-1].At(j));
01265 totqt+=strip->GetTPos()*strip->GetCharge();
01266 totq+=strip->GetCharge();
01267 }
01268 if(totq>0.0){
01269 opos=totqt/totq;
01270 for(j=0;j<1+fSliceStrpArr[i].GetLast();j++){
01271 CandStripHandle* strip = (CandStripHandle*)(fSliceStrpArr[i].At(j));
01272 upos=0.0; vpos=0.0; tpos=strip->GetTPos();
01273 if(strip->GetPlaneView()==PlaneView::kU){ upos=tpos; vpos=opos; }
01274 if(strip->GetPlaneView()==PlaneView::kV){ vpos=tpos; upos=opos; }
01275 xpos=0.7071*(upos-vpos); ypos=0.7071*(upos+vpos); rpos=4.0;
01276 up=4.0-upos; if(up<rpos) rpos=up;
01277 um=4.0+upos; if(um<rpos) rpos=um;
01278 vp=4.0-vpos; if(vp<rpos) rpos=vp;
01279 vm=4.0+vpos; if(vm<rpos) rpos=vm;
01280 xp=4.0-xpos; if(xp<rpos) rpos=xp;
01281 xm=4.0+xpos; if(xm<rpos) rpos=xm;
01282 yp=4.0-ypos; if(yp<rpos) rpos=yp;
01283 ym=4.0+ypos; if(ym<rpos) rpos=ym;
01284 if(rpos<begRdigits) begRdigits=rpos;
01285 if(strip->GetPlane()<begplndigits) begplndigits=strip->GetPlane();
01286 }
01287 }
01288 }
01289 }
01290
01291 endplndigits=epln;
01292 endRdigits=4.0;
01293 if(epln<249){
01294 plnP=endsm1; plnM=epln-3; if(plnM<begsm1) plnM=begsm1;
01295 }
01296 if(epln>249){
01297 plnP=endsm2; plnM=epln-3; if(plnM<begsm2) plnM=begsm2;
01298 }
01299 for(i=plnM;i<plnP;i++){
01300 if( i+1<500
01301 && 1+fSliceStrpArr[i].GetLast()>0 && 1+fSliceStrpArr[i+1].GetLast()>0 ){
01302 totqt=0.0; totq=0.0;
01303 for(j=0;j<1+fSliceStrpArr[i+1].GetLast();j++){
01304 CandStripHandle* strip = (CandStripHandle*)(fSliceStrpArr[i+1].At(j));
01305 totqt+=strip->GetTPos()*strip->GetCharge();
01306 totq+=strip->GetCharge();
01307 }
01308 if(totq>0.0){
01309 opos=totqt/totq;
01310 for(j=0;j<1+fSliceStrpArr[i].GetLast();j++){
01311 CandStripHandle* strip = (CandStripHandle*)(fSliceStrpArr[i].At(j));
01312 upos=0.0; vpos=0.0; tpos=strip->GetTPos();
01313 if(strip->GetPlaneView()==PlaneView::kU){ upos=tpos; vpos=opos; }
01314 if(strip->GetPlaneView()==PlaneView::kV){ vpos=tpos; upos=opos; }
01315 xpos=0.7071*(upos-vpos); ypos=0.7071*(upos+vpos); rpos=4.0;
01316 up=4.0-upos; if(up<rpos) rpos=up;
01317 um=4.0+upos; if(um<rpos) rpos=um;
01318 vp=4.0-vpos; if(vp<rpos) rpos=vp;
01319 vm=4.0+vpos; if(vm<rpos) rpos=vm;
01320 xp=4.0-xpos; if(xp<rpos) rpos=xp;
01321 xm=4.0+xpos; if(xm<rpos) rpos=xm;
01322 yp=4.0-ypos; if(yp<rpos) rpos=yp;
01323 ym=4.0+ypos; if(ym<rpos) rpos=ym;
01324 if(rpos<endRdigits) endRdigits=rpos;
01325 if(strip->GetPlane()>endplndigits) endplndigits=strip->GetPlane();
01326 }
01327 }
01328 }
01329 }
01330
01331 MSG("AlgTrackAtNu",Msg::kDebug) << " beg vtx topology : Umean=" << begUmean << " Uwidth=" << begUwidth << " Vmean=" << begVmean << " Vwidth=" << begVwidth << " Rmax=" << begRmax << " Qmax=" << begQmax << endl;
01332
01333 MSG("AlgTrackAtNu",Msg::kDebug) << " end vtx topology : Umean=" << endUmean << " Uwidth=" << endUwidth << " Vmean=" << endVmean << " Vwidth=" << endVwidth << " Rmax=" << endRmax << " Qmax=" << endQmax << endl;
01334
01335 }
01336
01337 for(i=0;i<500;i++){
01338 fSliceStrpArr[i].Clear();
01339 fTrackStrpArr[i].Clear();
01340 }
01341
01342
01343
01344
01345
01346 MSG("AlgTrackAtNu",Msg::kDebug) << " *** calculate momentum *** " << endl;
01347
01348
01349 PlexStripEndId plexseid;
01350 FloatErr sigcorr,sigmip,sigmap;
01351 FloatErr lpos;
01352 for(pln=0;pln<npln;pln++){
01353 for(i=0;i<1+fHitArr[pln].GetLast();i++){
01354 HitAtNu* hit = (HitAtNu*)(fHitArr[pln].At(i));
01355 CandStripHandle* strip = (CandStripHandle*)(hit->GetCandStripHandle());
01356
01357 lpos = 0.0;
01358 if( strip->GetPlaneView()==PlaneView::kU
01359 || strip->GetPlaneView()==PlaneView::kX ) lpos = V[pln];
01360 if( strip->GetPlaneView()==PlaneView::kV
01361 || strip->GetPlaneView()==PlaneView::kY ) lpos = U[pln];
01362 lpos.SetError(0.1);
01363
01364 if(strip->GetNDigit(StripEnd::kPositive)>0){
01365 plexseid=strip->GetStripEndId(StripEnd::kPositive);
01366 sigcorr=strip->GetCharge(StripEnd::kPositive,CalDigitType::kSigCorr);
01367 sigmap=cal.GetAttenCorrectedTpos(sigcorr,lpos,plexseid);
01368 sigmip=cal.GetMIP(sigmap,plexseid);
01369 track.CalibrateSigMapped(plexseid.GetEncoded(),sigmap);
01370 track.CalibrateMIP(plexseid.GetEncoded(),sigmip);
01371 }
01372
01373 if(strip->GetNDigit(StripEnd::kNegative)>0){
01374 plexseid=strip->GetStripEndId(StripEnd::kNegative);
01375 sigcorr=strip->GetCharge(StripEnd::kNegative,CalDigitType::kSigCorr);
01376 sigmap=cal.GetAttenCorrectedTpos(sigcorr,lpos,plexseid);
01377 sigmip=cal.GetMIP(sigmap,plexseid);
01378 track.CalibrateSigMapped(plexseid.GetEncoded(),sigmap);
01379 track.CalibrateMIP(plexseid.GetEncoded(),sigmip);
01380 }
01381 }
01382 }
01383
01384
01385 dir=0.0;
01386 if(time_score>0.0) dir=1.0; if(time_score<0.0) dir=-1.0;
01387
01388
01389 Double_t range_thru_steel_xtra=0.0;
01390 Double_t range_thru_detector_xtra=0.0;
01391 if(range_thru_detector>0.0){
01392 range_thru_steel_xtra = 0.0127*(1.0/bdirz+1.0/edirz);
01393 range_thru_detector_xtra = 0.0296*(1.0/bdirz+1.0/edirz);
01394 range_thru_steel+=range_thru_steel_xtra;
01395 range_thru_detector+=range_thru_detector_xtra;
01396 }
01397
01398
01399
01400
01401
01402
01403
01404
01405 Double_t dens=840.0;
01406 for(i=0;i<npln;i++){
01407 pln=bpln+i;
01408 if(plnvuw[i]>-1){
01409 if(dir>=0.0){
01410 track.SetdS(pln,range_thru_detector-dS[i]-0.5*range_thru_detector_xtra);
01411 track.SetRange(pln,dens*(range_thru_steel-dSsteel[i]-0.5*range_thru_steel_xtra));
01412 }
01413 else{
01414 track.SetdS(pln,dS[i]+0.5*range_thru_detector_xtra);
01415 track.SetRange(pln,dens*(dSsteel[i]+0.5*range_thru_steel_xtra));
01416 }
01417 }
01418 }
01419
01420
01421 Double_t momentum=0.0,momentum_err=0.0;
01422 if(range_thru_steel>0.0){
01423 momentum = 1.52*range_thru_steel;
01424
01425
01426
01427
01428 momentum_err = 0.0;
01429 }
01430
01431
01432
01433
01434
01435
01436 MSG("AlgTrackAtNu",Msg::kDebug) << " *** set track variables *** " << endl;
01437
01438
01439 Bool_t reseedflag=0;
01440 if(trku->GetReseedFlag()&&trkv->GetReseedFlag()) reseedflag=1;
01441
01442
01443 track.SetMinPlane(bpln);
01444 track.SetMaxPlane(epln);
01445 track.SetTimeSlope((ctimeslope)/(3.0e8));
01446 track.SetTimeOffset((ctimeoffset)/(3.0e8));
01447 track.SetDirTimeSlope((dir)/(3.0e8));
01448 track.SetDirTimeOffset((mytimeoffset)/(3.0e8));
01449 track.SetDirTimeScatter(ctimescatter);
01450 track.SetDirTimeScore(time_score);
01451 track.SetTimeForwardFitRMS(chisqp);
01452 track.SetTimeForwardFitNDOF(chisqndfp);
01453 track.SetTimeBackwardFitRMS(chisqm);
01454 track.SetTimeBackwardFitNDOF(chisqndfm);
01455 track.SetTrkPH(trkph);
01456 track.SetShwPH(shwph);
01457 track.SetAssocTrkPH(assocph);
01458 track.SetAssocTrkPHfrac(assocph/totph);
01459 track.SetTrackLikePlanes(ntrkplns);
01460 track.SetLinearDirCosU(dir*lindiru);
01461 track.SetLinearDirCosV(dir*lindirv);
01462 track.SetLinearDirCosZ(dir*lindirz);
01463 track.SetLinearDirFitChisq(lindirerr);
01464 track.SetLinearDirFitNdf(linndf);
01465 track.SetRangeThruSteel(range_thru_steel);
01466 track.SetRangeThruDetector(range_thru_detector);
01467 track.SetMomentum(momentum);
01468 track.SetMomentumErr(momentum_err);
01469 track.SetReseedFlag(reseedflag);
01470 if(dir>=0.0){
01471 track.SetVtxU(bvtxu);
01472 track.SetVtxV(bvtxv);
01473 track.SetVtxZ(bvtxz);
01474 track.SetVtxR(begr);
01475 track.SetVtxT((mytimeoffset)/(3.0e8));
01476 track.SetVtxPlane(bpln+plnm);
01477 track.SetVtxPlaneDigits(begplndigits);
01478 track.SetVtxTrace(begtrace);
01479 track.SetVtxTraceZ(begtraceZ);
01480 track.SetVtxRdigits(begRdigits);
01481 track.SetVtxUwidth(begUwidth);
01482 track.SetVtxUmean(begUmean);
01483 track.SetVtxVwidth(begVwidth);
01484 track.SetVtxVmean(begVmean);
01485 track.SetVtxRmax(begRmax);
01486 track.SetVtxQmax(begQmax);
01487 track.SetDirCosU(bdiru);
01488 track.SetDirCosV(bdirv);
01489 track.SetDirCosZ(bdirz);
01490 track.SetEndU(evtxu);
01491 track.SetEndV(evtxv);
01492 track.SetEndZ(evtxz);
01493 track.SetEndR(endr);
01494 track.SetEndT((mytimeoffset+range_thru_detector)/(3.0e8));
01495 track.SetEndPlane(bpln+plnp);
01496 track.SetEndPlaneDigits(endplndigits);
01497 track.SetEndTrace(endtrace);
01498 track.SetEndTraceZ(endtraceZ);
01499 track.SetEndRdigits(endRdigits);
01500 track.SetEndUwidth(endUwidth);
01501 track.SetEndUmean(endUmean);
01502 track.SetEndVwidth(endVwidth);
01503 track.SetEndVmean(endVmean);
01504 track.SetEndRmax(endRmax);
01505 track.SetEndQmax(endQmax);
01506 track.SetEndDirCosU(ediru);
01507 track.SetEndDirCosV(edirv);
01508 track.SetEndDirCosZ(edirz);
01509 track.SetVtxShw(begvtxshw);
01510 track.SetVtxShwStrips(begvtxshwstrips);
01511 track.SetVtxShwReseedFlag(begvtxreseeded);
01512 track.SetEndShw(endvtxshw);
01513 track.SetEndShwStrips(endvtxshwstrips);
01514 track.SetEndShwReseedFlag(endvtxreseeded);
01515 }
01516 else{
01517 track.SetVtxU(evtxu);
01518 track.SetVtxV(evtxv);
01519 track.SetVtxZ(evtxz);
01520 track.SetVtxR(endr);
01521 track.SetVtxT((mytimeoffset-range_thru_detector)/(3.0e8));
01522 track.SetVtxPlane(bpln+plnp);
01523 track.SetVtxPlaneDigits(endplndigits);
01524 track.SetVtxTrace(endtrace);
01525 track.SetVtxTraceZ(endtraceZ);
01526 track.SetVtxRdigits(endRdigits);
01527 track.SetVtxUwidth(endUwidth);
01528 track.SetVtxUmean(endUmean);
01529 track.SetVtxVwidth(endVwidth);
01530 track.SetVtxVmean(endVmean);
01531 track.SetVtxRmax(endRmax);
01532 track.SetVtxQmax(endQmax);
01533 track.SetDirCosU(-ediru);
01534 track.SetDirCosV(-edirv);
01535 track.SetDirCosZ(-edirz);
01536 track.SetEndU(bvtxu);
01537 track.SetEndV(bvtxv);
01538 track.SetEndZ(bvtxz);
01539 track.SetEndR(begr);
01540 track.SetEndT((mytimeoffset)/(3.0e8));
01541 track.SetEndPlane(bpln+plnm);
01542 track.SetEndPlaneDigits(begplndigits);
01543 track.SetEndTrace(begtrace);
01544 track.SetEndTraceZ(begtraceZ);
01545 track.SetEndRdigits(begRdigits);
01546 track.SetEndUmean(begUmean);
01547 track.SetEndUwidth(begUwidth);
01548 track.SetEndVmean(begVmean);
01549 track.SetEndVwidth(begVwidth);
01550 track.SetEndRmax(begRmax);
01551 track.SetEndQmax(begQmax);
01552 track.SetEndDirCosU(-bdiru);
01553 track.SetEndDirCosV(-bdirv);
01554 track.SetEndDirCosZ(-bdirz);
01555 track.SetVtxShw(endvtxshw);
01556 track.SetVtxShwStrips(endvtxshwstrips);
01557 track.SetVtxShwReseedFlag(endvtxreseeded);
01558 track.SetEndShw(begvtxshw);
01559 track.SetEndShwStrips(begvtxshwstrips);
01560 track.SetEndShwReseedFlag(begvtxreseeded);
01561 }
01562
01563 MSG("AlgTrackAtNu", Msg::kDebug)
01564 << " *** CREATING NEW TRACK *** " << endl
01565 << " VtxPln = " << track.GetVtxPlane() << endl
01566 << " VtxPlnDigits = " << track.GetVtxPlaneDigits() << endl
01567 << " VtxU = " << track.GetVtxU() << endl
01568 << " VtxV = " << track.GetVtxV() << endl
01569 << " VtxZ = " << track.GetVtxZ() << endl
01570 << " VtxT = " << track.GetVtxT() << endl
01571 << " VtxR = " << track.GetVtxR() << endl
01572 << " VtxTrace = " << track.GetVtxTrace() << endl
01573 << " VtxTraceZ = " << track.GetVtxTraceZ() << endl
01574 << " VtxRdigits = " << track.GetVtxRdigits() << endl
01575 << " VtxUwidth = " << track.GetVtxUwidth() << endl
01576 << " VtxUmean = " << track.GetVtxUmean() << endl
01577 << " VtxVwidth = " << track.GetVtxVwidth() << endl
01578 << " VtxVmean = " << track.GetVtxVmean() << endl
01579 << " VtxRmax = " << track.GetVtxRmax() << endl
01580 << " VtxQmax = " << track.GetVtxQmax() << endl
01581 << " DirCosU = " << track.GetDirCosU() << endl
01582 << " DirCosV = " << track.GetDirCosV() << endl
01583 << " DirCosZ = " << track.GetDirCosZ() << endl
01584 << " EndPln = " << track.GetEndPlane() << endl
01585 << " EndPlnDigits = " << track.GetEndPlaneDigits() << endl
01586 << " EndU = " << track.GetEndU() << endl
01587 << " EndV = " << track.GetEndV() << endl
01588 << " EndZ = " << track.GetEndZ() << endl
01589 << " EndT = " << track.GetEndT() << endl
01590 << " EndR = " << track.GetEndR() << endl
01591 << " EndTrace = " << track.GetEndTrace() << endl
01592 << " EndTraceZ = " << track.GetEndTraceZ() << endl
01593 << " EndRdigits = " << track.GetEndRdigits() << endl
01594 << " EndUwidth = " << track.GetEndUwidth() << endl
01595 << " EndUmean = " << track.GetEndUmean() << endl
01596 << " EndVwidth = " << track.GetEndVwidth() << endl
01597 << " EndVmean = " << track.GetEndVmean() << endl
01598 << " EndRmax = " << track.GetEndRmax() << endl
01599 << " EndQmax = " << track.GetEndQmax() << endl
01600 << " EndDirCosU = " << track.GetEndDirCosU() << endl
01601 << " EndDirCosV = " << track.GetEndDirCosV() << endl
01602 << " EndDirCosZ = " << track.GetEndDirCosZ() << endl
01603 << " MinPlane = " << track.GetMinPlane() << endl
01604 << " MaxPlane = " << track.GetMaxPlane() << endl
01605 << " TrackLikePlanes = " << track.GetTrackLikePlanes() << endl
01606 << " TrkPH = " << track.GetTrkPH() << endl
01607 << " ShwPH = " << track.GetShwPH() << endl
01608 << " AssocPH = " << track.GetAssocTrkPH() << endl
01609 << " AssocPHfrac = " << track.GetAssocTrkPHfrac() << endl
01610 << " TimeSlope = " << track.GetTimeSlope() << endl
01611 << " TimeOffset = " << track.GetTimeOffset() << endl
01612 << " DirTimeSlope = " << track.GetDirTimeSlope() << endl
01613 << " DirTimeOffset = " << track.GetDirTimeOffset() << endl
01614 << " DirTimeScatter = " << track.GetDirTimeScatter() << endl
01615 << " DirTimeScore = " << track.GetDirTimeScore()<< endl
01616 << " TimeForwardFitRms = " << track.GetTimeForwardFitRMS() << endl
01617 << " TimeForwardFitNdf = " << track.GetTimeForwardFitNDOF() << endl
01618 << " TimeBackwardFitRms = " << track.GetTimeBackwardFitRMS() << endl
01619 << " TimeBackwardFitNdf = " << track.GetTimeBackwardFitNDOF() << endl
01620 << " VtxShw = " << track.GetVtxShw() << endl
01621 << " VtxShwStrips = " << track.GetVtxShwStrips() << endl
01622 << " VtxShwReseedFlag = " << track.GetVtxShwReseedFlag() << endl
01623 << " EndShw = " << track.GetEndShw() << endl
01624 << " EndShwStrips = " << track.GetEndShwStrips() << endl
01625 << " EndShwReseedFlag = " << track.GetEndShwReseedFlag() << endl
01626 << " LinearDirCosU = " << track.GetLinearDirCosU() << endl
01627 << " LinearDirCosV = " << track.GetLinearDirCosV() << endl
01628 << " LinearDirCosZ = " << track.GetLinearDirCosZ() << endl
01629 << " LinearDirFitChisq = " << track.GetLinearDirFitChisq() << endl
01630 << " LinearDirFitNdf = " << track.GetLinearDirFitNdf() << endl
01631 << " RangeThruSteel = " << track.GetRangeThruSteel() << endl
01632 << " Momentum = " << track.GetMomentum() << endl
01633 << " MomentumErr = " << track.GetMomentumErr() << endl
01634 << " ReseedFlag = " << track.GetReseedFlag() << endl;
01635
01636
01637 for(i=0;i<npln;i++){
01638 fHitArr[i].Clear();
01639 }
01640
01641 delete [] U;
01642 delete [] V;
01643 delete [] Z;
01644 delete [] Q;
01645 delete [] dS;
01646 delete [] dSsteel;
01647 delete [] Qm;
01648 delete [] Qp;
01649 delete [] CTm;
01650 delete [] CTp;
01651 delete [] Wm;
01652 delete [] Wp;
01653 delete [] plnvuw;
01654 delete [] plnnum;
01655
01656 return;
01657
01658 }
01659
01660 void AlgTrackAtNu::Trace(const char * ) const
01661 {
01662
01663 }
01664
01665
01666
