AlgShowerCam Class Reference

#include <AlgShowerCam.h>

Inheritance diagram for AlgShowerCam:

AlgBase List of all members.

Public Member Functions

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

Private Member Functions

void ExtractHitProperties (const ShowerCamAtNu *shw, const int &bpln, CandShowerHandle &shower)
void FindTransversePositions ()
void SetupTimingInfo (UgliGeomHandle *ugh, const double &fFibreIndex)
void SetShowerCoordinates (const int &bpln, CandShowerAtNuHandle &shower)
void FindShowerVertex (ShowerCamAtNu *shwu, ShowerCamAtNu *shwv, double &trkscr, const int &bpln, const int &epln)
void DetermineDirection (ShowerCamAtNu *shwu, ShowerCamAtNu *shwv, double &offset, double &trkscr, double &shwscr, double &dirscr)
void CalculateShowerEnergy (CandShowerHandle &shower, const double &totph, double &energy)

Private Attributes

vector< PlaneInfo > fPlaneInfo
vector< HitCamAtNu * > fHitArr [500]
double vtxu
double vtxv
double vtxx
double vtxy
double vtxz
int vtxpln
bool vtxshw
double diru
double dirv
double dirz
double dirs
double direrru
double direrrv

Detailed Description

Definition at line 38 of file AlgShowerCam.h.

Constructor & Destructor Documentation

AlgShowerCam::AlgShowerCam (  ) 

Definition at line 39 of file AlgShowerCam.cxx.

00040 {
00041 }

AlgShowerCam::~AlgShowerCam (  ) 

Definition at line 43 of file AlgShowerCam.cxx.

00044 {
00045 }

Member Function Documentation

void AlgShowerCam::CalculateShowerEnergy ( CandShowerHandle &  shower,
const double &  totph,
double &  energy 
) [private]

Definition at line 691 of file AlgShowerCam.cxx.

References CandRecoHandle::CalibrateMIP(), CandRecoHandle::CalibrateSigMapped(), fHitArr, fPlaneInfo, Calibrator::GetAttenCorrectedTpos(), CandStripHandle::GetCharge(), PlexStripEndId::GetEncoded(), Calibrator::GetMIP(), CandStripHandle::GetNDigit(), CandStripHandle::GetPlaneView(), CandStripHandle::GetStripEndId(), Calibrator::Instance(), Msg::kDebug, StripEnd::kNegative, StripEnd::kPositive, CalDigitType::kSigCorr, PlaneView::kU, PlaneView::kV, PlaneView::kX, PlaneView::kY, MSG, and ValueErr< T >::SetError().

Referenced by RunAlg().

00692 {
00693   Calibrator& cal = Calibrator::Instance();
00694   MSG("AlgShowerCam", Msg::kDebug) << " *** calculate energy *** " << endl;
00695   // do the SIGMAP/SIGMIP calibration
00696   PlexStripEndId plexseid;
00697   FloatErr sigcorr,sigmip,sigmap;
00698   FloatErr lpos;
00699   unsigned int npln = fPlaneInfo.size();
00700   for(unsigned int pln=0;pln<npln;pln++)
00701     {
00702       for(unsigned int i=0;i<fHitArr[pln].size();i++)
00703         {
00704           HitCamAtNu* hit = fHitArr[pln][i];
00705           CandStripHandle* strip = hit->GetCandStripHandle();
00706           
00707           lpos = 0.0;
00708           if( strip->GetPlaneView()==PlaneView::kU || strip->GetPlaneView()==PlaneView::kX ) lpos = fPlaneInfo[pln].V;
00709           if( strip->GetPlaneView()==PlaneView::kV || strip->GetPlaneView()==PlaneView::kY ) lpos = fPlaneInfo[pln].U;
00710           lpos.SetError(1.0);
00711           
00712           if(strip->GetNDigit(StripEnd::kPositive)>0)
00713             {
00714               plexseid=strip->GetStripEndId(StripEnd::kPositive);
00715               sigcorr=strip->GetCharge(StripEnd::kPositive,CalDigitType::kSigCorr);
00716               sigmap=cal.GetAttenCorrectedTpos(sigcorr,lpos,plexseid);
00717               sigmip=cal.GetMIP(sigmap,plexseid);      
00718               shower.CalibrateSigMapped(plexseid.GetEncoded(),sigmap);
00719               shower.CalibrateMIP(plexseid.GetEncoded(),sigmip);
00720             }
00721           
00722           if(strip->GetNDigit(StripEnd::kNegative)>0)
00723             { 
00724               plexseid=strip->GetStripEndId(StripEnd::kNegative);
00725               sigcorr=strip->GetCharge(StripEnd::kNegative,CalDigitType::kSigCorr);
00726               sigmap=cal.GetAttenCorrectedTpos(sigcorr,lpos,plexseid);
00727               sigmip=cal.GetMIP(sigmap,plexseid);      
00728               shower.CalibrateSigMapped(plexseid.GetEncoded(),sigmap);
00729               shower.CalibrateMIP(plexseid.GetEncoded(),sigmip);
00730             }
00731         }
00732     }
00733   // calculate energy from pulse height
00734   if(totph>0.0)
00735     { 
00736       energy = totph; // NEEDS TUNING
00737       energy = 0.3 + totph/150.0;
00738     }
00739   return;
00740 }

void AlgShowerCam::DetermineDirection ( ShowerCamAtNu *  shwu,
ShowerCamAtNu *  shwv,
double &  offset,
double &  trkscr,
double &  shwscr,
double &  dirscr 
) [private]

Definition at line 572 of file AlgShowerCam.cxx.

References direrru, direrrv, dirs, diru, dirv, dirz, fHitArr, fPlaneInfo, ShowerCamAtNu::GetBegPlane(), ShowerCamAtNu::GetEndPlane(), Msg::kDebug, Msg::kVerbose, MSG, vtxshw, vtxu, vtxv, vtxz, and MinosMaterial::Z().

Referenced by RunAlg().

