Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

AlgShowerAtNu Class Reference

#include <AlgShowerAtNu.h>

Inheritance diagram for AlgShowerAtNu:

AlgBase List of all members.

Public Member Functions

 AlgShowerAtNu ()
 ~AlgShowerAtNu ()
void RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
void Trace (const char *c) const

Private Attributes

TObjArray fHitArr [500]

Constructor & Destructor Documentation

AlgShowerAtNu::AlgShowerAtNu  ) 
 

Definition at line 33 of file AlgShowerAtNu.cxx.

00034 {
00035 
00036 }

AlgShowerAtNu::~AlgShowerAtNu  ) 
 

Definition at line 38 of file AlgShowerAtNu.cxx.

00039 {
00040 
00041 }


Member Function Documentation

void AlgShowerAtNu::RunAlg AlgConfig ac,
CandHandle ch,
CandContext cx
[virtual]
 

Implements AlgBase.

Definition at line 43 of file AlgShowerAtNu.cxx.

References CandHandle::AddDaughterLink(), UgliStripHandle::ClearFiber(), fHitArr, ShowerAtNu::GetAssocTrackAt(), ObjTrackAtNu::GetBegPlane(), ObjShowerAtNu::GetBegPlane(), ShowerAtNu::GetBegVtxFlag(), CandContext::GetCandRecord(), ShowerAtNu::GetCandSliceHandle(), HitAtNu::GetCandStripHandle(), CandStripHandle::GetCharge(), HitAtNu::GetCharge(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), CandShowerAtNuHandle::GetDirCosErrU(), CandShowerAtNuHandle::GetDirCosErrV(), CandRecoHandle::GetDirCosU(), CandRecoHandle::GetDirCosV(), CandRecoHandle::GetDirCosZ(), CandShowerAtNuHandle::GetDirTimeScore(), CandShowerAtNuHandle::GetDirVtxShwScore(), Registry::GetDouble(), PlexSEIdAltL::GetEnd(), ObjTrackAtNu::GetEndPlane(), ObjShowerAtNu::GetEndPlane(), ShowerAtNu::GetEndVtxFlag(), CandShowerHandle::GetEnergy(), UgliStripHandle::GetHalfLength(), ObjAtNu::GetHitAt(), ObjAtNu::GetHitLast(), CandShowerAtNuHandle::GetMaxPlane(), CandShowerAtNuHandle::GetMinPlane(), CandShowerAtNuHandle::GetNPlanes(), CandShowerAtNuHandle::GetNStrips(), HitAtNu::GetPlane(), CandStripHandle::GetPlaneView(), HitAtNu::GetPlaneView(), CandDigitHandle::GetPlexSEIdAltL(), CandShowerAtNuHandle::GetReseedFlag(), ObjShowerAtNu::GetReseedFlag(), HitAtNu::GetShwFlag(), UgliGeomHandle::GetSteelPlnHandle(), CandStripHandle::GetStripEndId(), UgliGeomHandle::GetStripHandle(), CandStripHandle::GetTime(), CandRecoHandle::GetTimeOffset(), CandRecoHandle::GetTimeSlope(), HitAtNu::GetTPos(), RecMinos::GetVldContext(), CandRecoHandle::GetVtxPlane(), CandShowerAtNuHandle::GetVtxR(), CandShowerAtNuHandle::GetVtxShw(), CandRecoHandle::GetVtxT(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), CandRecoHandle::GetVtxZ(), UgliSteelPlnHandle::GetZ0(), HitAtNu::GetZPos(), UgliSteelPlnHandle::IsValid(), MSG, CandRecoHandle::SetCandSlice(), CandShowerAtNuHandle::SetDirCosErrU(), CandShowerAtNuHandle::SetDirCosErrV(), CandRecoHandle::SetDirCosU(), CandRecoHandle::SetDirCosV(), CandRecoHandle::SetDirCosZ(), CandShowerAtNuHandle::SetDirTimeScore(), CandShowerAtNuHandle::SetDirVtxShwScore(), CandShowerHandle::SetEnergy(), CandShowerAtNuHandle::SetMaxPlane(), CandShowerAtNuHandle::SetMinPlane(), CandShowerAtNuHandle::SetNPlanes(), CandShowerAtNuHandle::SetNStrips(), CandShowerAtNuHandle::SetReseedFlag(), CandRecoHandle::SetTimeOffset(), CandRecoHandle::SetTimeSlope(), CandShowerHandle::SetU(), CandShowerHandle::SetV(), CandRecoHandle::SetVtxPlane(), CandShowerAtNuHandle::SetVtxR(), CandShowerAtNuHandle::SetVtxShw(), CandRecoHandle::SetVtxT(), CandRecoHandle::SetVtxU(), CandRecoHandle::SetVtxV(), CandRecoHandle::SetVtxZ(), and UgliStripHandle::WlsPigtail().

