loon: ClusterCam.cxx Source File

00001 
00002 // Package: CandTrackCam
00003 //
00004 // CandTrackCam - ClusterCam.cxx
00005 //
00006 // marshall@hep.phy.cam.ac.uk
00008 #include <cmath>
00009 #include <cstdlib>
00010 
00011 #include "ClusterCam.h"
00012 #include "HitCam.h"
00013 
00014 ClassImp(ClusterCam)
00015 
00016   //CVSID("$Id: ClusterCam.cxx,v 1.5 2009/02/28 21:46:11 gmieg Exp $");
00017 
00018 
00019 ClusterCam::ClusterCam(HitCam* hit) :
00020   fPlane(-1),fBegStrip(-1), fEndStrip(-1), 
00021   fBegTime(0.), fEndTime(0.),
00022   fBegTPos(0.), fEndTPos(0.),
00023   fZPos(0.), fTPos(0.), fCharge(0.),
00024   fTrkFlag(0),fShwFlag(0),
00025   fTrkPlnFlag(0), fShwPlnFlag(0),
00026   fPlaneView(-1),fDigits(0),fNDFlag(1), 
00027   StripWidth(4.108e-2)
00028 {
00029   this->AddHit(hit);
00030 }
00032 
00033 
00035 ClusterCam::~ClusterCam()
00036 {
00037   HitsInCluster.clear();
00038 }
00040 
00041 
00043 void ClusterCam::AddHit(HitCam* hit)
00044 {
00045   if(HitsInCluster.size()==0) {
00046     HitsInCluster.push_back(hit);
00047 
00048     fPlane=hit->GetPlane();
00049     fBegStrip=hit->GetStrip();
00050     fEndStrip=hit->GetStrip();
00051     fBegTPos=hit->GetTPos();
00052     fEndTPos=hit->GetTPos();
00053     fBegTime=hit->GetTime();
00054     fEndTime=hit->GetTime();
00055     fZPos=hit->GetZPos();
00056     fPlaneView=hit->GetPlaneView();
00057   }
00058 
00059   else {
00060     if(this->ContainsHit(hit)==true) {return;}
00061     HitsInCluster.push_back(hit);   
00062 
00063     if(hit->GetStrip()<fBegStrip) fBegStrip=hit->GetStrip();
00064     if(hit->GetStrip()>fEndStrip) fEndStrip=hit->GetStrip();
00065     if(hit->GetTPos()<fBegTPos) fBegTPos=hit->GetTPos();
00066     if(hit->GetTPos()>fEndTPos) fEndTPos=hit->GetTPos();
00067     if(hit->GetTime()<fBegTime) fBegTime=hit->GetTime();
00068     if(hit->GetTime()>fEndTime) fEndTime=hit->GetTime();
00069   }
00070   fDigits += hit->GetCandStripHandle()->GetNDaughters();
00071   fTPos = (fTPos*fCharge+hit->GetTPos()*hit->GetCharge())/(fCharge+hit->GetCharge());
00072   fCharge += hit->GetCharge();
00073 
00074   return;
00075 
00076 }
00078 
00079 
00081 bool ClusterCam::ContainsHit(HitCam* hit)
00082 {
00083   for(unsigned int i=0; i<HitsInCluster.size(); ++i) {
00084     if(hit==HitsInCluster[i]) {return true;}
00085   }
00086 
00087   return false;
00088 }
00090 
00091 
00093 int ClusterCam::IsHitAssoc(HitCam* hit) const
00094 {
00095   double TimeWindow = 9999.9;
00096   if( hit->GetPlane()==fPlane 
00097       && hit->GetTPos()>(fBegTPos-2*StripWidth) && hit->GetTPos()<(fEndTPos+2*StripWidth)
00098       && hit->GetTime()>(fBegTime-TimeWindow) && hit->GetTime()<(fEndTime+TimeWindow) ) {
00099     return 1; 
00100   }
00101   else {return 0;}
00102 }
00104 
00105 
00107 int ClusterCam::IsShwAssoc(ClusterCam* clust) const
00108 {
00109   double TimeWindow = 99.9; int ShwAssocNum = 0;
00110 
00111   if( (fEndTime-clust->GetBegTime())>-TimeWindow || (clust->GetEndTime()-fBegTime)>-TimeWindow ) {
00112 
00113     if( abs(clust->GetPlane()-fPlane)<5 
00114         && (clust->GetEndTPos()-fBegTPos)>-6*StripWidth && (fEndTPos-clust->GetBegTPos())>-6*StripWidth ) {
00115       
00116       if( ( abs(clust->GetPlane()-fPlane)<3 
00117             && (clust->GetEndTPos()-fBegTPos)>-StripWidth && (fEndTPos-clust->GetBegTPos())>-StripWidth )
00118 
00119           || ( clust->GetPlane()==fPlane
00120                && (clust->GetEndTPos()-fBegTPos)>-3*StripWidth && (fEndTPos-clust->GetBegTPos())>-3*StripWidth ) )
00121         
00122         {ShwAssocNum=2;}
00123       
00124       else {ShwAssocNum=1;}
00125     }
00126 
00127   }
00128 
00129   if(ShwAssocNum==2 && this->GetEntries()<3 && clust->GetEntries()<2) {ShwAssocNum=1;}
00130  
00131   return ShwAssocNum;
00132 }
00134 
00135 
00136 
00138 int ClusterCam::IsTrkAssoc(ClusterCam* clustm, ClusterCam* clustp, int PlaneGap) const
00139 {
00140   double TimeWindow = 99.9; int TrkAssocNum = 0;
00141   double min=0.2; double max=0.8;
00142   double NDScale=1;
00143 
00144 
00145   // Configure for correct detector instrumentation
00146   if(PlaneGap==0) {
00147     if(this->GetNDFlag()==2) {PlaneGap=10;}
00148     else {PlaneGap=2;}
00149   }
00150 
00151   if(PlaneGap==10) {NDScale=2;}
00152 
00153 
00154 
00155   // Check timing proximity
00156   if(( (fEndTime-clustm->GetBegTime())>-TimeWindow && (clustm->GetEndTime()-fBegTime)>-TimeWindow)
00157      && ( (fEndTime-clustp->GetBegTime())>-TimeWindow && (clustp->GetEndTime()-fBegTime)>-TimeWindow) ) {
00158 
00159 
00160     // If more than two planes away, scale back width of cluster
00161     // and then treat as if only two planes away
00163     double mBegTPos=clustm->GetBegTPos();
00164     double mEndTPos=clustm->GetEndTPos();
00165     
00166     if((fPlane-clustm->GetPlane())>PlaneGap) {
00167       double mScale = double(PlaneGap)/double(fPlane-clustm->GetPlane());
00168     
00169       mBegTPos=fBegTPos+mScale*(mBegTPos-fBegTPos);
00170       mEndTPos=fEndTPos+mScale*(mEndTPos-fEndTPos);
00171     }
00172     
00173     double pBegTPos=clustp->GetBegTPos();
00174     double pEndTPos=clustp->GetEndTPos();
00175     
00176     if((clustp->GetPlane()-fPlane)>PlaneGap) {
00177       double pScale = double(PlaneGap)/double(clustp->GetPlane()-fPlane);
00178       
00179       pBegTPos=fBegTPos+pScale*(pBegTPos-fBegTPos);
00180       pEndTPos=fEndTPos+pScale*(pEndTPos-fEndTPos);
00181     } 
00183     
00184     
00185     // Scale tolerance for matching clusters for cases where there are gaps
00186     double k0 = 0.5*(clustp->GetPlane()-clustm->GetPlane()-(2*PlaneGap));
00187     // Make sure we do the same thing for the ND +/- 10 plane sections
00188     if(PlaneGap==10) {k0/=5;}
00189 
00190     min = min + 0.1*k0; 
00191     max = max - 0.1*k0; 
00192         
00193 
00194     // Determine associations
00196     // For tracks with +ve dtpos/dz
00197     if( (fEndTPos-clustm->GetBegTPos())>-(NDScale*1.1*StripWidth) && (clustp->GetEndTPos()-fBegTPos)>-(NDScale*1.1*StripWidth) ) {
00198 
00199       // Clusters don't overlap
00200       if( (fBegTPos-clustm->GetEndTPos())>-(NDScale*0.1*StripWidth) || (clustp->GetBegTPos()-fEndTPos)>-(NDScale*0.1*StripWidth) ) {
00201         if( abs( (clustp->GetBegTPos()-fEndTPos)-(fBegTPos-clustm->GetEndTPos()) )<(NDScale*2.1*StripWidth)
00202             || ( ((min*mEndTPos)+(max*pEndTPos))>(fBegTPos-(NDScale*0.5*StripWidth)) && ((max*mBegTPos)+(min*pBegTPos))<(fEndTPos+(NDScale*0.5*StripWidth)) ) )
00203           {TrkAssocNum=2;} 
00204       }
00205       
00206       // Overlapping clusters
00207       if( (clustm->GetEndTPos()-fBegTPos)>-(NDScale*1.1*StripWidth) && (fEndTPos-clustp->GetBegTPos())>-(NDScale*1.1*StripWidth) ) {
00208         if(TrkAssocNum<1) TrkAssocNum=1;
00209       }
00210     }
00211     
00212     
00213     // For tracks with -ve dtpos/dz
00214     if( (fBegTPos-clustm->GetEndTPos())<(NDScale*1.1*StripWidth) && (clustp->GetBegTPos()-fEndTPos)<(NDScale*1.1*StripWidth) ) {
00215       
00216       // Clusters don't overlap
00217       if( (fEndTPos-clustm->GetBegTPos())<(NDScale*0.1*StripWidth) || (clustp->GetEndTPos()-fBegTPos)<(NDScale*0.1*StripWidth) ) {
00218         if( abs( (clustp->GetEndTPos()-fBegTPos)-(fEndTPos-clustm->GetBegTPos()) )<(NDScale*2.1*StripWidth)
00219             || ( ((min*pEndTPos)+(max*mEndTPos))>(fBegTPos-(NDScale*0.5*StripWidth)) && ((max*pBegTPos)+(min*mBegTPos))<(fEndTPos+(NDScale*0.5*StripWidth)) ) )
00220           {TrkAssocNum=2;}
00221       }
00222       
00223       // Overlapping clusters
00224       if( (clustm->GetBegTPos()-fEndTPos)<(NDScale*1.1*StripWidth) && (fBegTPos-clustp->GetEndTPos())<(NDScale*1.1*StripWidth) ) {
00225         if(TrkAssocNum<1) TrkAssocNum=1;
00226       }
00227     }
00229     
00230   }
00231   
00232   return TrkAssocNum;
00233 }
00235 
00236 
00238 HitCam* ClusterCam::GetHit(unsigned int i) const
00239 {
00240   if(i<HitsInCluster.size()) {return HitsInCluster[i];}
00241   else {return 0;}
00242 }