00573 {
00574   MSG("AlgShowerCam", Msg::kDebug) << " *** determine direction *** " << endl;
00575 
00576   // measure overall timeslope
00577   double ctimeslope=0.0,ctimeoffset=0.0,ctimeaverage=0.0,ctimescatter=0.0;
00578   double sw,swz,swt,swzz,swzt,swtt;
00579   double du,dv,dz;
00580   double z,t,w;
00581 
00582   int flag,npts;
00583 
00584   npts=0;
00585   sw=0.0; swz=0.0; swzz=0.0; swt=0.0; swtt=0.0; swzt=0.0;
00586   unsigned int npln = fPlaneInfo.size();
00587   for(unsigned int i=0;i<npln;i++)
00588     {
00589       if( fPlaneInfo[i].plnvuw>-1 ){
00590         z=fPlaneInfo[i].Z; t=fPlaneInfo[i].CTm; w=fPlaneInfo[i].Qm;
00591         swz+=w*z; swzz+=w*z*z; swt+=w*t; swtt+=w*t*t; swzt+=w*z*t;
00592         sw+=w; 
00593         t=fPlaneInfo[i].CTp; w=fPlaneInfo[i].Qp;
00594         swz+=w*z; swzz+=w*z*z; swt+=w*t; swtt+=w*t*t; swzt+=w*z*t;
00595         sw+=w; 
00596         npts++;
00597       }
00598     }
00599   if(npts>1)
00600     {
00601       ctimeslope=(sw*swzt-swz*swt)/(sw*swzz-swz*swz); 
00602       ctimeoffset=(swt*swzz-swz*swzt)/(sw*swzz-swz*swz);
00603       ctimeaverage=(swt/sw);
00604       if( (swtt/sw)-(swt/sw)*(swt/sw)>0.0 && (swzz/sw)-(swz/sw)*(swz/sw)>0.0 )
00605         {
00606           ctimescatter=((swzt/sw)-(swz/sw)*(swt/sw))/((sqrt((swzz/sw)-(swz/sw)*(swz/sw)))*(sqrt((swtt/sw)-(swt/sw)*(swt/sw))));
00607         }
00608     }
00609   offset=ctimeaverage;
00610   shwscr=ctimescatter;
00611 
00612   if(vtxshw) dirscr=trkscr; else dirscr=shwscr;
00613 
00614   npts=0;
00615   diru=0.0; direrru=0.0;
00616   if(1+shwu->GetEndPlane()-shwu->GetBegPlane()>1)
00617     {
00618       sw=0.0; swz=0.0; swt=0.0; 
00619       swzz=0.0; swzt=0.0; swtt=0.0;
00620       for(unsigned int i=0;i<npln;i++)
00621         {
00622           flag=0;
00623           for(unsigned int j=0;j<fHitArr[i].size();j++)
00624             {
00625               HitCamAtNu* hit = fHitArr[i][j];
00626               if(hit->GetShwFlag()>1 && hit->GetPlaneView()==0)
00627                 {
00628                   du=hit->GetTPos()-vtxu; dz=hit->GetZPos()-vtxz; w=hit->GetCharge();
00629                   swz+=w*dz; swt+=w*du; 
00630                   swzz+=w*dz*dz; swzt+=w*dz*du; swtt+=w*du*du;
00631                   sw+=w; flag=1;
00632                 }           
00633             }
00634           if(flag) npts++;
00635         }
00636       if(npts>1)
00637         {
00638           diru=(sw*swzt-swz*swt)/(sw*swzz-swz*swz);
00639           if( (swtt/sw)-(swt/sw)*(swt/sw)>0.0 && (swzz/sw)-(swz/sw)*(swz/sw)>0.0 )
00640             {
00641               direrru=((swzt/sw)-(swz/sw)*(swt/sw))/((sqrt((swzz/sw)-(swz/sw)*(swz/sw)))*(sqrt((swtt/sw)-(swt/sw)*(swt/sw))));
00642             }
00643         }
00644     }
00645   
00646   npts=0;
00647   dirv=0.0; direrrv=0.0;
00648   if(1+shwv->GetEndPlane()-shwv->GetBegPlane()>1)
00649     {
00650       sw=0.0; swz=0.0; swt=0.0; 
00651       swzz=0.0; swzt=0.0; swtt=0.0;
00652       for(unsigned int i=0;i<npln;i++)
00653         {
00654           flag=0;
00655           for(unsigned int j=0;j<fHitArr[i].size();j++)
00656             {
00657               HitCamAtNu* hit = fHitArr[i][j];
00658               if(hit->GetShwFlag()>1 && hit->GetPlaneView()==1)
00659                 {
00660                   dv=hit->GetTPos()-vtxv; dz=hit->GetZPos()-vtxz; w=hit->GetCharge();
00661                   swz+=w*dz; swt+=w*dv; 
00662                   swzz+=w*dz*dz; swzt+=w*dz*dv; swtt+=w*dv*dv;
00663                   sw+=w; flag=1;
00664                 }
00665             }
00666           if(flag) npts++;
00667         }
00668       if(npts>1)
00669         {
00670           dirv=(sw*swzt-swz*swt)/(sw*swzz-swz*swz);
00671           if( (swtt/sw)-(swt/sw)*(swt/sw)>0.0 && (swzz/sw)-(swz/sw)*(swz/sw)>0.0 )
00672             {
00673               direrrv=((swzt/sw)-(swz/sw)*(swt/sw))/((sqrt((swzz/sw)-(swz/sw)*(swz/sw)))*(sqrt((swtt/sw)-(swt/sw)*(swt/sw))));
00674             }
00675         }
00676     }
00677   
00678   dirs=sqrt(diru*diru+dirv*dirv+1.0);
00679   diru=diru/dirs; dirv=dirv/dirs; dirz=1.0/dirs;
00680   
00681   MSG("AlgShowerCam", Msg::kVerbose) 
00682     << "   ... direction: "
00683     << "  (" << diru << "," << dirv << "," << dirz << ")" << endl;
00684   return;
00685 }

