/////////////////////////////////////////////////////////////////////////////// // Initialization of STNTUPLE electron block // P.Murat /////////////////////////////////////////////////////////////////////////////// #include "TROOT.h" #include "TFolder.h" #include "AbsEnv/AbsEnv.hh" #include "TrackingObjects/Storable/CdfTrackColl.hh" #include "TrackingObjects/Storable/CdfTrackView.hh" #include "ShowerMax/CesAlignCorr.hh" #include "ElectronObjects/CdfEmObject.hh" #include "Electron/emobj_alg.hh" #include "Electron/CemCorrAlg.hh" #include "Electron/PemCorrAlg.hh" #include "Electron/IsoCorrAlg.hh" #include "ElectronObjects/CdfEmObjectColl.hh" #include "ElectronObjects/CdfEmObjectView.hh" #include "ElectronObjects/SimpleExtrapolatedTrack.hh" // For Pes RDeltaPhi #include "Calor/calor_alg.hh" // and DeltaR #include "ElectronUser/StdCEMElectron.hh" #include "ElectronUser/StdPEMElectron.hh" #include "HighLevelObjects/PhotonVariables.hh" #include "HighLevelObjects/TrackingVariables.hh" #include #include #include #include #include #include #include #include using namespace PhotonVariables; namespace { TFolder* fgFolder = NULL; TStnErrorLogger* fgLogger = NULL; CT_ReFit* fgRefit = NULL; StntupleCotReFit* fgScrf = NULL; } //_____________________________________________________________________________ Int_t StntupleInitElectron(TStnElectron* ele, const CdfEmObject* bob) { // initialize STNTUPLE electron starting from CDF RunII Bob's electron CdfTrack* bobs_track; CesAlignCorr cesalign; double leak, cemcorr; float pemcorr=1.; StdCEMElectron* std_bob=NULL; StdPEMElectron* pem_bob=NULL; // char* report; AbsEnv* env = AbsEnv::instance(); //----------------------------------------------------------------------------- // backward pointer to CDF electron is used for internal bookkeeping only // during the STNTUPLE filling stage //----------------------------------------------------------------------------- ele->fCdfEmObject = (CdfEmObject*) bob; ele->fTrack = NULL; CdfEmObject::EmCluster_link cluster = bob->matchingEmCluster(); //----------------------------------------------------------------------------- // pack detector code //----------------------------------------------------------------------------- int barrel = cluster->side(); int wedge = cluster->seedPhiAddr(); ele->fDetCode = ((int) cluster->region() & 0xff); ele->fDetCode |= (cluster->seedEtaAddr () & 0xff) << 8; ele->fDetCode |= (wedge & 0xff) << 16; ele->fDetCode |= (barrel & 0xff) << 24; //----------------------------------------------------------------------------- // Store cluster size and set track number to "undefined" //----------------------------------------------------------------------------- ele->fTrind = ((cluster->nTowers() & 0xff) << 16)+0xffff; ele->fEtcor = 1; ele->fIdwrd = 0; ele->fIswrd = 0; ele->fConwrd = 0; //EMOBJECT 4-momentum (changed from EMCLUSTER 4-momentum 02-14-02) // Valid choices for string are // DEFAULT --- Energy from calorimeter // angles from track, shower max, or calorimeter centroid // tried in that order // CALONLY --- Energy and angles from calorimeter. Uses current z vertex // TRKONLY --- Energy = track total momentum. Angles from track. // SHWRMAX --- Energy from calorimeter. Angles from each view of shower // max. Uses current z vertex. // Why? Since fourMomentum() makes a new object. // A new object should be returned by value // HepLorentzVector why_it_is_returned_by_value = bob->fourMomentum(); // use 2-tower sum to be consistent with the ETF //... Bob, Bob, who forced you to return 4-vectors of doubles by value? HepLorentzVector why_it_is_returned_by_value = bob->fourMomentum(); HepLorentzVector* mom = &why_it_is_returned_by_value; //----------------------------------------------------------------------------- // fill TStnElectron's momentum: 1. default - 4-momentum //----------------------------------------------------------------------------- TLorentzVector* ele_mom = (TLorentzVector*) ele->Momentum(); ele_mom->SetX(mom->px()); ele_mom->SetY(mom->py()); ele_mom->SetZ(mom->pz()); ele_mom->SetT(mom->e()); //----------------------------------------------------------------------------- // try to improve accuracy in momentum for the central electrons //----------------------------------------------------------------------------- if (cluster->region() == 0) { why_it_is_returned_by_value = emobj_alg::CEMFourMomentum::TwoTower(bob,true); if (mom->e() == 0) { // no track, why_it_is_returned_by_value = emobj_alg::CEMFourMomentum::TwoTower(bob,false); } if (mom->e() != 0) { // one of the attempts succeeded ele_mom->SetX(mom->px()); ele_mom->SetY(mom->py()); ele_mom->SetZ(mom->pz()); ele_mom->SetT(mom->e()); } ele->fBcE = ele_mom->E(); //----------------------------------------------------------------------------- // store RAW CEM face correction in fEtcor (as name allows), mark those // electrons, for which correction failed // corrected momentum should be multiplied by fEtcor //----------------------------------------------------------------------------- if (CemCorrAlg::GetCemCorr(*bob,"USETRACK",cemcorr) == CemCorrAlg::OK) { ele->fEtcor = 1./cemcorr; } else { ele->fEtcor = -1.; } } //----------------------------------------------------------------------------- // common part again //----------------------------------------------------------------------------- ele->fEt = cluster->totalEt(); // BEWARE! This is cluster emEt+hadEt ele->fDteta = cluster->emEtaDetector(); ele->fEveta = cluster->emEtaEvent(); ele->fPhi = cluster->emPhi(); ele->fZv = cluster->vertexZ(); // Z vertex assigned to EM cluster; ele->fHadem = cluster->hadEm(); ele->fHadEmT = cluster->hadEmTriggerTower(); ele->fIso = cluster->totalIsolationEt4(); ele->fIso7 = cluster->totalIsolationEt7(); ele->fEmIso4 = cluster->emIsolationEt4(); ele->fHadIso4 = cluster->hadIsolationEt4(); ele->fLshr = bob->lshr(); ele->fLshr2 = bob->lshr(2); //----------------------------------------------------------------------------- // TRACK QUANTITIES //----------------------------------------------------------------------------- ele->fBcPt = TStnTypes::kUndefined; ele->fBcZ0 = TStnTypes::kUndefined; ele->fBcPhi0 = TStnTypes::kUndefined; ele->fBcLam0 = TStnTypes::kUndefined; ele->fCharge = TStnTypes::kUndefined; ele->fPt = TStnTypes::kUndefined; ele->fD0 = TStnTypes::kUndefined; ele->fZ0 = TStnTypes::kUndefined; ele->fEp = TStnTypes::kUndefined; ele->fXtrk = TStnTypes::kUndefined; ele->fZtrk = TStnTypes::kUndefined; ele->fTiso = TStnTypes::kUndefined; ele->fNTracks = bob->numberTracks(); bobs_track = (CdfTrack*) bob->maxPtTrack().operator->(); if ((! bobs_track) && (ele->fNTracks != 0)) { //----------------------------------------------------------------------------- // we are in deep trouble, force fNTracks to be 0 //----------------------------------------------------------------------------- // report = Form("%s with %i tracks but without maxPtTrack ", // "StntupleInitElectronBlock: Bobs electron",ele->fNTracks); // fgLogger->Report(-102,report); ele->fNTracks = 0; } //----------------------------------------------------------------------------- // These quantities must be undefined in case there is no matching cluster //----------------------------------------------------------------------------- ele->fDelX = TStnTypes::kUndefined; ele->fBcDelX = TStnTypes::kUndefined; ele->fXces = TStnTypes::kUndefined; ele->fChiw = TStnTypes::kUndefined; ele->fDelZ = TStnTypes::kUndefined; ele->fBcDelZ = TStnTypes::kUndefined; ele->fZces = TStnTypes::kUndefined; ele->fChis = TStnTypes::kUndefined; ele->fCesEnergy[0] = -1; ele->fCesEnergy[1] = -1; ele->fNCesClusters[0] = 0; ele->fNCesClusters[1] = 0; ele->fNCprClusters[0] = 0; ele->fNCprClusters[1] = 0; ele->fCprEnergy = 0; ele->fPprEnergy = 0.; ele->fPem3x3Eta = -8.; ele->fPem3x3Phi = -8.; ele->fChi3 = TStnTypes::kUndefined; ele->fPem3x3FitTower = -1; ele->fNPesClusters = 0; ele->fPesEta = TStnTypes::kUndefined; ele->fPesPhi = TStnTypes::kUndefined; ele->fPesEnergy[0] = -1; ele->fPesEnergy[1] = -1; ele->fPes5x9 [0] = TStnTypes::kUndefined; ele->fPes5x9 [1] = TStnTypes::kUndefined; ele->fPesXYZ [0] = TStnTypes::kUndefined; ele->fPesXYZ [1] = TStnTypes::kUndefined; ele->fPesXYZ [2] = TStnTypes::kUndefined; ele->fPlugHadEmCut = -1; ele->fPesRDeltaPhi = TStnTypes::kUndefined; ele->fPesDeltaR = TStnTypes::kUndefined; ele->fPesPemDeltaR = TStnTypes::kUndefined; ele->fPem3x3E = TStnTypes::kUndefined; if (ele->IsCentral()) { //----------------------------------------------------------------------------- // central electron //----------------------------------------------------------------------------- if (bobs_track) { std_bob = new StdCEMElectron(bob,*fgRefit); HepLorentzVector std_mob = std_bob->fourMomentum(); ele->fBcE = std_mob.e(); Hep3Vector bobs_mom = bobs_track->momentum(); ele->fCharge = bobs_track->charge(); ele->fPt = bobs_track->pt(); ele->fD0 = bobs_track->d0(); ele->fZ0 = bobs_track->z0(); // Bob says E/P=EM_E/P, not TOT_E/P ele->fEp = std_bob->eOverP(); //----------------------------------------------------------------------------- // These are X- and Z-coordinates of the track at CES, supposedly rotated into // the "local" coordinate of the wedge, but, as far as I understand, not // corrected for the CES alignment // COT-only beam-constrained track from std_bob is used to calculate residuals //----------------------------------------------------------------------------- Hep3Vector bobs_coord; bobs_coord = bob->maxPtTrackLocalCoord(); ele->fXtrk = bobs_coord.x(); ele->fZtrk = bobs_coord.z(); // Ray Culbertson 02-18-02 ele->fTiso = TrackingVariables::trkTrackIsolation(*bobs_track); } ele->fNCesClusters[0] = 0; ele->fNCesClusters[1] = 0; const IndexViewVector& ces_list=bob->matchingCesClusters(); vector< ValueVectorLink >::const_iterator cesCluster; for (cesCluster = ces_list.contents().begin(); cesCluster != ces_list.contents().end (); ++cesCluster){ if ((*cesCluster)->view() == CesCluster::WIRE ) ele->fNCesClusters[1]++; if ((*cesCluster)->view() == CesCluster::STRIP) ele->fNCesClusters[0]++; } ValueVectorLink strip_cluster; ValueVectorLink wire_cluster; wire_cluster = bob->bestMatchingCesCluster(CesCluster::WIRE); if (wire_cluster.is_nonnull()) { //----------------------------------------------------------------------------- // for some reason accounting for misalignments even for MC is necessary //----------------------------------------------------------------------------- ele->fXces = wire_cluster->fitted_z_strip(); ele->fChiw = wire_cluster->chi2(); ele->fCesEnergy[1] = wire_cluster->raw_energy(); if (bobs_track) { ele->fDelX = cesalign.xtrack_corr(ele->fXtrk,barrel,wedge)-ele->fXces; ele->fBcDelX = std_bob->qDeltaX()*ele->Charge(); } } strip_cluster = bob->bestMatchingCesCluster(CesCluster::STRIP); if (strip_cluster.is_nonnull()) { //----------------------------------------------------------------------------- // coordinate in the chamber ref. syst. //----------------------------------------------------------------------------- ele->fZces = strip_cluster->fitted_z_strip(); ele->fChis = strip_cluster->chi2(); ele->fCesEnergy[0] = strip_cluster->raw_energy(); if (bobs_track) { //----------------------------------------------------------------------------- // looks like alignment corrections assume track coordinates to be transformed // into the local wedge system, direction of Z-axis for which depends on the // barrel number... //----------------------------------------------------------------------------- if (barrel == 0) { ele->fDelZ = cesalign.ztrack_corr(-ele->fZtrk,barrel,wedge)-ele->fZces; } else { ele->fDelZ = cesalign.ztrack_corr(ele->fZtrk,barrel,wedge)-ele->fZces; } ele->fBcDelZ = std_bob->deltaZ(); } } //----------------------------------------------------------------------------- // here goes the CPR: N(clusters) and energy of the highest one //----------------------------------------------------------------------------- const CdfEmObject::Cpr_CollType& cpr_coll = bob->matchingCprClusters(); vector< ValueVectorLink >::const_iterator cprCluster; ele->fNCprClusters[0] = cpr_coll.contents().size(); for (cprCluster = cpr_coll.contents().begin(); cprCluster != cpr_coll.contents().end (); ++cprCluster) { if (ele->fCprEnergy < (*cprCluster)->energy()) { ele->fCprEnergy = (*cprCluster)->energy(); } } CdfEmObject::Cpr2EList elist = bob->cpr2PadEnergies(); CdfEmObject::Cpr2PadNoList ilist = bob->cpr2PadNumbers(); int ncpr2 = elist.size(); if(ilist[0]==0 && ilist[1]==0 && ilist[2]==0 && ilist[3]==0) ncpr2=0; ele->fCpr2I = (ncpr2<<24); ele->fCpr2E0 = ele->fCpr2E1 = 0; if(ncpr2>0) { ele->fCpr2E0 |= (elist[0]&0xFFFF); ele->fCpr2I |= (ilist[0]); } if(ncpr2>1) { ele->fCpr2E0 |= (elist[1]<<16); ele->fCpr2I |= (ilist[1]<<6); } if(ncpr2>2) { ele->fCpr2E1 |= (elist[2]&0xFFFF); ele->fCpr2I |= (ilist[2]<<12); } if(ncpr2>3) { ele->fCpr2E1 |= (elist[3]<<16); ele->fCpr2I |= (ilist[3]<<18); } //----------------------------------------------------------------------------- // central electron fiducial //----------------------------------------------------------------------------- int fid_smx = fidEle(*bob,"SHWRMAX"); int fid_trk = fidEle(*bob,"TRACK"); ele->fFid = ((fid_smx & 0xff) << 8) | (fid_trk & 0xff); } else if (ele->IsPlug()) { //----------------------------------------------------------------------------- // PLUG electron //----------------------------------------------------------------------------- pem_bob = new StdPEMElectron(bob); bobs_track = (CdfTrack*) pem_bob->theTrack(); HepLorentzVector pem_mob = pem_bob->fourMomentum(); ele->fBcE = pem_mob.e(); if (PemCorrAlg::GetPemCorr(*bob,pemcorr) == PemCorrAlg::OK) { ele->fEtcor = 1./pemcorr; } else { ele->fEtcor = -1.; } if (bobs_track) { ele->fCharge = pem_bob->charge(); ele->fPt = pem_bob->pt(); ele->fD0 = bobs_track->d0(); ele->fZ0 = pem_bob->trackZ0(); ele->fEp = pem_bob->eOverP(); //----------------------------------------------------------------------------- // These are X- and Z-coordinates of the track at CES, supposedly rotated into // the "local" coordinate of the wedge, but, as far as I understand, not // corrected for the CES alignment //----------------------------------------------------------------------------- Hep3Vector bobs_coord; bobs_coord = bob->maxPtTrackLocalCoord(); ele->fXtrk = bobs_coord.x(); ele->fZtrk = bobs_coord.z(); // Ray Culbertson 02-18-02 ele->fTiso = TrackingVariables::trkTrackIsolation(*bobs_track); } ele->fPem3x3FitTower = cluster->pem3x3FitTowers(); ele->fChi3 = cluster->pem3x3Chisq(); //Ray Culbertson 02-18-02 ele->fPem3x3Eta = cluster->pem3x3FitDetEta(); ele->fPem3x3Phi = cluster->pem3x3FitPhi(); ele->fPem3x3E = cluster->pem3x3EmEnergy(); ele->fPlugHadEmCut = pem_bob->hadEmCut(); // Alon Attal 3-24-04 ele->fNPesClusters = bob->matchingPes2dClusters().size(); ValueVectorLink pes_cluster; pes_cluster = bob->bestMatchingPes2dCluster(); if (pes_cluster.is_nonnull()) { ele->fPesEta = pes_cluster->plug2dClusterEta(); ele->fPesPhi = pes_cluster->plug2dClusterPhi(); //----------------------------------------------------------------------------- // 0: U-strips, 1: V-strips //----------------------------------------------------------------------------- ele->fPesEnergy[0] = pes_cluster->plug2dUenergy(); ele->fPesEnergy[1] = pes_cluster->plug2dVenergy(); ele->fPes5x9[0] = pes_cluster->pesCluster0().pesProfileRatio5by9(); ele->fPes5x9[1] = pes_cluster->pesCluster1().pesProfileRatio5by9(); ele->fPesXYZ[0] = pes_cluster->plug2dClusterX(); ele->fPesXYZ[1] = pes_cluster->plug2dClusterY(); ele->fPesXYZ[2] = pes_cluster->plug2dClusterZ(); ele->fPesPemDeltaR = pem_bob->pesPemDeltaR(); if (bobs_track) { // Determine RDeltaPhi and DeltaR between PES and max Pt track. // These variables are the most relevant due to the cylindrical // nature of the plug. Alon Attal 3-24-04 SimpleExtrapolatedTrack exttrk; calor_alg::ExtrapolateToShowerMax::extrapolate(bobs_track, exttrk); float trackX, trackY, trackR, trackPhi, pesX, pesY, pesR, pesPhi, avgR; trackX = exttrk.currentSpacePoint().x(); trackY = exttrk.currentSpacePoint().y(); trackR = sqrt(trackX*trackX + trackY*trackY); trackPhi = atan2(trackY,trackX); trackPhi = trackPhi < 0 ? trackPhi += 2*M_PI : trackPhi; pesX = pes_cluster->plug2dClusterX(); pesY = pes_cluster->plug2dClusterY(); pesR = sqrt(pesX*pesX + pesY*pesY); pesPhi = atan2(pesY,pesX); pesPhi = pesPhi < 0 ? pesPhi += 2*M_PI : pesPhi; avgR = .5 * (pesR + trackR); ele->fPesRDeltaPhi = avgR * (pesPhi - trackPhi); ele->fPesDeltaR = pesR - trackR; } } ele->fPprEnergy = pem_bob->pprEnergy(); //----------------------------------------------------------------------------- // plug electron fiducial //----------------------------------------------------------------------------- ele->fFid = (pem_bob->fiducial() & 0xff); } ele->fNasl = -1; ele->fNssl = -1; if (bobs_track) { //----------------------------------------------------------------------------- // loop over the COT SL's and figure the number of segments; // SL 0,2,4,6 are stereo, the rest - axial //----------------------------------------------------------------------------- ele->fNasl = 0; ele->fNssl = 0; for (int sl=0; sl<8; sl++) { if (bobs_track->numCTHits(sl) >= TStnElectron::kMinCotHits[sl]) { if (sl%2 == 0) ele->fNssl++; else ele->fNasl++; } } } ele->fTime = phoTime(*bob); //----------------------------------------------------------------------------- // relative correction for the leakage //----------------------------------------------------------------------------- leak = 0; IsoCorrAlg::GetIsoCorr(*bob,leak); ele->fLeakCorr = leak/cluster->emEt(); if (std_bob) delete std_bob; if (pem_bob) delete pem_bob; return 0; } //_____________________________________________________________________________ Int_t StntupleInitElectronBlock(TStnDataBlock* block, AbsEvent* event, int mode){ int ev_number, rn_number; char process[100], description[100], coll_type[100]; int inum[100]; char* report; ev_number = AbsEnv::instance()->trigNumber(); rn_number = AbsEnv::instance()->runNumber(); if (block->Initialized(ev_number,rn_number)) return 0; TStnElectronBlock* data = (TStnElectronBlock*) block; data->Clear(); if (! fgFolder) { //----------------------------------------------------------------------------- // initialize local static variables //----------------------------------------------------------------------------- fgFolder = (TFolder*) gROOT->GetRootFolder()->FindObject("StntupleFolder"); fgLogger = (TStnErrorLogger* ) fgFolder->FindObject("StntupleErrorLogger"); fgScrf = (StntupleCotReFit*) fgFolder->FindObject("CotReFit"); fgRefit = fgScrf->GetReFit(); } //----------------------------------------------------------------------------- // Find event info by collection name ... //----------------------------------------------------------------------------- StntupleGetProcessName(data->CollName()->Data(), process, description, coll_type); EventRecord::ConstIterator ci; int nele; const CdfEmObject *e1, *e2; vector elist; if (strcmp(coll_type,"VIEW") == 0) { ci = StntupleGetIterator(event,"CdfEmObjectView",process,description); if (! ci.is_valid()) return -2; CdfEmObjectView::const_handle handle(ci); const CdfEmObjectView::collection_type& bobs_list = handle->contents(); nele = bobs_list.size(); for (int i=0; i(); elist.push_back((CdfEmObject*) e); } } else { ci = StntupleGetIterator(event,"CdfEmObjectColl",process,description); if (! ci.is_valid()) return -1; CdfEmObjectColl::const_handle handle(ci); const vector& bobs_list = handle->contents(); nele = bobs_list.size(); for (int i=0; i 100) { //----------------------------------------------------------------------------- // protect case of too many EM objects //----------------------------------------------------------------------------- report = Form("StntupleInitElectronBlock: too many EM objects %i", nele); fgLogger->Report(-101,report); nele = 100; } // order electrons in Et for (int i=0; igetEmCluster()->totalEt() < e2->getEmCluster()->totalEt()) { int itmp = inum[i1]; inum[i1] = inum[i2]; inum[i2] = itmp; e1 = e2; } } } for (int i=0; iNewElectron(); StntupleInitElectron(ele,bob); } //----------------------------------------------------------------------------- // mark event and run as initialized //----------------------------------------------------------------------------- data->f_RunNumber = rn_number; data->f_EventNumber = ev_number; return 0; } //_____________________________________________________________________________ Int_t StntupleElectronBlockLinks(TStnDataBlock* block,AbsEvent* event,int mode) { // set links between the electrons and the objects in the other blocks // returns: 0 if everything is OK // -1 if block itself has not been initialized int ev_number, rn_number; TStnTrackBlock* track_block; TStnElectronBlock* ele_block = (TStnElectronBlock*) block; TStnElectron* ele; TStnTrack* trk; const CdfEmObject* bob; const CdfTrack* bobs_track; Int_t nele, it, rc; // assume that the block exists, it // is supposed to be initialized ev_number = AbsEnv::instance()->trigNumber(); rn_number = AbsEnv::instance()->runNumber(); if (! block->Initialized(ev_number,rn_number)) return -1; // do not do initialize links 2nd time if (block->LinksInitialized()) return 0; //----------------------------------------------------------------------------- // pointers to the default track block and to the list of HPT leptons //----------------------------------------------------------------------------- TStnEvent* ev = ele_block->GetEvent(); track_block = (TStnTrackBlock*) ev->GetDataBlock("TrackBlock"); TObjArray* hptl = ev->GetListOfHptl(); rc = 0; nele = ele_block->NElectrons(); for (int iele=0; ieleElectron(iele); bob = ele->GetCdfEmObject(); bobs_track = bob->maxPtTrack().operator->(); if (bobs_track) { it = track_block->TrackNumber(bobs_track); if (it != -1) { trk = track_block->Track(it); ele->SetTrackNumber(it); ele->SetTIso(trk->Iso4()); ele->SetBcPt(fabs(TStntuple::TrackPt(trk->BcC0()))); //----------------------------------------------------------------------------- // define coordinates of the seed track extrapolated to CES //----------------------------------------------------------------------------- // ele->SetGenpNumber( track_block->Track(it)->GenpNumber() ); // ele->SetObspNumber( track_block->Track(it)->ObspNumber() ); } } //----------------------------------------------------------------------------- // fill tower link block //----------------------------------------------------------------------------- const std::vector& towlist = bob->getEmCluster()->towerIds(); int ntowers = towlist.size(); for (int itower=0; itowerTowerLinkList()->Add(iele,tower_key); } //----------------------------------------------------------------------------- // update list of good HPT leptons to calculate MET (if any...) - // we do have isolation at this point // hptl may not exist... //----------------------------------------------------------------------------- if (hptl && TStntuple::GoodHptElectron(ele)) { hptl->Add(ele); } } //----------------------------------------------------------------------------- // mark links as initialized //----------------------------------------------------------------------------- ele_block->fLinksInitialized = 1; return rc; }