//_____________________________________________________________________________ // H1Algo algorithm to combine tracks and clusters for more precise // determination of the hadronic final state //_____________________________________________________________________________ #include "Stntuple/jet/TH1Algo.hh" #include "Stntuple/loop/TStnAna.hh" #include "Stntuple/obj/TStnEvent.hh" #include "TLorentzVector.h" #include "Stntuple/geom/TTrajectoryPoint.hh" using namespace std; ClassImp(TH1Algo) // Average energy per tower from minimum bias events. // Use to correct for multiple interactions const float multipleInteractionEnergy[52]={0, 0, 0, 0, 1.17101, 0.697716, 0.398352, 0.251982, 0.0840322, 0.0604389, 0.0465943, 0.0362292, 0.0289948, 0.023739, 0.0472224, 0.0247379, 0.0169786, 0.0186056, 0.0243106, 0.024291, 0.0241344, 0.0251232, 0.0249255, 0.0255015, 0.0253386, 0.0211731, 0.0211904, 0.0251327, 0.0251705, 0.024785, 0.0245487, 0.0236294, 0.0237527, 0.0239289, 0.0175915, 0.0155199, 0.0232503, 0.0446146, 0.0226859, 0.0273827, 0.0345807, 0.0451293, 0.0584592, 0.0818644, 0.244477, 0.386263, 0.679165, 1.13689, 0, 0, 0, 0}; //Parameters are for a simple polynomial par0+par1*x+par2*x*x+par3*x*x*x //where x=log10(pt) // //Parameters for H1 corrections cone size 0.4 in 28 eta bins const float par0h1cone0[28]={1.12904, 1.08522, 1.09614, 1.14226, 1.13969, 1.10205, 1.04654, 1.17784, 1.08757, 1.22537, 1.2402, 1.1822, 1.11348, 1.07837, 1.197, 0.934583, 0.855385, 0.704081, 0.834574, 0.871973, 0.76866, 0.866306, 0.732271, 0.801908, 0.734397, 0.855577, 0.985105, 1.06361}; const float par1h1cone0[28]={-0.238226, -0.171702, -0.187206, -0.292032, -0.29133, -0.210006, -0.122367, -0.392999, -0.231779, -0.520036, -0.608374, -0.58913, -0.517142, -0.532814, -0.837773, -0.442636, -0.343009, -0.0766073, -0.352898, -0.406942, -0.204635, -0.383296, -0.122097, -0.25181, -0.0833368, -0.317949, -0.574084, -0.782453}; const float par2h1cone0[28]={0.0694756, 0.0457471, 0.0578375, 0.12887, 0.129845, 0.0746982, 0.0301097, 0.197912, 0.0951751, 0.263543, 0.319135, 0.32922, 0.304322, 0.359081, 0.590177, 0.403336, 0.364528, 0.211403, 0.389162, 0.412154, 0.284883, 0.392434, 0.239906, 0.326253, 0.209557, 0.363342, 0.53694, 0.694098}; const float par3h1cone0[28]={0.00175895, 0.00276848, -0.000200591, -0.0146188, -0.0143435, -0.0016142, 0.00645396, -0.0257823, -0.00393989, -0.0338138, -0.0416447, -0.0449447, -0.0424114, -0.0617433, -0.114272, -0.0863601, -0.0825614, -0.0545295, -0.0915425, -0.0950123, -0.0694529, -0.0904055, -0.0622019, -0.0806304, -0.0548023, -0.086432, -0.124035, -0.159061}; //Parameters for H1 corrections cone size 0.7 in 28 eta bins const float par0h1cone1[28]={1.0697, 1.06333, 1.08457, 1.15353, 1.06602, 1.09471, 1.12413, 1.11394, 1.1333, 1.07299, 1.12905, 1.14139, 1.10066, 1.0761, 1.24302, 1.19137, 1.16548, 1.02731, 1.09106, 0.99677, 1.01333, 1.08204, 1.05767, 1.03275, 1.07314, 1.06165, 1.03222, 1.12234}; const float par1h1cone1[28]={-0.0809223, -0.0797861, -0.14536, -0.278051, -0.126696, -0.176014, -0.244573, -0.255957, -0.300614, -0.210207, -0.363534, -0.472478, -0.508617, -0.564873, -0.990927, -0.974299, -0.980912, -0.746204, -0.88189, -0.683827, -0.699305, -0.857857, -0.754937, -0.697351, -0.749071, -0.674091, -0.5942, -0.758647}; const float par2h1cone1[28]={-0.0215294, -0.0116542, 0.0434479, 0.122683, 0.04044, 0.0649521, 0.109876, 0.126096, 0.145224, 0.0825464, 0.169251, 0.251225, 0.313251, 0.391841, 0.70258, 0.733685, 0.762884, 0.635201, 0.725169, 0.598567, 0.602276, 0.729839, 0.640274, 0.608815, 0.638045, 0.562338, 0.517472, 0.638334}; const float par3h1cone1[28]={0.0163186, 0.0124725, -0.000461225, -0.0148843, 6.96668e-05, -0.00272131, -0.0111072, -0.0151572, -0.0169678, -0.00230503, -0.015726, -0.0321498, -0.0492311, -0.0708905, -0.138696, -0.151076, -0.16109, -0.139274, -0.159052, -0.134133, -0.134563, -0.167437, -0.145036, -0.140046, -0.146734, -0.125853, -0.120111, -0.151427}; //Parameters for H1 corrections cone size 1.0 in 28 eta bins //w/o silicon standalone track const float par0h1cone2[28]={1.02321, 1.06001, 1.03588, 1.07222, 1.00664, 0.993627, 1.05339, 1.08524, 1.00147, 1.04812, 1.01346, 1.07818, 1.08718, 1.18601, 1.26728, 1.32749, 1.27865, 1.11992, 1.17405, 1.05201, 1.03297, 1.07646, 1.07066, 1.03982, 0.861052, 0.991925, 0.98374, 0.998515}; const float par1h1cone2[28]={0.000340644, -0.0868955, -0.0566876, -0.139339, -0.0316727, -0.0127186, -0.151048, -0.243651, -0.105208, -0.208924, -0.190043, -0.377278, -0.494274, -0.783894, -1.05062, -1.25672, -1.21147, -0.957884, -1.07991, -0.8624, -0.808685, -0.882157, -0.829103, -0.747275, -0.332926, -0.558913, -0.523679, -0.535761}; const float par2h1cone2[28]={-0.0625178, -0.000841379, -0.00694706, 0.0483186, -0.00881702, -0.0223343, 0.0673014, 0.131517, 0.042479, 0.0918668, 0.0778581, 0.198645, 0.302757, 0.519355, 0.729896, 0.897996, 0.891989, 0.759843, 0.842868, 0.719177, 0.68009, 0.729833, 0.677354, 0.627197, 0.341388, 0.483651, 0.473072, 0.486857}; const float par3h1cone2[28]={0.0215234, 0.00858735, 0.0081102, -0.00272808, 0.00760535, 0.0115503, -0.00592988, -0.0192625, -0.000341113, -0.0069588, -0.00212012, -0.0250815, -0.0489383, -0.0960288, -0.142671, -0.181503, -0.18405, -0.162218, -0.180672, -0.158556, -0.150364, -0.162358, -0.148402, -0.139731, -0.0778546, -0.108532, -0.111564, -0.116571}; //Parameters for H1 corrections kt cone size 0.7 in 28 eta bins const float par0h1kt[28]={0.890553, 0.906432, 0.927289, 0.914774, 0.828871, 0.874047, 0.868488, 0.917702, 0.770858, 0.751653, 0.741766, 0.661683, 0.612315, 0.637236, 0.784618, 0.75515, 0.674387, 0.646504, 0.588591, 0.568239, 0.629466, 0.631187, 0.701155, 0.758762, 0.510231, 0.676013, 0.588541, 0.797555}; const float par1h1kt[28]={0.214661, 0.177329, 0.118089, 0.12774, 0.277217, 0.18389, 0.189287, 0.0736355, 0.325971, 0.334429, 0.301708, 0.354751, 0.352619, 0.196709, -0.174249, -0.184366, -0.063758, -0.0250013, 0.0699472, 0.139763, 0.0412342, 0.0535046, -0.0708959, -0.124181, 0.442159, 0.115294, 0.407385, -0.014388}; const float par2h1kt[28]={-0.172771, -0.141691, -0.0937847, -0.092704, -0.174695, -0.11913, -0.119343, -0.0459569, -0.197067, -0.208706, -0.193679, -0.203597, -0.170633, -0.0256794, 0.246301, 0.287977, 0.23257, 0.21526, 0.164602, 0.111817, 0.16558, 0.1592, 0.245624, 0.257431, -0.134129, 0.0953007, -0.145862, 0.149934}; const float par3h1kt[28]={0.0413288, 0.0336952, 0.0227743, 0.0221209, 0.0371918, 0.0278548, 0.0280528, 0.0140172, 0.0438761, 0.0483942, 0.0488429, 0.0493759, 0.0394927, 0.00356778, -0.0561517, -0.0700508, -0.0627596, -0.0609499, -0.0528145, -0.041954, -0.0520452, -0.0519195, -0.0724655, -0.071548, 0.0142716, -0.0393562, 0.0195147, -0.0490906}; //_____________________________________________________________________________ TH1Algo::TH1Algo() { //Array of locked tracks (will expand if more than 100 tracks) fLockedTrack.Set(100); fHad4Vector.SetPxPyPzE(0,0,0,0); fHad4VectorESysUp.SetPxPyPzE(0,0,0,0); fHad4VectorESysDown.SetPxPyPzE(0,0,0,0); fExtrap = new TSimpleExtrapolator(); fGeom = new TStnCalGeometry(); fUseMiniPlug=kFALSE; //Mini plug is not used bu default fUseTracks=kTRUE; //Use tracks. If Unset algorithm is calo only //Systematics fSysMultiInter=0.25; fSysTrack=0.03; //Track quality cuts fTrMinPt=0.5; //minimum Pt fTrMaxPt=15.0; //maximum Pt fTrMinCotHits=25; //minimum COT hits in central region //Swith off silicon stand alone tracks // fTrEtaCotBoundary=0.8; // Boundary between COT and Si stand alone regions fTrEtaCotBoundary=99; // Boundary between COT and Si stand alone regions fTrZ0Max=60; // maximum Z0 fTrZ0MinusZVetexMax=2.5; // maximum |Z0-zvertex| //Calo quality cuts fCaloMinEt=0.0; //minimum Et //Extra calibration on top of everything else. //If level >4 is used these should be left to default as calibration is //done after these scales are applied fExtraCaloFactor=1.12; fExtraTrackFactor=1.12; // fOldCombinationMethod=kFALSE; fMaxTrackTowerEtaDist=0.1; fMaxTrackTowerPhiDist=0.2; fNvtx=0; fZv=0; fConeSize[0]=0.4; //cone fConeSize[1]=0.7; //cone fConeSize[2]=1.0; //cone fSysUseMultiInter=kTRUE; fSysUseRel=kTRUE; fSysUseAbs=kTRUE; fSysUseTrack=kTRUE; fSysUseLevel6=kTRUE; fSysUseLevel7=kTRUE; ResetLockedTowers(); ResetLockedTracks(); fPrintLevel = 0; } //_____________________________________________________________________________ TH1Algo::~TH1Algo() { if(fGeom) delete fGeom; if(fExtrap) delete fExtrap; } //_____________________________________________________________________________ Int_t TH1Algo::MakeH1Algo(TStnEvent* event, const char* jetName, Int_t level, Float_t zvtx, Int_t nvtx) { //Make H1 algo jets and hadronic 4 vector for this event //Level specifies the level of jet corrections to be applied // if(_level==0) no corrections // if(_level>=1) relativeEnergyCorrections // if(_level>=4) multipleInteractionEnergyCorrections // if(_level>=5) absoluteEnergyCorrections // if(_level>=6) underlyingEnergyCorrections // if(_level>=7) outOfConeEnergyCorrections //level 5 is hadron level //level 7 is parton level int rc = UnpackBlocks(event,jetName); if(rc!=0) { std::cout << "Error: could not unpack event" << endl; return -1; } float rcone = fJetBlock->ConeSize(); fCone = -1; if(fabs(rcone-0.4)<0.01) fCone = 0; if(fabs(rcone-0.7)<0.01) fCone = 1; if(fabs(rcone-1.0)<0.01) fCone = 2; if(fCone<0) { std::cout << "Error: could not find cone size" << endl; return -1; } fJetType = -1; if(strstr(jetName,"JetClu") || strstr(jetName,"JetBlock")) { fJetType = fCone; } else if(strstr(jetName,"MidPoint")) { if(fCone==0) fJetType = 3; if(fCone==1) fJetType = 4; } else if(strstr(jetName,"KtCLus")) { if(fCone==1) fJetType = 5; } else { std::cout << "Error: could not parse jet algorithm" << endl; return -1; } fLevel=level; if(zvtx<-998.0) { MakeZVertexQuantities(); } else { fZv=zvtx; fNvtx=nvtx; } bool mc=fHeaderBlock->McFlag(); static const Int_t version=5; ResetArrarys() ; fJetcor=new JetEnergyCorrections(7,fNvtx,fCone,version,0,fHeaderBlock->fRunNumber,!mc); JetLoop() ; TowerLoop(); if(fPrintLevel>9) { for(Int_t i=0; i<52; i++){ for(Int_t j=0; j<48; j++){ if(fTowerTrack[i][j]>0 || fTowerJet[i][j]>0) { printf("TH1Algo::Make JetLoop Dump %2d %2d %2d %2d\n", i,j,fTowerTrack[i][j],fTowerJet[i][j]); } } } } TrackLoop(); if(fPrintLevel>9) { for(Int_t i=0; i<52; i++){ for(Int_t j=0; j<48; j++){ if(fTowerTrack[i][j]>0 || fTowerJet[i][j]>0) { printf("TH1Algo::Make TrkLoop Dump %2d %2d %2d %2d\n", i,j,fTowerTrack[i][j],fTowerJet[i][j]); } } } } MakeH1Jets(); if(fPrintLevel>0) Print(); if(fJetcor) delete fJetcor; return 0; } //_____________________________________________________________________________ void TH1Algo::MakeZVertexQuantities(){ fNvtx=0; fZv=0; float maxpt=0; for(Int_t ivtx=0; ivtxNVertices(); ivtx++) { TStnVertex &vert = *(fZVertexBlock->Vertex(ivtx)); if (vert.VClass()>=12) { fNvtx++; if (vert.SumPt()>maxpt) { maxpt=vert.SumPt(); fZv=vert.Z(); } } } if (fNvtx<1) fNvtx=1; } //_____________________________________________________________________________ void TH1Algo::JetLoop() { //loop over jets if(fJetBlock) fNjets=fJetBlock->NJets(); else fNjets=0; if(fNjets>500) fNjets=500; for(Int_t ijet=0; ijetJet(ijet); Double_t jDeta = jet->DetEta(); Double_t jPhi = jet->Momentum()->Phi(); Double_t jPt = jet->Momentum()->Pt(); TStnLinkBlock* towlink = fJetBlock->TowerLinkList(); Int_t nl = towlink->NLinks(ijet); //#of related towers // Float_t multie=0; if(fPrintLevel>=2) printf("TH1Algo Jet %2d has %3d towers, Pt=%7.2f\n", ijet,nl,jPt); for (Int_t j=0; jIndex(ijet,j); //word containing ie & ip Int_t iphi = (icl & 0x3f); Int_t ieta = ((icl>>8) & 0x3f); if(fPrintLevel>=5) printf("TH1Algo Jet %2d iphi=%3d ieta=%3d\n",ijet,iphi,ieta); //If we want any info from the towers the link is in the comment below // TCalTower* t = fCalDataBlock->Tower(ieta,iphi); if(!fUseMiniPlug&&(ieta<=3||ieta>=48)) continue; //Do not look in Miniplug if(fPrintLevel>9) { printf("TH1Algo::JetLoop:Assoc %2d %2d %2d %2d\n", ijet,ieta,iphi,fLockedTower[ieta][iphi]); } if(!fLockedTower[ieta][iphi]) { fTowerJet[ieta][iphi] = ijet+2; // Flag shows to which jet energy is associated fTowerjetup[ieta][iphi] = fTowerJet[ieta][iphi] ; fTowerjetdo[ieta][iphi] = fTowerJet[ieta][iphi] ; fTowerjettrup[ieta][iphi] = fTowerJet[ieta][iphi] ; fTowerjettrdo[ieta][iphi] = fTowerJet[ieta][iphi] ; //Store the Et of the jet (used for systematic errors) fPtJet[ieta][iphi]=jPt; } } } } //_____________________________________________________________________________ void TH1Algo::TowerLoop() { // Collect Calorimeter energies Int_t ntowers = fCalDataBlock->NTowers(); for (Int_t i=0; iTower(i); Int_t ieta = to->IEta(); Int_t iphi = to->IPhi(); if(!fUseMiniPlug&&(ieta<=3||ieta>=48)) continue; //Do not look in Miniplug if(fLockedTower[ieta][iphi]) continue; if(to->Energy()<=0.0001) continue; Double_t thprime = fGeom->Theta(ieta); if(to->Energy()*sin(thprime)EmEnergy(); Double_t Ehad=to->HadEnergy(); Double_t Etot=(Ehad+Eem); //Upward and downward systematic error Double_t EtotUp=Etot; Double_t EtotDo=Etot; //Here we should correct for number of extra interactions if(fLevel>=4&&fNvtx>1) { if(fPtJet[ieta][iphi]>1){ Float_t ptafter=fPtJet[ieta][iphi]; fJetcor->multipleInteractionEnergyCorrections(ptafter); Float_t cor=ptafter/fPtJet[ieta][iphi]; Etot *= cor; EtotUp *=cor+MISYSV5*fSysUseMultiInter*(cor-1); EtotDo *=cor-MISYSV5*fSysUseMultiInter*(cor-1); } } if(Etot<0) Etot=0; if(EtotDo<0) EtotDo=0; if(EtotUp<0) continue; //This means all Etot and EtotDo is also <0 so we //shouldn't continue //Here we perform relative (eta dependent) corrections on the calo cluster Float_t pt=25.0; if(fLevel>=1) { Float_t eta= fGeom->Eta(ieta); fJetcor->relativeEnergyCorrections(pt,eta); Double_t corrfactor=pt/25.0; //correction factor Etot *= corrfactor; EtotUp *=corrfactor; EtotDo *=corrfactor; //Get relative corrections error add in quad to running total shift Float_t sysrel=fJetcor->sysRelativeCor(eta,pt); if(fSysUseRel) { EtotUp =Etot+sqrt(pow(EtotUp-Etot,2)+pow(sysrel*Etot,2)); EtotDo =Etot-sqrt(pow(Etot-EtotDo,2)+pow(sysrel*Etot,2)); } } if(fLevel>=5) { //Get absolute corrections error (at jet scale) if(fSysUseAbs&&fPtJet[ieta][iphi]>1){ Float_t sysabs= fJetcor->Cal2HadError(fPtJet[ieta][iphi]); EtotUp =Etot+sqrt(pow(EtotUp-Etot,2)+pow(sysabs*Etot,2));; EtotDo =Etot-sqrt(pow(Etot-EtotDo,2)+pow(sysabs*Etot,2)); } } fETot[ieta][iphi] = Etot; fETottrdo[ieta][iphi] = Etot; // track systematic fETottrup[ieta][iphi] = Etot; // track systematic fETotup[ieta][iphi] = EtotUp; //calo systematic fETotdo[ieta][iphi] = EtotDo; //calo systematic if(fPrintLevel>=5) printf("TH1Algo Tower E's %2d %2d %7.4f %7.4f %7.4f \n", ieta,iphi,Etot,EtotUp,EtotDo); if(fTowerJet[ieta][iphi]==0) { //If tower is not in a jet flag to show there is energy in this tower fTowerJet[ieta][iphi]=1; fTowerjetup[ieta][iphi] = fTowerJet[ieta][iphi]; fTowerjetdo[ieta][iphi] = fTowerJet[ieta][iphi]; fTowerjettrup[ieta][iphi] = fTowerJet[ieta][iphi]; fTowerjettrdo[ieta][iphi] = fTowerJet[ieta][iphi]; } } } //_____________________________________________________________________________ void TH1Algo::TrackLoop() { // Collect tracks and perform matching Int_t ntr; // We ignore tracks if fUseTracks is unset if(fUseTracks) ntr = fTrackBlock->NTracks(); else ntr=0; for (Int_t i=0; iTrack(i); if(!TrackSelection(tr)) continue; //select quality tracks TTrajectoryPoint* pt=GetPointAtCalo(tr); //extrapolate to calo // Get eta and phi at calo Double_t r=sqrt(pow(pt->X(),2)+pow(pt->Y(),2)+pow(pt->Z(),2)); Double_t exetatrack=0; if(r>0) { Double_t exthetatrack=acos(pt->Z()/r); exetatrack=-log(tan(exthetatrack/2)); } else { std::cout <<"TH1Algo Something very wrong with track extrapolation- Track not taken"<Y(),pt->X()); //Get tower indices Int_t ieta=fGeom->IEta(pt->X(),pt->Y(),pt->Z()); Int_t iphi = fGeom->IPhi(ieta,exphitrack); Double_t Etrack= sqrt(pow(tr->Momentum()->P(),2)+(0.14*0.14)); Double_t EtrackUp =Etrack*(1+fSysTrack*fSysUseTrack) ; Double_t EtrackDo =Etrack*(1-fSysTrack*fSysUseTrack) ; //Get sorted list of nearest towers with energy>0 NearestTowers(exetatrack, exphitrack, ieta, iphi); if(fPrintLevel>=7) printf("Track extrap %7.4f %7.4f %3d %3d\n", exetatrack, exphitrack, ieta, iphi); Double_t EKilled=0; //Energy killed in cal for(Int_t itower=0; itowerfETot[itoweta][itowphi]) { //Here the tower energy is less than that left to kill so //we kill the tower fTowerJet[itoweta][itowphi]=0; //kill tower EKilled += fETot[itoweta][itowphi]; } else { //Here the tower energy is more than that left to kill so //we rescale the cluster fETot[itoweta][itowphi] -= Etrack-EKilled; break; // we now have subtracted all the energy of the track // so we shouldn't continue } } //Now systematics //CaloUp EKilled=0; //Energy killed in cal for(Int_t itower=0; itowerfETotup[itoweta][itowphi]) { fTowerjetup[itoweta][itowphi]=0; //kill tower EKilled += fETotup[itoweta][itowphi]; } else { fETotup[itoweta][itowphi] -= Etrack-EKilled; break; } } //Calo Down EKilled=0; //Energy killed in cal for(Int_t itower=0; itowerfETotdo[itoweta][itowphi]) { fTowerjetdo[itoweta][itowphi]=0; //kill tower EKilled += fETotdo[itoweta][itowphi]; } else { fETotdo[itoweta][itowphi] -= Etrack-EKilled; break; } } //Track Up EKilled=0; //Energy killed in cal for(Int_t itower=0; itowerfETottrup[itoweta][itowphi]) { fTowerjettrup[itoweta][itowphi]=0; //kill tower EKilled += fETottrup[itoweta][itowphi]; } else { fETottrup[itoweta][itowphi] -= EtrackUp-EKilled; break; } } //Track Do EKilled=0; //Energy killed in cal for(Int_t itower=0; itowerfETottrdo[itoweta][itowphi]) { fTowerjettrdo[itoweta][itowphi]=0; //kill tower EKilled += fETottrdo[itoweta][itowphi]; } else { fETottrdo[itoweta][itowphi] -= EtrackDo-EKilled; break; } } fNtra[ieta][iphi]++; //Here we associate the tracks to jets if(!fTowerTrack[ieta][iphi]) fTowerTrack[ieta][iphi]=2+NearJet(tr->Phi0(),tr->Eta()); //Store the total pt fTrpt[ieta][iphi] += tr->Pt(); fTrpx[ieta][iphi] += tr->Momentum()->Px(); fTrpy[ieta][iphi] += tr->Momentum()->Py(); fTrpz[ieta][iphi] += tr->Momentum()->Pz(); fTre[ieta][iphi] += Etrack; if(fPrintLevel>9) { printf("TH1Algo::TrackLoop adding track %7.4f %7.4f\n", fTrpt[ieta][iphi],fTre[ieta][iphi]); } } } //_____________________________________________________________________________ void TH1Algo::MakeH1Jets(){ //====================================== // Here we make the new jet 4 vector // In the first run through energies not in a jet are taken Double_t trackEtot=0; //total track energy summed over all hadrons Double_t caloEtot=0; //total calo energy summed over all hadrons (uncalibrated) if(fPrintLevel>9) { for(Int_t i=0; i<52; i++){ for(Int_t j=0; j<48; j++){ if(fTowerTrack[i][j]>0 || fTowerJet[i][j]>0) { printf("TH1Algo::Make Dump %2d %2d %2d %2d\n", i,j,fTowerTrack[i][j],fTowerJet[i][j]); } } } } for(Int_t ijet=-1; ijet=5) printf("TH1Algo:MakeJets starting %2d\n",ijet); Double_t jcorpx=0; Double_t jcorpy=0; Double_t jcorpz=0; Double_t jcorE=0; Double_t trackE=0; Double_t jcorEtrup=0; Double_t jcorEtrdo=0; Double_t jcorEup=0; Double_t jcorEdo=0; Int_t ijetplus2=ijet+2; Bool_t TakeTrack=kFALSE; Bool_t TakeCalo=kFALSE; if(ijet==-1) fNCombTracks=0; if(ijet==-1) fNCombTowers=0; for(Int_t i=0; i<52; i++){ for(Int_t j=0; j<48; j++){ TakeTrack=kFALSE; TakeCalo=kFALSE; if(ijet==-1) { //This is energy that is not in jets if(fTowerTrack[i][j]==1) TakeTrack=kTRUE; if(fTowerJet[i][j]==1) TakeCalo=kTRUE; // else continue; } else { //Make decision for each jet if(fTowerTrack[i][j]==ijetplus2) TakeTrack=kTRUE; if(fTowerJet[i][j]==ijetplus2) TakeCalo=kTRUE; // else continue; } if(TakeTrack){ jcorpx += fTrpx[i][j]*fExtraTrackFactor; jcorpy += fTrpy[i][j]*fExtraTrackFactor; jcorpz += fTrpz[i][j]*fExtraTrackFactor; jcorE += fTre[i][j]*fExtraTrackFactor; // TakeCalo=kFALSE; trackE +=fTre[i][j]*fExtraTrackFactor; if(ijet==-1) fNCombTracks++; //Systematics if(ijet>-1) { jcorEtrup += fTre[i][j]*fExtraTrackFactor*(1+fSysTrack*fSysUseTrack); jcorEtrdo += fTre[i][j]*fExtraTrackFactor*(1-fSysTrack*fSysUseTrack); jcorEup += fTre[i][j]*fExtraTrackFactor; jcorEdo += fTre[i][j]*fExtraTrackFactor; } } if(TakeCalo) { Double_t thprime = fGeom->Theta(i); Double_t thcentr = CorrectedTheta(thprime,fZv); //Zv correted Double_t phcentr = fGeom->Phi(i,j); jcorpx += fExtraCaloFactor*fETot[i][j]*sin(thcentr)*cos(phcentr); jcorpy += fExtraCaloFactor*fETot[i][j]*sin(thcentr)*sin(phcentr); jcorpz += fExtraCaloFactor*fETot[i][j]*cos(thcentr); jcorE += fExtraCaloFactor*fETot[i][j]; // TakeTrack=kFALSE; fNCombTowers++; } if(fPrintLevel>=5 && (TakeCalo||TakeTrack)) { printf("TH1Algo:MakeJets cal %1d tr %1d ie %2d ip %2d E %7.4f\n", TakeCalo,TakeTrack,i,j,jcorE); } //Systematics if(fTowerjettrup[i][j]==ijetplus2) jcorEtrup += fExtraCaloFactor*fETottrup[i][j]; if(fTowerjettrdo[i][j]==ijetplus2) jcorEtrdo += fExtraCaloFactor*fETottrdo[i][j]; if(fTowerjetup[i][j]==ijetplus2) jcorEup += fExtraCaloFactor*fETotup[i][j]; if(fTowerjetdo[i][j]==ijetplus2) jcorEdo += fExtraCaloFactor*fETotdo[i][j]; } } //Now we form the total shift which is the individual contributions added //in quadrature jcorEup= jcorE+sqrt(pow(jcorEup-jcorE,2)+pow(jcorEtrup-jcorE,2)); jcorEdo= jcorE-sqrt(pow(jcorE-jcorEdo,2)+pow(jcorE-jcorEtrdo,2)); if(ijet>=0) { fJetESysUp[ijet] =jcorEup; fJetESysDown[ijet] =jcorEdo; } Double_t trackfrac; if(jcorE>0) trackfrac=trackE/jcorE; else trackfrac=1; trackEtot +=trackE; caloEtot +=jcorE; if(ijet==-1) { //These are total hadronic sums for event fHad4Vector.SetPxPyPzE(jcorpx,jcorpy,jcorpz,jcorE); fHad4VectorESysUp.SetPxPyPzE(jcorpx,jcorpy,jcorpz,jcorE); fHad4VectorESysDown.SetPxPyPzE(jcorpx,jcorpy,jcorpz,jcorE); if(jcorE>0) { fHad4VectorESysUp *=jcorEup/jcorE; fHad4VectorESysDown *=jcorEdo/jcorE; } //Needs updating fTrackFraction=trackfrac; } else { //New method just scale energy and keep angles the same //Get origianl jet energy TStnJet* jet = fJetBlock->Jet(ijet); Float_t jcorEo = jet->Momentum()->E(); fJet4Vector[ijet]=*jet->Momentum(); fJet4Vector[ijet] *=jcorE/jcorEo; //Take recalculated jet 4 vector // fJet4Vector[ijet].SetPxPyPzE(jcorpx,jcorpy,jcorpz,jcorE); //Perform absolute corrections if(fLevel>=5&&jcorE>0) { Float_t h1abscor=AbsoluteEnergyCorrections(fJet4Vector[ijet].Pt(), fJet4Vector[ijet].Eta(),fJetType); fJet4Vector[ijet] *= h1abscor; fJetESysUp[ijet]*=h1abscor; fJetESysDown[ijet]*=h1abscor; } //Now corrections to take us to parton level if(fLevel>=6) Level6Corrections(ijet); if(fLevel>=7) Level7Corrections(ijet); fHad4Vector +=fJet4Vector[ijet]; fHad4VectorESysUp +=fJet4Vector[ijet]*(fJetESysUp[ijet]/fJet4Vector[ijet].E()); fHad4VectorESysDown +=fJet4Vector[ijet]*(fJetESysDown[ijet]/fJet4Vector[ijet].E()); fH1trackfraction[ijet]=trackfrac; } } fTrackFraction=trackEtot/caloEtot; } //_____________________________________________________________________________ void TH1Algo::Level6Corrections(Int_t ijet){ Float_t ptunc=fJet4Vector[ijet].Pt(); if(ptunc<0.1) return; //underlying event Float_t eunc=fJet4Vector[ijet].E(); Float_t pt=ptunc; fJetcor->underlyingEnergyCorrections(pt); //Shift up fJetcor->setSysUnderlyingEventCorrection(+1); Float_t ptup=fJet4Vector[ijet].Pt(); fJetcor->underlyingEnergyCorrections(ptup); //Shift down fJetcor->setSysUnderlyingEventCorrection(-1); Float_t ptdo=fJet4Vector[ijet].Pt(); fJetcor->underlyingEnergyCorrections(ptdo); fJetcor->setSysUnderlyingEventCorrection(0); if(fSysUseLevel6) { fJetESysUp[ijet] =eunc+sqrt(pow(fJetESysUp[ijet]-eunc,2)+pow(eunc*(ptup/pt-1),2)); fJetESysDown[ijet] =eunc-sqrt(pow(eunc-fJetESysDown[ijet],2)+pow(eunc*(1-ptdo/pt),2)); } fJet4Vector[ijet] *=pt/ptunc; fJetESysUp[ijet] *=pt/ptunc; fJetESysDown[ijet] *=pt/ptunc; } //_____________________________________________________________________________ void TH1Algo::Level7Corrections(Int_t ijet){ Float_t ptunc=fJet4Vector[ijet].Pt(); if(ptunc<0.1) return; Float_t eunc=fJet4Vector[ijet].E(); //Get out of cone corrections error //Default out of cone corrections Float_t pt=ptunc; fJetcor->outOfConeEnergyCorrections(pt); //Shift up fJetcor->setSysOutOfConeCorrection(+1); Float_t ptup=fJet4Vector[ijet].Pt(); fJetcor->outOfConeEnergyCorrections(ptup); //Shift down fJetcor->setSysOutOfConeCorrection(-1); Float_t ptdo=fJet4Vector[ijet].Pt(); fJetcor->outOfConeEnergyCorrections(ptdo); fJetcor->setSysOutOfConeCorrection(0); if(fSysUseLevel7) { fJetESysUp[ijet] =eunc+sqrt(pow(fJetESysUp[ijet]-eunc,2)+pow(eunc*(ptup/pt-1),2)); fJetESysDown[ijet] =eunc-sqrt(pow(eunc-fJetESysDown[ijet],2)+pow(eunc*(1-ptdo/pt),2)); } fJet4Vector[ijet] *=pt/ptunc; fJetESysUp[ijet] *=pt/ptunc; fJetESysDown[ijet] *=pt/ptunc; } //_____________________________________________________________________________ float TH1Algo::AbsoluteEnergyCorrections(const double pt, const double eta, int cone) { //Perform absolute corrections for H1 Algo jets to bring //them close to the hadron level ie level 5 corrections float lptuse=log10(pt); if(lptuse<0.85) lptuse=0.85; //Get the eta bins below and above the true eta int ietam1; int ietap1; if(fabs(eta)<0.05) { //first bin we can't average ietam1=0; ietap1=0; } else{ ietam1=(int)((fabs(eta)-0.05)*10); ietap1=ietam1+1; } float etam1=ietam1/10.0+0.05; float etap1=ietap1/10.0+0.05; //Prevent going outside of range of validity of calibrated kinematic region. if(ietam1>27) ietam1=27; if(ietap1>27) ietap1=27; if(ietap1>20&&lptuse>1.95) lptuse=1.95; else if(ietap1>15&&lptuse>2.3) lptuse=2.3; else if(ietap1>10&&lptuse>2.45) lptuse=2.45; else if(ietap1>5&&lptuse>2.75) lptuse=2.75; else if(lptuse>2.8) lptuse=2.8; float corm1=-1; float corp1=-1; static Bool_t err=kTRUE; if(fJetType==0) { corm1=par0h1cone0[ietam1]+par1h1cone0[ietam1]*lptuse+par2h1cone0[ietam1]*pow(lptuse,2)+par3h1cone0[ietam1]*pow(lptuse,3); corp1=par0h1cone0[ietap1]+par1h1cone0[ietap1]*lptuse+par2h1cone0[ietap1]*pow(lptuse,2)+par3h1cone0[ietap1]*pow(lptuse,3); } else if(fJetType==1) { corm1=par0h1cone1[ietam1]+par1h1cone1[ietam1]*lptuse+par2h1cone1[ietam1]*pow(lptuse,2)+par3h1cone1[ietam1]*pow(lptuse,3); corp1=par0h1cone1[ietap1]+par1h1cone1[ietap1]*lptuse+par2h1cone1[ietap1]*pow(lptuse,2)+par3h1cone1[ietap1]*pow(lptuse,3); } else if(fJetType==2) { corm1=par0h1cone2[ietam1]+par1h1cone2[ietam1]*lptuse+par2h1cone2[ietam1]*pow(lptuse,2)+par3h1cone2[ietam1]*pow(lptuse,3); corp1=par0h1cone2[ietap1]+par1h1cone2[ietap1]*lptuse+par2h1cone2[ietap1]*pow(lptuse,2)+par3h1cone2[ietap1]*pow(lptuse,3); } else if(fJetType==5) { corm1=par0h1kt[ietam1]+par1h1kt[ietam1]*lptuse+par2h1kt[ietam1]*pow(lptuse,2)+par3h1kt[ietam1]*pow(lptuse,3); corp1=par0h1kt[ietap1]+par1h1kt[ietap1]*lptuse+par2h1kt[ietap1]*pow(lptuse,2)+par3h1kt[ietap1]*pow(lptuse,3); } else { if(err) { std::cout <<"TH1Algo::AbsoluteCorrections: We don't yet have corrections implemented for JetType " < TMath::Pi()/2.) my_Z_PES = -1.* Z_PES; //Double_t r=tan(thetauncor)*Z_PES; Double_t r=tan(thetauncor)*my_Z_PES; return fabs(atan2(r,my_Z_PES)); } } //_____________________________________________________________________________ void TH1Algo::NearestTowers(Float_t etatrack, Float_t phitrack,Int_t ietatrack,Int_t iphitrack){ //Here we loop over the Towers and // see which are nearest in eta-phi space // Inputs should be eta and phi coords at calo after extrapolating // to calo without Z-vertex corrections // Output is sorted list of nearest towers fTowerDistance, //and eta, phi coords fNearTowersIeta, fNearTowersIphi static const Int_t etatol=4; //+- # towers that will be searched for energy // static const Int_t phitol=4; // in eta-phi grid fOldCombinationMethod=kFALSE; fNtow=0; if(fOldCombinationMethod){ //Here we only take the nearest tower in eta/phi and ignore all //possible energy deposits in other towers // Float_t delphi=fabs(fGeom->Phi(ietatrack,iphitrack)-phitrack); //if(delphi > M_PI) delphi = 2*M_PI-delphi; Float_t delphi=acos(cos(fGeom->Phi(ietatrack,iphitrack)-phitrack)); Float_t Dist2=pow(fGeom->Eta(ietatrack)-etatrack,2)+pow(delphi,2); //square of distance of track to tower in eta-phi space fTowerDistance[fNtow]=sqrt(Dist2); fNearTowersIeta[fNtow]=ietatrack; fNearTowersIphi[fNtow]=iphitrack; fNtow=1; } else{ for(Int_t ieta=ietatrack-etatol; ieta<=ietatrack+etatol ; ieta++) { if(ieta<0) continue; //Protect against out of bounds error if(ieta>51) continue; for(Int_t iphi=0; iphi<48 ; iphi++){ if (fTowerJet[ieta][iphi]<=0) continue; Float_t deleta=fabs(fGeom->Eta(ieta)-etatrack); // Float_t delphi=fabs(fGeom->Phi(ieta,iphi)-phitrack); Float_t delphi=acos(cos(fGeom->Phi(ieta,iphi)-phitrack)); if(fPrintLevel>9) { printf("TH1Algo nearest tow %2d %2d %7.4f %7.4f %7.4f %7.4f \n", ieta,iphi,fGeom->Eta(ieta),etatrack, fGeom->Phi(ieta,iphi),phitrack); } if((deletaEta(ieta)-etatrack,2)+pow(delphi,2); //square of distance of track to tower in eta-phi space fTowerDistance[fNtow]=sqrt(Dist2); fNearTowersIeta[fNtow]=ieta; fNearTowersIphi[fNtow]=iphi; fNtow++; } } } } //Sort the towers in increasing distance TMath::Sort(fNtow,fTowerDistance,fNearTowersIndex,kFALSE); if(fPrintLevel>=5) { printf("TH1Algo Found %3d nearest towers\n",fNtow); if(fPrintLevel>=7) { for(int k=0; kJet(ijet); Double_t jeta = jet->Momentum()->Eta(); Double_t jPhi = jet->Momentum()->Phi(); // Double_t dphi = fabs(PhiTrack - jPhi); // if(dphi > TMath::Pi()) dphi = 2*TMath::Pi()-dphi; Double_t dphi=acos(cos(PhiTrack - jPhi)); Double_t deta = fabs(EtaTrack - jeta); Double_t drr = sqrt(dphi*dphi + deta*deta); if(drr < fConeSize[fCone]) return ijet; } return -1; } //_____________________________________________________________________________ Bool_t TH1Algo::TrackSelection( TStnTrack* tr){ // return kFALSE; if(fPrintLevel>7) printf("TH1Algo Track Pt=%7.4f Eta=%7.4f NCotH=%3d Z=%7.4f\n", tr->Pt(),tr->Eta(),tr->NCotHitsTot(),tr->Z0()); if( tr->Pt() < fTrMinPt) return kFALSE; if( tr->Pt() > fTrMaxPt) return kFALSE; if( tr->Eta()<-9 ) return kFALSE; //rejects silly tracks if(tr->NCotHitsTot() < fTrMinCotHits && fabs(tr->Eta())Eta())>0.8) return kFALSE; if(fabs(tr->Z0()) > fTrZ0Max ) return kFALSE; if(fabs(tr->Z0()-fZv)> fTrZ0MinusZVetexMax ) return kFALSE; if(fPrintLevel>7) printf("TH1Algo Track selected\n"); return kTRUE; } //_____________________________________________________________________________ TTrajectoryPoint* TH1Algo::GetPointAtCalo(TStnTrack* track) { //Swims track to calorimeter and returns trajectory at that //point static TTrajectoryPoint p0; double xyz[8]; xyz[0] = -track->D0() * TMath::Sin(track->Phi0()); xyz[1] = track->D0() * TMath::Cos(track->Phi0()); xyz[2] = track->Z0(); xyz[3] = track->Momentum()->Px()/track->Momentum()->P(); xyz[4] = track->Momentum()->Py()/track->Momentum()->P(); xyz[5] = track->Momentum()->Pz()/track->Momentum()->P(); xyz[6] = 0; xyz[7] = track->Momentum()->P(); p0.SetPoint(xyz); Int_t intcharge=0; if(track->Charge()>0) intcharge=1; else if(track->Charge()<0) intcharge=-1; int status = fExtrap->SwimToCal(&p0,intcharge); //if(status!=0) return -1; return fExtrap->GetTrajectoryPoint(); } //_____________________________________________________________________________ void TH1Algo::LockTowersInJet(Int_t ijet) { //Exclude towers corresponding to jet with index ijet from FScomb algorithm //Must be same jets as defined in constructor TStnLinkBlock* towlink = fJetBlock->TowerLinkList(); Int_t nl = towlink->NLinks(ijet); //#of related towers for (Int_t j=0; jIndex(ijet,j); //word containing ie & ip Int_t iphi = (icl & 0x3f); Int_t ieta = ((icl>>8) & 0x3f); LockTower(ieta, iphi); } } //_____________________________________________________________________________ void TH1Algo::LockTrack(TStnTrack* track) { // Tracks locked in similar way to towers (slower than index method) // Tracks must be pointers to tracks in TStnTrackBlock. Copies will not work // Int_t ntr = fTrackBlock->NTracks(); for (Int_t i=0; iTrack(i); if(tracktest==track) { LockTrack(i); return; } } std::cout <<"TH1Algo::LockTrack: Track not found in TStnTrackBlock"<GetCurrentTreeEntry(); fHeaderBlock = (TStnHeaderBlock*) event->GetDataBlock("HeaderBlock"); fCalDataBlock = (TCalDataBlock*) event->GetDataBlock("CalDataBlock"); if(!fCalDataBlock) { printf("TH1Algo: Error CalDataBlock not registered\n"); return 1; } fCalDataBlock->GetEntry(ientry); fTrackBlock = (TStnTrackBlock*) event->GetDataBlock("TrackBlock"); if(!fTrackBlock) { printf("TH1Algo: Error TrackBlock not registered\n"); return 1; } fTrackBlock->GetEntry(ientry); fZVertexBlock = (TStnVertexBlock*) event->GetDataBlock("ZVertexBlock"); if(!fZVertexBlock) { printf("TH1Algo: Error ZVertexBlock not registered\n"); return 1; } fZVertexBlock->GetEntry(ientry); fJetBlock = (TStnJetBlock*) event->GetDataBlock(jetName); if(!fJetBlock) { printf("TH1Algo: Error JetBlock not registered\n"); return 1; } fJetBlock->GetEntry(ientry); return 0; } //_____________________________________________________________________________ void TH1Algo::Print() { if(!fJetBlock) return; printf(" i Et H1Et Eta H1Eta EtUp EtDo TrFrac\n"); for(Int_t ijet=0; ijet< fJetBlock->NJets(); ijet++) { TStnJet* jet = fJetBlock->Jet(ijet); TLorentzVector* ojet= jet->Momentum(); Float_t jetEta = ojet->PseudoRapidity(); Float_t jetEt = ojet->E()*sin(ojet->Theta()); TLorentzVector* h1jet= GetJet4Vector(ijet); Float_t njetEta = h1jet->PseudoRapidity(); Float_t njetEt = h1jet->E()*sin(h1jet->Theta()); //Jet pt after upward systematic shift Float_t jetEtUp =jetEt*GetJetSysUp(ijet); //Jet pt after downward systematic shift Float_t jetEtDo =jetEt*GetJetSysDown(ijet); Float_t trackfraction= GetTrackFraction(ijet); printf("%2i %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f \n",ijet, jetEt,njetEt,jetEta,njetEta,jetEtUp,jetEtDo,trackfraction); } }