#include "evt/Event.hh" #include #include "AbsEnv/AbsEnv.hh" #include "TrackingObjects/Storable/CdfTrackColl.hh" #include "TrackingObjects/Storable/CdfTrackView.hh" #include "ShowerMax/CesAlignCorr.hh" #include #include "Electron/emobj_alg.hh" #include "Electron/CemCorrAlg.hh" #include "Electron/PemCorrAlg.hh" #include "Electron/IsoCorrAlg.hh" #include #include #include "ElectronObjects/SimpleExtrapolatedTrack.hh" // For Pes RDeltaPhi #include "Calor/calor_alg.hh" // and DeltaR #include "ElectronUser/StdPEMElectron.hh" #include #include #include #include #include #include #include #include #include #include using namespace PhotonVariables; namespace { // TStnErrorLogger* fgLogger; } //_____________________________________________________________________________ Int_t StntupleInitPhoenixElectron(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.; AbsEnv* env = AbsEnv::instance(); //----------------------------------------------------------------------------- // backward pointer to CDF electron is used for internal bookkeeping only // during the STNTUPLE filling stage //----------------------------------------------------------------------------- Ele->fCdfEmObject = (CdfEmObject*) bob; 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; //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 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()); } } //----------------------------------------------------------------------------- // store RAW CEM face correction in fEtcor (as name allows) // also mark those electrons, for which correction failed // corrected momentum should be multiplied by fEtcor //----------------------------------------------------------------------------- Ele->fEtcor = 1; if (cluster->region() == 0) { // central electron if (CemCorrAlg::GetCemCorr(*bob,"USETRACK",cemcorr) == CemCorrAlg::OK) { Ele->fEtcor = 1./cemcorr; } else { Ele->fEtcor = -1.; } } if (cluster->region() == 1) { // plug electron if (PemCorrAlg::GetPemCorr(*bob,pemcorr) == PemCorrAlg::OK) { Ele->fEtcor = 1./pemcorr; } else { Ele->fEtcor = -1.; } } 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->fEmIso4 = cluster->emIsolationEt4(); Ele->fHadIso4 = cluster->hadIsolationEt4(); Ele->fLshr = bob->lshr(); Ele->fLshr2 = bob->lshr(2); //----------------------------------------------------------------------------- // TRACK QUANTITIES //----------------------------------------------------------------------------- Ele->fNTracks = bob->numberTracks(); bobs_track = (CdfTrack*) bob->maxPtTrack().operator->(); if ((! bobs_track) && (Ele->fNTracks != 0)) { //----------------------------------------------------------------------------- // we are in deep trouble //----------------------------------------------------------------------------- // fgLogger->Report(-102, // Form("StntupleInitPhoenixElectron: bad Bobs electron")); Ele->fNTracks = 0; } if (bobs_track) { 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 is EM_E/P, not TOT_E/P Ele->fEp = mom->e()/bobs_mom.mag(); //----------------------------------------------------------------------------- // 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); } else { 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; } //----------------------------------------------------------------------------- // 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()) { 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(); Ele->fDelX = cesalign.xtrack_corr(Ele->fXtrk,barrel,wedge)-Ele->fXces; } 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; } } } //----------------------------------------------------------------------------- // 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(); } } } else if (Ele->IsPlug()) { //----------------------------------------------------------------------------- // PLUG electron //----------------------------------------------------------------------------- StdPEMElectron* pem_bob = new StdPEMElectron(bob); 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(); Ele->fPesEnergy[0] = pes_cluster->plug2dUenergy(); Ele->fPesEnergy[1] = pes_cluster->plug2dVenergy(); Ele->fPes5x9[0] = pes_cluster->pesCluster0().pesProfileRatio5by9(); //5by9 shower profile ration for u strips Ele->fPes5x9[1] = pes_cluster->pesCluster1().pesProfileRatio5by9(); //5by9 shower profile ration for v strips 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 geometry // of the PES strips. 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; } //here we will the PPR energy using bob's function. Ele->fPprEnergy=emobj_alg::PPREnergy(bob); if (pem_bob) delete pem_bob; } if (bobs_track) { // loop over the COT SL's and figure the number of segments Ele->fNasl = 0; Ele->fNssl = 0; for (int sl=0; sl<8; sl++) { if (bobs_track->numCTHits(sl) >= TStnElectron::kMinCotHits[sl]) { // SL 0,2,4,6 are stereo, the rest - axial if (sl%2 == 0) { Ele->fNssl++; }else { Ele->fNasl++; } } } } else { Ele->fNasl = -1; Ele->fNssl = -1; } Ele->fIdwrd = 0; Ele->fIswrd = 0; Ele->fConwrd = 0; int fid_smx = fidEle(*bob,"SHWRMAX"); int fid_trk = fidEle(*bob,"TRACK"); Ele->fFid = ((fid_smx & 0xff) << 8) | (fid_trk & 0xff); Ele->fTime = phoTime(*bob); // relative correction for the leakage leak = 0; IsoCorrAlg::GetIsoCorr(*bob,leak); Ele->fLeakCorr = leak/cluster->emEt(); return 0; } //_____________________________________________________________________________ Int_t StntupleInitPhoenixElectronBlock(TStnDataBlock* block, AbsEvent* event, int mode){ int ev_number, rn_number; char process[100], description[100]; int inum[100]; ev_number = AbsEnv::instance()->trigNumber(); rn_number = AbsEnv::instance()->runNumber(); if (block->Initialized(ev_number,rn_number)) return 0; TPhoenixElectronBlock* data = (TPhoenixElectronBlock*) block; data->Clear(); TStnEvent* ev = data->GetEvent(); // fgLogger = ev->GetErrorLogger(); // Find event info by collection name ... StntupleGetProcessName(data->CollName()->Data(), process, description); EventRecord::ConstIterator ci; int nele; const CdfEmObject *e1, *e2; vector elist; 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; igetEmCluster()->totalEt() < e2->getEmCluster()->totalEt()) { // swap the tracks int itmp = inum[i1]; inum[i1] = inum[i2]; inum[i2] = itmp; e1 = e2; } } } for (int i=0; iNewElectron(); StntupleInitPhoenixElectron(ele,bob); } //----------------------------------------------------------------------------- // mark event and run as initialized //----------------------------------------------------------------------------- data->f_RunNumber = rn_number; data->f_EventNumber = ev_number; return 0; } //_____________________________________________________________________________ Int_t StntuplePhoenixElectronBlockLinks(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; char process[100], description[100]; // TStnTrackBlock* track_block; TStnTrackBlock* phx_seed_block; TStnTrackBlock* phx_trk_block; TStnElectronBlock* stn_ele_block; TStnElectron* ele; const CdfEmObject* phoenix_emobject; const CdfEmObject* bob; TPhoenixElectronBlock* phx_ele_block = (TPhoenixElectronBlock*) block; StntupleGetProcessName(phx_ele_block->CollName()->Data(), process, description); char name[100]; Int_t nele, nele_original,rc; ev_number = AbsEnv::instance()->trigNumber(); rn_number = AbsEnv::instance()->runNumber(); //----------------------------------------------------------------------------- // assume that the block exists, it is supposed to be initialized // do not initialize the links 2nd time //----------------------------------------------------------------------------- if (! block->Initialized(ev_number,rn_number)) return -1; if (block->LinksInitialized()) return 0; //----------------------------------------------------------------------------- // pointers to the default track block and to the list of HPT leptons //----------------------------------------------------------------------------- TStnEvent* ev = phx_ele_block->GetEvent(); // track_block = (TStnTrackBlock*) ev->GetDataBlock("TrackBlock"); strcpy(name,process); strcat(name,"@Phoenix_Tracking"); phx_seed_block = (TStnTrackBlock*) ev->GetDataBlock(name); strcpy(name,process); strcat(name,"@PhoenixSI_Tracking"); phx_trk_block = (TStnTrackBlock*) ev->GetDataBlock(name); stn_ele_block = (TStnElectronBlock*) ev->GetDataBlock("ElectronBlock"); TObjArray* hptl = ev->GetListOfHptl(); rc = 0; nele = phx_ele_block->NElectrons(); nele_original = stn_ele_block->NElectrons(); //----------------------------------------------------------------------------- // loop over the Phoenix Electrons //----------------------------------------------------------------------------- for (int iele=0; ieleElectron(iele); phoenix_emobject = ele->GetCdfEmObject(); ele->SetPhoenixSeedID(-1,-1); // maxPt track not defined for PhoenixElectrons. ele->SetTrackNumber(-1); // for each Phoenix EmObject find Bob's EmObject for (int ori = 0; ori < nele_original; ori++) { bob = stn_ele_block->Electron(ori)->GetCdfEmObject(); if (bob->matchingEmCluster() == phoenix_emobject->matchingEmCluster()) { phx_ele_block->fOriginalElectron->Add(iele,ori); break; } } // end of loop over the bob's electrons //----------------------------------------------------------------------------- // Phoenix Stuff //----------------------------------------------------------------------------- const CdfTrack* seed_track[2]; const CdfTrackView &emtracks = phoenix_emobject->matchingTracks(); const CdfTrackView::collection_type &tracks = emtracks.contents(); int iseed = 0; for (CdfTrackView::const_iterator track_iter = tracks.begin(); track_iter != tracks.end(); ++track_iter) { const CdfTrack* track = &(**track_iter); if ((track->numCTHits() == 0) && (track->numSIHits() == 0)) { //----------------------------------------------------------------------------- // no hits - assume that this is a Phoenix track //----------------------------------------------------------------------------- seed_track[iseed] = &(**track_iter); int sicount=0; for (CdfTrackView::const_iterator tracksi_iter = tracks.begin(); tracksi_iter != tracks.end(); ++tracksi_iter) { const CdfTrack* tracksi = &(**tracksi_iter); const CdfTrack* parentcig=&(*tracksi->parent()); if (parentcig == seed_track[iseed]) { int itrk = phx_trk_block->TrackNumber(tracksi); phx_ele_block->fPhoenixSiTrkLink[iseed]->Add(iele,itrk); sicount++; } } //----------------------------------------------------------------------------- // if seed track has no si track set sitrklink to -1 //----------------------------------------------------------------------------- if (sicount == 0) phx_ele_block->fPhoenixSiTrkLink[iseed]->Add(iele,-1); iseed++; } } ele->SetPhoenixSeedID((phx_seed_block->TrackNumber(seed_track[0])), (phx_seed_block->TrackNumber(seed_track[1]))); //----------------------------------------------------------------------------- // fill tower link block //----------------------------------------------------------------------------- const std::vector& towlist = phoenix_emobject->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 //----------------------------------------------------------------------------- phx_ele_block->fLinksInitialized = 1; return rc; }