void AlgShowerCam::ExtractHitProperties ( const ShowerCamAtNu *  shw,
const int &  bpln,
CandShowerHandle &  shower 
) [private]

Definition at line 251 of file AlgShowerCam.cxx.

References CandHandle::AddDaughterLink(), fHitArr, fPlaneInfo, HitCamAtNu::GetCandStripHandle(), ShowerCamAtNu::GetHit(), and CandStripHandle::GetPlane().

Referenced by RunAlg().

00252 {
00253   //Loop over all the HitCamAtNu objects in this ShowerCamAtNu object
00254   //Extract U/V/Z postion info as will as charge from each one
00255   //Add daughter links for the strips they represent to the CandShowerHandle.
00256   unsigned int nhits = shw-> GetNHits();
00257   int pln;
00258   for(unsigned int i=0;i<nhits;++i)
00259     {
00260       HitCamAtNu* hit = shw->GetHit(i);
00261       CandStripHandle* strip = hit->GetCandStripHandle();
00262       pln = hit->GetPlane()-bpln;
00263       fHitArr[pln].push_back(hit);
00264       if(fPlaneInfo[pln].plnvuw<0)
00265         {
00266           fPlaneInfo[pln].plnvuw = hit->GetPlaneView();
00267           fPlaneInfo[pln].Z = hit->GetZPos();
00268         }
00269       if(fPlaneInfo[pln].plnvuw==0) fPlaneInfo[pln].U += hit->GetCharge()*hit->GetTPos();
00270       if(fPlaneInfo[pln].plnvuw==1) fPlaneInfo[pln].V += hit->GetCharge()*hit->GetTPos();
00271       fPlaneInfo[pln].Q += hit->GetCharge();
00272       shower.AddDaughterLink(*strip);
00273       //shower.SetInShower(strip,hit->GetShwFlag());
00274       (fPlaneInfo[pln].plnnum)++; 
00275     }
00276 }

void AlgShowerCam::FindShowerVertex ( ShowerCamAtNu *  shwu,
ShowerCamAtNu *  shwv,
double &  trkscr,
const int &  bpln,
const int &  epln 
) [private]

Definition at line 476 of file AlgShowerCam.cxx.

References fHitArr, fPlaneInfo, ShowerCamAtNu::GetAssocTrackAt(), TrackCamAtNu::GetBegPlane(), ShowerCamAtNu::GetBegVtxFlag(), TrackCamAtNu::GetEndPlane(), ShowerCamAtNu::GetEndVtxFlag(), VHS::GetPlane(), Msg::kDebug, Munits::km, Msg::kVerbose, MSG, vtxpln, vtxshw, vtxu, vtxv, vtxx, vtxy, and vtxz.

Referenced by RunAlg().

00477 {
00478 
00479   MSG("AlgShowerCam", Msg::kDebug) << " *** determine vertex  *** " << endl;
00480   int trkbpln=-1,trkepln=-1,shwbpln=-1,shwepln=-1;
00481   unsigned int npln =  fPlaneInfo.size();
00482   unsigned int NHitsInPlane;
00483   // determine vertex plane
00484   double totchgpln=0.0,totchg=0.0;
00485   int pln,k,km,kp;
00486   for(unsigned int i=0;i<npln;i++)
00487     {
00488       NHitsInPlane= fHitArr[i].size();
00489       for(unsigned int j=0;j<NHitsInPlane;j++)
00490         {
00491           HitCamAtNu* hit = fHitArr[i][j];
00492           if(hit->GetShwFlag()>1)
00493             {
00494               totchg+=hit->GetCharge();
00495               totchgpln+=hit->GetCharge()*hit->GetPlane(); 
00496               if(shwbpln<0||hit->GetPlane()<shwbpln) shwbpln=hit->GetPlane();
00497               if(shwepln<0||hit->GetPlane()>shwepln) shwepln=hit->GetPlane(); 
00498             } 
00499         }
00500     }
00501   if(totchg>0.0)
00502     {
00503       vtxpln=(int)(totchgpln/totchg);
00504     }
00505   else
00506     {
00507       vtxpln=(int)(0.5*(bpln+epln));
00508       shwbpln=bpln; shwepln=epln;
00509     }
00510 
00511   const TrackCamAtNu* vtxtrkU = shwu->GetAssocTrackAt(0);
00512   const TrackCamAtNu* vtxtrkV = shwv->GetAssocTrackAt(0);
00513   if(vtxtrkU && vtxtrkV)
00514     {
00515       vtxshw=1;
00516       if(vtxtrkU->GetBegPlane()<vtxtrkV->GetBegPlane()) trkbpln=vtxtrkU->GetBegPlane(); else trkbpln=vtxtrkV->GetBegPlane();
00517       if(vtxtrkU->GetEndPlane()>vtxtrkV->GetEndPlane()) trkepln=vtxtrkU->GetEndPlane(); else trkepln=vtxtrkV->GetEndPlane();
00518       trkscr = (trkbpln+trkepln-2.0*vtxpln)/(1.0+trkepln-trkbpln);
00519       if(trkscr>0.0){ trkscr = (shwbpln+shwepln-2.0*trkbpln)/(1.0+shwepln-shwbpln); }
00520       if(trkscr<0.0){ trkscr = (shwbpln+shwepln-2.0*trkepln)/(1.0+shwepln-shwbpln); }
00521       if( ( shwu->GetBegVtxFlag() && shwv->GetBegVtxFlag())
00522           && !( shwu->GetEndVtxFlag() && shwv->GetEndVtxFlag()) ) vtxpln=trkbpln;
00523       if( ( shwu->GetEndVtxFlag() && shwv->GetEndVtxFlag())
00524           && !( shwu->GetBegVtxFlag() && shwv->GetBegVtxFlag()) ) vtxpln=trkepln;
00525     }
00526   
00527   // determine vertex position
00528   pln=vtxpln-bpln;
00529   if( pln>-1 && pln<(int)npln && fPlaneInfo[pln].plnvuw>-1 )
00530     {
00531       vtxu=fPlaneInfo[pln].U; vtxv=fPlaneInfo[pln].V; vtxz=fPlaneInfo[pln].Z;
00532     }
00533   else
00534     {
00535       k=0; km=-1;
00536       while(pln-k>0 && km<0)
00537         {
00538           k++; if(pln-k<(int)npln && fPlaneInfo[pln-k].plnvuw>-1) km=k;
00539         }
00540       k=0; kp=-1;
00541       while(pln+k<(int)npln-1 && kp<0)
00542         {
00543           k++; if(pln+k>-1 && fPlaneInfo[pln+k].plnvuw>-1) kp=k; 
00544     }
00545       if(km>0 && kp>0)
00546         {
00547           vtxu=(km*fPlaneInfo[pln+kp].U+kp*fPlaneInfo[pln-km].U)/(km+kp); 
00548           vtxv=(km*fPlaneInfo[pln+kp].V+kp*fPlaneInfo[pln-km].V)/(km+kp); 
00549           vtxz=(km*fPlaneInfo[pln+kp].Z+kp*fPlaneInfo[pln-km].Z)/(km+kp); 
00550         }
00551       else
00552         {
00553           if(km>0){ vtxu=fPlaneInfo[pln-km].U; vtxv=fPlaneInfo[pln-km].V; vtxz=fPlaneInfo[pln-km].Z; } 
00554           if(kp>0){ vtxu=fPlaneInfo[pln+kp].U; vtxv=fPlaneInfo[pln+kp].V; vtxz=fPlaneInfo[pln+kp].Z; }
00555         }
00556     }
00557   
00558   vtxx=0.7071*(vtxu-vtxv); vtxy=0.7071*(vtxu+vtxv); 
00559   
00560   MSG("AlgShowerCam", Msg::kVerbose) 
00561     << "   ... vertex: "
00562     << " plane = " << vtxpln
00563     << " position = (" << vtxu << "," << vtxv << "," << vtxz << ")" << endl;
00564   
00565   return;
00566 }