00044 {
00045 
00046   MSG("AlgShowerAtNu", Msg::kDebug) << "AlgShowerAtNu::RunAlg(...)" << endl;
00047 
00048   // Get CandShowerAtNuHandle
00049   CandShowerAtNuHandle& shower = dynamic_cast<CandShowerAtNuHandle&>(ch);
00050 
00051   // Unpack AlgConfig
00052   Double_t fFibreIndex;
00053   fFibreIndex = ac.GetDouble("FibreIndex");
00054   MSG("AlgShowerAtNu", Msg::kDebug) << " FibreIndex=" << fFibreIndex << endl;
00055 
00056   // Unpack CandContext
00057   const TObjArray *arr = dynamic_cast<const TObjArray*>(cx.GetDataIn());
00058   ShowerAtNu* shwu = (ShowerAtNu*)(arr->At(0));
00059   ShowerAtNu* shwv = (ShowerAtNu*)(arr->At(1)); 
00060   shower.SetCandSlice(shwu->GetCandSliceHandle());
00061 
00062   CandRecord* candrec = (CandRecord*)(cx.GetCandRecord());
00063   VldContext *vldc = (VldContext*)(candrec->GetVldContext());
00064   UgliGeomHandle ugh(*vldc);
00065 
00066   // Get Geometry
00067   Int_t atmosflag,modtype,modnum,modpln,modstr;
00068   switch(candrec->GetVldContext()->GetDetector()){
00069     case DetectorType::kFar:
00070       atmosflag=1; modtype=1; modnum=2; modpln=248; modstr=192; break;
00071     case DetectorType::kNear:
00072       atmosflag=0; modtype=2; modnum=1; modpln=282; modstr=96; break;
00073     case DetectorType::kCalib:
00074       atmosflag=0; modtype=3; modnum=1; modpln=60; modstr=24; break;
00075     default:
00076       atmosflag=0; modtype=-1; modnum=0; modpln=0; modstr=0; break;
00077   }
00078 
00079   Int_t i,j,k;
00080   Int_t vuw;
00081   Int_t pln,bpln,epln,npln;
00082   Double_t begZsm1=-9999.9,endZsm1=-9999.9;
00083   Double_t begZsm2=-9999.9,endZsm2=-9999.9;
00084   Int_t begsm1=-1,endsm1=-1,begsm2=-1,endsm2=-1;
00085 
00086   for(j=0;j<486;j++){
00087     PlexPlaneId myplaneid(vldc->GetDetector(),j,1);
00088     UgliSteelPlnHandle myplanehandle = ugh.GetSteelPlnHandle(myplaneid);
00089     if(myplanehandle.IsValid()){
00090       if(j<=modpln && begZsm1<0.0){ begsm1=j; begZsm1=myplanehandle.GetZ0(); }
00091       if(j<=modpln && myplanehandle.GetZ0()>endZsm1){ endsm1=j; endZsm1=myplanehandle.GetZ0(); }
00092       if(j>modpln && begZsm2<0.0){ begsm2=j; begZsm2=myplanehandle.GetZ0(); }
00093       if(j>modpln && myplanehandle.GetZ0()>endZsm2){ endsm2=j; endZsm2=myplanehandle.GetZ0(); }
00094     }
00095   }
00096 
00097   MSG("AlgShowerAtNu", Msg::kDebug) << " *** detector geometry *** " << endl
00098                                    << "  SM1 :  pln " << begsm1 << " -> " << endsm1
00099                                    << "    Z " << begZsm1 << " -> " << endZsm1 << endl
00100                                    << "  SM2 :  pln " << begsm2 << " -> " << endsm2 
00101                                    << "    Z " << begZsm2 << " -> " << endZsm2 << endl;
00102 
00103 
00104   if(shwu->GetBegPlane()<shwv->GetBegPlane()) bpln=shwu->GetBegPlane(); else bpln=shwv->GetBegPlane();
00105   if(shwu->GetEndPlane()>shwv->GetEndPlane()) epln=shwu->GetEndPlane(); else epln=shwv->GetEndPlane();
00106 
00107   // ***********************************************
00108   // * C A L C U L A T E   C O - O R D I N A T E S *
00109   // ***********************************************
00110 
00111   MSG("AlgShowerAtNu", Msg::kDebug) << " *** calculate co-ordinates *** " << endl;
00112 
00113   npln = 1+epln-bpln;  
00114   Double_t* Z = new Double_t[npln];
00115   Double_t* U = new Double_t[npln];
00116   Double_t* V = new Double_t[npln];
00117   Double_t* CT = new Double_t[npln];
00118   Double_t* Q = new Double_t[npln];
00119   Int_t* plnvuw = new Int_t[npln]; 
00120   Int_t* plnnum = new Int_t[npln];
00121   for(i=0;i<npln;i++){ 
00122     plnvuw[i]=-1; plnnum[i]=0;
00123     U[i]=0.0; V[i]=0.0; Z[i]=0.0; 
00124     Q[i]=0.0; CT[i]=0.0; 
00125   }
00126 
00127   // loop over hits and fill arrays
00128   for(i=0;i<1+shwu->GetHitLast();i++){
00129     HitAtNu* hit = (HitAtNu*)(shwu->GetHitAt(i));
00130     CandStripHandle* strip = (CandStripHandle*)(hit->GetCandStripHandle());
00131     pln = hit->GetPlane()-bpln;
00132     fHitArr[pln].Add(hit);
00133     if(plnvuw[pln]<0){
00134       plnvuw[pln]=hit->GetPlaneView();
00135       Z[pln]=hit->GetZPos();
00136     }
00137     U[pln]+=hit->GetTPos()*hit->GetCharge();
00138     Q[pln]+=hit->GetCharge();
00139     shower.AddDaughterLink(*strip);
00140     plnnum[pln]++;
00141   }
00142 
00143   for(i=0;i<1+shwv->GetHitLast();i++){
00144     HitAtNu* hit = (HitAtNu*)(shwv->GetHitAt(i));
00145     CandStripHandle* strip = (CandStripHandle*)(hit->GetCandStripHandle());
00146     pln = hit->GetPlane()-bpln;
00147     fHitArr[pln].Add(hit);
00148     if(plnvuw[pln]<0){
00149       plnvuw[pln]=hit->GetPlaneView();
00150       Z[pln]=hit->GetZPos();
00151     }
00152     V[pln]+=hit->GetTPos()*hit->GetCharge();
00153     Q[pln]+=hit->GetCharge();
00154     shower.AddDaughterLink(*strip);
00155     plnnum[pln]++;
00156   }
00157 
00158   Int_t nplanes=0,nstrips=0;
00159   for(i=0;i<npln;i++){
00160     if(plnvuw[i]>-1){
00161       if(Q[i]>0.0){
00162         U[i]=U[i]/Q[i];
00163         V[i]=V[i]/Q[i];
00164         nstrips+=plnnum[i];
00165         nplanes++;
00166       }
00167     }
00168   }
00169 
00170   // estimate transverse co-ordinates
00171   Int_t km,kp;
00172   Double_t q,sq,sqt;
00173   for(pln=0;pln<npln;pln++){
00174     if(plnvuw[pln]>-1){
00175       vuw=plnvuw[pln];
00176       sq=0.0; sqt=0.0; 
00177       k=0; km=-1;
00178       while(pln-k>0 && km<0){
00179         k++; if(plnvuw[pln-k]>-1 && plnvuw[pln-k]!=vuw) km=k;
00180       }
00181       k=0; kp=-1;
00182       while(pln+k<npln-1 && kp<0){
00183         k++; if(plnvuw[pln+k]>-1 && plnvuw[pln+k]!=vuw) kp=k; 
00184       }
00185 
00186       if(km>0){
00187         if(vuw==0){ q=1.0/km; sq+=q*Q[pln-km]; sqt+=q*Q[pln-km]*V[pln-km]; }
00188         if(vuw==1){ q=1.0/km; sq+=q*Q[pln-km]; sqt+=q*Q[pln-km]*U[pln-km]; }
00189       }
00190       if(kp>0){
00191         if(vuw==0){ q=1.0/kp; sq+=q*Q[pln+kp]; sqt+=q*Q[pln+kp]*V[pln+kp]; }
00192         if(vuw==1){ q=1.0/kp; sq+=q*Q[pln+kp]; sqt+=q*Q[pln+kp]*U[pln+kp]; }
00193       }
00194 
00195       if(sq>0.0){
00196         if(vuw==0){ V[pln]=sqt/sq; } if(vuw==1){ U[pln]=sqt/sq; }
00197       }
00198       else{
00199         if(vuw==0){ V[pln]=0.0; } if(vuw==1){ U[pln]=0.0; }
00200       }
00201     }
00202   }
00203 
00204   // get timing information
00205   Double_t tpos,fibre,corr;
00206   Double_t tot_evtchg=0.0,tot_evtchg_fibrecorr=0.0;
00207   Double_t digchg,tot_digchg;
00208   Double_t ctime,tot_ctime;
00209   Double_t indx=fFibreIndex;
00210   for(pln=0;pln<npln;pln++){
00211     tot_ctime=0.0; tot_digchg=0.0;
00212     for(i=0;i<1+fHitArr[pln].GetLast();i++){
00213       HitAtNu* hit = (HitAtNu*)(fHitArr[pln].At(i));
00214 
00215       // timing/charge
00216       CandStripHandle* strip = (CandStripHandle*)(hit->GetCandStripHandle());
00217       if(strip->GetPlaneView()==PlaneView::kU) tpos = -V[pln];
00218       if(strip->GetPlaneView()==PlaneView::kV) tpos = U[pln];
00219       PlexStripEndId stripid = strip->GetStripEndId();
00220       UgliStripHandle striphandle = ugh.GetStripHandle(stripid);
00221       TIter digitr(strip->GetDaughterIterator());
00222       while(CandDigitHandle* digit = dynamic_cast<CandDigitHandle*>(digitr())){
00223         ctime=0.0; fibre=0.0; digchg=0.0;
00224         if(digit->GetPlexSEIdAltL().GetEnd()==StripEnd::kPositive){ fibre = striphandle.ClearFiber(StripEnd::kPositive)+striphandle.WlsPigtail(StripEnd::kPositive)+striphandle.GetHalfLength()-tpos; ctime = 3.0e8*strip->GetTime(StripEnd::kPositive); digchg=strip->GetCharge(StripEnd::kPositive); }
00225         if(digit->GetPlexSEIdAltL().GetEnd()==StripEnd::kNegative){ fibre = striphandle.ClearFiber(StripEnd::kNegative)+striphandle.WlsPigtail(StripEnd::kNegative)+striphandle.GetHalfLength()+tpos; ctime = 3.0e8*strip->GetTime(StripEnd::kNegative); digchg=strip->GetCharge(StripEnd::kNegative); }
00226         ctime = ctime-indx*fibre; corr = exp(-(8.0-fibre)/7.05); 
00227         tot_ctime += digchg*ctime; tot_digchg += digchg;
00228         tot_evtchg += digchg; tot_evtchg_fibrecorr += digchg*corr;   
00229       }
00230     }
00231     if(tot_digchg>0.0){
00232       CT[pln] = tot_ctime/tot_digchg; 
00233     }
00234   }
00235   
00236   // set map values
00237   for(i=0;i<npln;i++){
00238     pln=bpln+i;
00239     if(plnvuw[i]>-1){
00240       shower.SetU(pln,U[i]);
00241       shower.SetV(pln,V[i]);    
00242     }
00243   }
00244 
00245   MSG("AlgShowerAtNu", Msg::kDebug) << " *** shower hits *** " << endl;
00246   for(pln=0;pln<npln;pln++){
00247     if( plnvuw[pln]>-1 ){
00248       MSG("AlgShowerAtNu",Msg::kDebug)
00249         << " pln=" << pln << " vuw=" << plnvuw[pln]
00250         << " z=" << Z[pln] << " u=" << U[pln] << " v=" << V[pln] << endl;
00251     }
00252   }
00253 
00254 
00255   // ***********************************
00256   // * D E T E R M I N E   V E R T E X * 
00257   // ***********************************
00258 
00259   MSG("AlgShowerAtNu", Msg::kDebug) << " *** determine vertex  *** " << endl;
00260 
00261   Bool_t vtxshw=0;
00262   Double_t trkscr=0.0,shwscr=0.0,dirscr=0.0;
00263   Double_t vtxu=0.0,vtxv=0.0,vtxx=0.0,vtxy=0.0,vtxz=0.0; 
00264   Int_t vtxpln=-1,trkbpln=-1,trkepln=-1,shwbpln=-1,shwepln=-1;
00265 
00266   // determine vertex plane
00267   Double_t totchgpln=0.0,totchg=0.0;
00268   for(i=0;i<npln;i++){
00269     for(j=0;j<1+fHitArr[i].GetLast();j++){
00270       HitAtNu* hit = (HitAtNu*)(fHitArr[i].At(j));
00271       if(hit->GetShwFlag()>1){
00272         totchg+=hit->GetCharge();
00273         totchgpln+=hit->GetCharge()*hit->GetPlane(); 
00274         if(shwbpln<0||hit->GetPlane()<shwbpln) shwbpln=hit->GetPlane();
00275         if(shwepln<0||hit->GetPlane()>shwepln) shwepln=hit->GetPlane(); 
00276       } 
00277     }
00278   }
00279   if(totchg>0.0){
00280     vtxpln=(Int_t)(totchgpln/totchg);
00281   }
00282   else{
00283     vtxpln=(Int_t)(0.5*(bpln+epln));
00284     shwbpln=bpln; shwepln=epln;
00285   }
00286 
00287   TrackAtNu* vtxtrkU = (TrackAtNu*)(shwu->GetAssocTrackAt(0));
00288   TrackAtNu* vtxtrkV = (TrackAtNu*)(shwv->GetAssocTrackAt(0));
00289   if(vtxtrkU && vtxtrkV){
00290     vtxshw=1;
00291     if(vtxtrkU->GetBegPlane()<vtxtrkV->GetBegPlane()) trkbpln=vtxtrkU->GetBegPlane(); else trkbpln=vtxtrkV->GetBegPlane();
00292     if(vtxtrkU->GetEndPlane()>vtxtrkV->GetEndPlane()) trkepln=vtxtrkU->GetEndPlane(); else trkepln=vtxtrkV->GetEndPlane();
00293     trkscr = (trkbpln+trkepln-2.0*vtxpln)/(1.0+trkepln-trkbpln);
00294     if(trkscr>0.0){ trkscr = (shwbpln+shwepln-2.0*trkbpln)/(1.0+shwepln-shwbpln); }
00295     if(trkscr<0.0){ trkscr = (shwbpln+shwepln-2.0*trkepln)/(1.0+shwepln-shwbpln); }
00296     if( ( shwu->GetBegVtxFlag() && shwv->GetBegVtxFlag())
00297      && !( shwu->GetEndVtxFlag() && shwv->GetEndVtxFlag()) ) vtxpln=trkbpln;
00298     if( ( shwu->GetEndVtxFlag() && shwv->GetEndVtxFlag())
00299      && !( shwu->GetBegVtxFlag() && shwv->GetBegVtxFlag()) ) vtxpln=trkepln;
00300   }
00301 
00302   // determine vertex position
00303   pln=vtxpln-bpln;
00304   if( pln>-1 && pln<npln
00305    && plnvuw[pln]>-1 ){
00306     vtxu=U[pln]; vtxv=V[pln]; vtxz=Z[pln];
00307   }
00308   else{
00309     k=0; km=-1;
00310     while(pln-k>0 && km<0){
00311       k++; if(pln-k<npln && plnvuw[pln-k]>-1) km=k;
00312     }
00313     k=0; kp=-1;
00314     while(pln+k<npln-1 && kp<0){
00315       k++; if(pln+k>-1 && plnvuw[pln+k]>-1) kp=k; 
00316     }
00317     if(km>0 && kp>0){
00318       vtxu=(km*U[pln+kp]+kp*U[pln-km])/(km+kp); 
00319       vtxv=(km*V[pln+kp]+kp*V[pln-km])/(km+kp); 
00320       vtxz=(km*Z[pln+kp]+kp*Z[pln-km])/(km+kp); 
00321     }
00322     else{
00323       if(km>0){ vtxu=U[pln-km]; vtxv=V[pln-km]; vtxz=Z[pln-km]; } 
00324       if(kp>0){ vtxu=U[pln+kp]; vtxv=V[pln+kp]; vtxz=Z[pln+kp]; }
00325     }
00326   }
00327 
00328   vtxx=0.7071*(vtxu-vtxv); vtxy=0.7071*(vtxu+vtxv); 
00329 
00330   MSG("AlgShowerAtNu", Msg::kVerbose) 
00331     << "   ... vertex: "
00332     << " plane = " << vtxpln
00333     << " position = (" << vtxu << "," << vtxv << "," << vtxz << ")" << endl;
00334 
00335 
00336   // *************************************
00337   // * C O N T A I N M E N T   S T U F F *
00338   // *************************************
00339 
00340   MSG("AlgShowerAtNu", Msg::kDebug) << " *** containment stuff *** " << endl;
00341 
00342   Double_t vtxr=0.0;
00343 
00344   if(atmosflag){
00345 
00346     // closest distance to edge
00347     Double_t xm,xp,ym,yp,um,up,vm,vp;
00348     Double_t rmin;
00349 
00350     rmin=4.0; 
00351     up=4.0-vtxu; if(up<rmin) rmin=up;
00352     um=4.0+vtxu; if(um<rmin) rmin=um;
00353     vp=4.0-vtxv; if(vp<rmin) rmin=vp;
00354     vm=4.0+vtxv; if(vm<rmin) rmin=vm;
00355     xp=4.0-vtxx; if(xp<rmin) rmin=xp;
00356     xm=4.0+vtxx; if(xm<rmin) rmin=xm;
00357     yp=4.0-vtxy; if(yp<rmin) rmin=yp;
00358     ym=4.0+vtxy; if(ym<rmin) rmin=ym;
00359     vtxr=rmin;
00360 
00361   }
00362 
00363   //*****************************************
00364   //* D E T E R M I N E   D I R E C T I O N *
00365   //*****************************************
00366 
00367   MSG("AlgShowerAtNu", Msg::kDebug) << " *** determine direction *** " << endl;
00368 
00369   // measure overall timeslope
00370   Double_t ctimeslope=0.0,ctimeoffset=0.0,ctimeaverage=0.0,ctimescatter=0.0;
00371   Double_t sw,swz,swt,swzz,swzt,swtt;
00372   Double_t du,dv,dz;
00373   Double_t z,t,w;
00374   Double_t diru,dirv,dirz,dirs;
00375   Double_t direrru,direrrv;
00376   Double_t offset;
00377   Int_t flag,npts;
00378 
00379   npts=0;
00380   sw=0.0; swz=0.0; swzz=0.0; swt=0.0; swtt=0.0; swzt=0.0;
00381   for(i=0;i<npln;i++){
00382     if( plnvuw[i]>-1 ){
00383       z=Z[i]; t=CT[i]; w=Q[i];
00384       swz+=w*z; swzz+=w*z*z; swt+=w*t; swtt+=w*t*t; swzt+=w*z*t;
00385       sw+=w; npts++;
00386     }
00387   }
00388   if(npts>1){
00389     ctimeslope=(sw*swzt-swz*swt)/(sw*swzz-swz*swz); 
00390     ctimeoffset=(swt*swzz-swz*swzt)/(sw*swzz-swz*swz);
00391     ctimeaverage=(swt/sw);
00392     if( (swtt/sw)-(swt/sw)*(swt/sw)>0.0 && (swzz/sw)-(swz/sw)*(swz/sw)>0.0 ){
00393       ctimescatter=((swzt/sw)-(swz/sw)*(swt/sw))/((sqrt((swzz/sw)-(swz/sw)*(swz/sw)))*(sqrt((swtt/sw)-(swt/sw)*(swt/sw))));
00394     }
00395   }
00396   offset=ctimeaverage;
00397   shwscr=ctimescatter;
00398 
00399   if(vtxshw) dirscr=trkscr; else dirscr=shwscr;
00400 
00401   npts=0;
00402   diru=0.0; direrru=0.0;
00403   if(1+shwu->GetEndPlane()-shwu->GetBegPlane()>1){
00404     sw=0.0; swz=0.0; swt=0.0; 
00405     swzz=0.0; swzt=0.0; swtt=0.0;
00406     for(i=0;i<npln;i++){
00407       flag=0;
00408       for(j=0;j<1+fHitArr[i].GetLast();j++){
00409         HitAtNu* hit = (HitAtNu*)(fHitArr[i].At(j));
00410         if(hit->GetShwFlag()>1 && hit->GetPlaneView()==0){
00411           du=hit->GetTPos()-vtxu; dz=hit->GetZPos()-vtxz; w=hit->GetCharge();
00412           swz+=w*dz; swt+=w*du; 
00413           swzz+=w*dz*dz; swzt+=w*dz*du; swtt+=w*du*du;
00414           sw+=w; flag=1;
00415         }           
00416       }
00417       if(flag) npts++;
00418     }
00419     if(npts>1){
00420       diru=(sw*swzt-swz*swt)/(sw*swzz-swz*swz);
00421       if( (swtt/sw)-(swt/sw)*(swt/sw)>0.0 && (swzz/sw)-(swz/sw)*(swz/sw)>0.0 ){
00422         direrru=((swzt/sw)-(swz/sw)*(swt/sw))/((sqrt((swzz/sw)-(swz/sw)*(swz/sw)))*(sqrt((swtt/sw)-(swt/sw)*(swt/sw))));
00423       }
00424     }
00425   }
00426 
00427   npts=0;
00428   dirv=0.0; direrrv=0.0;
00429   if(1+shwv->GetEndPlane()-shwv->GetBegPlane()>1){
00430     sw=0.0; swz=0.0; swt=0.0; 
00431     swzz=0.0; swzt=0.0; swtt=0.0;
00432     for(i=0;i<npln;i++){
00433       flag=0;
00434       for(j=0;j<1+fHitArr[i].GetLast();j++){
00435         HitAtNu* hit = (HitAtNu*)(fHitArr[i].At(j));
00436         if(hit->GetShwFlag()>1 && hit->GetPlaneView()==1){
00437           dv=hit->GetTPos()-vtxv; dz=hit->GetZPos()-vtxz; w=hit->GetCharge();
00438           swz+=w*dz; swt+=w*dv; 
00439           swzz+=w*dz*dz; swzt+=w*dz*dv; swtt+=w*dv*dv;
00440           sw+=w; flag=1;
00441         }
00442       }
00443       if(flag) npts++;
00444     }
00445     if(npts>1){
00446       dirv=(sw*swzt-swz*swt)/(sw*swzz-swz*swz);
00447       if( (swtt/sw)-(swt/sw)*(swt/sw)>0.0 && (swzz/sw)-(swz/sw)*(swz/sw)>0.0 ){
00448         direrrv=((swzt/sw)-(swz/sw)*(swt/sw))/((sqrt((swzz/sw)-(swz/sw)*(swz/sw)))*(sqrt((swtt/sw)-(swt/sw)*(swt/sw))));
00449       }
00450     }
00451   }
00452     
00453   dirs=sqrt(diru*diru+dirv*dirv+1.0);
00454   diru=diru/dirs; dirv=dirv/dirs; dirz=1.0/dirs;
00455 
00456   MSG("AlgShowerAtNu", Msg::kVerbose) 
00457     << "   ... direction: "
00458     << "  (" << diru << "," << dirv << "," << dirz << ")" << endl;
00459 
00460 
00461   // ***********************************
00462   // * C A L C U L A T E   E N E R G Y *
00463   // ***********************************
00464 
00465   MSG("AlgShowerAtNu", Msg::kDebug) << " *** calculate energy *** " << endl;
00466 
00467   // calculate energy from pulse height
00468   Double_t energy;
00469   Double_t tot_ph;
00470   if(tot_evtchg>0.0){
00471     tot_ph = tot_evtchg; 
00472     energy = tot_ph; // NEEDS TUNING
00473     energy = 0.3 + tot_ph/150.0;
00474   }
00475 
00476   // *******************************************
00477   // * S E T   S H O W E R   V A R I A B L E S *
00478   // *******************************************
00479 
00480   // set reseed flag
00481   Bool_t reseedflag=0;
00482   if(shwu->GetReseedFlag()||shwv->GetReseedFlag()) reseedflag=1; 
00483    
00484   // determine overall direction 
00485   Double_t dir=1.0;
00486   if(dirscr>0.0) dir=1.0; if(dirscr<0.0) dir=-1.0;  
00487 
00488   // set shower variables
00489   shower.SetMinPlane(bpln);
00490   shower.SetMaxPlane(epln);
00491   shower.SetNPlanes(nplanes);
00492   shower.SetNStrips(nstrips);
00493   shower.SetTimeSlope((dir)/(3.0e8)); 
00494   shower.SetTimeOffset((offset)/(3.0e8)); 
00495   shower.SetDirTimeScore(shwscr);
00496   shower.SetDirVtxShwScore(dirscr);
00497   shower.SetVtxShw(vtxshw); 
00498   shower.SetVtxU(vtxu);
00499   shower.SetVtxV(vtxv);
00500   shower.SetVtxZ(vtxz);
00501   shower.SetVtxR(vtxr);
00502   shower.SetVtxT((offset)/(3.0e8));
00503   shower.SetVtxPlane(vtxpln);
00504   shower.SetDirCosU(dir*diru);
00505   shower.SetDirCosV(dir*dirv);
00506   shower.SetDirCosZ(dir*dirz);
00507   shower.SetDirCosErrU(dir*direrru);
00508   shower.SetDirCosErrV(dir*direrrv);
00509   shower.SetEnergy(energy);
00510   shower.SetReseedFlag(reseedflag);
00511 
00512   MSG("AlgShowerAtNu", Msg::kDebug) 
00513     << " *** CREATING NEW SHOWER *** " << endl
00514     << "  MinPlane = " << shower.GetMinPlane() << endl
00515     << "  MaxPlane = " << shower.GetMaxPlane() << endl
00516     << "  NPlanes = " << shower.GetNPlanes() << endl
00517     << "  NStrips = " << shower.GetNStrips() << endl
00518     << "  TimeSlope = " << shower.GetTimeSlope() << endl
00519     << "  TimeOffset = " << shower.GetTimeOffset() << endl
00520     << "  DirTimeScore = " << shower.GetDirTimeScore() << endl
00521     << "  DirVtxShwScore= " << shower.GetDirVtxShwScore() << endl
00522     << "  VtxShw = " << shower.GetVtxShw() << endl
00523     << "  VtxU = " << shower.GetVtxU() << endl
00524     << "  VtxV = " << shower.GetVtxV() << endl
00525     << "  VtxZ = " << shower.GetVtxZ() << endl
00526     << "  VtxR = " << shower.GetVtxR() << endl
00527     << "  VtxT = " << shower.GetVtxT() << endl
00528     << "  VtxPln = " << shower.GetVtxPlane() << endl
00529     << "  DirCosU = " << shower.GetDirCosU() << endl
00530     << "  DirCosV = " << shower.GetDirCosV() << endl
00531     << "  DirCosZ = " << shower.GetDirCosZ() << endl
00532     << "  DirCosErrU = " << shower.GetDirCosErrU() << endl
00533     << "  DirCosErrV = " << shower.GetDirCosErrV() << endl
00534     << "  Energy = " << shower.GetEnergy() << endl
00535     << "  ReseedFlag = " << shower.GetReseedFlag() << endl;
00536 
00537 
00538   for(i=0;i<npln;i++){
00539     fHitArr[i].Clear();
00540   }
00541 
00542   delete [] Z;
00543   delete [] U;
00544   delete [] V; 
00545   delete [] CT;
00546   delete [] Q;
00547   delete [] plnvuw;
00548   delete [] plnnum;
00549 
00550   return;
00551 }

void AlgShowerAtNu::Trace const char *  c  )  const [virtual]
 

Reimplemented from AlgBase.

Definition at line 553 of file AlgShowerAtNu.cxx.

00554 {
00555 
00556 }


Member Data Documentation

TObjArray AlgShowerAtNu::fHitArr[500] [private]
 

Definition at line 18 of file AlgShowerAtNu.h.

Referenced by RunAlg().


The documentation for this class was generated from the following files:
Generated on Thu Nov 1 15:55:27 2007 for loon by  doxygen 1.3.9.1