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

Go to the documentation of this file.
00001 // ----------------------------------------------------------------------------
00002 // $Id: LCDJetFinder.cxx,v 1.3 2001/10/08 17:53:00 toshi Exp $
00003 // ----------------------------------------------------------------------------
00004 //
00005 // Revision 1.1  2001/05/10 19:32:06  toshi
00006 // C++ file name convention has been changed from .C to .cxx .
00007 //
00008 //
00009 // Jet finder base class routines.
00010 //
00011 // V0.0 Mar 01 99 : R. Shanks  Derived from Java routines written by G.Bower. 
00012 // V0.1 Jul 02 99 : M.Iwasaki  Make some Modification on masscut and ymass 
00013 //                             calculation (in JadeEJetfonder.cxx 
00014 //                             or DurhamJetFinder.cxx) on doFindJets().
00015 //                             Adding Jade cluster (yclus).
00016 // V0.2 Jul 07 99 : M.Iwasaki  Change Mod to Mag function in TVector3
00017 //                             so as to use root 2.22.x
00018 // V0.3 Jul 21 99 : M.Iwasaki  Make temporary 4-vectors copied from input 
00019 //                             particle 4-vectors in doFindJets() 
00020 //                             to protect the initial 4-vectors.
00021 // V0.4 Sep 23 99 : M.Iwasaki  Fix setEvent memory leak.
00022 //
00023 // V1.0 May 15 00 : M.Iwasaki  Make necessary modifications, and change classes
00024 //                             Merge JadeEJetFinder, JadeJetFinder, 
00025 //                             DurhamJetfinder, and JetFinder into one JetFinder.
00026 // V1.1 Aug 15 00 : M.Iwasaki  Change calcinvmass. Remove DurhamJetFinder, 
00027 //                             JadeJetFinder, and JadeEJetFinder classes.
00028 //                             Add SetDURHAM, SetJADE, or SetJADEE. 
00029 //                             Now we can change algorithms by one of them.
00030 //                             Change Class name from JetFinder to LCDJetFinder.
00031 
00032 #include "LCDJetFinder.h"
00033 
00034 //_____________________________________________________________________
00035 //  ------------------
00036 //   JetFinder Class
00037 //  ------------------
00038 //
00039 ClassImp(LCDJetFinder)
00040 
00041 const Int_t LCDJetFinder::UNASSOC = -999;
00042 
00043 const Int_t LCDJetFinder::DURHAM = 1;
00044 const Int_t LCDJetFinder::JADE   = 2;
00045 const Int_t LCDJetFinder::JADEE  = 3;
00046 
00047 LCDJetFinder::LCDJetFinder(): 
00048     m_ymin(1e20),m_evis(0.),m_dycut(0.1),
00049     m_algorithm(LCDJetFinder::DURHAM){ // Default constructor  uses Durham
00050     m_jet  = new TClonesArray("TLorentzVector", 20);
00051     m_part = new TClonesArray("TLorentzVector",200);
00052 }
00053 
00054 LCDJetFinder::LCDJetFinder(Double_t ycut): 
00055     m_ymin(0.0),m_evis(0.),m_dycut(ycut),
00056     m_algorithm(LCDJetFinder::DURHAM){ // Default constructor  uses Durham
00057     m_jet  = new TClonesArray("TLorentzVector", 20);
00058     m_part = new TClonesArray("TLorentzVector",200);
00059 }
00060 //______________________________________________________
00061 
00062 LCDJetFinder::~LCDJetFinder() {
00063   m_jet->Delete();  delete m_jet;
00064   m_part->Delete(); delete m_part;
00065 }
00066 //______________________________________________________
00067 
00068 TLorentzVector* LCDJetFinder::jet4vec(Int_t index) {
00069   return (TLorentzVector*)m_jet->At(index);
00070 }
00071 //_______________________________________________________
00072 
00073 Int_t LCDJetFinder::nParticlesPerJet(Int_t index) {
00074   return m_inparts_per_jet[index];
00075 }
00076 //_______________________________________________________
00077 
00078 void LCDJetFinder::SetYCut(Double_t ycut) {
00079   m_dycut = ycut;
00080 }
00081 //_______________________________________________________
00082 
00083 // Input the particle 4(3)-vector list
00084 // e: 4-vector  TLorentzVector ..(px,py,pz,E) or
00085 //    3-vector  TVector3       ..(px,py,pz) 
00086 // If you input TVector3, the energy of particle
00087 // will be E = sqrt(px**2 + py**2 + pz**2)
00088 void LCDJetFinder::SetPartList(TObjArray* e) {
00089 
00090   m_evis = 0;
00091   m_4vec = e;
00092 
00093   Int_t ne = e->GetEntries();
00094   TObject* o;
00095   for(Int_t i=0;i<ne;i++) {
00096     o = m_4vec->At(i);
00097     TString nam(o->IsA()->GetName());
00098     if (nam.Contains("TLorentzVector")) {
00099       m_evis += ((TLorentzVector *) o)->T();
00100     } else if (nam.Contains("TVector3")) {
00101       m_evis += ((TVector3 *) o)->Mag();
00102     } else {
00103       printf("LCDJetFinder::SetEvent input is not a TVector3 or a TLorentzVector\n");
00104     }
00105   }
00106 
00107 }
00108 //_______________________________________________________
00109 
00110 void LCDJetFinder::doFindJets(Double_t ycut){
00111   SetYCut(ycut);
00112   doFindJets();
00113 }
00114 //______________________________________________________
00115 
00116 void LCDJetFinder::doFindJets(){
00117 
00118   Int_t np = m_4vec->GetEntries();      
00119   if (np < 2) return;
00120 
00121   TObject* o;
00122   Int_t ipart;
00123   for(ipart=0; ipart < np ; ipart++) {
00124     o = m_4vec->At(ipart);
00125     TString nam(o->IsA()->GetName());
00126     if (nam.Contains("TLorentzVector")) {
00127       new((*m_part)[ipart]) TLorentzVector(*((TLorentzVector*)o));
00128     } else if (nam.Contains("TVector3")) {
00129       new((*m_part)[ipart]) TLorentzVector(((TVector3 *) o)->X(),
00130                                            ((TVector3 *) o)->Y(),
00131                                            ((TVector3 *) o)->Z(),
00132                                            ((TVector3 *) o)->Mag());
00133     } else {
00134       printf("LCDJetFinder::SetEvent input is not a TVector3 or a TLorentzVector\n");
00135     }
00136   }
00137 
00138   m_ipart_jet_assoc.Reset();
00139   m_ipart_jet_assoc.Set(np);
00140   for (Int_t m=0; m<np; m++) m_ipart_jet_assoc[m] = UNASSOC;
00141 
00142   m_njets = 0;
00143 
00144   //
00145   // create invariant mass pair array.
00146   //
00147   TMatrix ymass = TMatrix(np,np);
00148   
00149   for (Int_t i1 = 0; i1 < np - 1; i1++ ) {
00150     for (Int_t i2 = i1 + 1 ; i2 < np ; i2++ ) {
00151       ymass(i1,i2) = calcinvmass(*(TLorentzVector*)m_part->At(i1),
00152                                  *(TLorentzVector*)m_part->At(i2));
00153     }
00154   }
00155   
00156   Double_t masscut = m_dycut * m_evis * m_evis ;
00157 
00158   Int_t im,jm;
00159   Double_t minmass;
00160   Int_t i,j;
00161   Int_t imin;
00162   Int_t imax;
00163   for (;;) {
00164     im = -1;
00165     jm = -1;
00166     minmass = 100000000.;
00167     if ( minmass <= masscut ) minmass = masscut * 10.;
00168     //
00169     // find least invariant mass pair.
00170     //
00171     for(i = 0 ; i < np - 1 ; i++ ) {
00172       for(j = i + 1 ; j < np ; j++ ) {
00173         if (m_ipart_jet_assoc[i] != LCDJetFinder::UNASSOC) continue;
00174         if (m_ipart_jet_assoc[j] != LCDJetFinder::UNASSOC) continue;
00175         if (ymass(i,j) > minmass) continue;
00176           
00177         minmass = ymass(i,j);
00178         im = i;  jm = j;
00179       }
00180     }
00181 
00182     m_ymin=minmass/m_evis/m_evis;
00183 
00184     if (minmass > masscut) break;
00185     //
00186     // combine particles im and jm.
00187     //
00188     *((TLorentzVector*)m_part->At(im)) += *((TLorentzVector*)m_part->At(jm));
00189 
00190     for(ipart = 0; ipart < np; ipart++ ){
00191       if(m_ipart_jet_assoc[ipart] == jm) m_ipart_jet_assoc[ipart] = im;
00192     }
00193     //
00194     // Recalculate invariant masses for newly combined particle
00195     //
00196     m_ipart_jet_assoc[jm] = im;
00197     for (ipart = 0; ipart < np ; ipart++) {
00198       if (ipart == im) continue;
00199       if (m_ipart_jet_assoc[ipart] != UNASSOC ) continue;
00200       
00201       imin = TMath::Min(ipart,im);
00202       imax = TMath::Max(ipart,im);
00203       
00204       ymass(imin,imax) = calcinvmass(*((TLorentzVector*)m_part->At(imin)),
00205                                      *((TLorentzVector*)m_part->At(imax)));
00206     }
00207     
00208   }
00209   //
00210   // finish up by filling jet array.
00211   //
00212   
00213   for (ipart = 0 ; ipart < np ; ipart++) {
00214     if (m_ipart_jet_assoc[ipart] == UNASSOC) m_njets++;                 
00215   }
00216   
00217   m_jet->Delete();
00218   m_inparts_per_jet.Reset();
00219   m_inparts_per_jet.Set(m_njets);
00220 
00221   Int_t nj = 0;
00222   Int_t npart;
00223   m_ifewest_parts = 100000; // Starting min value
00224   for(i = 0 ; i < np ; i++ ){
00225     if (m_ipart_jet_assoc[i] != UNASSOC) continue;
00226 
00227     new((*m_jet)[nj]) TLorentzVector(*((TLorentzVector*)m_part->At(i)));
00228 
00229     npart = 1;
00230     for (j = 0 ; j < np ; j++) {
00231       if(m_ipart_jet_assoc[j] == i) {
00232         m_ipart_jet_assoc[j] = nj;
00233         npart++;
00234       }
00235     }
00236     m_ipart_jet_assoc[i] = nj;
00237     m_inparts_per_jet[nj] = npart;
00238     if( npart < m_ifewest_parts) m_ifewest_parts = npart;
00239     nj++;
00240   }  
00241   m_part->Delete();
00242 
00243 }
00244 //______________________________________________________
00245 
00246 Int_t LCDJetFinder::doFindJets(Int_t njet){
00247 
00248   Int_t np = m_4vec->GetEntries();      
00249   if (np < njet) return -1;
00250 
00251   TObject* o;
00252   Int_t ipart;
00253   for(ipart=0; ipart < np ; ipart++) {
00254     o = m_4vec->At(ipart);
00255     TString nam(o->IsA()->GetName());
00256     if (nam.Contains("TLorentzVector")) {
00257       new((*m_part)[ipart]) TLorentzVector(*((TLorentzVector*)o));
00258     } else if (nam.Contains("TVector3")) {
00259       new((*m_part)[ipart]) TLorentzVector(((TVector3 *) o)->X(),
00260                                            ((TVector3 *) o)->Y(),
00261                                            ((TVector3 *) o)->Z(),
00262                                            ((TVector3 *) o)->Mag());
00263     } else {
00264       printf("LCDJetFinder::SetEvent input is not a TVector3 or a TLorentzVector\n");
00265     }
00266   }
00267 
00268   m_ipart_jet_assoc.Reset();
00269   m_ipart_jet_assoc.Set(np);
00270   for (Int_t m=0; m<np; m++) m_ipart_jet_assoc[m] = UNASSOC;
00271 
00272   m_njets = 0;
00273 
00274   //
00275   // create invariant mass pair array.
00276   //
00277   TMatrix ymass = TMatrix(np,np);
00278   
00279   for (Int_t i1 = 0; i1 < np - 1; i1++ ) {
00280     for (Int_t i2 = i1 + 1 ; i2 < np ; i2++ ) {
00281       ymass(i1,i2) = calcinvmass(*(TLorentzVector*)m_part->At(i1),
00282                                  *(TLorentzVector*)m_part->At(i2));
00283     }
00284   }
00285   
00286   Int_t im,jm;
00287   Int_t i,j;
00288   Int_t imin;
00289   Int_t imax;
00290   Int_t nj=np;
00291   for (;;) {
00292     im = -1;
00293     jm = -1;
00294     m_ymin=1E20;
00295     //
00296     // find least invariant mass pair.
00297     //
00298     for(i = 0 ; i < np - 1 ; i++ ) {
00299       for(j = i + 1 ; j < np ; j++ ) {
00300         if (m_ipart_jet_assoc[i] != LCDJetFinder::UNASSOC) continue;
00301         if (m_ipart_jet_assoc[j] != LCDJetFinder::UNASSOC) continue;
00302         if (ymass(i,j) > m_ymin) continue;
00303           
00304         m_ymin = ymass(i,j);
00305         im = i;  jm = j;
00306       }
00307     }
00308 
00309     //
00310     // combine particles im and jm.
00311     //
00312     *((TLorentzVector*)m_part->At(im)) += *((TLorentzVector*)m_part->At(jm));
00313 
00314     for(ipart = 0; ipart < np; ipart++ ){
00315       if(m_ipart_jet_assoc[ipart] == jm) m_ipart_jet_assoc[ipart] = im;
00316     }
00317     m_ipart_jet_assoc[jm] = im;
00318     m_ymin/= m_evis;
00319     m_ymin/= m_evis;
00320     if (--nj == njet) break;
00321 
00322     //
00323     // Recalculate invariant masses for newly combined particle
00324     //
00325     for (ipart = 0; ipart < np ; ipart++) {
00326       if (ipart == im) continue;
00327       if (m_ipart_jet_assoc[ipart] != UNASSOC ) continue;
00328       
00329       imin = TMath::Min(ipart,im);
00330       imax = TMath::Max(ipart,im);
00331       
00332       ymass(imin,imax) = calcinvmass(*((TLorentzVector*)m_part->At(imin)),
00333                                      *((TLorentzVector*)m_part->At(imax)));
00334     }
00335     
00336   }
00337   //
00338   // finish up by filling jet array.
00339   //
00340   
00341   for (ipart = 0 ; ipart < np ; ipart++) {
00342     if (m_ipart_jet_assoc[ipart] == UNASSOC) m_njets++;                 
00343   }
00344 
00345   //For debug...
00346   if (m_njets != njet) {
00347     fprintf(stderr,"LCDJetFinder::doFindJets error m_njets=%d njet=%d\n",
00348             m_njets,njet);
00349     return -1;
00350   }
00351   
00352   m_jet->Delete();
00353   m_inparts_per_jet.Reset();
00354   m_inparts_per_jet.Set(m_njets);
00355 
00356   nj = 0;
00357   Int_t npart;
00358   m_ifewest_parts = 100000; // Starting min value
00359   for(i = 0 ; i < np ; i++ ){
00360     if (m_ipart_jet_assoc[i] != UNASSOC) continue;
00361 
00362     new((*m_jet)[nj]) TLorentzVector(*((TLorentzVector*)m_part->At(i)));
00363 
00364     npart = 1;
00365     for (j = 0 ; j < np ; j++) {
00366       if(m_ipart_jet_assoc[j] == i) {
00367         m_ipart_jet_assoc[j] = nj;
00368         npart++;
00369       }
00370     }
00371     m_ipart_jet_assoc[i] = nj;
00372     m_inparts_per_jet[nj] = npart;
00373     if( npart < m_ifewest_parts) m_ifewest_parts = npart;
00374     nj++;
00375   }
00376   m_part->Delete();
00377 
00378   return 0;
00379 }
00380 //______________________________________________________
00381 
00382 void LCDJetFinder::SetDURHAM(){
00383   m_algorithm = LCDJetFinder::DURHAM; 
00384 }
00385 void LCDJetFinder::SetJADE(){
00386   m_algorithm = LCDJetFinder::JADE;
00387 }
00388 void LCDJetFinder::SetJADEE(){
00389   m_algorithm = LCDJetFinder::JADEE;
00390 }
00391 
00392 Double_t LCDJetFinder::calcinvmass(const TLorentzVector& jet1,
00393                                    const TLorentzVector& jet2){
00394   TVector3 P_jet1 = jet1.Vect();
00395   TVector3 P_jet2 = jet2.Vect();
00396   Double_t costh = (P_jet1 * P_jet2)/P_jet1.Mag()/P_jet2.Mag();
00397 
00398   if     (m_algorithm == LCDJetFinder::DURHAM) {  // DURHAM
00399     Double_t minT = TMath::Min(jet1.E(),jet2.E());
00400     return 2. * minT*minT * (1.-costh);
00401   } else if (m_algorithm == LCDJetFinder::JADE)   {  // JADE    
00402     return 2. * (jet1.E())*(jet2.E()) * (1.-costh) ; 
00403   } else if (m_algorithm == LCDJetFinder::JADEE)  {  // JADE E
00404     return (jet1 + jet2).M2();
00405   }
00406 
00407   printf(" Strange Algorithm!!!! \n");
00408   return 0.;
00409 
00410 };

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