void AlgShowerCam::FindTransversePositions (  )  [private]

Definition at line 281 of file AlgShowerCam.cxx.

References fPlaneInfo, Munits::km2, and MinosMaterial::Z().

Referenced by RunAlg().

00282 {
00283   const int npln = (int)fPlaneInfo.size();
00284   // interpolate between views for U,V
00285   int flag,vuw;
00286   int pln;
00287   int k,km1,kp1,km2,kp2;
00288   for(pln=0;pln<npln;pln++)
00289     {
00290       if(fPlaneInfo[pln].plnvuw>-1)
00291         {
00292           flag=0;
00293           vuw=fPlaneInfo[pln].plnvuw;
00294           k=0; km1=-1; km2=-1;
00295           while(pln-k>0 && km2<0)
00296             {
00297               k++; 
00298               if(fPlaneInfo[pln-k].plnvuw>-1 && fPlaneInfo[pln-k].plnvuw!=vuw)
00299                 {
00300                   if(km1<0) km1=k; else km2=k;
00301                 }
00302             }
00303           k=0; kp1=-1; kp2=-1;
00304           while(pln+k<npln-1 && kp2<0)
00305             {
00306               k++;
00307               if(fPlaneInfo[pln+k].plnvuw>-1 && fPlaneInfo[pln+k].plnvuw!=vuw)
00308                 {
00309                   if(kp1<0) kp1=k; else kp2=k;
00310                 }
00311             }
00312           if(km1>0 && kp1>0)
00313             { 
00314               if(vuw==0) fPlaneInfo[pln].V=(double(kp1)*fPlaneInfo[pln-km1].V+double(km1)*fPlaneInfo[pln+kp1].V)/double(km1+kp1);
00315               if(vuw==1) fPlaneInfo[pln].U=(double(kp1)*fPlaneInfo[pln-km1].U+double(km1)*fPlaneInfo[pln+kp1].U)/double(km1+kp1);
00316               flag=1;
00317             }
00318           else
00319             {
00320               if(km1>0)
00321                 {
00322                   if(vuw==0) 
00323                     {
00324                       if(km2>0) fPlaneInfo[pln].V=fPlaneInfo[pln-km1].V+(fPlaneInfo[pln].Z-fPlaneInfo[pln-km1].Z)*(fPlaneInfo[pln-km1].V-fPlaneInfo[pln-km2].V)/(fPlaneInfo[pln-km1].Z-fPlaneInfo[pln-km2].Z); 
00325                       else fPlaneInfo[pln].V=fPlaneInfo[pln-km1].V;
00326                     }
00327                   if(vuw==1) 
00328                     {
00329                       if(km2>0) fPlaneInfo[pln].U=fPlaneInfo[pln-km1].U+(fPlaneInfo[pln].Z-fPlaneInfo[pln-km1].Z)*(fPlaneInfo[pln-km1].U-fPlaneInfo[pln-km2].U)/(fPlaneInfo[pln-km1].Z-fPlaneInfo[pln-km2].Z); 
00330                       else fPlaneInfo[pln].U=fPlaneInfo[pln-km1].U;
00331                     }
00332                   flag=1;
00333                 }
00334               if(kp1>0)
00335                 {
00336                   if(vuw==0) 
00337                     {
00338                       if(kp2>0) fPlaneInfo[pln].V=fPlaneInfo[pln+kp1].V+(fPlaneInfo[pln].Z-fPlaneInfo[pln+kp1].Z)*(fPlaneInfo[pln+kp1].V-fPlaneInfo[pln+kp2].V)/(fPlaneInfo[pln+kp1].Z-fPlaneInfo[pln+kp2].Z); 
00339                       else fPlaneInfo[pln].V=fPlaneInfo[pln+kp1].V;
00340                     }
00341                   if(vuw==1) 
00342                     {
00343                       if(kp2>0) fPlaneInfo[pln].U=fPlaneInfo[pln+kp1].U+(fPlaneInfo[pln].Z-fPlaneInfo[pln+kp1].Z)*(fPlaneInfo[pln+kp1].U-fPlaneInfo[pln+kp2].U)/(fPlaneInfo[pln+kp1].Z-fPlaneInfo[pln+kp2].Z); 
00344                       else fPlaneInfo[pln].U=fPlaneInfo[pln+kp1].U;
00345                     }
00346                   flag=1;
00347                 }
00348             }
00349           if(flag==0){ fPlaneInfo[pln].plnvuw=-1; }
00350         }
00351     }
00352   return;
00353 //   // estimate transverse co-ordinates
00354 //   int km,kp;
00355 //   double q,sq,sqt;
00356 //   for(pln=0;pln<npln;pln++){
00357 //     if(plnvuw[pln]>-1){
00358 //       vuw=plnvuw[pln];
00359 //       sq=0.0; sqt=0.0; 
00360 //       k=0; km=-1;
00361 //       while(pln-k>0 && km<0){
00362 //         k++; if(plnvuw[pln-k]>-1 && plnvuw[pln-k]!=vuw) km=k;
00363 //       }
00364 //       k=0; kp=-1;
00365 //       while(pln+k<npln-1 && kp<0){
00366 //         k++; if(plnvuw[pln+k]>-1 && plnvuw[pln+k]!=vuw) kp=k; 
00367 //       }
00368 
00369 //       if(km>0){
00370 //         if(vuw==0){ q=1.0/km; sq+=q*Q[pln-km]; sqt+=q*Q[pln-km]*V[pln-km]; }
00371 //         if(vuw==1){ q=1.0/km; sq+=q*Q[pln-km]; sqt+=q*Q[pln-km]*U[pln-km]; }
00372 //       }
00373 //       if(kp>0){
00374 //         if(vuw==0){ q=1.0/kp; sq+=q*Q[pln+kp]; sqt+=q*Q[pln+kp]*V[pln+kp]; }
00375 //         if(vuw==1){ q=1.0/kp; sq+=q*Q[pln+kp]; sqt+=q*Q[pln+kp]*U[pln+kp]; }
00376 //       }
00377 
00378 //       if(sq>0.0){
00379 //         if(vuw==0){ V[pln]=sqt/sq; } if(vuw==1){ U[pln]=sqt/sq; }
00380 //       }
00381 //       else{
00382 //         if(vuw==0){ V[pln]=0.0; } if(vuw==1){ U[pln]=0.0; }
00383 //       }
00384 //     }
00385 //   }
00386 }

