/Home/aooliver/hepsoft/lcdroot/LCDRootApps/PhUtil/src/LCDTrjPart.cxx

Go to the documentation of this file.
00001 // ----------------------------------------------------------------------------
00002 // $Id: LCDTrjPart.cxx,v 1.4 2001/11/25 07:44:03 toshi Exp $
00003 // ----------------------------------------------------------------------------
00004 
00005 #include "LCDTrjPart.h"
00006 #include "LCDTrack.h"
00007 #include <iostream.h>
00008 
00009 ClassImp(LCDTrjPart)   
00010 
00011 //______________________________________________________________________
00012 //                                                                       
00013 // LCDTrjPart
00014 //                                                                
00015 LCDTrjPart::LCDTrjPart() {
00016   Init();
00017 }
00018 
00019 LCDTrjPart::LCDTrjPart(LCDGetParameters* gp)  { 
00020   // Constructor.   Primary job is to assemble various parameters
00021   // describing detector geometry and resolution.
00022   Init();
00023   SetDetectorParameters(gp);
00024 }
00025 
00026 void LCDTrjPart::Init() {
00027   nVolumes     = 0;
00028   pInnerVolume = 0;
00029   pVolumes     = 0;
00030   f_dblevel    = 0;
00031 }
00032 
00033 void LCDTrjPart::SetDetectorParameters(LCDGetParameters* gp)  { 
00034   // Constructor.   Primary job is to assemble various parameters
00035   // describing detector geometry and resolution.
00036   nVolumes     = gp->GetNVolumes();
00037   if (nVolumes == 0) {
00038     Int_t ok = gp->SetDetectorVolumes();
00039     if (ok) {
00040       nVolumes     = gp->GetNVolumes();
00041       pInnerVolume = gp->GetInnerVolume();
00042       pVolumes     = gp->GetDetectorVolumes();
00043     } else  { // make it clear we have no detector
00044       nVolumes     = 0;
00045       pInnerVolume = 0;
00046       pVolumes     = 0;
00047       printf(" Error in TrjPart :: Cannot define detector volumes");
00048     }
00049   }else{
00050     pInnerVolume = gp->GetInnerVolume();
00051     pVolumes     = gp->GetDetectorVolumes();
00052   }
00053   
00054 }
00055 
00056 void LCDTrjPart::Swim(Int_t charge,
00057                       const TVector3& momentum, 
00058                       const TVector3& position, 
00059                       Int_t barend, Double_t tRZ) {  
00060 
00061   LCDSwimTraj* pTraj=&traj;
00062   traj.SetParticleParams(position, momentum, charge);
00063     
00064   LCDSwimTraj::PropagateRet status;
00065   LCDDetectorVolume*  pCur = (LCDDetectorVolume *) pInnerVolume;
00066 
00067   Double_t decrS;              // measured in cm
00068   
00069   //printf("pInnerVolume=%d pVolumes=%d\n",pInnerVolume,pVolumes);
00070   // traj has been initialized with position, momentum, charge
00071   
00072   if (f_dblevel > 10)
00073     cout<<"LCDTrjPart::Swim Start!!!"<<endl;
00074 
00075   while(!(pCur->isInside(position.X(), position.Y(), position.Z()))) {
00076     if (pCur->pNextZ == 0) {
00077 
00078       if (f_dblevel > 10) {
00079         cout<<"   "
00080             <<"pCur->pNextZ==0!!! @ step 0"<<" "
00081             <<"Z="<<pCur->outerZ<<endl;
00082       }
00083 
00084       return;
00085     }
00086     pCur = (LCDDetectorVolume *) pCur->pNextZ;
00087   }
00088 
00089   // First swim through empty central volume
00090   if ((pCur->field != 0.0) && (charge != 0)) { // helix
00091     status = pTraj->HelixToCyl(pCur, &decrS);
00092   } else {  // straight line
00093     status = pTraj->LineToCyl(pCur, &decrS);
00094   }
00095 
00096   if (f_dblevel > 10) {
00097     cout<<"   "
00098         <<"First swim status="<<(Int_t)status<<" "
00099         <<"z="<<pTraj->GetNewPosPtr()->Z()<<" "
00100         <<endl;
00101   }
00102 
00103   if ((status == LCDSwimTraj::INTERSECT_NONE) ||
00104       (status == LCDSwimTraj::INTERSECT_ERROR)   ) return;
00105   
00106   // Check that we didn't exit to the air volume next to the beampipe.
00107   TVector3 pPos(*(pTraj->GetNewPosPtr()));
00108   if ((status == LCDSwimTraj::INTERSECT_ZPLUS) || 
00109       (status == LCDSwimTraj::INTERSECT_ZMINUS)) {
00110     if (pCur->pNextZ == 0) {
00111 
00112       if (f_dblevel > 10) {
00113         cout<<"   "
00114             <<"pCur->pNextZ==0!!! @ step1"<<endl;
00115       }
00116 
00117       return;
00118     }
00119     LCDDetectorVolume* pNext = (LCDDetectorVolume *) pCur->pNextZ;
00120     if ( pPos.Perp() < pNext->outerR ) {
00121 
00122       if (f_dblevel > 10) {
00123         cout<<"   "
00124             <<"pPos.Perp() < pNextVolume->outerR!!!"<<" "
00125             <<"pPos.Perp="<<pPos.Perp()<<" "
00126             <<"outerR="<<pNext->outerR<<" "
00127             <<endl;
00128       }
00129 
00130       return;
00131     }
00132   }
00133 
00134   pTraj->Update();
00135 
00136   while (kTRUE) {
00137 
00138     if (f_dblevel > 10) {
00139       cout<<"   "
00140           <<"In while loop"<<" "
00141           <<"status="<<(Int_t)status<<" "
00142           <<endl;
00143     }
00144 
00145     switch(status) {   // find new volume
00146     case LCDSwimTraj::INTERSECT_ZMINUS:
00147     case LCDSwimTraj::INTERSECT_ZPLUS:
00148       pCur = (LCDDetectorVolume *) pCur->pNextZ;
00149       if (!pCur) {
00150 
00151         if (f_dblevel > 10) {
00152           cout<<"   "
00153               <<"pCur==0!!! in while loop"<<" "
00154               <<"status="<<(Int_t)status<<" "
00155               <<endl;
00156         }
00157 
00158         break;
00159       }
00160       // In the special case in which we're coming from an endcap air
00161       // volume headed in a z-direction, there will typically be more 
00162       // than one adjacent volume.  pNextZ pointer gives us innermost.
00163       // Check this against our current radius.  If necessary, follow
00164       // pNextR pointers until we find an acceptable volume
00165       while ( pTraj->GetNewPosPtr()->Perp() > pCur->outerR ) {
00166         pCur = (LCDDetectorVolume *) pCur->pNextR;
00167         if (!pCur) {
00168 
00169           if (f_dblevel > 10) {
00170             cout<<"   "
00171                 <<"pCur==0!!! in while loop"<<" "
00172                 <<"status="<<(Int_t)status<<" "
00173                 <<endl;
00174           }
00175 
00176           break;
00177         }
00178       }
00179       break;
00180     case LCDSwimTraj::INTERSECT_CYLIN:
00181       // One extra thing to check here as well.  In case we go inwards, 
00182       // there may be more than one component next to the original. 
00183       // pPrevR points
00184       // to the one with smallest z.  If this isn't right (based on
00185       // new position), follow its pNextZ pointer.
00186       pCur = (LCDDetectorVolume *) pCur->pPrevR;
00187       if (!pCur) {
00188 
00189         if (f_dblevel > 10) {
00190           cout<<"   "
00191               <<"pCur==0!!! in while loop"<<" "
00192               <<"status="<<(Int_t)status<<" "
00193               <<endl;
00194         }
00195 
00196         break;
00197       }
00198       while (TMath::Abs(pTraj->GetNewPosPtr()->Z()) > pCur->outerZ) {
00199         pCur = (LCDDetectorVolume *) pCur->pNextZ;
00200         if (!pCur) {
00201 
00202           if (f_dblevel > 10) {
00203             cout<<"   "
00204                 <<"pCur==0!!! in while loop"<<" "
00205                 <<"status="<<(Int_t)status<<" "
00206                 <<endl;
00207           }
00208 
00209           break;
00210         }
00211       }
00212       break;
00213     case LCDSwimTraj::INTERSECT_CYLOUT:
00214       pCur = (LCDDetectorVolume *) pCur->pNextR;
00215       break;
00216     default:
00217       return;
00218     }
00219     
00220     if (!pCur) return;  // Went outside detector
00221     
00222     // Propagate to boundary of volume
00223     Bool_t helix = ((pCur->field != 0.0) && (charge != 0));
00224     
00225     Int_t targetvol=1;
00226     Double_t tmprz=0;
00227     if (barend == 0) { // barrel
00228       if (tRZ < pCur->outerR) {
00229         tmprz = pCur->outerR;
00230         pCur->outerR = tRZ;
00231         targetvol=0;
00232       }
00233     } else { // endcap
00234       if (TMath::Abs(tRZ) < pCur->outerZ ) {
00235         tmprz = pCur->outerZ;
00236         pCur->outerZ = TMath::Abs(tRZ);
00237         targetvol=0;
00238       }
00239     }
00240     if (helix) {//helix
00241       status = pTraj->HelixToCyl(pCur, &decrS);
00242     } else {  // straight line
00243       status = pTraj->LineToCyl(pCur, &decrS);
00244     }
00245     if (!targetvol) {
00246       if (barend==0) {
00247         pCur->outerR=tmprz;
00248       } else {
00249         pCur->outerZ=tmprz;
00250       }
00251     }
00252     if ((status == LCDSwimTraj::INTERSECT_NONE) ||
00253         (status == LCDSwimTraj::INTERSECT_ERROR)  ) return;
00254     
00255     if (targetvol) { // keep going
00256       pTraj->Update();
00257     } else {  // We'll finish in this volume
00258       if (helix) {
00259         status = pTraj->HelixPath(pCur->field, decrS);
00260       } else {
00261         status = pTraj->LinePath(decrS); 
00262       }
00263 
00264       if (f_dblevel > 10) {
00265         cout<<"   "
00266             <<"Swim end"<<" "
00267             <<"status"<<(Int_t)status<<" "
00268             <<"xposZ="<<traj.GetNewPosPtr()->Z()<<" "
00269             <<"helix="<<(Int_t)helix<<" "
00270             <<"tRZ="<<tRZ<<" "
00271             <<endl;
00272       }
00273       
00274       pTraj->Update();
00275 
00276       if (barend == 0 || 
00277           (barend != 0 && status == LCDSwimTraj::INTERSECT_ZPLUS) ||
00278           (barend != 0 && status == LCDSwimTraj::INTERSECT_ZMINUS)) {
00279         return;
00280       }
00281 
00282     }
00283   }
00284   return;
00285 }
00286 
00287 void LCDTrjPart::Swim(LCDEvent* event,Int_t itrk,Int_t barend, Double_t tRZ) {
00288   LCDTrack*  trk=(LCDTrack*) (event->Tracks()->At(itrk));
00289   LCDMcPart* mcp=(LCDMcPart*)(event->MCparticles()->At(trk->GetParticle()));
00290   Swim(trk,mcp,barend,tRZ);
00291 }
00292 
00293 void LCDTrjPart::Swim(LCDTrack* trk,LCDMcPart* pmc,Int_t barend, Double_t tRZ){
00294   TVector3 ptrk(trk->GetMomentumVector(0.0));
00295   TVector3 xtrk(trk->GetPositionVector(0.0));
00296 
00297   if (f_dblevel > 10) {
00298     cout<<"LCDTrjPart::Swim Original Track Swim"<<endl;
00299   }
00300 
00301   Int_t chg=trk->GetCharge() > 0 ? 1 : -1 ;
00302   Swim(chg,ptrk,xtrk,barend,tRZ);
00303   if (barend != 0) { //endcap cluster...
00304     Double_t aomega=trk->GetTrackParameter(2);//=1/rho
00305     Double_t tanl  =trk->GetTrackParameter(4);//=tan lambda
00306     Double_t cosl  =1.0/TMath::Sqrt(1.0 + tanl*tanl);
00307     Double_t phi   =traj.GetS()*cosl*aomega;
00308     if (phi > 2*TMath::Pi()) {
00309       TVector3 xold(*(traj.GetNewPosPtr()));
00310 
00311       TVector3 p_orig(pmc->Get4MomentumPtr()->Vect());
00312       TVector3 x_orig(*(pmc->GetPositionPtr()));
00313       TVector3 p_poca;
00314       TVector3 x_poca;
00315       CalcPOCA(pmc->GetCharge(),trk->GetMagneticField(),
00316                p_orig,x_orig,p_poca,x_poca);
00317 
00318       if (f_dblevel > 10) {
00319         cout<<"LCDTrjPart::Swim MC track Swim"<<endl;
00320       }
00321 
00322       Int_t chgmc=pmc->GetCharge() > 0 ? 1 : -1 ;
00323       Swim(chgmc,p_poca,x_poca,barend,tRZ);
00324       TVector3 xtru(*(traj.GetNewPosPtr()));
00325 
00326       Double_t omega0  =-pmc->GetCharge()*trk->GetMagneticField()
00327         /333.567/p_orig.Pt();
00328       Double_t z00     = x_poca.Z();
00329       Double_t tanl0   = p_orig.Pz()/p_orig.Pt();
00330       Double_t tmpphi  = (tRZ - z00)/tanl0*omega0;
00331       Int_t    n2pi    = (Int_t)(TMath::Abs(tmpphi)/TMath::Pi()/2.0);
00332       Double_t l       = TMath::Abs(n2pi*2*TMath::Pi()/omega0);
00333       Double_t zoffset = tanl0*l;
00334       xtrk.SetZ(xtrk.Z()+zoffset);
00335 
00336       if (f_dblevel > 10) {
00337         cout<<"LCDTrjPart::Swim new track Swim"<<endl;
00338       }
00339 
00340       Int_t chgtk=trk->GetCharge() > 0 ? +1 : -1 ;
00341       Swim(chgtk,ptrk,xtrk,barend,tRZ);
00342 
00343       TVector3 xnew(*(traj.GetNewPosPtr()));
00344 
00345       if (f_dblevel > 0) {
00346         cout<<"---"
00347             <<"xyz mcp="<<xtru.X()<<","<<xtru.Y()<<","<<xtru.Z()<<" "
00348             <<"tRZ="<<tRZ
00349             <<endl;
00350         cout<<"   "
00351             <<"xyz old="<<xold.X()<<","<<xold.Y()<<","<<xold.Z()<<" "
00352             <<"dist="<<(xold-xtru).Perp()
00353             <<endl;
00354         cout<<"   "
00355             <<"xyz new="<<xnew.X()<<","<<xnew.Y()<<","<<xnew.Z()<<" "
00356             <<"dist="<<(xnew-xtru).Perp()
00357             <<endl;
00358       }
00359 
00360     }
00361   }
00362 }
00363 
00364 void LCDTrjPart::CalcPOCA(Double_t  qMC,
00365                           Double_t  bfld_z,
00366                           const TVector3&  p_orig, 
00367                           const TVector3&  x_orig, 
00368                           TVector3& p_poca, 
00369                           TVector3& x_poca) {
00370   Double_t pT   = p_orig.Pt();
00371   Double_t tanL = p_orig.Z()/pT;
00372   Double_t arho = 333.567/bfld_z*pT;
00373   
00374   TVector3 nzv(0.0,0.0,qMC/TMath::Abs(qMC));
00375   TVector3 xxc = p_orig.Cross(nzv);
00376   xxc.SetZ(0.0);
00377   TVector3 nxxc = xxc.Unit();
00378   TVector3 xc   = x_orig + arho*nxxc;
00379   xc.SetZ(0.0);
00380   TVector3 nxc  = xc.Unit();
00381   
00382   TVector3 catMC=nzv.Cross(nxc);
00383   catMC.SetZ(0.0);
00384   TVector3 ncatMC=catMC.Unit();
00385 
00386   p_poca    = pT*ncatMC;
00387   p_poca.SetZ(p_orig.Z());
00388   
00389   Double_t dphi     = TMath::ASin(nxxc.Cross(nxc).Z());
00390   Double_t dl       =-dphi*arho*qMC;
00391   x_poca            = xc - arho*nxc;
00392   x_poca.SetZ(x_orig.Z() + dl*tanL);
00393   
00394 }

Generated on Tue Jul 18 18:33:59 2006 for LCDROOT by  doxygen 1.4.6