#include "RhoEventShapes.h" #include "RhoCandidate.h" ClassImp(RhoEventShapes); RhoEventShapes::RhoEventShapes(RhoCandList &l, TLorentzVector cms) : fnChrg(0), fnNeut(0), fN(0), fpmaxlab(0.), fpmaxcms(0.), fptmax(0.), fptsumlab(0.),fneutetsumlab(0.),fneutesumlab(0.),fchrgptsumlab(0.),fchrgpsumlab(0.), fptsumcms(0.),fneutetsumcms(0.),fneutesumcms(0.),fchrgptsumcms(0.),fchrgpsumcms(0.), fsph(-1.), fapl(-1.), fpla(-1.), fFWready(false), fthr(-1.) { int i; // initialize more complex things fThrVect.SetXYZ(0.,0.,0.); fBoost=cms.BoostVector(); for (i=0;i<=FWMAX;++i) fFWmom[i]=0.; double pmax=0., ptmax=0., pmaxcms=0.; fLabList.clear(); fCmsList.clear(); for (i=0;iP4()); int chrg(cand->Charge()); fN++; // cache multiplicities if (chrg==0) fnNeut++; else fnChrg++; // cache unboosted 4-vectors fLabList.push_back(lv); // cache charges fCharge.push_back(chrg); // sum momentum variables (lab) fptsumlab += lv.Pt(); if (chrg==0) { fneutetsumlab += lv.Pt(); fneutesumlab += lv.E(); } else { fchrgptsumlab += lv.Pt(); fchrgpsumlab += lv.P(); } // cache maximum momenta in lab if (lv.P()>pmax) pmax=lv.P(); if (lv.Pt()>ptmax) ptmax=lv.Pt(); // cache boosted vectors lv.Boost(-fBoost); fCmsList.push_back(lv); // sum momentum variables (cms) fptsumcms += lv.Pt(); if (chrg==0) { fneutetsumcms += lv.Pt(); fneutesumcms += lv.E(); } else { fchrgptsumcms += lv.Pt(); fchrgpsumcms += lv.P(); } // cache maximum momenta in cms if (lv.P()>pmaxcms) pmaxcms=lv.P(); } fpmaxlab = pmax; fptmax = ptmax; fpmaxcms = pmaxcms; } // ---------------------------- void RhoEventShapes::ComputeSphericity() { if( fN==0 ) return; double stot=0, sxx=0, sxy=0, sxz=0, syy=0, syz=0, szz=0; int i; for (i=0;i0 double RhoEventShapes::Thrust() { // did we already compute? if (fthr>-1.) return fthr; TVector3 n0(0,0,0); if( fN==0 ) return -1.; int i,j; double pmax=0; //find starting vector as maximum momentum vector for (i=0;ipmax) { n0=fCmsList[i].Vect(); pmax=fCmsList[i].Vect().Mag(); } } TVector3 nnew(0,0,0); // find thrust axis (10 iterations) for (i=0;i<10;++i) { for (j=0;j1.) return -999.; if (l==0) return 1.; double pmm = 1.; double pmmp1 = x; if (l>1) { for(int ll=2; ll<=l; ll++) { double pll = (x * (2 * ll - 1) * pmmp1 - (ll - 1) * pmm) / (ll); pmm = pmmp1; pmmp1 = pll; } } return pmmp1; } // --------------------------------------- double RhoEventShapes::FoxWolfMomH(int order) { if (order>FWMAX) return -1.; // already computed FW moments if (fFWready) return fFWmom[order]; if( fN==0 ) return -1.; double s = 0.; int i,j,l; for (i=0; ipmin) cnt++; return cnt; } // --------------------------------------- int RhoEventShapes::MultPmaxLab(double pmax) { int cnt=0; for (int i=0;ipmin) cnt++; return cnt; } // --------------------------------------- int RhoEventShapes::MultPmaxCms(double pmax) { int cnt=0; for (int i=0;i