//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // // Component: PhysicsTowerDataMaker.cc // // Purpose: This is the class that makes the PhysicsTowerData object. // // Created: 13/09/1999 Pierre Savard // // History: // 29/08/2000 Matthias Toennesmann - Added FourVectorCalculator in clump function. No changes in StandardCalculator, // because it is used in JetCluModule. // 30/08/2000 Matthias Toennesmann - In function clump: description string needed some diversity -> must now be given // when calling the function. // 17/05/2001 Matthias Toennesmann - Added a call to PhysicsTowerData::find() at the beginning of methods create and // clump. If a PTD object with the right description string and parameters already // exists, just grab it from the event. // 11/02/2002 Matthias Toennesmann - Always fill totEta and totPhi in clump(). Needed in JetClu's pre-clustering. // Remove use of the Locations object (_locs) in clump(). Use the data stored in the // towers that get clumped instead. // 15/08/2002 Jean-Francois Arguin - Activate making of PhysicsTowerData from different sources than CalData: // now CdfTrackView and HEPG also available. // "Clump" function changed for the new sources. It just stores towers above seed // threshold. // 25/11/2002 Matthias Toennesmann - Add SourceHEPGHerwigPartons and SourceHEPGHerwigHadrons. // Add clump function for CdfTracks and HEPG particles. // 01/05/2003 Dan MacQueen - Changes made to allow use of three-vertex corrections to Locations object. // // // Notes: // - This maker class uses a data input source and a tower calculator object to make physics towers (see SourceHEPG.hh // and SourceCalData.hh for input source examples and StdCalculator and FourVectorCalculator for calculator examples). // The default constructor of this class uses SourceCalData and StdCalculator, i.e. the Run1b default. // // - The description string and the parameter set are used to find a PhysicsTowerData object in an event. // We will move to RCPID when available. // // 01/05/2003 Ben Cooper - Added SourceHEPGHadronsWJets //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //============================================================================= // Edm //============================================================================= #include "BaBar/Cdf.hh" #include "Edm/EventRecord.hh" #include "Edm/GenericConstHandle.hh" #include "Edm/ConstHandle.hh" #include "Edm/Handle.hh" //============================================================================= // The Class header //============================================================================= #include "ErrorLogger_i/gERRLOG.hh" #include "Calor/PhysicsTowerDataMaker.hh" #include "Calor/SourceCalData.hh" #include "Calor/SourceHEPG.hh" #include "Calor/SourceCdfTrackView.hh" #include "Calor/SourceTrackCal.hh" #include "Calor/SourceHEPGHerwigPartons.hh" #include "Calor/SourceHEPGHerwigHadrons.hh" #include "Calor/SourceHEPGPythiaPartons.hh" #include "Calor/SourceHEPGPythiaPartons.hh" #include "Calor/SourceHEPGHerwigHadronsWjets.hh" #include "Calor/SourceHEPGPythiaHadronsWjets.hh" #include "Calor/StandardCalculator.hh" #include "Calor/FourVectorCalculator.hh" #include "Calor/MetCalculator.hh" //============================================================================= // Standard Library //============================================================================= #include PhysicsTowerDataMaker::PhysicsTowerDataMaker() { _source = 0; _locs = new Locations(0.0); _calc = 0; } PhysicsTowerDataMaker::~PhysicsTowerDataMaker() { if (_source) delete _source; if (_locs) delete _locs ; if (_calc) delete _calc; } PhysicsTowerData_ch PhysicsTowerDataMaker::create(AbsEvent* anEvent, PhysicsTowerParams& params, const std::string& description) { const TowerCorrectionColl *TowCorr = NULL; TowCorr=new TowerCorrectionColl(); PhysicsTowerData_ch theHandle; // TowerCorrectionColl* towercorr=new TowerCorrectionColl(); theHandle = create(anEvent, params, description, *TowCorr); // towercorr->destroy(); return theHandle; } // pass here the handle instead? PhysicsTowerData_ch PhysicsTowerDataMaker::create(AbsEvent* anEvent, PhysicsTowerParams& params, const std::string& description, const TowerCorrectionColl& towercorr) { PhysicsTowerData_ch theHandle; float correm[TENETA][TENPHI]; float corrhad[TENETA][TENPHI]; // Define the source: CalData, HEPG bank or CdfTrackView if(params.inputType() == 0){ _source = new SourceCalData(); memset(correm,0,TENETA*TENPHI*sizeof(float)); memset(corrhad,0,TENETA*TENPHI*sizeof(float)); for (int ieta=0;ietainitEvent(anEvent); params.setInputID(foundSource); } else{ foundSource = _source->initEvent(anEvent,params.inputID()); } if(foundSource == 0){ ERRLOG(ELerror,"PhysicsTowerDataMaker::create") << "Could not create towers." << endmsg; theHandle.set_null(); } else{ // ID _source->getID(); // get CalData identifier // int sourceID = 0; // temp until rcpid is ready PhysicsTowerData::Error result = PhysicsTowerData::find(theHandle,¶ms,description,anEvent); if(!result){ if(params.calculatorType() != 0 && params.calculatorType() != 1 && params.calculatorType() != 2){ ERRLOG(ELerror,"PhysicsTowerDataMaker::create") << "Could not find requested calculator." << endmsg; theHandle.set_null(); } else{ switch(params.calculatorType()){ case 0: _calc = new StandardCalculator(params.calcThreshold()); break; case 1: _calc = new FourVectorCalculator(params.calcThreshold()); break; case 2: _calc = new MetCalculator(params.calcThreshold()); } float vtx[3]; float slope[2]; vtx[0] = params.vertex(0); vtx[1] = params.vertex(1); vtx[2] = params.vertex(2); slope[0] = params.slope(0); slope[1] = params.slope(1); // Make a new Locations object if any component of the 3-vertex has changed. if(vtx[0] != _locs->threeVertex(0) || vtx[1] != _locs->threeVertex(1) || vtx[2] != _locs->threeVertex(2)){ delete _locs; _locs = new Locations(vtx, slope); } // std::cout << *_locs; _pData = new PhysicsTowerData(params); while(!(_source->done())){ EnergyData* eData = _source->next(); if(eData){ bool miniplug=false; if ((params.inputType() == 0 || params.inputType() == 6) && eData->iEta()<4 || eData->iEta()>47) miniplug=true; if (!miniplug) { PhysicsTower* pTower; // Input from CalData if(params.inputType() == 0){ int ie=eData->iEta(); int ip=eData->iPhi(); if (correm[ie][ip]==0) correm[ie][ip]=1; if (corrhad[ie][ip]==0) corrhad[ie][ip]=1; pTower = _calc->makeCorrPhysicsTower(eData,_locs,correm[ie][ip],corrhad[ie][ip]); } // Input from CdfTrack or HEPG else{ pTower = _calc->makePhysicsTower(eData,_locs); } if (pTower) { bool ok = _pData->add(pTower); if(!ok) ERRLOG(ELerror,"PhysicsTowerDataMaker::create") << "Tried to add already existent tower." << endmsg; } } delete eData; } } _pData->set_description(description); GenericConstHandle writeHandle(anEvent->append(_pData)); result = PhysicsTowerData::find(theHandle,¶ms,description,anEvent); if(!result){ ERRLOG(ELerror,"PhysicsTowerDataMaker::create") << "Could not retrieve PhysicsTowerData." << endmsg; theHandle.set_null(); } // cleanup delete _calc; _calc = 0; } } } delete _source; _source = 0; return theHandle; } PhysicsTowerData_h PhysicsTowerDataMaker::create_h(AbsEvent* anEvent, PhysicsTowerParams& params, const std::string& description) { PhysicsTowerData_h theHandle; TowerCorrectionColl* towercorr=new TowerCorrectionColl(); theHandle = create_h(anEvent, params, description, *towercorr); towercorr->destroy(); return theHandle; } //New method created by Mircea: PhysicsTowerData_h PhysicsTowerDataMaker::create_h(AbsEvent* anEvent, PhysicsTowerParams& params, const std::string& description, const TowerCorrectionColl& towercorr) { PhysicsTowerData_h theHandle; float correm[TENETA][TENPHI]; float corrhad[TENETA][TENPHI]; memset(correm,0,TENETA*TENPHI*sizeof(float)); memset(corrhad,0,TENETA*TENPHI*sizeof(float)); for(TowerCorrectionColl::const_iterator towc = towercorr.contents().begin(); towc != towercorr.contents().end();++towc) { int ieta=int(towc->iEta()); int iphi=int(towc->iPhi()); correm[ieta][iphi]=1.0; corrhad[ieta][iphi]=1.0; } // if(params.inputType() == 1){ // // _input = new SourceHEPG(); // ERRLOG(ELerror,"PhysicsTowerDataMaker::create_h") << "HEPG source not implemented yet!" << endmsg; // theHandle.set_null(); // } if(params.inputType() == 0){ _source = new SourceCalData(); } else if(params.inputType() == 1){ _source = new SourceHEPG(); } else if(params.inputType() == 2){ _source = new SourceCdfTrackView(); } else if(params.inputType() == 3){ _source = new SourceHEPGHerwigPartons(); } else if(params.inputType() == 4){ _source = new SourceHEPGHerwigHadrons(); } else if(params.inputType() == 5){ _source = new SourceHEPGPythiaPartons(); } else if(params.inputType() == 6){ _source = new SourceTrackCal(); } else { ERRLOG(ELerror,"PhysicsTowerDataMaker::create") << "Unrecognized input type" << endmsg; theHandle.set_null(); return theHandle; } // Tell source to fetch new object in the event. Id foundSource = 0; if(params.inputID() == 0){ foundSource = _source->initEvent(anEvent); params.setInputID(foundSource); } else{ foundSource = _source->initEvent(anEvent,params.inputID()); } if(foundSource == 0){ ERRLOG(ELerror,"PhysicsTowerDataMaker::create") << "Could not create towers." << endmsg; theHandle.set_null(); } else{//1 // ID _source->getID(); // get CalData identifier // int sourceID = 0; // temp until rcpid is ready // PhysicsTowerData::Error result = PhysicsTowerData::find(theHandle,¶ms,description,anEvent); bool result = 0; if(!result){ //ura ! if(params.calculatorType() != 0 && params.calculatorType() != 1 && params.calculatorType() != 2){ ERRLOG(ELerror,"PhysicsTowerDataMaker::create") << "Could not find requested calculator." << endmsg; theHandle.set_null(); } else{//3 switch(params.calculatorType()){ case 0: _calc = new StandardCalculator(params.calcThreshold()); break; case 1: _calc = new FourVectorCalculator(params.calcThreshold()); break; case 2: _calc = new MetCalculator(params.calcThreshold()); } float vtx[3]; float slope[2]; vtx[0] = params.vertex(0); vtx[1] = params.vertex(1); vtx[2] = params.vertex(2); slope[0] = params.slope(0); slope[1] = params.slope(1); // Make a Locations object if the three-vertex has changed. if(vtx[0] != _locs->threeVertex(0) || vtx[1] != _locs->threeVertex(1) || vtx[2] != _locs->threeVertex(2)){ delete _locs; _locs = new Locations(vtx, slope); } theHandle = new PhysicsTowerData(params); while(!(_source->done())){ EnergyData* eData = _source->next(); bool miniplug=false; if ((params.inputType() == 0 || params.inputType() == 6)&& eData->iEta()<4 || eData->iEta()>47) miniplug=true; if (!miniplug) { PhysicsTower* pTower; // Input from CalData if(params.inputType() == 0){ int ie=eData->iEta(); int ip=eData->iPhi(); if (correm[ie][ip]==0) correm[ie][ip]=1; if (corrhad[ie][ip]==0) corrhad[ie][ip]=1; pTower = _calc->makeCorrPhysicsTower(eData,_locs,correm[ie][ip],corrhad[ie][ip]); } // Input from CdfTrack or HEPG else{ pTower = _calc->makePhysicsTower(eData,_locs); } if(pTower){ // std::cout << " theHandle ready to have a tower added "<< std::endl; bool ok = theHandle->add(pTower); if(!ok) ERRLOG(ELerror,"PhysicsTowerDataMaker::create") << "Tried to add already existent tower." << endmsg; } } delete eData; } theHandle->set_description(description); //Append the new collection to the event record // GenericConstHandle writeHandle(anEvent->append(_pData)); //result = PhysicsTowerData::find(theHandle,¶ms,description,anEvent); bool result = 1 ; if(!result){ ERRLOG(ELerror,"PhysicsTowerDataMaker::create") << "Could not retrieve PhysicsTowerData." << endmsg; theHandle.set_null(); } // cleanup delete _calc; _calc = 0; }//3 } //ura! }//1 delete _source; _source = 0; return theHandle; } PhysicsTowerData_ch PhysicsTowerDataMaker::clump(AbsEvent* anEvent, PhysicsTowerData_ch& theData, float etThreshold, const std::string& description) { // I will simplify this method in the future ... // Warning: Some of the PhysicsTowers made with this method will have no emEta, hadEta, emPhi, hadPhi, emEt, hadEt, 4-vector // information. They may be set to 0 instead. Note that JetClu (which is currently the only module using this method) // doesn't need those numbers. // ID _source->getID(); // get CalData identifier PhysicsTowerParams params = theData->parameters(); params.setJetSeedThreshold(etThreshold); PhysicsTowerData_ch theHandle; PhysicsTowerData::Error result = PhysicsTowerData::find(theHandle,¶ms,description,anEvent); if(result) return theHandle; _pData = new PhysicsTowerData(params); // For CalData input if(params.inputType() ==0){ // for(int iEta = 4; iEta != TOWER_NETA-4; ++iEta){ for(int iEta = 0; iEta != TOWER_NETA; ++iEta){ for(int iPhi = 0; iPhi != TOWER_PHI_SEG[TOWER_TYPE[iEta]]; ++iPhi){ const PhysicsTower* tower1 = theData->tower(iEta,iPhi); if(tower1){ if(TOWER_TYPE[iEta] == 5){ if(iPhi%2 == 0){ // an iPhi even tower const PhysicsTower* tower2 = theData->tower(iEta,iPhi + 1); if(tower2){ // we have two towers to deal with if(params.calculatorType() != 1){ size_t iEta = tower1->iEta(); size_t iPhi = tower1->iPhi()/2; float emEnergy = tower1->emEnergy() + tower2->emEnergy(); float hadEnergy = tower1->hadEnergy() + tower2->hadEnergy(); float emEt = tower1->emEt() + tower2->emEt(); float hadEt = tower1->hadEt() + tower2->hadEt(); float totEt = emEt + hadEt; float emEta = emEnergy ? (tower1->emEta()*tower1->emEnergy() + tower2->emEta()*tower2->emEnergy())/emEnergy : 0; float hadEta = hadEnergy ? (tower1->hadEta()*tower1->hadEnergy() + tower2->hadEta()*tower2->hadEnergy())/hadEnergy : 0; float totEta = (emEta*emEt + hadEta*hadEt)/totEt; float emPhi = emEnergy ? (tower1->emPhi()*tower1->emEnergy() + tower2->emPhi()*tower2->emEnergy())/emEnergy : 0; float hadPhi = hadEnergy ? (tower1->hadPhi()*tower1->hadEnergy() + tower2->hadPhi()*tower2->hadEnergy())/hadEnergy : 0; float totPhi = (emPhi*emEt + hadPhi*hadEt)/totEt; if(totEt >= etThreshold){ PhysicsTower* newTower = new PhysicsTower(totEta,emEta,hadEta, totPhi,emPhi,hadPhi, emEnergy,hadEnergy, HepLorentzVector(), emEt,hadEt,totEt, iEta,iPhi); bool ok = _pData->add(newTower); if(!ok) ERRLOG(ELerror,"PhysicsTowerDataMaker::clump") << "Tried to add already existent tower:" << "iEta =" << iEta << "iPhi =" << iPhi << endmsg; } } else{ // FourVectorCalculator - add four-momenta of towers. HepLorentzVector fourMomentum = tower1->fourMomentum() + tower2->fourMomentum(); float pt = fourMomentum.vect().perp(); if(pt >= etThreshold){ size_t iEta = tower1->iEta(); size_t iPhi = tower1->iPhi()/2; float emEnergy = tower1->emEnergy() + tower2->emEnergy(); float hadEnergy = tower1->hadEnergy() + tower2->hadEnergy(); float totP = fourMomentum.rho(); float totE = fourMomentum.e(); float totEt = totE/totP*pt; float totEta = fourMomentum.pseudoRapidity(); float totPhi = fourMomentum.phi(); if(totPhi < 0) totPhi += TWOPI; PhysicsTower* newTower = new PhysicsTower(totEta,0,0, totPhi,0,0, emEnergy,hadEnergy, fourMomentum, 0,0,totEt, iEta,iPhi); bool ok = _pData->add(newTower); if(!ok) ERRLOG(ELerror,"PhysicsTowerDataMaker::clump") << "Tried to add already existent tower:" << "iEta =" << iEta << "iPhi =" << iPhi << endmsg; } } } else{ // we only have one tower to deal with if(params.calculatorType() != 1){ // StandardCalculator - no changes here because this is used in JetClu (Run I)! size_t iEta = tower1->iEta(); size_t iPhi = tower1->iPhi()/2; float emEnergy = tower1->emEnergy(); float hadEnergy = tower1->hadEnergy(); float emEt = tower1->emEt(); float hadEt = tower1->hadEt(); float totEt = emEt + hadEt; float emEta = tower1->emEta(); float hadEta = tower1->hadEta(); float totEta = (emEta*emEt + hadEta*hadEt)/totEt; float emPhi = tower1->emPhi(); float hadPhi = tower1->hadPhi(); float totPhi = (emPhi*emEt + hadPhi*hadEt)/totEt; if(totEt >= etThreshold){ PhysicsTower* newTower = new PhysicsTower(totEta,emEta,hadEta, totPhi,emPhi,hadPhi, emEnergy,hadEnergy, HepLorentzVector(), emEt,hadEt,totEt, iEta,iPhi); bool ok = _pData->add(newTower); if(!ok) ERRLOG(ELerror,"PhysicsTowerDataMaker::clump") << "Tried to add already existent tower:" << "iEta =" << iEta << "iPhi =" << iPhi << endmsg; } } else{ // FourVectorCalculator - change only iPhi; rest is copied. if(tower1->fourMomentum().vect().perp() >= etThreshold){ PhysicsTower* newTower = new PhysicsTower(tower1->totEta(),tower1->emEta(),tower1->hadEta(), tower1->totPhi(),tower1->emPhi(),tower1->hadPhi(), tower1->emEnergy(),tower1->hadEnergy(), tower1->fourMomentum(), tower1->emEt(),tower1->hadEt(),tower1->totEt(), tower1->iEta(),tower1->iPhi()/2); bool ok = _pData->add(newTower); if(!ok) ERRLOG(ELerror,"PhysicsTowerDataMaker::clump") << "Tried to add already existent tower:" << "iEta =" << tower1->iEta() << "iPhi =" << tower1->iPhi()/2 << endmsg; } } } } else{ // make sure we added tower before and if not, add it. const PhysicsTower* evenTower = theData->tower(iEta,iPhi - 1); if(!evenTower){ // add this tower ... if(params.calculatorType() != 1){ // StandardCalculator - no changes here because this is used in JetClu (Run I)! size_t iEta = tower1->iEta(); size_t iPhi = (tower1->iPhi() - 1)/2; float emEnergy = tower1->emEnergy(); float hadEnergy = tower1->hadEnergy(); float emEt = tower1->emEt(); float hadEt = tower1->hadEt(); float totEt = emEt + hadEt; float emEta = tower1->emEta(); float hadEta = tower1->hadEta(); float totEta = (emEta*emEt + hadEta*hadEt)/totEt; float emPhi = tower1->emPhi(); float hadPhi = tower1->hadPhi(); float totPhi = (emPhi*emEt + hadPhi*hadEt)/totEt; if(totEt >= etThreshold){ PhysicsTower* newTower = new PhysicsTower(totEta,emEta,hadEta, totPhi,emPhi,hadPhi, emEnergy,hadEnergy, HepLorentzVector(), emEt,hadEt,totEt, iEta,iPhi); bool ok = _pData->add(newTower); if(!ok) ERRLOG(ELerror,"PhysicsTowerDataMaker::clump") << "Tried to add already existent tower:" << "iEta =" << iEta << "iPhi =" << iPhi << endmsg; } } else{ // FourVectorCalculator - change only iPhi; rest is copied. if(tower1->fourMomentum().vect().perp() >= etThreshold){ PhysicsTower* newTower = new PhysicsTower(tower1->totEta(),tower1->emEta(),tower1->hadEta(), tower1->totPhi(),tower1->emPhi(),tower1->hadPhi(), tower1->emEnergy(),tower1->hadEnergy(), tower1->fourMomentum(), tower1->emEt(),tower1->hadEt(),tower1->totEt(), tower1->iEta(),(tower1->iPhi() - 1)/2); bool ok = _pData->add(newTower); if(!ok) ERRLOG(ELerror,"PhysicsTowerDataMaker::clump") << "Tried to add already existent tower:" << "iEta =" << tower1->iEta() << "iPhi =" << (tower1->iPhi() - 1)/2 << endmsg; } } } } } else{ // not a type 5 tower so just copy the pointer if(params.calculatorType() != 1){ // StandardCalculator - no changes here because this is used in JetClu (Run I)! if(tower1->totEt() >= etThreshold){ PhysicsTower* newTower = new PhysicsTower(*tower1); bool ok = _pData->add(newTower); if(!ok) ERRLOG(ELerror,"PhysicsTowerDataMaker::clump") << "Tried to add already existent tower:" << "iEta =" << tower1->iEta() << "iPhi =" << tower1->iPhi() << endmsg; } } else{ // FourVectorCalculator - check Pt, not Et! if(tower1->fourMomentum().vect().perp() >= etThreshold){ PhysicsTower* newTower = new PhysicsTower(*tower1); bool ok = _pData->add(newTower); if(!ok) ERRLOG(ELerror,"PhysicsTowerDataMaker::clump") << "Tried to add already existent tower:" << "iEta =" << tower1->iEta() << "iPhi =" << tower1->iPhi() << endmsg; } } } } } } } // End of CalData // If inputType is CdfTracks or HEPG ... else{ if(params.calculatorType() != 1){ // First collect tracks/particles pointing to the same calorimeter tower. vector iClumpEta[TOWER_NETA][TOWER_MAX_NPHI],iClumpPhi[TOWER_NETA][TOWER_MAX_NPHI]; for(int iEta = 0; iEta < TOWER_NETA; iEta++){ for(int iPhi = 0; iPhi < TOWER_PHI_SEG[TOWER_TYPE[iEta]]; iPhi++){ const PhysicsTower* tower = theData->tower(iEta,iPhi); if(tower){ float totEta = tower->totEta(); if(abs(totEta) < -log(tan(TOWER_THETA[0]*DEGRAD/2))){ int iRealEta; if(totEta <= 0){ for(size_t i = 0; i < TOWER_NETA/2; i++) if(totEta < -log(tan((180-TOWER_THETA[i + 1])*DEGRAD/2))){ iRealEta = i; break; } } else{ for(size_t i = 0; i < TOWER_NETA/2; i++) if(-totEta < -log(tan((180-TOWER_THETA[i + 1])*DEGRAD/2))){ iRealEta = TOWER_NETA - 1 - i; break; } } float totPhi = tower->totPhi(); if(totPhi < 0) totPhi += TWOPI; int iRealPhi = int(totPhi/TWOPI * TEPSEG[TETTYP[iRealEta]]); iClumpEta[iRealEta][iRealPhi].push_back(iEta); iClumpPhi[iRealEta][iRealPhi].push_back(iPhi); } } } } // Now clump the tracks/particles which point to the same calorimeter tower. // Ignore tracks/particles pointing to the miniplug region; the adjacency check in JetClu doesn't work here // because of the phi segmentation of the miniplug (12,18,24,30). for(int iRealEta = 4; iRealEta < TOWER_NETA - 4; iRealEta++){ for(int iRealPhi = 0; iRealPhi < TOWER_PHI_SEG[TOWER_TYPE[iRealEta]]; iRealPhi++){ // Don't create a new clumped tower for phi segments 1, 3, 5, ... if tower type is 5 (48 phi segments): if(TOWER_TYPE[iRealEta] != 5 || iRealPhi%2 == 0){ float emEnergy = 0, hadEnergy = 0; float emEt = 0, hadEt = 0; float emEta = 0, hadEta = 0; float emPhi0 = -100, hadPhi0 = -100; float emPhi = 0, hadPhi = 0; for(int iParticle = 0; iParticle < iClumpEta[iRealEta][iRealPhi].size(); iParticle++){ const PhysicsTower* tower = theData->tower(iClumpEta[iRealEta][iRealPhi][iParticle], iClumpPhi[iRealEta][iRealPhi][iParticle]); emEnergy += tower->emEnergy(); hadEnergy += tower->hadEnergy(); emEt += tower->emEt(); hadEt += tower->hadEt(); emEta += tower->emEt()*tower->emEta(); hadEta += tower->hadEt()*tower->hadEta(); if(emPhi0 == -100) emPhi0 = tower->emPhi(); else{ float dEmPhi = tower->emPhi() - emPhi0; if(dEmPhi > PI) dEmPhi -= TWOPI; else if(dEmPhi < -PI) dEmPhi += TWOPI; emPhi += tower->emEt()*dEmPhi; } if(hadPhi0 == -100) hadPhi0 = tower->hadPhi(); else{ float dHadPhi = tower->hadPhi() - hadPhi0; if(dHadPhi > PI) dHadPhi -= TWOPI; else if(dHadPhi < -PI) dHadPhi += TWOPI; hadPhi += tower->hadEt()*dHadPhi; } } if(TOWER_TYPE[iRealEta] == 5){ // Add tracks/particles pointing to the next tower in phi. for(int iParticle = 0; iParticle < iClumpEta[iRealEta][iRealPhi + 1].size(); iParticle++){ const PhysicsTower* tower = theData->tower(iClumpEta[iRealEta][iRealPhi + 1][iParticle], iClumpPhi[iRealEta][iRealPhi + 1][iParticle]); emEnergy += tower->emEnergy(); hadEnergy += tower->hadEnergy(); emEt += tower->emEt(); hadEt += tower->hadEt(); emEta += tower->emEt()*tower->emEta(); hadEta += tower->hadEt()*tower->hadEta(); if(emPhi0 == -100) emPhi0 = tower->emPhi(); else{ float dEmPhi = tower->emPhi() - emPhi0; if(dEmPhi > PI) dEmPhi -= TWOPI; else if(dEmPhi < -PI) dEmPhi += TWOPI; emPhi += tower->emEt()*dEmPhi; } if(hadPhi0 == -100) hadPhi0 = tower->hadPhi(); else{ float dHadPhi = tower->hadPhi() - hadPhi0; if(dHadPhi > PI) dHadPhi -= TWOPI; else if(dHadPhi < -PI) dHadPhi += TWOPI; hadPhi += tower->hadEt()*dHadPhi; } } } float totEt = emEt + hadEt; if(totEt >= etThreshold){ emEta = emEt ? emEta/emEt : 0; hadEta = hadEt ? hadEta/hadEt : 0; emPhi = emEt ? emPhi0 + emPhi/emEt : 0; if(emPhi > TWOPI) emPhi -= TWOPI; else if(emPhi < 0) emPhi += TWOPI; hadPhi = hadEt ? hadPhi0 + hadPhi/hadEt : 0; if(hadPhi > TWOPI) hadPhi -= TWOPI; else if(hadPhi < 0) hadPhi += TWOPI; float totEta = (emEta*emEt + hadEta*hadEt)/totEt; float dTotPhi = hadPhi - emPhi; if(dTotPhi > PI) dTotPhi -= TWOPI; else if(dTotPhi < -PI) dTotPhi += TWOPI; float totPhi = emPhi + hadEt/totEt*dTotPhi; if(totPhi > TWOPI) totPhi -= TWOPI; else if(totPhi < 0) totPhi += TWOPI; int iNewRealPhi; if(TOWER_TYPE[iRealEta] != 5) iNewRealPhi = iRealPhi; else iNewRealPhi = iRealPhi/2; PhysicsTower* newTower = new PhysicsTower(totEta,emEta,hadEta, totPhi,emPhi,hadPhi, emEnergy,hadEnergy, HepLorentzVector(), emEt,hadEt,totEt, iRealEta,iNewRealPhi); bool ok = _pData->add(newTower); if(!ok) ERRLOG(ELerror,"PhysicsTowerDataMaker::clump") << "Tried to add already existent tower:" << "iEta =" << iRealEta << "iPhi =" << iNewRealPhi << endmsg; } } } } } else{ // This is used by the MidPoint algorithm only. It is sufficient to just copy the pointers. for(int iEta = 0; iEta != TOWER_NETA; ++iEta){ for(int iPhi = 0; iPhi != TOWER_PHI_SEG[TOWER_TYPE[iEta]]; ++iPhi){ const PhysicsTower* tower = theData->tower(iEta,iPhi); if(tower){ if(tower->fourMomentum().vect().perp() >= etThreshold){ PhysicsTower* newTower = new PhysicsTower(*tower); bool ok = _pData->add(newTower); if(!ok) ERRLOG(ELerror,"PhysicsTowerDataMaker::clump") << "Tried to add already existent tower:" << "iEta =" << tower->iEta() << "iPhi =" << tower->iPhi() << endmsg; } // End of four-vector } } // Iteration on iphi } // Iteration on ieta } } // End of HEPG or CdfTrackView part // Check that we are returning something: // PhysicsTowerData::Iterator beginIter = _pData->begin(); // PhysicsTowerData::Iterator endIter = _pData->end(); // if(beginIter == endIter) std::cerr << "Warning!!! clump collection is empty"; _pData->set_description(description); GenericConstHandle writeHandle(anEvent->append(_pData)); result = PhysicsTowerData::find(theHandle,¶ms,description,anEvent); return theHandle; }