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

Go to the documentation of this file.
00001 // ----------------------------------------------------------------------------
00002 // $Id: LCDEventShape.cxx,v 1.2 2001/10/08 17:52:54 toshi Exp $
00003 // ----------------------------------------------------------------------------
00004 //
00005 //
00006 // V0.0           :            Thrust finder utility translated from 
00007 //                             Java routines written by G.Bower. 
00008 // V0.2 Jul 07/99 : M.Iwasaki  Change Mod to Mag function in TVector3
00009 //                             so as to use root 2.2.x
00010 // V0.3 Sep 23/99 : M.Iwasaki  Fix setEvent memory leak, 
00011 //                             apply nessesary modifications in Thrust, Major, 
00012 //                             Minor axis, and Thrust, and add ~EventSape().
00013 // V0.4 May 12/00 : M.Iwasaki  Make necessary modifications.
00014 //                                        
00015 // V1.0 Aug   /00 : M.Iwasaki  Change Class name from EventShape to 
00016 //                             LCDEventShape.
00017 //
00018 // V1.1 Aug 29/00 : M.Iwasaki  Remove Thrust(), and add GetThrust(), GetMajor(),
00019 //                             and GetMiner(), so as to provide thrust, Major, 
00020 //                             and Miner values, respectively.
00021 //                             change class names:: thrustAxis -> GetThrustAxis etc..
00022 
00023 #include "LCDEventShape.h"
00024 
00025 ClassImp(LCDEventShape)
00026 
00027 Int_t LCDEventShape::m_maxpart = 1000;
00028 
00029 LCDEventShape::LCDEventShape():
00030   m_dSphMomPower(2.0),m_dDeltaThPower(0),
00031   m_iFast(4),m_dConv(0.0001),m_iGood(2)
00032 {
00033   m_dAxes.ResizeTo(4,4);
00034 }
00035 //______________________________________________________________
00036 
00037 LCDEventShape::~LCDEventShape(){
00038 }
00039 //______________________________________________________________
00040 
00041 // Input the particle 3(4)-vector list
00042 // e: 3-vector  TVector3       ..(px,py,pz) or
00043 //    4-vector  TLorentzVector ..(px,py,pz,E) 
00044 // Even input the TLorentzVector, we don't use Energy 
00045 // 
00046 void LCDEventShape::SetPartList(TObjArray* e)
00047 {       
00048   //To make this look like normal physics notation the
00049   //zeroth element of each array, mom[i][0], will be ignored
00050   //and operations will be on elements 1,2,3...
00051   TMatrixD mom(m_maxpart,6);
00052   Double_t tmax = 0;
00053   Double_t phi = 0.;
00054   Double_t the = 0.;
00055   Double_t sgn;
00056   TMatrixD fast(m_iFast + 1,6);
00057   TMatrixD work(11,6);
00058   Double_t tdi[] = {0.,0.,0.,0.};
00059   Double_t tds;
00060   Double_t tpr[] = {0.,0.,0.,0.};
00061   Double_t thp;
00062   Double_t thps;
00063   TMatrixD temp(3,5);
00064   Int_t np = 0;
00065   Int_t numElements = e->GetEntries();
00066   TObject* o;
00067   for(Int_t elem=0;elem<numElements;elem++) {
00068     o = e->At(elem);
00069     
00070     if (np >= m_maxpart) { 
00071       printf("Too many particles input to LCDEventShape");
00072       return;
00073     }
00074 
00075     TString nam(o->IsA()->GetName());
00076     if (nam.Contains("TVector3")) {
00077       mom(np,1) = ((TVector3 *) o)->X();
00078       mom(np,2) = ((TVector3 *) o)->Y();
00079       mom(np,3) = ((TVector3 *) o)->Z();
00080       mom(np,4) = TMath::Sqrt(mom(np,1)*mom(np,1) + mom(np,2)*mom(np,2)
00081                               + mom(np,3)*mom(np,3));
00082     } else if (nam.Contains("TLorentzVector")) {
00083       mom(np,1) = ((TLorentzVector *) o)->X();
00084       mom(np,2) = ((TLorentzVector *) o)->Y();
00085       mom(np,3) = ((TLorentzVector *) o)->Z();
00086       mom(np,4) = TMath::Sqrt(mom(np,1)*mom(np,1) + mom(np,2)*mom(np,2)
00087                               + mom(np,3)*mom(np,3));
00088     } else {
00089       printf("LCDEventShape::SetEvent input is not a TVector3 or a TLorentzVector\n");
00090       continue;
00091     }
00092 
00093     if ( TMath::Abs( m_dDeltaThPower ) <= 0.001 ) {
00094       mom(np,5) = 1.0;
00095     } else {
00096       mom(np,5) = TMath::Power(mom(np,4),m_dDeltaThPower);
00097     }
00098     tmax = tmax + mom(np,4)*mom(np,5);
00099     np++;
00100   }
00101   if ( np < 2 ) {
00102     m_dThrust[1] = -1.0;
00103     m_dOblateness = -1.0;
00104     return;
00105   }
00106   // for pass = 1: find thrust axis.
00107   // for pass = 2: find major axis.
00108   for ( Int_t pass=1; pass < 3; pass++) {
00109     if ( pass == 2 ) {
00110       phi = ulAngle(m_dAxes(1,1), m_dAxes(1,2));
00111       ludbrb( &mom, 0, -phi, 0., 0., 0. );
00112       for ( Int_t i = 0; i < 3; i++ ) {
00113         for ( Int_t j = 1; j < 4; j++ ) {
00114           temp(i,j) = m_dAxes(i+1,j);
00115         }
00116         temp(i,4) = 0;
00117       }
00118       ludbrb(&temp,0.,-phi,0.,0.,0.);
00119       for ( Int_t ib = 0; ib < 3; ib++ ) {
00120         for ( Int_t j = 1; j < 4; j++ ) {
00121           m_dAxes(ib+1,j) = temp(ib,j);
00122         }
00123       }
00124       the = ulAngle( m_dAxes(1,3), m_dAxes(1,1) );
00125       ludbrb( &mom, -the, 0., 0., 0., 0. );
00126       for ( Int_t ic = 0; ic < 3; ic++ ) {
00127         for ( Int_t j = 1; j < 4; j++ ) {
00128           temp(ic,j) = m_dAxes(ic+1,j);
00129         }
00130         temp(ic,4) = 0;
00131       }
00132       ludbrb(&temp,-the,0.,0.,0.,0.);
00133       for ( Int_t id = 0; id < 3; id++ ) {      
00134         for ( Int_t j = 1; j < 4; j++ ) {
00135           m_dAxes(id+1,j) = temp(id,j);
00136         }
00137       }
00138     }
00139     for ( Int_t ifas = 0; ifas < m_iFast + 1 ; ifas++ ) {
00140       fast(ifas,4) = 0.;
00141     }
00142     // Find the m_iFast highest momentum particles and
00143     // put the highest in fast[0], next in fast[1],....fast[m_iFast-1].
00144     // fast[m_iFast] is just a workspace.
00145     for ( Int_t i = 0; i < np; i++ ) {
00146       if ( pass == 2 ) {
00147         mom(i,4) = TMath::Sqrt( mom(i,1)*mom(i,1) 
00148                                 + mom(i,2)*mom(i,2) ); 
00149       }
00150       for ( Int_t ifas = m_iFast - 1; ifas > -1; ifas-- ) {
00151         if ( mom(i,4) > fast(ifas,4) ) {
00152           for ( Int_t j = 1; j < 6; j++ ) {
00153             fast(ifas+1,j) = fast(ifas,j);
00154             if ( ifas == 0 ) fast(ifas,j) = mom(i,j);       
00155           }
00156         } else {
00157           for ( Int_t j = 1; j < 6; j++ ) {
00158             fast(ifas+1,j) = mom(i,j);
00159           }
00160           break;
00161         }
00162       }
00163     }
00164     // Find axis with highest thrust (case 1)/ highest major (case 2).
00165     for ( Int_t ie = 0; ie < work.GetNrows(); ie++ ) {
00166       work(ie,4) = 0.;
00167     }
00168     Int_t p = TMath::Min( m_iFast, np ) - 1;
00169     // Don't trust Math.pow to give right answer always.
00170     // Want nc = 2**p.
00171     Int_t nc = iPow(2,p); 
00172     for ( Int_t n = 0; n < nc; n++ ) {
00173       for ( Int_t j = 1; j < 4; j++ ) {
00174         tdi[j] = 0.;
00175       }
00176       for ( Int_t i = 0; i < TMath::Min(m_iFast,n); i++ ) {
00177         sgn = fast(i,5);
00178         if (iPow(2,(i+1))*((n+iPow(2,i))/iPow(2,(i+1))) >= i+1){
00179           sgn = -sgn;
00180         }
00181         for ( Int_t j = 1; j < 5-pass; j++ ) {
00182           tdi[j] = tdi[j] + sgn*fast(i,j);
00183         }
00184       }
00185       tds = tdi[1]*tdi[1] + tdi[2]*tdi[2] + tdi[3]*tdi[3];
00186       for ( Int_t iw = TMath::Min(n,9); iw > -1; iw--) {
00187         if( tds > work(iw,4) ) {
00188           for ( Int_t j = 1; j < 5; j++ ) {
00189             work(iw+1,j) = work(iw,j);
00190             if ( iw == 0 ) {
00191               if ( j < 4 ) {
00192                 work(iw,j) = tdi[j];
00193               } else {
00194                 work(iw,j) = tds;
00195               }
00196             }
00197           }
00198         } else {
00199           for ( Int_t j = 1; j < 4; j++ ) {
00200             work(iw+1,j) = tdi[j];
00201           }
00202           work(iw+1,4) = tds;
00203         }
00204       }
00205     }
00206     // Iterate direction of axis until stable maximum.
00207     m_dThrust[pass] = 0;
00208     thp = -99999.;
00209     Int_t nagree = 0;
00210     for ( Int_t iw = 0; 
00211           iw < TMath::Min(nc,10) && nagree < m_iGood; iw++ ){
00212       thp = 0.;
00213       thps = -99999.;
00214       while ( thp > thps + m_dConv ) {
00215         thps = thp;
00216         for ( Int_t j = 1; j < 4; j++ ) {
00217           if ( thp <= 1E-10 ) {
00218             tdi[j] = work(iw,j);
00219           } else {
00220             tdi[j] = tpr[j];
00221             tpr[j] = 0;
00222           }
00223         }
00224         for ( Int_t i = 0; i < np; i++ ) {
00225           sgn = sign(mom(i,5), 
00226                      tdi[1]*mom(i,1) + 
00227                      tdi[2]*mom(i,2) + 
00228                      tdi[3]*mom(i,3));
00229           for ( Int_t j = 1; j < 5 - pass; j++ ){
00230             tpr[j] = tpr[j] + sgn*mom(i,j);
00231           }
00232         }
00233         thp = TMath::Sqrt(tpr[1]*tpr[1] 
00234                           + tpr[2]*tpr[2] 
00235                           + tpr[3]*tpr[3])/tmax;
00236       }
00237       // Save good axis. Try new initial axis until enough
00238       // tries agree.
00239       if ( thp < m_dThrust[pass] - m_dConv ) {
00240         break;
00241       }
00242       if ( thp > m_dThrust[pass] + m_dConv ) {
00243         nagree = 0;
00244         sgn = iPow( -1, (Int_t)TMath::Nint(m_random.Rndm()) );
00245         for ( Int_t j = 1; j < 4; j++ ) {
00246           m_dAxes(pass,j) = sgn*tpr[j]/(tmax*thp);
00247         }
00248         m_dThrust[pass] = thp;
00249       }
00250       nagree = nagree + 1;
00251     }
00252   }
00253   // Find minor axis and value by orthogonality.
00254   sgn = iPow( -1, (Int_t)TMath::Nint(m_random.Rndm()));
00255   m_dAxes(3,1) = -sgn*m_dAxes(2,2);
00256   m_dAxes(3,2) = sgn*m_dAxes(2,1);
00257   m_dAxes(3,3) = 0.;
00258   thp = 0.;
00259   for ( Int_t i = 0; i < np; i++ ) {
00260     thp += mom(i,5)*TMath::Abs(m_dAxes(3,1)*mom(i,1) + 
00261                                m_dAxes(3,2)*mom(i,2));
00262   }
00263   m_dThrust[3] = thp/tmax;
00264   // Rotate back to original coordinate system.
00265   for ( Int_t i6 = 0; i6 < 3; i6++ ) {
00266     for ( Int_t j = 1; j < 4; j++ ) {
00267       temp(i6,j) = m_dAxes(i6+1,j);
00268     }
00269     temp(i6,4) = 0;
00270   }
00271   ludbrb(&temp,the,phi,0.,0.,0.);
00272   for ( Int_t i7 = 0; i7 < 3; i7++ ) {
00273     for ( Int_t j = 1; j < 4; j++ ) {
00274       m_dAxes(i7+1,j) = temp(i7,j);
00275     }
00276   }
00277   m_dOblateness = m_dThrust[2] - m_dThrust[3];
00278   
00279 }
00280 //______________________________________________________________
00281 
00282 // Setting and getting parameters.
00283 
00284 void LCDEventShape::SetThMomPower(Double_t tp)
00285 {
00286   // Error if sp not positive.
00287   if ( tp > 0. ) m_dDeltaThPower = tp - 1.0;
00288   return;
00289 }
00290 //______________________________________________________________
00291 
00292 Double_t LCDEventShape::GetThMomPower()
00293 {
00294   return 1.0 + m_dDeltaThPower;
00295 }
00296 //______________________________________________________________
00297 
00298 void LCDEventShape::SetFast(Int_t nf)
00299 {
00300   // Error if sp not positive.
00301   if ( nf > 3 ) m_iFast = nf;
00302   return;
00303 }
00304 //______________________________________________________________
00305 
00306 Int_t LCDEventShape::GetFast()
00307 {
00308   return m_iFast;
00309 }
00310 //______________________________________________________________
00311 
00312 // Returning results
00313 
00314 TVector3 LCDEventShape::GetThrustAxis() {
00315   return TVector3(m_dAxes(1,1),m_dAxes(1,2),m_dAxes(1,3));
00316 }
00317 //______________________________________________________________
00318 
00319 TVector3 LCDEventShape::GetMajorAxis() {
00320   return TVector3(m_dAxes(2,1),m_dAxes(2,2),m_dAxes(2,3));
00321 }
00322 //______________________________________________________________
00323 
00324 TVector3 LCDEventShape::GetMinorAxis() {
00325   return TVector3(m_dAxes(3,1),m_dAxes(3,2),m_dAxes(3,3));
00326 }
00327 //______________________________________________________________
00328 
00329 Double_t LCDEventShape::oblateness() {
00330   return m_dOblateness;
00331 }
00332 //______________________________________________________________
00333 
00334 // utilities(from Jetset):
00335 Double_t LCDEventShape::ulAngle(Double_t x, Double_t y)
00336 {
00337   Double_t ulangl = 0;
00338   Double_t r = TMath::Sqrt(x*x + y*y);
00339   if ( r < 1.0E-20 ) {
00340     return ulangl; 
00341   }
00342   if ( TMath::Abs(x)/r < 0.8 ) {
00343     ulangl = sign(TMath::ACos(x/r),y);
00344   } else {
00345     ulangl = TMath::ASin(y/r);
00346     if ( x < 0. && ulangl >= 0. ) {
00347       ulangl = TMath::Pi() - ulangl;
00348     } else if ( x < 0. ) {
00349       ulangl = -TMath::Pi() - ulangl;
00350     }
00351   }
00352   return ulangl;
00353 }
00354 //______________________________________________________________
00355 
00356 Double_t LCDEventShape::sign(Double_t a, Double_t b) {
00357   if ( b < 0 ) {
00358     return -TMath::Abs(a);
00359   } else {
00360     return TMath::Abs(a);
00361   }
00362 }
00363 //______________________________________________________________
00364 
00365 void LCDEventShape::ludbrb(TMatrixD* mom, 
00366                         Double_t the, 
00367                         Double_t phi, 
00368                         Double_t bx, 
00369                         Double_t by,
00370                         Double_t bz)
00371 {
00372   // Ignore "zeroth" elements in rot,pr,dp.
00373   // Trying to use physics-like notation.
00374   TMatrixD rot(4,4);
00375   Double_t pr[4];
00376   Double_t dp[5];
00377   Int_t np = mom->GetNrows();
00378   if ( the*the + phi*phi > 1.0E-20 )
00379     {
00380       rot(1,1) = TMath::Cos(the)*TMath::Cos(phi);
00381       rot(1,2) = -TMath::Sin(phi);
00382       rot(1,3) = TMath::Sin(the)*TMath::Cos(phi);
00383       rot(2,1) = TMath::Cos(the)*TMath::Sin(phi);
00384       rot(2,2) = TMath::Cos(phi);
00385       rot(2,3) = TMath::Sin(the)*TMath::Sin(phi);
00386       rot(3,1) = -TMath::Sin(the);
00387       rot(3,2) = 0.0;
00388       rot(3,3) = TMath::Cos(the);
00389       for ( Int_t i = 0; i < np; i++ ) {
00390         for ( Int_t j = 1; j < 4; j++ ) {
00391           pr[j] = (*mom)(i,j);
00392           (*mom)(i,j) = 0;
00393         }
00394         for ( Int_t jb = 1; jb < 4; jb++) {
00395           for ( Int_t k = 1; k < 4; k++) {
00396             (*mom)(i,jb) = (*mom)(i,jb) + rot(jb,k)*pr[k];
00397           }
00398         }
00399       }
00400       Double_t beta = TMath::Sqrt( bx*bx + by*by + bz*bz );
00401       if ( beta*beta > 1.0E-20 ) {
00402         if ( beta >  0.99999999 ) {
00403           //send message: boost too large, resetting to <~1.0.
00404           bx = bx*(0.99999999/beta);
00405           by = by*(0.99999999/beta);
00406           bz = bz*(0.99999999/beta);
00407           beta =   0.99999999;
00408         }
00409         Double_t gamma = 1.0/TMath::Sqrt(1.0 - beta*beta);
00410         for ( Int_t i = 0; i < np; i++ ) {
00411           for ( Int_t j = 1; j < 5; j++ ) {
00412             dp[j] = (*mom)(i,j);
00413           }
00414           Double_t bp = bx*dp[1] + by*dp[2] + bz*dp[3];
00415           Double_t gbp = gamma*(gamma*bp/(1.0 + gamma) + dp[4]);
00416           (*mom)(i,1) = dp[1] + gbp*bx;
00417           (*mom)(i,2) = dp[2] + gbp*by;
00418           (*mom)(i,3) = dp[3] + gbp*bz;
00419           (*mom)(i,4) = gamma*(dp[4] + bp);
00420         }
00421       }
00422     }
00423   return;
00424 }
00425 //______________________________________________________________
00426 
00427 Int_t LCDEventShape::iPow(Int_t man, Int_t exp)
00428 {
00429   Int_t ans = 1;
00430   for( Int_t k = 0; k < exp; k++) {
00431     ans = ans*man;
00432   }
00433   return ans;
00434 }
00435 
00436 

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