#include "hparticletool.h" #include "hphysicsconstants.h" #include "hcategorymanager.h" #include "hgeomvector.h" #include "hlinearcategory.h" //--------category definitions--------- #include "hparticledef.h" #include "haddef.h" #include "hmdcdef.h" #include "hmdctrackddef.h" #include "hmdctrackgdef.h" #include "richdef.h" #include "rpcdef.h" #include "showerdef.h" #include "tofdef.h" #include "walldef.h" //------------------------------------- //-------objects----------------------- #include "hgeantkine.h" #include "hgeanttof.h" #include "hgeantshower.h" #include "hgeantrpc.h" #include "hparticlecand.h" #include "hparticlecandsim.h" #include "hparticlepairmaker.h" #include "hmetamatch2.h" #include "hrichhit.h" #include "hrichhitsim.h" #include "htofhit.h" #include "htofhitsim.h" #include "htofcluster.h" #include "htofclustersim.h" #include "hrpccluster.h" #include "hrpcclustersim.h" #include "hshowerhit.h" #include "hshowerhitsim.h" #include "hmdctrkcand.h" #include "hmdcseg.h" #include "hmdcsegsim.h" #include "hmdchit.h" #include "hmdcclusinf.h" #include "hmdcclusfit.h" #include "hmdcwirefit.h" #include "hmdccal1.h" #include "hmdcclus.h" //------------------------------------- #include "TMath.h" #include "TAxis.h" #include "TVector3.h" //_HADES_CLASS_DESCRIPTION /////////////////////////////////////////////////////////// // HParticleTool // // library of static functions // look here for auxiliary functions for: // physics analysis, simulation analysis, ... /////////////////////////////////////////////////////////// ClassImp(HParticleTool) HParticleTool::HParticleTool() { } HParticleTool::~HParticleTool() { } Float_t HParticleTool::phiSecToLabDeg(Int_t sec, Float_t phiRad) { // Convert phi (rad) angle from sec to Lab phi in deg phiRad *= TMath::RadToDeg(); if(sec < 5) return phiRad += 60.0f * sec; else return phiRad -= 60.0f; } Float_t HParticleTool::thetaToLabDeg(Float_t thetaRad) { // Convert theta angle (rad) to the coordinate lab system return TMath::RadToDeg() * thetaRad; } Float_t HParticleTool::phiLabToPhiSecDeg(Float_t phiLabDeg) { // Convert phi angle (deg,lab: 0-360) to the coordinate sector system (60,120) return fmod(phiLabDeg,60.F) + 60; } Int_t HParticleTool::phiLabToSec(Float_t phiLabDeg) { // Convert phi angle (deg,lab: 0-360) to the sector number Int_t sec; if(phiLabDeg >= 60) sec = ((Int_t) phiLabDeg/60) - 1; else sec = 5; return sec; } Float_t HParticleTool::calcRichQA(HMdcSeg* seg, HRichHit* richHit) { // return -1 if fails. if(!seg || ! richHit) return -1; Float_t rphi = fmod(richHit->getPhi(),60.F) + 60; Float_t dTheta = richHit->getTheta() - seg->getTheta()*TMath::RadToDeg(); Float_t dPhi = ( rphi - seg->getPhi()*TMath::RadToDeg() ) * TMath::Sin(seg->getTheta()); return sqrt(dPhi*dPhi + dTheta*dTheta); } Float_t HParticleTool::calcRichQA(HMdcSeg* seg,Float_t richTheta,Float_t richPhi) { // return -1 if fails. richTheta and richPhi in deg in lab system // if(!seg ) return -1; Float_t rphi = fmod(richPhi,60.F) + 60; Float_t dTheta = richTheta - seg->getTheta()*TMath::RadToDeg(); Float_t dPhi = ( rphi - seg->getPhi()*TMath::RadToDeg() ) * TMath::Sin(seg->getTheta()); return sqrt(dPhi*dPhi + dTheta*dTheta); } Float_t HParticleTool::getOpeningAngle(Float_t phi1,Float_t theta1,Float_t phi2,Float_t theta2) { // phi and theta angles of particle 1 and 2 in lab coordinates [deg] // returns opening angle [deg] //Temporary vector objects for scalar product computation TVector3 vec1; TVector3 vec2; // convert angles to radians phi1 *= TMath::DegToRad() ; theta1 *= TMath::DegToRad() ; phi2 *= TMath::DegToRad() ; theta2 *= TMath::DegToRad() ; //Above all angles are in degree! We need them in radians! vec1.SetMagThetaPhi(1.0,theta1,phi1); vec2.SetMagThetaPhi(1.0,theta2,phi2); return TMath::RadToDeg() * vec1.Angle(vec2); } Float_t HParticleTool::getOpeningAngle(TLorentzVector& vec1,TLorentzVector& vec2) { // returns opening angle [deg] return vec1.Angle(vec2.Vect())*TMath::RadToDeg(); } Float_t HParticleTool::getOpeningAngle(HParticleCand* cand1,HParticleCand* cand2) { // returns opening angle [deg] return HParticleTool::getOpeningAngle(cand1->getPhi(),cand1->getTheta(),cand2->getPhi(),cand2->getTheta()); } Float_t HParticleTool::getOpeningAngle(HGeantKine* kine1,HGeantKine* kine2) { // returns opening angle [deg] Float_t theta1,theta2,phi1,phi2; kine1->getPhiThetaDeg(theta1,phi1); kine2->getPhiThetaDeg(theta2,phi2); return HParticleTool::getOpeningAngle(phi1,theta1,phi2,theta2); } void HParticleTool::getTLorentzVector(HGeantKine* kine, TLorentzVector& vec,Int_t pid) { // fills TLorentzVector vec from HGeantKine (if no id is provide id from // kine will be used). // CAUTION : vec.Phi() will be -pi to pi (-180 to +180 deg) // To be compatibel with the HADES lab system phi (0-360 deg) // one has to use // HParticleTool::getLabPhiDeg(TLorentzVector& vec) if(kine == NULL) { vec.SetXYZM(-1000,-1000,-1000,-1); return; } Int_t id = ( pid != -1 ) ? pid : kine->getID(); Float_t mass = HPhysicsConstants::mass(id); Float_t xmom,ymom,zmom; kine->getMomentum(xmom,ymom,zmom); vec.SetPxPyPzE(xmom,ymom,zmom,TMath::Sqrt(mass*mass + xmom*xmom + ymom*ymom + zmom*zmom)); } void HParticleTool::fillTLorentzVector(TLorentzVector& v,HParticleCand* cand,Float_t mass) { // fill v from cand. cand will be unmodified. // user has to provide the mass. if the candidate has // no reconstructed momentum the vector will be filled // by angles only (enough for opening angle calculation, but invMass // for such vectors will be wrong). FOR PARTICLES WITH CHARGE != 1 // USE fillTLorentzVector(TLorentzVector& v,HParticleCand* cand,Int_t pid,Bool_t correctMom=kTRUE) if(cand->getMomentum() != -1){ v.SetXYZM( cand->getMomentum() * TMath::Sin( TMath::DegToRad() * cand->getTheta()) * TMath::Cos( TMath::DegToRad() * cand->getPhi()), cand->getMomentum() * TMath::Sin( TMath::DegToRad() * cand->getTheta() ) * TMath::Sin( TMath::DegToRad() * cand->getPhi() ), cand->getMomentum() * TMath::Cos( TMath::DegToRad() * cand->getTheta() ), mass ); } else { TVector3 v1; v1.SetMagThetaPhi(1.0,cand->getTheta() * TMath::DegToRad(),cand->getPhi() * TMath::DegToRad()); v.SetVect(v1); } } void HParticleTool::fillTLorentzVector(TLorentzVector& v,HParticleCand* cand,Int_t pid,Bool_t correctMom) { // fill v from cand. cand will be unmodified. // user has to provide the PID. if the candidate has // no reconstructed momentum the vector will be filled // by angles only (enough for opening angle calculation, but invMass // for such vectors will be wrong). correctMom will // switch to correction of momentum (HEnergylossCorrPar needed) // The function takes care about particle charge != 1 when // fill the momenta. Float_t mass = HPhysicsConstants::mass(pid); Float_t mom = cand->getMomentumPID(pid); if(correctMom) mom = cand->getCorrectedMomentumPID(pid); if(mom != -1){ v.SetXYZM( mom * TMath::Sin( TMath::DegToRad() * cand->getTheta()) * TMath::Cos( TMath::DegToRad() * cand->getPhi()), mom * TMath::Sin( TMath::DegToRad() * cand->getTheta() ) * TMath::Sin( TMath::DegToRad() * cand->getPhi() ), mom * TMath::Cos( TMath::DegToRad() * cand->getTheta() ), mass ); } else { TVector3 v1; v1.SetMagThetaPhi(1.0,cand->getTheta() * TMath::DegToRad(),cand->getPhi() * TMath::DegToRad()); v.SetVect(v1); } } Float_t HParticleTool::getLabPhiDeg(TLorentzVector& vec) { // CAUTION : vec.Phi() will be -pi to pi (-180 to +180 deg) // To be compatibel with the HADES lab system phi (0-360 deg) // one has to use shift negative values by 360deg return vec.Phi() < 0 ? vec.Phi()*TMath::RadToDeg()+360 : vec.Phi()*TMath::RadToDeg(); } HGeomVector HParticleTool::getGlobalVertex(Int_t v,Bool_t warn) { // return the global vertex position // v can be Particle::kVertexCluster, Particle::kVertexSegment, // Particle::kVertexParticle. In case of Particle::kVertexSegment // and Particle::kVertexParticle the chi2 of the vertex fit is // checked. If the the fit Particle::kVertexSegment failed the position of // Particle::kVertexCluster is returned instead. If Particle::kVertexParticle // fails Particle::kVertexSegment is used if possible otherwise fallback // is Particle::kVertexCluster. if(!gHades) { HGeomVector vec; cout<<"HParticleTool::getGlobalVertex() : No Hades object found !"<getCurrentEvent()) { HGeomVector vec; cout<<"HParticleTool::getGlobalVertex() : No Event structure found !"<getCurrentEvent()->getHeader()->getVertexCluster().getPos(); } else if (v == Particle::kVertexSegment){ if(gHades->getCurrentEvent()->getHeader()->getVertex().getChi2() < 0){ if(warn) cout<<"HParticleTool::getGlobalVertex() : Vertex option v = "<getCurrentEvent()->getHeader()->getVertexCluster().getPos(); } return gHades->getCurrentEvent()->getHeader()->getVertex().getPos(); } else if (v == Particle::kVertexParticle ){ if(gHades->getCurrentEvent()->getHeader()->getVertexReco().getChi2() < 0){ if(gHades->getCurrentEvent()->getHeader()->getVertex().getChi2() < 0){ if(warn) cout<<"HParticleTool::getGlobalVertex() : Vertex option v = "<getCurrentEvent()->getHeader()->getVertexCluster().getPos(); } if(warn) cout<<"HParticleTool::getGlobalVertex() : Vertex option v = "<getCurrentEvent()->getHeader()->getVertex().getPos(); } return gHades->getCurrentEvent()->getHeader()->getVertexReco().getPos(); } else { cout<<"HParticleTool::getGlobalVertex() : Vertex option v = "<getCurrentEvent()->getHeader()->getVertex().getPos(); } }; //------------------------------------------------------------------------------- // Functions from GeomFunct.h (A.Schmah) void HParticleTool::calcSegVector(Double_t z, Double_t rho, Double_t phi, Double_t theta, HGeomVector &base, HGeomVector &dir) { base.setXYZ(rho*cos(phi+TMath::PiOver2()), rho*sin(phi+TMath::PiOver2()), z); dir.setXYZ(sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)); } Double_t HParticleTool::calcRMS(const Double_t* valArr, Double_t Mean,Int_t valNum) { // Calculates RMS of valNum numbers in valArr using the Mean of the values Double_t RMS = 0.0; Double_t sum = 0.0; for(Int_t i = 0; i < valNum; i++) { sum += pow(valArr[i]-Mean,2.0); } RMS = sqrt(sum/valNum); return RMS; } Double_t HParticleTool::calcDeterminant(HGeomVector& v1, HGeomVector& v2, HGeomVector& v3) { // calculating the Determinant of a 3 x 3 Matrix // with the column vectors [v1, v2, v3] // using the RULE of SARRUS // // | v1(0) v2(0) v3(0) | | v1(0) v2(0) v3(0)|v1(0) v2(0) . // | | | \\ \\ X | / / . // | | | \\ \\ / \\ | / / . // | | | \\ \\ / \\ |/ / . // | | | \\ X \\/ / . // | | | \\ / \\ /\\ / . // | | | \\ / \\ / |\\ / . // | v1(1) v2(1) v3(1) | = | v1(1) v2(1) v3(1)|v1(1) v2(1) . // | | | / \\ /\\ | /\\ . // | | | / \\ / \\ |/ \\ . // | | | / \\/ \\/ \\ . // | | | / /\\ /\\ \\ . // | | | / / \\ / |\\ \\ . // | | | / / \\/ | \\ \\ . // | v1(2) v2(2) v3(2) | | v1(2) v2(2) v3(2)| v1(2) v2(2) . // / / / \\ \\ \\ . // // - - - + + + . return ( v1(0) * v2(1) * v3(2) + v2(0) * v3(1) * v1(2) + v3(0) * v1(1) * v2(2) - v3(0) * v2(1) * v1(2) - v1(0) * v3(1) * v2(2) - v2(0) * v1(1) * v3(2)); } Double_t HParticleTool::calculateMinimumDistanceStraightToPoint(HGeomVector &base, HGeomVector &dir,HGeomVector &point) { //calculates the minimum distance of a point to a straight //given as parametric straight x = base + n * dir if (!(dir.length()>0)) { return -1000000.; } HGeomVector diff = base-point; HGeomVector cross = dir.vectorProduct(diff); return cross.length()/dir.length(); } Double_t HParticleTool::calculateMinimumDistance(HGeomVector &base1, HGeomVector &dir1, HGeomVector &base2, HGeomVector &dir2) { // calculates the minimum distance of two tracks given // as parametric straights x = base + n * dir HGeomVector cross = dir1.vectorProduct(dir2); HGeomVector ab = base1 - base2; if ( !( fabs(cross.length())>0.)) // dir1 || dir2 { return HParticleTool::calculateMinimumDistanceStraightToPoint(base1, dir1, base2); } return fabs(ab.scalarProduct(cross)/cross.length()); } HGeomVector HParticleTool::calculatePointOfClosestApproach(HGeomVector &base1, HGeomVector &dir1, HGeomVector &base2, HGeomVector &dir2) { // calculating point of closest approach // // from the equations of the straight lines of g and h // g: x1 = base1 + l * dir1 // h: x2 = base2 + m * dir2 // // you can construct the following planes: // // E1: e1 = base1 + a * dir1 + b * (dir1 x dir2) // E2: e2 = base2 + s * dir2 + t * (dir1 x dir2) // // now the intersection point of E1 with g2 = {P1} // and the intersection point of E2 with g1 = {P2} // // form the base points of the perpendicular to both straight lines. // // The point of closest approach is the middle point between P1 and P2: // // vertex = (p2 - p1)/2 // // E1 ^ g2: // // e1 = x2 // --> base1 + a * dir1 + b * (dir1 x dir2) = base2 + m * dir2 // --> base1 - base2 = m * dir2 - a * dir1 - b * (dir1 x dir2) // (m) // --> [ dir2, -dir1, -(dir1 x dir2)] (a) = base1 - base2 // (b) // // using CRAMER's RULE you can find the solution for m (a,b, not used) // // using the rules for converting determinants: // // D12 = det [dir2, -dir1, -(dir1 x dir2)] // = det [dir2, dir1, (dir1 x dir2)] // // Dm = det [base1 - base2, -dir1, -(dir1 x dir2)] // = det [base1 - base2, dir1, (dir1 x dir2)] // // m = Dm/D12 // // P1: p1 = x2(m) // = base2 + Dm/D12 * dir2 // // E2 ^ g1: // // e2 = x1 // --> base2 + s * dir2 + t * (dir1 x dir2) = base1 + l * dir1 // --> base2 - base1 = l * dir1 - s * dir2 - t * (dir1 x dir2) // (l) // --> [ dir1, -dir2, -(dir1 x dir2)] (s) = base2 - base1 // (t) // // again using CRAMER's RULE you can find the solution for m (a,b, not used) // // using the rules for converting determinants: // // D21 = det [dir1, -dir2, -(dir1 x dir2)] // = det [dir1, dir2, (dir1 x dir2)] // = -det [dir2, dir1, (dir1 x dir2)] // = -D12 // // Dl = det [base2 - base1, -dir2, -(dir1 x dir2)] // = det [base2 - base1, dir1, (dir1 x dir2)] // = -det [base1 - base2, dir1, (dir1 x dir2)] // // l = Dl/D21 // = - Dl/D12 // // P2: p2 = x1(m) // = base1 - Dl/D12 * dir1 // // // vertex = p1 + 1/2 * (p2 - p1) // = 1/2 * (p2 + p1) // = 1/2 *( (base1 + base2) + 1/D12 * ( Dm * dir2 - Dl * dir1) ) // HGeomVector cross = dir1.vectorProduct(dir2); // cross product: dir1 x dir2 // straight lines are either skew or have a cross point HGeomVector diff = base1; diff-=base2; // Difference of two base vectors base1 - base2 Double_t D; D = calcDeterminant(dir2, dir1 ,cross); if (!(fabs(D) > 0.)) { ::Warning(":calculatePointOfClosestApproach","Dirs and cross-product are lin. dependent: returning default Vertex (-20000,-20000,-20000)"); return HGeomVector(-20000.,-20000.,-20000.); } Double_t Dm = calcDeterminant(diff , dir1, cross); Double_t Dl = -calcDeterminant(diff , dir2, cross); HGeomVector vertex; HGeomVector dm; HGeomVector dl; dm = dir2; dm *= Dm; dl = dir1; dl *= Dl; vertex = dm - dl; vertex *= ((1.)/D); vertex+=base1; vertex+=base2; vertex*=0.5; return HGeomVector(vertex); } HGeomVector HParticleTool::calculateCrossPoint(HGeomVector &base1, HGeomVector &dir1, HGeomVector &base2, HGeomVector &dir2) { // calculating cross point // taking all three equations into account solving the overdetermined set of lin. equations // of // base1 + l * dir2 = base1 + m * dir2 // // set of lin. equations: // // base1(0) + l * dir1(0) = base2(0) + m * dir2(0) // base1(1) + l * dir1(1) = base2(1) + m * dir2(1) // base1(2) + l * dir1(2) = base2(2) + m * dir2(2) this line is ignored // // written in matrix form // // l \\ . // M * | | = base2 - base1 // \\ m / // // M is a 3x2 matrix // // to solve multiply the equation by the transposed Matrix of M from the left: M // // T / l \\ . // M * M * | | = M * (base2 - base1) // \\ -m / // MIND THE '-' of m // // / dir1(0) dir2(0) \\ . // | | T / dir1(0) dir1(1) dir1(2) \\ . // M = | dir1(1) dir2(1) |, M = | | // | | \\ dir2(0) dir2(1) dir2(2) / . // \\ dir1(2) dir2(2) / // // T / (dir1(0)*dir1(0) + dir1(1)*dir1(1) + dir1(2)*dir1(2)) (dir1(0)*dir2(0) + dir1(1)*dir2(1) + dir1(2)*dir2(2)) \\ . // // M * M = | | // // \\ (dir1(0)*dir2(0) + dir1(1)*dir2(1) + dir1(2)*dir2(2)) (dir2(0)*dir2(0) + dir2(1)*dir2(1) + dir2(2)*dir2(2)) / // // // T / d1d1 d1d2 \\ . // M * M = | | // \\ d1d2 d2d2 / // // diff = base2 - base1 // // T / (dir1(0)*diff(0) + dir1(1)*diff(1) + dir1(2)*diff(2)) \\ . // M * diff = | | // \\ (dir2(0)*diff(0) + dir2(1)*diff(1) + dir2(2)*diff(2)) / // // T / d1diff \\ . // M * diff = | | // \\ d2diff / // // now the new Matrix set is to be solved by CRAMER'S Rule: // // / d1d1 d1d2 \\ / l \\ / d1diff \\ . // | | * | | = | | // \\ d1d2 d2d2 / \\ -m / \\ d2diff / // // | d1d1 d1d2 | // D = | | = d1d1*d2d2 - d1d2*d1d2; // | d1d2 d2d2 | // // | d1diff d1d2 | // Dl= | | = d1diff*d2d2 - d1d2*d2diff; // | d2diff d2d2 | // // l = Dl/D = l_cross // // vertex = base1 + l_cross * dir1 // Double_t d1d1 = dir1(0)*dir1(0) + dir1(1)*dir1(1) + dir1(2)*dir1(2); Double_t d2d2 = dir2(0)*dir2(0) + dir2(1)*dir2(1) + dir2(2)*dir2(2); Double_t d1d2 = dir1(0)*dir2(0) + dir1(1)*dir2(1) + dir1(2)*dir2(2); Double_t D = d1d1*d2d2 - (d1d2*d1d2); if (!(fabs(D) > 0.)) { ::Warning("calculateCrossPoint","Error while calculating cross point ... eqns are lin. dependent:returning default Vertex (-20000,-20000,-20000)"); return HGeomVector(-20000.,-20000.,-20000.); } Double_t d1diff = dir1(0)*(base2(0)-base1(0))+dir1(1)*(base2(1)-base1(1))+dir1(2)*(base2(2)-base1(2)); Double_t d2diff = dir2(0)*(base2(0)-base1(0))+dir2(1)*(base2(1)-base1(1))+dir2(2)*(base2(2)-base1(2)); Double_t Dlambda = d1diff*d2d2-d1d2*d2diff; Double_t lambda = Dlambda/D; HGeomVector vertex; vertex += dir1; vertex *= lambda; vertex += base1; return HGeomVector(vertex); } HGeomVector HParticleTool::calcVertexAnalytical(HGeomVector &base1, HGeomVector &dir1, HGeomVector &base2, HGeomVector &dir2) { // Calculates the Vertex of two straight lines defined by the vectors base and dir // // g: x1 = base1 + l * dir1 // h: x2 = base2 + m * dir2 , where l,m are real numbers // h // 1. are g and h // parallel / identical, i.e. are dir1 and dir2 linear dependent? // // /- // | // | = 0 linear dependent, no unique solution, returning dummy // => cross product : dir1 x dir2 = -| // | != 0 linear independent // | // \\- // // 2. are g and h // skew or do they have a crossing point, i.e are dir1, dir2 and (base1 - base2) linear dependent ? // // /- // | // | = 0 linear dependent // | g and h are intersecting // | calculating vertex as point of intersection // | // => determinant: det[ dir1, dir2, base1-base2]= -| // | != 0 linear independent // | g and h are skew // | calulating vertex as point of closest approach // | // \\- // // 3. // (a) calculating intersection point // (b) calculating point of closest approach // 1. exists a unique solution ? if ((dir1.vectorProduct(dir2)).length()> 0.) // dir1 and dir2 linear independent { // straight lines are either skew or have a cross point HGeomVector diff = base1; diff-=base2; // Difference of two base vectors base1 - base2 // 2. skew or intersecting ? if (fabs(HParticleTool::calcDeterminant(dir2, dir1 ,diff))>0.) { // 3. (b) skew return HGeomVector(HParticleTool::calculatePointOfClosestApproach(base1, dir1, base2, dir2)); } else { // 3. (a) intersection return HGeomVector(HParticleTool::calculateCrossPoint(base1 ,dir1, base2 ,dir2)); } } else { // dir1 and dir2 linear dependent -> g1 and g2 identical or parallel return HGeomVector(-10000000.,-10000000.,-10000000.); } return HGeomVector(-10000000.,-10000000.,-10000000.); } // end functions from GeomFunct.h //------------------------------------------------------------------------------- Int_t HParticleTool::findFirstHitInTof(Int_t trackID,Int_t modeTrack) { //------------------------------------------------- // find the first track ID entering the TOF // Used to suppress the secondaries created in the // TOF itself. // 0 = realistic (secondaries included) // 1 primary particle is stored // 2 (default) the first track number entering the tof in SAME SECTOR is stored // 3 as 2 but condition on SAME SECTOR && MOD // 4 as 2 but SAME SECTOR && MOD && ROD Int_t numTrack = trackID; if(numTrack <= 0) return numTrack; // nothing to do for negative track numbers //------------------------------------------ // getting categories HLinearCategory* fGeantKineCat = (HLinearCategory*)HCategoryManager::getCategory(catGeantKine, kFALSE); if(!fGeantKineCat){ return trackID; } HLinearCategory* fGeantCat = (HLinearCategory*) HCategoryManager::getCategory(catTofGeantRaw,kFALSE); if(!fGeantCat){ return trackID; } //------------------------------------------ HGeantTof *poldTof; Int_t first = 0; HGeantKine* kine = 0; kine = HCategoryManager::getObject(kine,fGeantKineCat,numTrack - 1); if(kine){ first=kine->getFirstTofHit(); if(first != -1){ poldTof = (HGeantTof*)fGeantCat->getObject(first); } else { ::Error("findFirstHitInTof()","No first tof hit!"); return numTrack; } } else { ::Error("findFirstHitInTof()","Received Zero pointer for kine index = %i !",numTrack); return numTrack; } if(numTrack != poldTof->getTrack()){ ::Error("findFirstHitInTof()","First tof hit not same trackID!"); return numTrack; } //-------------------------------------------------------- // return the track number for // the selected option storeFirstTrack //-------------------------------------- // case 0 if(modeTrack == 0) return numTrack; //-------------------------------------- // case 1 if(modeTrack == 1) { // return track number of primary particle // of the given track kine = (HGeantKine*)fGeantKineCat->getObject(numTrack-1); kine = HGeantKine::getPrimary(numTrack,fGeantKineCat); return kine->getTrack(); } //-------------------------------------- // case 2 and 3 kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1); if(kine->getParentTrack() == 0){return numTrack;} // nothing to do Int_t s,m,c; s = poldTof->getSector(); // GEANT 15-22 (TOF), 23-26 (TOFINO) // numbering reverse ! m = 21-poldTof->getModule(); c = 7 -poldTof->getCell(); first = 0; Int_t tempTrack = numTrack; while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0) { first = kine->getFirstTofHit(); if(first != -1) { // track is still in TOF // now we have to check if it is in TOF or TOFINO HGeantTof* gtof = (HGeantTof*)fGeantCat->getObject(first); Int_t s1,m1,c1; s1 = m1 = c1 = 0; m1 = 21 - gtof->getModule(); if(m1 >= 0) { // inside TOF s1 = gtof->getSector(); c1 = 7-gtof->getCell(); if(modeTrack == 2 && s == s1) { // case 2 :: check only sector tempTrack = kine->getTrack(); } if(modeTrack == 3 && s == s1 && m == m1) { // case 3 :: check only sector,module tempTrack = kine->getTrack(); } else if(modeTrack == 4 && s == s1 && m == m1 && c == c1) { // case 4 :: check for same rod tempTrack = kine->getTrack(); } else { // track has no TOF hit any longer // which fulfills the condition break; } } else { // track is in TOFINO break; } } else { // track has no TOF hit any longer, // so the previous object was the one we // searched for break; } } return tempTrack; } Int_t HParticleTool::findFirstHitShowerInTofino(Int_t trackID, Int_t modeTrack) { //------------------------------------------------- // Used to suppress the secondaries created in the // SHOWER itself. // 0 = realistic (secondaries included) // 1 primary particle is stored // 2 the first track number entering the SHOWER in SAME SECTOR is stored // 3 the first track number entering the TOFINO in SAME SECTOR is stored // or the primary track if no TOFINO was found // 4 (default) the first track number entering the TOFINO in SAME SECTOR is stored // or the first track in SHOWER if no TOFINO was found Int_t numTrack = trackID; if(numTrack <= 0) return numTrack; // nothing to do for negative track numbers //-------------------------------------- // case 0 : in = out if(modeTrack == 0) { return numTrack; } HGeantShower *poldShower; Int_t first = 0; Int_t parent= 0; //------------------------------------------ // getting categories HLinearCategory* fGeantKineCat = (HLinearCategory*)HCategoryManager::getCategory(catGeantKine, kFALSE); if(!fGeantKineCat){ return trackID; } HLinearCategory* fGeantShowerCat = (HLinearCategory*) HCategoryManager::getCategory(catShowerGeantRaw,kFALSE); if(!fGeantShowerCat){ return trackID; } HLinearCategory* fGeantTofCat = (HLinearCategory*) HCategoryManager::getCategory(catTofGeantRaw,kFALSE); if(!fGeantTofCat){ return trackID; } //------------------------------------------ HGeantKine* kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1); if(kine){ parent = kine->getParentTrack(); if( parent == 0) { return numTrack; } // nothing todo first = kine->getFirstShowerHit(); if(first != -1){ poldShower = (HGeantShower*)fGeantShowerCat->getObject(first); } else { ::Error("findFirstHitShowerInTofino()","No first shower hit!"); return numTrack; } } else { ::Error("findFirstHitShowerInTofino()","Received Zero pointer for kine!"); return numTrack; } if(numTrack != poldShower->getTrack()){ ::Error("findFirstHitShowerInTofino()","First shower hit not same trackID!"); return numTrack; } //-------------------------------------------------------- // return the track number for // the selected option modeTrack Int_t s = poldShower->getSector(); //-------------------------------------- // case 1 : in -> primary if(modeTrack == 1) { // return track number of primary particle // of the given track kine = HGeantKine::getPrimary(numTrack,fGeantKineCat); return kine->getTrack(); } //-------------------------------------- // case 2 and 3 : entering SHOWER/TOF Int_t tempTrack = numTrack; first = 0; //-------------------------------------- // case 2 : entering SHOWER if(modeTrack == 2) { while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0) { first = kine->getFirstShowerHit(); if(first != -1) { // track is still in SHOWER HGeantShower* gshower = (HGeantShower*)fGeantShowerCat->getObject(first); Int_t s1 = gshower->getSector(); if(s == s1) { // check only sector tempTrack = kine->getTrack(); } else { // track has no SHOWER hit any longer // which fulfills the condition break; } } else { // track has no SHOWER hit any longer, // so the previous object was the one we // searched for break; } } return tempTrack; } //-------------------------------------- //-------------------------------------- // case 3 : entering TOFINO if(modeTrack >= 3) { Bool_t foundTof = kFALSE; Int_t tempTrack2 = tempTrack; do { first = kine->getFirstTofHit(); if(first != -1) { // we are in TOF HGeantTof* gtof = (HGeantTof*)fGeantTofCat->getObject(first); Int_t s1 = gtof->getSector(); Int_t m = gtof->getModule(); if(s == s1 && m > 21 ) { // check only sector + TOFINO foundTof = kTRUE; tempTrack = tempTrack2; } } tempTrack2 = kine->getParentTrack(); } while( tempTrack2 > 0 && (kine = (HGeantKine*)fGeantKineCat->getObject(tempTrack2 - 1)) != 0); if(foundTof) { tempTrack += 10000; } else { kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1); if( modeTrack == 3 ){ // store primaries if no TOFino was found kine = HGeantKine::getPrimary(numTrack,fGeantKineCat); tempTrack = kine->getTrack() + 100000000; } else if (modeTrack == 4){ //-------------------------------------- // recover first particle entering SHOWER while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0) { first = kine->getFirstShowerHit(); if(first != -1) { // track is still in SHOWER HGeantShower* gshower = (HGeantShower*)fGeantShowerCat->getObject(first); Int_t s1 = gshower->getSector(); if(s == s1) { // check only sector tempTrack = kine->getTrack(); } else { // track has no SHOWER hit any longer // which fulfills the condition break; } } else { // track has no SHOWER hit any longer, // so the previous object was the one we // searched for break; } } tempTrack += 100000000; //-------------------------------------- } } } //-------------------------------------- return tempTrack; } Int_t HParticleTool::findFirstHitShowerInRpc(Int_t trackID,Int_t modeTrack) { //------------------------------------------------- // Used to suppress the secondaries created in the // SHOWER itself. // 0 = realistic (secondaries included) // 1 primary particle is stored // 2 the first track number entering the SHOWER in SAME SECTOR is stored // 3 the first track number entering the RPC in SAME SECTOR is stored // or the primary track if no RPC was found // 4 (default) the first track number entering the RPC in SAME SECTOR is stored // or the first track in SHOWER if no RPC was found Int_t numTrack = trackID; if(numTrack <= 0) return numTrack; // nothing to do for negative track numbers //-------------------------------------- // case 0 : in = out if(modeTrack == 0) { return numTrack; } HGeantShower *poldShower; Int_t first = 0; Int_t parent= 0; //------------------------------------------ // getting categories HLinearCategory* fGeantKineCat = (HLinearCategory*)HCategoryManager::getCategory(catGeantKine, kFALSE); if(!fGeantKineCat){ return trackID; } HLinearCategory* fGeantShowerCat = (HLinearCategory*) HCategoryManager::getCategory(catShowerGeantRaw,kFALSE); if(!fGeantShowerCat){ return trackID; } HLinearCategory* fGeantRpcCat = (HLinearCategory*) HCategoryManager::getCategory(catRpcGeantRaw,kFALSE); if(!fGeantRpcCat){ return trackID; } //------------------------------------------ HGeantKine* kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1); if(kine){ parent = kine->getParentTrack(); if( parent == 0) { return numTrack; } // nothing todo first = kine->getFirstShowerHit(); if(first != -1){ poldShower = (HGeantShower*)fGeantShowerCat->getObject(first); } else { ::Error("findFirstHitShowerInRpc()","No first shower hit!"); return numTrack; } } else { ::Error("findFirstHitShowerInRpc()","Received Zero pointer for kine!"); return numTrack; } if(numTrack != poldShower->getTrack()){ ::Error("findFirstHitShowerInRpc()","First shower hit not same trackID!"); return numTrack; } //-------------------------------------------------------- // return the track number for // the selected option modeTrack Int_t s = poldShower->getSector(); //-------------------------------------- // case 1 : in -> primary if(modeTrack == 1) { // return track number of primary particle // of the given track kine = HGeantKine::getPrimary(numTrack,fGeantKineCat); return kine->getTrack(); } //-------------------------------------- // case 2 and 3 : entering SHOWER/RPC Int_t tempTrack = numTrack; first = 0; //-------------------------------------- // case 2 : entering SHOWER if(modeTrack == 2) { while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0) { first = kine->getFirstShowerHit(); if(first != -1) { // track is still in SHOWER HGeantShower* gshower = (HGeantShower*)fGeantShowerCat->getObject(first); Int_t s1 = gshower->getSector(); if(s == s1) { // check only sector tempTrack = kine->getTrack(); } else { // track has no SHOWER hit any longer // which fulfills the condition break; } } else { // track has no SHOWER hit any longer, // so the previous object was the one we // searched for break; } } return tempTrack; } //-------------------------------------- //-------------------------------------- // case 3 : entering RPC if(modeTrack >= 3) { Bool_t foundRPC = kFALSE; Int_t tempTrack2 = tempTrack; do { first = kine->getFirstRpcHit(); if(first != -1) { // we are in TOF HGeantRpc* grpc = (HGeantRpc*)fGeantRpcCat->getObject(first); Int_t s1 = grpc->getSector(); if(s == s1) { // check only sector + RPC foundRPC = kTRUE; tempTrack = tempTrack2; } } tempTrack2 = kine->getParentTrack(); } while( tempTrack2 > 0 && (kine = (HGeantKine*)fGeantKineCat->getObject(tempTrack2 - 1)) != 0); if(foundRPC) { tempTrack += 10000; } else { kine = (HGeantKine*)fGeantKineCat->getObject(numTrack - 1); if( modeTrack == 3 ){ // store primaries if no RPC was found kine = HGeantKine::getPrimary(numTrack,fGeantKineCat); tempTrack = kine->getTrack() + 100000000; } else if (modeTrack == 4){ //-------------------------------------- // recover first particle entering SHOWER while((kine = kine->getParent(tempTrack,fGeantKineCat)) != 0) { first = kine->getFirstShowerHit(); if(first != -1) { // track is still in SHOWER HGeantShower* gshower = (HGeantShower*)fGeantShowerCat->getObject(first); Int_t s1 = gshower->getSector(); if(s == s1) { // check only sector tempTrack = kine->getTrack(); } else { // track has no SHOWER hit any longer // which fulfills the condition break; } } else { // track has no SHOWER hit any longer, // so the previous object was the one we // searched for break; } } tempTrack += 100000000; //-------------------------------------- } } } //-------------------------------------- return tempTrack; } Float_t HParticleTool::getInterpolatedValue(TH1* h, Float_t xVal, Bool_t warn) { // retrieve content of 1-dim Hist corresponding to x val. // The values will be linear interpolated. If warn == kTRUE // warning will be emitted if the xval is out of the range // of the hist. In case the value is out of range the value // of first or last bin will be returned respectively. if(!h) return 0; TAxis *pA = h->GetXaxis(); Int_t binX = pA->FindBin(xVal); // We need to start interpolation from the next-lower bin boundary if(binX > 1 && binX <= pA->GetNbins()) { // only here it make sense to interpolate Int_t binXLower = binX - 1; Float_t startVal = h->GetBinContent(binXLower); Float_t endVal = h->GetBinContent(binX); Float_t startx = pA->GetBinUpEdge(binXLower); Float_t slope = (endVal - startVal)/pA->GetBinWidth(binX); Float_t calc = startVal + slope * (xVal - startx); return calc; } else { // noting to interpolate.... if(binX == 1) { return h->GetBinContent(1); } else if(binX < 1) { // lower boundary violation ...return lowest bin if(warn) ::Warning("HParticleTool::getInterpolatedValue()","Non-intercepted bin-out-of-range-situation (bin too low)"); return h->GetBinContent(1); } else { // upper boundary violation ...return highest bin if(warn) ::Warning("HParticleTool::getInterpolatedValue()","Non-intercepted bin-out-of-range-situation (bin too high)"); return h->GetBinContent(pA->GetNbins()); } } } Stat_t HParticleTool::getValue(TH1* h,Float_t xVal, Float_t yVal, Float_t zVal) { // retrieve value from Histogram (1,2 and 3 dim) corresponding // to x val, y val and z val. If the values are out of range of // the histogram axis the lowest/highest bin of the axis will be // used. No interpolation performed. if(!h) return 0; TAxis *pA = h->GetXaxis(); Int_t binX, binY, binZ; binY = binZ = 0; pA = h->GetXaxis(); binX = pA->FindBin(xVal); if(binX <= 0) { binX = 1; } else { if(binX > pA->GetNbins()) binX = pA->GetNbins(); } pA = h ->GetYaxis(); binY = pA->FindBin(yVal); if(binY <= 0) { binY = 1; } else { if(binY > pA->GetNbins()) binY = pA->GetNbins(); } pA = h ->GetZaxis(); binZ = pA->FindBin(zVal); if(binZ <= 0) { binZ = 1; } else { if(binZ > pA->GetNbins()) binZ = pA->GetNbins(); } return h->GetBinContent(h->GetBin(binX, binY, binZ)); } HRichHit* HParticleTool::getRichHit(Int_t richind) { // return HRichHit pointer from the object index // in the category. In case of no success returns NULL HRichHit* richhit = 0; richhit = HCategoryManager::getObject(richhit,catRichHit,richind,kTRUE); return richhit; } HTofHit* HParticleTool::getTofHit(Int_t tofind) { // return HTofHit pointer from the object index // in the category. In case of no success returns NULL HTofHit* tofhit=0; tofhit = HCategoryManager::getObject(tofhit,catTofHit,tofind,kTRUE); return tofhit; } HTofCluster* HParticleTool::getTofCluster(Int_t tofind) { // return HTofCluster pointer from the object index // in the category. In case of no success returns NULL HTofCluster* tofclst=0; tofclst = HCategoryManager::getObject(tofclst,catTofCluster,tofind,kTRUE); return tofclst; } HRpcCluster* HParticleTool::getRpcCluster(Int_t rpcind) { // return HRpcCluster pointer from the object index // in the category. In case of no success returns NULL HRpcCluster* rpcclst=0; rpcclst = HCategoryManager::getObject(rpcclst,catRpcCluster,rpcind,kTRUE); return rpcclst; } HShowerHit* HParticleTool::getShowerHit(Int_t showerind) { // return HShowerHit pointer from the object index // in the category. In case of no success returns NULL HShowerHit* showerhit=0; showerhit = HCategoryManager::getObject(showerhit,catShowerHit,showerind,kTRUE); return showerhit; } HMetaMatch2* HParticleTool::getMetaMatch(Int_t metaind) { // return HMetaMatch2 pointer from the object index // in the category. In case of no success returns NULL HMetaMatch2* meta = 0; meta = HCategoryManager::getObject(meta,catMetaMatch,metaind,kTRUE); return meta; } HMdcTrkCand* HParticleTool::getMdcTrkCand(Int_t metaind) { // return HMdcTrkCand pointer from the object index of MetaMatch2 // in the category (indMeta->HMetaMatch2->HMdcTrkCand). In case of // no success returns NULL HMetaMatch2* meta = 0; HMdcTrkCand* trk = 0; meta = HCategoryManager::getObject(meta,catMetaMatch,metaind); if(meta){ Int_t mdctrkInd = meta->getTrkCandInd(); if(mdctrkInd >= 0){ trk = HCategoryManager::getObject(trk,catMdcTrkCand,mdctrkInd,kTRUE); } } return trk; } HMdcSeg* HParticleTool::getMdcSeg(Int_t segind) { // return HMdcSeg pointer from the object index // in the category. In case of no success returns NULL HMdcSeg* seg = 0; seg = HCategoryManager::getObject(seg,catMdcSeg,segind,kTRUE); return seg; } HMdcHit* HParticleTool::getMdcHit(Int_t segind,Int_t nhit ) { // return HMdcHit pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcHit). nhit (0 or 1) // is the number of the hit inside the segment which should // be retrived. In case of no success returns NULL HMdcSeg* seg = 0; HMdcHit* hit = 0; seg = HCategoryManager::getObject(seg,catMdcSeg,segind,kTRUE); if(seg){ Int_t hitind = seg->getHitInd(nhit); if(hitind != -1) { hit = HCategoryManager::getObject(hit,catMdcHit,hitind,kTRUE); } } return hit; } HMdcClusInf* HParticleTool::getMdcClusInf(Int_t segind,Int_t nhit) { // return HMdcClusInf pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcClusInf). nhit (0 or 1) // is the number of the hit inside the segment which should // be retrived. In case of no success returns NULL HMdcSeg* seg = 0; HMdcClusInf* clusinf = 0; seg = HCategoryManager::getObject(seg,catMdcSeg,segind,kTRUE); if(seg){ Int_t hitind = seg->getHitInd(nhit); if(hitind != -1) { clusinf = HCategoryManager::getObject(clusinf,catMdcClusInf,hitind,kTRUE); } } return clusinf; } HMdcClusFit* HParticleTool::getMdcClusFit(Int_t segind) { // return HMdcClusFit pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcClsuInf->HMdcClusFit). // In case of no success returns NULL HMdcSeg* seg = 0; HMdcClusInf* clusinf = 0; HMdcClusFit* clusfit = 0; seg = HCategoryManager::getObject(seg,catMdcSeg,segind,kTRUE); if(seg){ Int_t hitind = seg->getHitInd(0); if(hitind != -1) { clusinf = HCategoryManager::getObject(clusinf,catMdcClusInf,hitind,kTRUE); } else { hitind = seg->getHitInd(1); clusinf = HCategoryManager::getObject(clusinf,catMdcClusInf,hitind,kTRUE); } if(clusinf){ Int_t clusfitind = clusinf ->getClusFitIndex(); if(clusfitind !=-1){ clusfit = HCategoryManager::getObject(clusfit,catMdcClusFit,clusfitind,kTRUE); } } } return clusfit; } HMdcClus* HParticleTool::getMdcClus(Int_t segind) { // return HMdcClus pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcClus). // In case of no success returns NULL HMdcSeg* seg = 0; HMdcClus* clus = 0; seg = HCategoryManager::getObject(seg,catMdcSeg,segind,kTRUE); if(seg){ Int_t clusind = seg->getClusIndex(); if(clusind != -1) { clus = HCategoryManager::getObject(clus,catMdcClus,clusind,kTRUE); } } return clus; } TObjArray* HParticleTool::getMdcWireFitSeg(Int_t segind) { // return TObjArray pointer to an array of HMdcWireFit // pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcClusInf->HMdcClusFit->HMdcWireFit). // The user has to delete the TObjArray object // (not the content) by himself. TObjArray* ar = new TObjArray(); HMdcClusFit* clusfit = HParticleTool::getMdcClusFit(segind); if(clusfit){ Int_t first = clusfit->getFirstWireFitInd(); Int_t last = clusfit->getLastWireFitInd(); HCategory *wirefitcat = HCategoryManager::getCategory(catMdcWireFit,2);// silent if(wirefitcat){ for(Int_t i = first; i<= last; i++){ HMdcWireFit* wf = 0; wf = HCategoryManager::getObject(wf,wirefitcat,i,kTRUE); if(wf) ar->Add(wf); } } } return ar; } Int_t HParticleTool::getMdcWireFitSeg(Int_t segind, vector& v, Bool_t clear) { // fill vector of HMdcWireFit // pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcClusInf->HMdcClusFit->HMdcWireFit). // The user has to delete the TObjArray object // (not the content) by himself. Returns the number of // collected objects .the vector will be cleared // if clear=kTRUE if(clear) v.clear(); HMdcClusFit* clusfit = HParticleTool::getMdcClusFit(segind); if(clusfit){ Int_t first = clusfit->getFirstWireFitInd(); Int_t last = clusfit->getLastWireFitInd(); HCategory *wirefitcat = HCategoryManager::getCategory(catMdcWireFit,2);// silent if(wirefitcat){ for(Int_t i = first; i<= last; i++){ HMdcWireFit* wf = 0; wf = HCategoryManager::getObject(wf,wirefitcat,i,kTRUE); v.push_back(wf); } } } return v.size(); } TObjArray* HParticleTool::getMdcCal1Seg(Int_t segind) { // return TObjArray pointer to an array of HMdcCal1 // pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcCal1). // user has to delete the TObjArray object // (not the content) by himself. TObjArray* ar = new TObjArray(); HMdcSeg* seg = HParticleTool::getMdcSeg(segind); HCategory *cal1cat = HCategoryManager::getCategory(catMdcCal1,2);// silent if(seg && cal1cat){ HLocation loccal; loccal.set(4,0,0,0,0); // t2-t1 from segment for(Int_t l = 0; l < 12; l ++) // loop over layers in segment { Int_t nCell = seg->getNCells(l); loccal[0] = seg->getSec(); Int_t io = seg->getIOSeg(); for(Int_t ic = 0; ic < nCell; ic ++) { //--------------------------------------------- // find cal1 coordinates // find module if(io == 0){ (l < 6)? loccal[1] = 0 : loccal[1] = 1; } else if(io == 1) { (l < 6)? loccal[1] = 2 : loccal[1] = 3; } // find layer (l < 6)? loccal[2] = l : loccal[2] = l - 6; // find cell number loccal[3] = seg->getCell(l,ic); //--------------------------------------------- HMdcCal1* cal1 = (HMdcCal1*)cal1cat->getObject(loccal); if(cal1) { ar->Add(cal1); } } } } return ar; } Int_t HParticleTool::getMdcCal1Seg(Int_t segind, vector& v, Bool_t clear) { // fill vector of HMdcCal1 // pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcCal1). // user has to delete the TObjArray object // (not the content) by himself. Returns number // of collected objects. the vector will be cleared // if clear=kTRUE if(clear) v. clear(); HMdcSeg* seg = HParticleTool::getMdcSeg(segind); HCategory *cal1cat = HCategoryManager::getCategory(catMdcCal1,2);// silent if(seg && cal1cat){ HLocation loccal; loccal.set(4,0,0,0,0); // t2-t1 from segment for(Int_t l = 0; l < 12; l ++) // loop over layers in segment { Int_t nCell = seg->getNCells(l); loccal[0] = seg->getSec(); Int_t io = seg->getIOSeg(); for(Int_t ic = 0; ic < nCell; ic ++) { //--------------------------------------------- // find cal1 coordinates // find module if(io == 0){ (l < 6)? loccal[1] = 0 : loccal[1] = 1; } else if(io == 1) { (l < 6)? loccal[1] = 2 : loccal[1] = 3; } // find layer (l < 6)? loccal[2] = l : loccal[2] = l - 6; // find cell number loccal[3] = seg->getCell(l,ic); //--------------------------------------------- HMdcCal1* cal1 = (HMdcCal1*)cal1cat->getObject(loccal); if(cal1) { v.push_back(cal1); } } } } return v.size(); } TObjArray* HParticleTool::getMdcCal1Cluster(Int_t segind) { // return TObjArray pointer to an array of HMdcCal1 // pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcClus->HMdcCal1). // user has to delete the TObjArray object // (not the content) by himself. TObjArray* ar = new TObjArray(); HMdcClus* clus = HParticleTool::getMdcClus(segind); HCategory *cal1cat = HCategoryManager::getCategory(catMdcCal1,2);// silent if(clus && cal1cat){ HLocation loccal; loccal.set(4,0,0,0,0); // t2-t1 from segment for(Int_t l = 0; l < 12; l ++) // loop over layers in segment { Int_t nCell = clus->getNCells(l); loccal[0] = clus->getSec(); Int_t io = clus->getIOSeg(); for(Int_t ic = 0; ic < nCell; ic ++) { //--------------------------------------------- // find cal1 coordinates // find module if(io == 0){ (l < 6)? loccal[1] = 0 : loccal[1] = 1; } else if(io == 1) { (l < 6)? loccal[1] = 2 : loccal[1] = 3; } // find layer (l < 6)? loccal[2] = l : loccal[2] = l - 6; // find cell number loccal[3] = clus->getCell(l,ic); //--------------------------------------------- HMdcCal1* cal1 = (HMdcCal1*)cal1cat->getObject(loccal); if(cal1) { ar->Add(cal1); } } } } return ar; } Int_t HParticleTool::getMdcCal1Cluster(Int_t segind, vector& v, Bool_t clear) { // fill vector of HMdcCal1 // pointer from the object index of HMdcSeg // in the category (indSeg->HMdcSeg->HMdcClus->HMdcCal1). // user has to delete the TObjArray object // (not the content) by himself. Returns the number // of collected objects. The vector will be cleared // if clear=kTRUE if(clear) v. clear(); HMdcClus* clus = HParticleTool::getMdcClus(segind); HCategory *cal1cat = HCategoryManager::getCategory(catMdcCal1,2);// silent if(clus && cal1cat){ HLocation loccal; loccal.set(4,0,0,0,0); // t2-t1 from segment for(Int_t l = 0; l < 12; l ++) // loop over layers in segment { Int_t nCell = clus->getNCells(l); loccal[0] = clus->getSec(); Int_t io = clus->getIOSeg(); for(Int_t ic = 0; ic < nCell; ic ++) { //--------------------------------------------- // find cal1 coordinates // find module if(io == 0){ (l < 6)? loccal[1] = 0 : loccal[1] = 1; } else if(io == 1) { (l < 6)? loccal[1] = 2 : loccal[1] = 3; } // find layer (l < 6)? loccal[2] = l : loccal[2] = l - 6; // find cell number loccal[3] = clus->getCell(l,ic); //--------------------------------------------- HMdcCal1* cal1 = (HMdcCal1*)cal1cat->getObject(loccal); if(cal1) { v.push_back(cal1); } } } } return v.size(); } Bool_t HParticleTool::getSimTracks(HParticleCandSim* cand, vector&tracksMeta, vector&tracksShower, vector&tracksRich, vector&weightRich, vector&tracksInnerMdc, vector&weightInnerMdc, vector&tracksOuterMdc, vector&weightOuterMdc, Bool_t print ){ tracksMeta.clear(); tracksShower.clear(); tracksRich.clear(); tracksInnerMdc.clear(); tracksOuterMdc.clear(); weightInnerMdc.clear(); weightOuterMdc.clear(); Bool_t ret =kTRUE; if(print)cout<<"SIM tracks----------------------------------------------"<getInnerSegInd()); seg [1]= (HMdcSegSim*) HParticleTool::getMdcSeg(cand->getOuterSegInd()); for(Int_t io=0;io<2;io++){ if(seg[io]==0) continue; for(Int_t tr=0;trgetNTracks(); tr++){ if(io==0){ tracksInnerMdc.push_back((Int_t)seg[io]->getTrack(tr)); weightInnerMdc.push_back((Int_t)seg[io]->getNTimes(tr)); } else { tracksOuterMdc.push_back(seg[io]->getTrack(tr)); weightOuterMdc.push_back(seg[io]->getNTimes(tr)); } } } if(print){ if(seg[0]){ cout<<"seg : "<<0<<" , ntr "<getRichInd()>-1){ HRichHitSim* richHit = (HRichHitSim*) HParticleTool::getRichHit(cand->getRichInd()); if(richHit){ if(richHit->track1>0) { tracksRich.push_back(richHit->track1); weightRich.push_back(richHit->weigTrack1); } if(richHit->track2>0) { tracksRich.push_back(richHit->track2); weightRich.push_back(richHit->weigTrack2); } if(richHit->track3>0) { tracksRich.push_back(richHit->track3); weightRich.push_back(richHit->weigTrack3); } if(print){ cout<<"Rich : "<track1 <<" ntimes "<weigTrack1<track2 <<" ntimes "<weigTrack2<track3 <<" ntimes "<weigTrack3<isTofHitUsed()) { tofhit =(HTofHitSim*) HParticleTool::getTofHit(cand->getTofHitInd()); if(tofhit){ tracksMeta.push_back(tofhit->getNTrack1()); tracksMeta.push_back(tofhit->getNTrack2()); if(print){ cout<<"TOF hit : "<getNTrack1()<<" "<getNTrack2() <isTofClstUsed()){ tofclst=(HTofClusterSim*)HParticleTool::getTofCluster(cand->getTofClstInd()); if(tofclst){ if(print){ cout<<"TOF clst : "<getNTrack1()<<" "<getNTrack2() <getNTrack1()); tracksMeta.push_back(tofclst->getNTrack2()); } else { if(print)cout<<"HParticleTool::getSimTracks() : retrieved zero pointer of tofclst!"<isRpcClstUsed()){ Int_t tracks[4]; rpcclst=(HRpcClusterSim*)HParticleTool::getRpcCluster(cand->getRpcInd()); if(rpcclst){ rpcclst->getTrackList(tracks); for(Int_t i=0;i<4;i++){ if(tracks[i]>0)tracksMeta.push_back(tracks[i]); } if(print){ cout<<"RPC clst : "<getShowerInd()>-1) { shower=(HShowerHitSim*)HParticleTool::getShowerHit(cand->getShowerInd()); if(shower){ for(Int_t i=0;igetNTracks();i++){ tracksShower.push_back(shower->getTrack(i)); if(cand->isShowerUsed()) tracksMeta.push_back(shower->getTrack(i)); } if(print) { cout<<"ShowerHit : "<getTrack(i)<<" "<getInnerSegInd()); seg [1]= (HMdcSegSim*) HParticleTool::getMdcSeg(cand->getOuterSegInd()); for(Int_t io=0;io<2;io++){ if(seg[io]==0) continue; cout<<"seg : "<getNTracks()<<" : "<getNTracks(); tr++){ cout<<" "<getTrack(tr)<<" ntimes "<getNTimes(tr)<getRichInd()>-1){ HRichHitSim* richHit = (HRichHitSim*) HParticleTool::getRichHit(cand->getRichInd()); if(richHit){ cout<<"Rich : "<track1 <<" ntimes "<weigTrack1<track2 <<" ntimes "<weigTrack2<track3 <<" ntimes "<weigTrack3<isTofHitUsed()) { tofhit =(HTofHitSim*) HParticleTool::getTofHit(cand->getTofHitInd()); if(tofhit){ cout<<"TOF hit : "<getNTrack1()<<" "<getNTrack2() <isTofClstUsed()){ tofclst=(HTofClusterSim*)HParticleTool::getTofCluster(cand->getTofClstInd()); if(tofclst){ cout<<"TOF clst : "<getNTrack1()<<" "<getNTrack2() <isRpcClstUsed()){ Int_t tracks[4]; rpcclst=(HRpcClusterSim*)HParticleTool::getRpcCluster(cand->getRpcInd()); if(rpcclst){ rpcclst->getTrackList(tracks); cout<<"RPC clst : "<getShowerInd()>-1){ shower = (HShowerHitSim*)HParticleTool::getShowerHit(cand->getShowerInd()); if(shower){ cout<<"ShowerHit : "<getNTracks();i++){ cout<getTrack(i)<<" "<getRichInd() >= 0 && cand2->getRichInd() >= 0) || (cand1->getRichInd() >= 0 && cand2->getRichInd() < 0)){ if( cand1->getRichInd() == cand2->getRichInd() ) flag = kSameRICH; else flag = kNoSameRICH; } if(cand1->getInnerSegInd() >= 0 && cand2->getInnerSegInd() >= 0){ if( cand1->getInnerSegInd() == cand2->getInnerSegInd() ) flag = flag|kSameInnerMDC; else flag = flag|kNoSameInnerMDC; } if(cand1->getOuterSegInd() >= 0 && cand2->getOuterSegInd() >= 0 ){ if( cand1->getOuterSegInd() == cand2->getOuterSegInd() ) flag = flag|kSameOuterMDC; else flag = flag|kNoSameOuterMDC; } if(cand1->getSystemUsed() != -1 && cand2->getSystemUsed() != -1 ){ if( cand1->isTofHitUsed() && cand2->isTofHitUsed()){ if( cand1->getTofHitInd() == cand2->getTofHitInd() ) flag = flag|kSameMETA; else flag = flag|kNoSameMETA; } else if(cand1->isTofClstUsed() && cand2->isTofClstUsed() ){ if(cand1->getTofClstInd() == cand2->getTofClstInd() ) flag = flag|kSameMETA; else flag = flag|kNoSameMETA; } else if(cand1->isRpcClstUsed() && cand2->isRpcClstUsed() ) { if (cand1->getRpcInd() == cand2->getRpcInd()) flag = flag|kSameMETA; else flag = flag|kNoSameMETA; } else if(cand1->isShowerUsed() && cand2->isShowerUsed()){ if(cand1->getShowerInd() == cand2->getShowerInd() ) flag = flag|kSameMETA; else flag = flag|kNoSameMETA; } else { flag = flag|kNoSameMETA;} } if(cand1->getCharge() != 0 && cand2->getCharge() != 0){ if (cand1->getCharge() > 0 && cand2->getCharge() > 0) flag = flag|kSamePosPolarity|kSamePolarity; else if(cand1->getCharge() < 0 && cand2->getCharge() < 0) flag = flag|kSameNegPolarity|kSamePolarity; else if(cand1->getCharge() < 0 && cand2->getCharge() > 0) flag = flag|kNoSamePosPolarity|kNoSameNegPolarity|kNoSamePolarity; else if(cand1->getCharge() > 0 && cand2->getCharge() < 0) flag = flag|kNoSamePosPolarity|kNoSameNegPolarity|kNoSamePolarity; } } //---------------------------------------------------------------------- if(cand2->getRichInd() >= 0) flag = flag|kRICH2; else flag = flag|kNoRICH2; if(cand2->getInnerSegmentChi2() >= 0) flag = flag|kFittedInnerMDC2; else flag = flag|kNoFittedInnerMDC2; if(cand2->getOuterSegmentChi2() >= 0) flag = flag|kFittedOuterMDC2; else flag = flag|kNoFittedOuterMDC2; if(cand2->getOuterSegInd() >= 0 ) flag = flag|kOuterMDC2; else flag = flag|kNoOuterMDC2; if(cand2->getChi2() >= 0 && cand2->getChi2() < 1e6) flag = flag|kRK2; else flag = flag|kNoRK2; if(cand2->getSystemUsed() != -1 ) flag = flag|kMETA2; else flag = flag|kNoMETA2; if(cand2->isFlagBit(kIsLepton)) flag = flag|kIsLepton2; else flag = flag|kNoIsLepton2; if(cand2->isFlagBit(kIsUsed)) flag = flag|kIsUsed2; else flag = flag|kNoIsUsed2; return kTRUE; } Bool_t HParticleTool::evalPairsFlags(UInt_t flag,UInt_t fl) { // evalutates the bitwise flags defined in flag (see eClosePairSelect + ePairCase (hparticledef.h) // against fl. return kTRUE if all set bit in flag are set in fl too. /* if( (flag&kSameRICH) == kSameRICH && (fl&kSameRICH) != kSameRICH ) return kFALSE; if( (flag&kSameInnerMDC) == kSameInnerMDC && (fl&kSameInnerMDC) != kSameInnerMDC) return kFALSE; if( (flag&kSameOuterMDC) == kSameOuterMDC && (fl&kSameOuterMDC) != kSameOuterMDC) return kFALSE; if( (flag&kSameMETA) == kSameMETA && (fl&kSameMETA) != kSameMETA ) return kFALSE; if( (flag&kSamePosPolarity) == kSamePosPolarity && (fl&kSamePosPolarity) != kSamePosPolarity ) return kFALSE; if( (flag&kSameNegPolarity) == kSameNegPolarity && (fl&kSameNegPolarity) != kSameNegPolarity ) return kFALSE; if( (flag&kSamePolarity) == kSamePolarity && (fl&kSamePolarity) != kSamePolarity ) return kFALSE; if( (flag&kRICH2) == kRICH2 && (fl&kRICH2) != kRICH2 ) return kFALSE; if( (flag&kFittedInnerMDC2) == kFittedInnerMDC2 && (fl&kFittedInnerMDC2) != kFittedInnerMDC2 ) return kFALSE; if( (flag&kFittedOuterMDC2) == kFittedOuterMDC2 && (fl&kFittedOuterMDC2) != kFittedOuterMDC2 ) return kFALSE; if( (flag&kOuterMDC2) == kOuterMDC2 && (fl&kOuterMDC2) != kOuterMDC2 ) return kFALSE; if( (flag&kRK2) == kRK2 && (fl&kRK2) != kRK2 ) return kFALSE; if( (flag&kMETA2) == kMETA2 && (fl&kMETA2) != kMETA2 ) return kFALSE; if( (flag&kIsLepton2) == kIsLepton2 && (fl&kIsLepton2) != kIsLepton2 ) return kFALSE; if( (flag&kIsUsed2) == kIsUsed2 && (fl&kIsUsed2) != kIsUsed2 ) return kFALSE; if( (flag&kNoSameRICH) == kNoSameRICH && (fl&kNoSameRICH) != kNoSameRICH ) return kFALSE; if( (flag&kNoSameInnerMDC) == kNoSameInnerMDC && (fl&kNoSameInnerMDC) != kNoSameInnerMDC) return kFALSE; if( (flag&kNoSameOuterMDC) == kNoSameOuterMDC && (fl&kNoSameOuterMDC) != kNoSameOuterMDC) return kFALSE; if( (flag&kNoSameMETA) == kNoSameMETA && (fl&kNoSameMETA) != kNoSameMETA ) return kFALSE; if( (flag&kNoSamePosPolarity) == kNoSamePosPolarity && (fl&kNoSamePosPolarity) != kNoSamePosPolarity ) return kFALSE; if( (flag&kNoSameNegPolarity) == kNoSameNegPolarity && (fl&kNoSameNegPolarity) != kNoSameNegPolarity ) return kFALSE; if( (flag&kNoSamePolarity) == kNoSamePolarity && (fl&kNoSamePolarity) != kNoSamePolarity ) return kFALSE; if( (flag&kNoRICH2) == kNoRICH2 && (fl&kNoRICH2) != kNoRICH2 ) return kFALSE; if( (flag&kNoFittedInnerMDC2) == kNoFittedInnerMDC2 && (fl&kNoFittedInnerMDC2) != kNoFittedInnerMDC2 ) return kFALSE; if( (flag&kNoFittedOuterMDC2) == kNoFittedOuterMDC2 && (fl&kNoFittedOuterMDC2) != kNoFittedOuterMDC2 ) return kFALSE; if( (flag&kNoOuterMDC2) == kNoOuterMDC2 && (fl&kNoOuterMDC2) != kNoOuterMDC2 ) return kFALSE; if( (flag&kNoRK2) == kNoRK2 && (fl&kNoRK2) != kNoRK2 ) return kFALSE; if( (flag&kNoMETA2) == kNoMETA2 && (fl&kNoMETA2) != kNoMETA2 ) return kFALSE; if( (flag&kNoIsLepton2) == kNoIsLepton2 && (fl&kNoIsLepton2) != kNoIsLepton2 ) return kFALSE; if( (flag&kNoIsUsed2) == kNoIsUsed2 && (fl&kNoIsUsed2) != kNoIsUsed2 ) return kFALSE; if( (flag&kNoUseRICH) == kNoUseRICH && (fl&kNoUseRICH) != kNoUseRICH ) return kFALSE; */ if ( (flag&fl) == flag ) return kTRUE; else return kFALSE; // return kTRUE; } Bool_t HParticleTool::isPairsFlagsBit(UInt_t flag,UInt_t fl) { // evalutates the bitwise flags defined fl (see eClosePairSelect + ePairCase (hparticledef.h)) // return kTRUE if all set bits in fl are set in flags too if ( (flag&fl) == fl ) return kTRUE; else return kFALSE; } Bool_t HParticleTool::evalPairsFlags(UInt_t flag,HParticleCand* cand1,HParticleCand* cand2) { // evalutates the bitwise flags defined in eClosePairSelect + ePairCase (hparticledef.h) // return kTRUE if all set bit are identical. cand1 as reference for Pair. if flag == 0 // result will be kTRUE if (flag == 0) return kTRUE; UInt_t fl = 0; HParticleTool::setPairFlags(fl,cand2,cand1); return evalPairsFlags(flag,fl); } Bool_t HParticleTool::evalPairsFlags(UInt_t flag,HParticlePair& pair) { // evalutates the bitwise flags defined in eClosePairSelect + ePairCase (hparticledef.h) // return kTRUE if all set bit are identical . if flag == 0 // result will be kTRUE if (flag == 0) return kTRUE; UInt_t fl = pair.getPairFlags(); return evalPairsFlags(flag,fl); return kTRUE; } Bool_t HParticleTool::evalPairsFlags(vector& flags,vector& results,HParticlePair& pair) { // evalutates the bitwise flags defined in eClosePairSelect + ePairCase (hparticledef.h) // return kTRUE if one of the flags in vector flags are identical. The result for // each individual flag is store in vector results (which will be automaticly cleared before). results.clear(); Bool_t val =kFALSE; Bool_t val2=kFALSE; for(UInt_t i = 0 ; i < flags.size(); i ++){ val = HParticleTool::evalPairsFlags(flags[i],pair); if(val) val2 = kTRUE; results.push_back(val); } return val2; } Bool_t HParticleTool::evalPairsFlags(vector& flags,HParticlePair& pair) { // evalutates the bitwise flags defined in eClosePairSelect + ePairCase (hparticledef.h) // return kTRUE if one of the flags in vector flags are identical. for(UInt_t i = 0 ; i < flags.size(); i ++){ if(HParticleTool::evalPairsFlags(flags[i],pair)) return kTRUE; } return kFALSE; }