void AlgShowerCam::RunAlg ( AlgConfig &  ac,
CandHandle &  ch,
CandContext &  cx 
) [virtual]

Implements AlgBase.

Definition at line 47 of file AlgShowerCam.cxx.

References CalculateShowerEnergy(), BeamSummary::clear(), PlaneInfo::CTm, PlaneInfo::CTp, DetermineDirection(), direrru, direrrv, dirs, diru, dirv, dirz, ExtractHitProperties(), fHitArr, FindShowerVertex(), FindTransversePositions(), fPlaneInfo, ShowerCamAtNu::GetBegPlane(), CandContext::GetCandRecord(), CandContext::GetDataIn(), VldContext::GetDetector(), CandRecoHandle::GetDirCosU(), CandRecoHandle::GetDirCosV(), CandRecoHandle::GetDirCosZ(), Registry::GetDouble(), ShowerCamAtNu::GetEndPlane(), CandShowerHandle::GetEnergy(), Registry::GetInt(), CandHandle::GetNDaughters(), UgliGeomHandle::GetSteelPlnHandle(), CandRecoHandle::GetTimeOffset(), CandRecoHandle::GetTimeSlope(), RecMinos::GetVldContext(), CandRecoHandle::GetVtxPlane(), CandRecoHandle::GetVtxT(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), CandRecoHandle::GetVtxZ(), UgliSteelPlnHandle::GetZ0(), Calibrator::Instance(), UgliSteelPlnHandle::IsValid(), Detector::kCalib, Msg::kDebug, Registry::KeyExists(), Detector::kFar, Detector::kNear, MSG, PlaneInfo::plnnum, PlaneInfo::plnvuw, PlaneInfo::Q, PlaneInfo::Qm, PlaneInfo::Qp, Calibrator::ReInitialise(), CandRecoHandle::SetCandSlice(), CandRecoHandle::SetDirCosU(), CandRecoHandle::SetDirCosV(), CandRecoHandle::SetDirCosZ(), CandShowerHandle::SetEnergy(), SetShowerCoordinates(), CandRecoHandle::SetTimeOffset(), CandRecoHandle::SetTimeSlope(), SetupTimingInfo(), CandRecoHandle::SetVtxPlane(), CandRecoHandle::SetVtxT(), CandRecoHandle::SetVtxU(), CandRecoHandle::SetVtxV(), CandRecoHandle::SetVtxZ(), PlaneInfo::U, PlaneInfo::V, vtxpln, vtxshw, vtxu, vtxv, vtxx, vtxy, vtxz, PlaneInfo::Wm, PlaneInfo::Wp, and PlaneInfo::Z.

00048 {
00049 
00050   MSG("AlgShowerCam", Msg::kDebug) << "AlgShowerCam::RunAlg(...)" << endl;
00051 
00052   // Get CandShowerAtNuHandle
00053   CandShowerAtNuHandle& shower = dynamic_cast<CandShowerAtNuHandle&>(ch);
00054 
00056   // Unpack AlgConfig
00057   double fFibreIndex=1.77;
00058   int fCamAnaOnOff=0;
00059   if( ac.KeyExists("FibreIndex") )  fFibreIndex  = ac.GetDouble("FibreIndex");
00060   if( ac.KeyExists("CamAnaOnOff") ) fCamAnaOnOff = ac.GetInt("CamAnaOnOff");
00061   MSG("AlgShowerCam", Msg::kDebug) << " FibreIndex=" << fFibreIndex << endl;
00062   MSG("AlgShowerCam", Msg::kDebug) << " CamAnaOnOff=" << fCamAnaOnOff << endl;
00063 
00065   // Unpack CandContext
00066   const TObjArray *arr = dynamic_cast<const TObjArray*>(cx.GetDataIn());
00067   ShowerCamAtNu* shwu = (ShowerCamAtNu*)(arr->At(0));
00068   ShowerCamAtNu* shwv = (ShowerCamAtNu*)(arr->At(1)); 
00069   shower.SetCandSlice(shwu->GetCandSliceHandle());
00070 
00071   CandRecord* candrec = (CandRecord*)(cx.GetCandRecord());
00072   VldContext *vldc = (VldContext*)(candrec->GetVldContext());
00073   UgliGeomHandle ugh(*vldc);
00074 
00076   // Reset Calibrator
00077   Calibrator& cal = Calibrator::Instance();
00078   cal.ReInitialise(*vldc);
00079 
00081   // Get Geometry
00082   int atmosflag,modtype,modnum,modpln,modstr;
00083   switch(candrec->GetVldContext()->GetDetector()){
00084     case Detector::kFar:
00085       atmosflag=1; modtype=1; modnum=2; modpln=248; modstr=192; break;
00086     case Detector::kNear:
00087       atmosflag=0; modtype=2; modnum=1; modpln=282; modstr=96; break;
00088     case Detector::kCalib:
00089       atmosflag=0; modtype=3; modnum=1; modpln=60; modstr=24; break;
00090     default:
00091       atmosflag=0; modtype=-1; modnum=0; modpln=0; modstr=0; break;
00092   }
00093 
00094   int i,j;//,k;
00095   //int vuw;
00096 
00097   double begZsm1=-9999.9,endZsm1=-9999.9;
00098   double begZsm2=-9999.9,endZsm2=-9999.9;
00099   int begsm1=-1,endsm1=-1,begsm2=-1,endsm2=-1;
00100 
00101   for(j=0;j<486;j++){
00102     PlexPlaneId myplaneid(vldc->GetDetector(),j,1);
00103     UgliSteelPlnHandle myplanehandle = ugh.GetSteelPlnHandle(myplaneid);
00104     if(myplanehandle.IsValid()){
00105       if(j<=modpln && begZsm1<0.0){ begsm1=j; begZsm1=myplanehandle.GetZ0(); }
00106       if(j<=modpln && myplanehandle.GetZ0()>endZsm1){ endsm1=j; endZsm1=myplanehandle.GetZ0(); }
00107       if(j>modpln && begZsm2<0.0){ begsm2=j; begZsm2=myplanehandle.GetZ0(); }
00108       if(j>modpln && myplanehandle.GetZ0()>endZsm2){ endsm2=j; endZsm2=myplanehandle.GetZ0(); }
00109     }
00110   }
00111 
00112   MSG("AlgShowerCam", Msg::kDebug) << " *** detector geometry *** " << endl
00113                                    << "  SM1 :  pln " << begsm1 << " -> " << endsm1
00114                                    << "    Z " << begZsm1 << " -> " << endZsm1 << endl
00115                                    << "  SM2 :  pln " << begsm2 << " -> " << endsm2 
00116                                    << "    Z " << begZsm2 << " -> " << endZsm2 << endl;
00117 
00119   // Find the extent of the shower in planes
00120   int bpln,epln;
00121   if(shwu->GetBegPlane()<shwv->GetBegPlane()) bpln=shwu->GetBegPlane(); else bpln=shwv->GetBegPlane();
00122   if(shwu->GetEndPlane()>shwv->GetEndPlane()) epln=shwu->GetEndPlane(); else epln=shwv->GetEndPlane();
00123   const int npln = 1+epln-bpln;  
00124 
00125 
00127   // Fill the vector with an empty TrkPlaneInfo object for each plane in the range spanned by the track
00128   PlaneInfo tmpplane;
00129   tmpplane.plnvuw=-1; tmpplane.plnnum=0;
00130   tmpplane.U=0.0; tmpplane.V=0.0; tmpplane.Z=0.0; tmpplane.Q=0.0; 
00131   tmpplane.Qm=0.0; tmpplane.Qp=0.0; tmpplane.CTm=0.0; tmpplane.CTp=0.0; 
00132   tmpplane.Wm=0.0; tmpplane.Wp=0.0;
00133   for(i=0;i<npln;i++) { fPlaneInfo.push_back(tmpplane); }
00134 
00136   // Calculate Track Co-ordinates
00137 
00138   MSG("AlgShowerCam", Msg::kDebug) << " *** calculate co-ordinates *** " << endl;
00139   ExtractHitProperties(shwu, bpln, shower);
00140   ExtractHitProperties(shwv, bpln, shower);
00141   int nplanes=0,nstrips=0;
00142   double totph=0.0;
00143   for(i=0;i<npln;i++)
00144     {
00145       if(fPlaneInfo[i].plnvuw>-1 && fPlaneInfo[i].Q>0.0)
00146         {
00147           totph+=fPlaneInfo[i].Q;
00148           fPlaneInfo[i].U/=fPlaneInfo[i].Q; fPlaneInfo[i].V/=fPlaneInfo[i].Q;
00149           nstrips+=fPlaneInfo[i].plnnum; ++nplanes;
00150         }
00151     }
00152 
00153   // Interpolate to get V for U planes and vice versa
00154   FindTransversePositions();
00155   // Calculate The Timing Info For Each Shower StripEnd
00156   SetupTimingInfo(&ugh, fFibreIndex);
00157   // Set U, V for each plane/stripend in the shower
00158   SetShowerCoordinates(bpln, shower);
00159 
00160 //   MSG("AlgShowerCam", Msg::kDebug) << " *** shower hits *** " << endl;
00161 //   for(pln=0;pln<npln;pln++){
00162 //     if( plnvuw[pln]>-1 ){
00163 //       MSG("AlgShowerCam",Msg::kDebug)
00164 //         << " pln=" << pln << " vuw=" << plnvuw[pln]
00165 //         << " z=" << Z[pln] << " u=" << U[pln] << " v=" << V[pln] << endl;
00166 //     }
00167 //   }
00168 
00169 
00170   // ***********************************
00171   // * D E T E R M I N E   V E R T E X * 
00172   // ***********************************
00173   vtxu=0.0,vtxv=0.0,vtxx=0.0,vtxy=0.0,vtxz=0.0; 
00174   vtxpln=-1; vtxshw=false;
00175   double trkscore=0.0,showerscore=0.0,dirscore=0.0;
00176   FindShowerVertex(shwu, shwv, trkscore, bpln, epln);
00177 
00178   //*****************************************
00179   //* D E T E R M I N E   D I R E C T I O N *
00180   //*****************************************
00181 
00182   diru=0.0; dirv=0.0; dirz=0.0; dirs=0.0; 
00183   direrru=0.0; direrrv=0.0; 
00184   double offset=0.0;
00185   DetermineDirection( shwu, shwv, offset, trkscore, showerscore, dirscore);
00186 
00187   // ***********************************
00188   // * C A L C U L A T E   E N E R G Y *
00189   // ***********************************
00190   double energy(0.0);
00191   CalculateShowerEnergy( shower, totph, energy );
00192 
00193   // determine overall direction 
00194   double dir=1.0;
00195   //if(dirscore>0.0) dir=1.0; 
00196   if(dirscore<0.0) dir=-1.0;  
00197 
00198   // *******************************************
00199   // * S E T   S H O W E R   V A R I A B L E S *
00200   // *******************************************
00201 
00202   // set shower variables
00203   shower.SetTimeSlope((dir)/(3.0e8)); 
00204   shower.SetTimeOffset((offset)/(3.0e8)); 
00205   shower.SetVtxU(vtxu);
00206   shower.SetVtxV(vtxv);
00207   shower.SetVtxZ(vtxz);
00208   //shower.SetVtxR(vtxr);
00209   shower.SetVtxT((offset)/(3.0e8));
00210   shower.SetVtxPlane(vtxpln);
00211   shower.SetDirCosU(dir*diru);
00212   shower.SetDirCosV(dir*dirv);
00213   shower.SetDirCosZ(dir*dirz);
00214   shower.SetEnergy(energy);
00215 
00216   MSG("AlgShowerCam", Msg::kDebug) 
00217     << " *** CREATING NEW SHOWER *** " << endl
00218     << "  NStrips = " << shower.GetNDaughters() << endl
00219     << "  TimeSlope = " << shower.GetTimeSlope() << endl
00220     << "  TimeOffset = " << shower.GetTimeOffset() << endl
00221     << "  VtxU = " << shower.GetVtxU() << endl
00222     << "  VtxV = " << shower.GetVtxV() << endl
00223     << "  VtxZ = " << shower.GetVtxZ() << endl
00224     << "  VtxT = " << shower.GetVtxT() << endl
00225     << "  VtxPln = " << shower.GetVtxPlane() << endl
00226     << "  DirCosU = " << shower.GetDirCosU() << endl
00227     << "  DirCosV = " << shower.GetDirCosV() << endl
00228     << "  DirCosZ = " << shower.GetDirCosZ() << endl
00229     << "  Energy = " << shower.GetEnergy() << endl;
00230 
00231   for(i=0;i<npln;i++)
00232     {
00233       fHitArr[i].clear();
00234     }
00235   fPlaneInfo.clear();
00236   
00237   return;
00238 }

void AlgShowerCam::SetShowerCoordinates ( const int &  bpln,
CandShowerAtNuHandle &  shower 
) [private]

Definition at line 431 of file AlgShowerCam.cxx.

References fPlaneInfo, CandShowerHandle::SetU(), and CandShowerHandle::SetV().

Referenced by RunAlg().

00432 {//int& plnm, int& plnp, 
00433   int pln;
00434   const int npln = (int)fPlaneInfo.size();
00435   for(int i=0;i<npln;i++)
00436     {
00437       pln=bpln+i;
00438       if(fPlaneInfo[i].plnvuw>-1)
00439         {
00440           shower.SetU(pln,fPlaneInfo[i].U); shower.SetV(pln,fPlaneInfo[i].V);
00441           if(fPlaneInfo[i].Qm>0.0)
00442             {
00443               fPlaneInfo[i].CTm/=fPlaneInfo[i].Qm;
00444             }
00445           if(fPlaneInfo[i].Qp>0.0)
00446             {
00447               fPlaneInfo[i].CTp/=fPlaneInfo[i].Qp;
00448             }
00449           //if(plnm<0||i<plnm) plnm=i; if(plnp<0||i>plnp) plnp=i;
00450         }
00451     }
00452   return;
00453 
00454 //   // set map values
00455 //   for(i=0;i<npln;i++){
00456 //     pln=bpln+i;
00457 //     if(plnvuw[i]>-1){
00458 //       shower.SetU(pln,U[i]);
00459 //       shower.SetV(pln,V[i]);    
00460 //       if(Qm[i]>0.0){
00461 //         CTm[i]=CTm[i]/Qm[i];
00462 //       }
00463 //       if(Qp[i]>0.0){
00464 //         CTp[i]=CTp[i]/Qp[i];
00465 //       }
00466 //       if(Qm[i]+Qp[i]>0){
00467 //         CT[i]=(Qm[i]*CTm[i]+Qp[i]*CTp[i])/(Qm[i]+Qp[i]);
00468 //       }
00469 //     }
00470 //   }
00471 }

void AlgShowerCam::SetupTimingInfo ( UgliGeomHandle *  ugh,
const double &  fFibreIndex 
) [private]

Definition at line 391 of file AlgShowerCam.cxx.

References UgliStripHandle::ClearFiber(), fHitArr, fPlaneInfo, Munits::g, HitCamAtNu::GetCandStripHandle(), CandStripHandle::GetCharge(), UgliStripHandle::GetHalfLength(), CandStripHandle::GetNDigit(), CandStripHandle::GetPlaneView(), CandStripHandle::GetStripEndId(), UgliGeomHandle::GetStripHandle(), CandStripHandle::GetTime(), StripEnd::kNegative, StripEnd::kPositive, PlaneView::kU, PlaneView::kV, PlaneView::kX, PlaneView::kY, and UgliStripHandle::WlsPigtail().

Referenced by RunAlg().

00392 {
00393   double tpos = 0., fibre = 0., ctime =0., digchg = 0.;
00394   HitCamAtNu* hit = 0; CandStripHandle* strip = 0;
00395   const unsigned int npln = fPlaneInfo.size();
00396   for(unsigned int pln=0;pln<npln;pln++)
00397     {
00398       for(unsigned int g=0;g<fHitArr[pln].size();++g)
00399         {
00400           hit = fHitArr[pln][g];
00401           strip = hit->GetCandStripHandle();  
00402           if( strip->GetPlaneView()==PlaneView::kU
00403               || strip->GetPlaneView()==PlaneView::kX ) tpos = -fPlaneInfo[pln].V;
00404           if( strip->GetPlaneView()==PlaneView::kV
00405               || strip->GetPlaneView()==PlaneView::kY ) tpos = fPlaneInfo[pln].U;
00406           PlexStripEndId stripid = strip->GetStripEndId();
00407           UgliStripHandle striphandle = ugh->GetStripHandle(stripid);
00408 
00409           if(strip->GetNDigit(StripEnd::kPositive)>0)
00410             { 
00411               fibre = striphandle.ClearFiber(StripEnd::kPositive)+striphandle.WlsPigtail(StripEnd::kPositive)+striphandle.GetHalfLength()-tpos; 
00412               ctime = 3.0e8*strip->GetTime(StripEnd::kPositive)-fFibreIndex*fibre; 
00413               digchg=strip->GetCharge(StripEnd::kPositive); 
00414               fPlaneInfo[pln].Qp+=digchg; fPlaneInfo[pln].CTp+=digchg*ctime;
00415             }
00416         
00417           if(strip->GetNDigit(StripEnd::kNegative)>0)
00418             { 
00419               fibre = striphandle.ClearFiber(StripEnd::kNegative)+striphandle.WlsPigtail(StripEnd::kNegative)+striphandle.GetHalfLength()+tpos; 
00420               ctime = 3.0e8*strip->GetTime(StripEnd::kNegative)-fFibreIndex*fibre; 
00421               digchg=strip->GetCharge(StripEnd::kNegative); 
00422               fPlaneInfo[pln].Qm+=digchg; fPlaneInfo[pln].CTm+=digchg*ctime;
00423             }   
00424         }
00425     }
00426 }

void AlgShowerCam::Trace ( const char *  c  )  const [virtual]

Reimplemented from AlgBase.

Definition at line 243 of file AlgShowerCam.cxx.

00244 {
00245 
00246 }

Member Data Documentation

double AlgShowerCam::direrru [private]

Definition at line 61 of file AlgShowerCam.h.

Referenced by DetermineDirection(), and RunAlg().

double AlgShowerCam::direrrv [private]

Definition at line 62 of file AlgShowerCam.h.

Referenced by DetermineDirection(), and RunAlg().

double AlgShowerCam::dirs [private]

Definition at line 60 of file AlgShowerCam.h.

Referenced by DetermineDirection(), and RunAlg().

double AlgShowerCam::diru [private]

Definition at line 57 of file AlgShowerCam.h.

Referenced by DetermineDirection(), and RunAlg().

double AlgShowerCam::dirv [private]

Definition at line 58 of file AlgShowerCam.h.

Referenced by DetermineDirection(), and RunAlg().

double AlgShowerCam::dirz [private]

Definition at line 59 of file AlgShowerCam.h.

Referenced by DetermineDirection(), and RunAlg().

vector<HitCamAtNu*> AlgShowerCam::fHitArr[500] [private]

Definition at line 49 of file AlgShowerCam.h.

Referenced by CalculateShowerEnergy(), DetermineDirection(), ExtractHitProperties(), FindShowerVertex(), RunAlg(), and SetupTimingInfo().

vector<PlaneInfo> AlgShowerCam::fPlaneInfo [private]

Definition at line 48 of file AlgShowerCam.h.

Referenced by CalculateShowerEnergy(), DetermineDirection(), ExtractHitProperties(), FindShowerVertex(), FindTransversePositions(), RunAlg(), SetShowerCoordinates(), and SetupTimingInfo().

int AlgShowerCam::vtxpln [private]

Definition at line 55 of file AlgShowerCam.h.

Referenced by FindShowerVertex(), and RunAlg().

bool AlgShowerCam::vtxshw [private]

Definition at line 56 of file AlgShowerCam.h.

Referenced by DetermineDirection(), FindShowerVertex(), and RunAlg().

double AlgShowerCam::vtxu [private]

Definition at line 50 of file AlgShowerCam.h.

Referenced by DetermineDirection(), FindShowerVertex(), and RunAlg().

double AlgShowerCam::vtxv [private]

Definition at line 51 of file AlgShowerCam.h.

Referenced by DetermineDirection(), FindShowerVertex(), and RunAlg().

double AlgShowerCam::vtxx [private]

Definition at line 52 of file AlgShowerCam.h.

Referenced by FindShowerVertex(), and RunAlg().

double AlgShowerCam::vtxy [private]

Definition at line 53 of file AlgShowerCam.h.

Referenced by FindShowerVertex(), and RunAlg().

double AlgShowerCam::vtxz [private]

Definition at line 54 of file AlgShowerCam.h.

Referenced by DetermineDirection(), FindShowerVertex(), and RunAlg().

The documentation for this class was generated from the following files:
Generated on Tue May 21 22:00:48 2013 for loon by  doxygen 1.4.7