XANADOO: writeESDForFAMOS.cpp Source File

00001 
00002 //we have to minimize includes in order not to have compile problems!!!!!
00004 #include <Utilities/Configuration/interface/Architecture.h>
00005 
00006 #include <XANADOO/domain/interface/writeESDForFAMOS.h>
00007 #include <XANADOO/domain/interface/XANAEsdEvent.h>
00008 #include <XANADOO/XANAClusters/interface/XANASuperCluster.h>
00009 #include <XANADOO/XANAClusters/interface/XANAEmCluster.h>
00010 #include <XANADOO/XANAClusters/interface/XANAHadCluster.h>
00011 #include <XANADOO/XANAClusters/interface/XANACaloRecHit.h>
00012 #include <XANADOO/XANATracks/interface/XANATrackHit.h>
00013 #include <XANADOO/XANAVertex/interface/XANAVertex.h>
00014 #include <XANADOO/XANAMcInfo/interface/XANAGeneratorEvent.h>
00015 #include <XANADOO/XANAMcInfo/interface/XANAGeneratorParticle.h>
00016 #include <XANADOO/XANAMcInfo/interface/XANAGeantEvent.h>
00017 #include <XANADOO/XANAMcInfo/interface/XANAGeantVertex.h>
00018 #include <XANADOO/XANAMcInfo/interface/XANAGeantTrack.h>
00019 #include <XANADOO/XANATriggerInfo/interface/XANATriggerInfo.h>
00020 #include <XANADOO/XANAUtilities/interface/XANAToolkit.h>
00021 #include <XANADOO/XANAUtilities/interface/XANAMemLog.h>
00022 #include <XANADOO/XANAElectronCandidate/interface/XANAElectronCandidate.h>
00023 #include <XANADOO/XANAElectronCandidate/interface/XANAElectronGSFTrack.h>
00024 #include <XANADOO/XANAElectronCandidate/interface/XANAElectronTools.h>
00025 #include <XANADOO/XANAElectronCandidate/interface/XANAElectronSeed.h>
00026 #include <XANADOO/XANAMuons/interface/XANAMuonCandidate.h>
00027 #include <XANADOO/XANAMuons/interface/XANAMuonTrack.h>
00028 #include <XANADOO/XANAJets/interface/XANAJet.h>
00029 #include <XANADOO/XANAMet/interface/XANAMet.h>
00030 
00031 #include <ElectronPhoton/ClusterTools/interface/EgammaCandidate.h>
00032 #include <ElectronPhoton/ElectronFromPixels/interface/EPHLTElectronSeed.h>
00033 #include <CommonReco/PatternTools/interface/TTrack.h>
00034 #include <ElectronPhoton/EgammaElectron/interface/ElectronCandidate.h>
00035 #include <ElectronPhoton/EgammaPhoton/interface/PhotonCandidate.h>
00036 #include <ElectronPhoton/EgammaAnalysis/interface/EgammaItr.h>
00037 #include <ElectronPhoton/EgammaClusters/interface/EgammaBasicCluster.h>
00038 #include <ElectronPhoton/EgammaClusters/interface/EgammaCluster.h>
00039 #include <ElectronPhoton/ClusterTools/interface/EgammaSuperCluster.h>
00040 #include <ElectronPhoton/EgammaPreshower/interface/EgammaEndcapCluster.h>
00041 #include <ElectronPhoton/EgammaAnalysisFactory/interface/EgammaAnalysisFactory.h>
00042 #include <ElectronPhoton/ClusterTools/interface/IsoCalculator.h>
00043 
00044 #include <CommonReco/PatternTools/interface/ImpactPointExtrapolator.h>
00045 #include <CommonReco/MaterialEffects/interface/MultipleScatteringUpdator.h>
00046 #include <CommonReco/GeomPropagators/interface/AnalyticalPropagator.h>
00047 #include <TrackerReco/TkSeedGenerator/interface/CombinatorialSeedGeneratorFromPixel.h>
00048 
00049 #include <CommonDet/BasicDet/interface/DetUnit.h>
00050 #include <CommonDet/BasicDet/interface/DetType.h>
00051 #include <CommonDet/BasicDet/interface/Readout.h>
00052 
00053 #include <GeneratorInterface/HepEvent/interface/RawHepEvent.h>
00054 #include <GeneratorInterface/EventProxyInterface/interface/HepEventG3EventProxyReader.h>
00055 #include <GeneratorInterface/Particle/interface/RawParticle.h>
00056 #include <GeneratorInterface/Particle/interface/RawHepEventParticle.h>
00057 #include <GeneratorInterface/HepPDT/interface/HepPDTable.h>
00058 #include <GeneratorInterface/HepPDT/interface/HepParticleData.h>
00059 
00060 #include <Utilities/Notification/interface/PackageInitializer.h>
00061 #include <Utilities/Notification/interface/Observer.h> 
00062 #include <Utilities/UI/interface/SimpleConfigurable.h>
00063 
00064 #include <CARF/G3Event/interface/G3EventProxy.h>
00065 #include <CARF/Reco/interface/RecCollection.h>
00066 #include <CARF/Reco/interface/RecQuery.h>
00067 
00068 
00069 #include <ORCAInterface/EgammaInterface/interface/FastBasicReco.h>
00070 #include <ORCAInterface/EgammaInterface/interface/doFamosElectronTracks.h>
00071 #include <FamosGeneric/FamosManager/interface/FamosEventMgr.h>
00072 #include <FamosGeneric/FamosSimEvent/interface/FSimEvent.h>
00073 #include <FamosGeneric/FamosSimEvent/interface/FSimTrack.h>
00074 #include <FamosGeneric/FamosSimEvent/interface/FSimVertex.h>
00075 #include <ElectronPhoton/EgammaAnalysis/interface/EGDummy.h>
00076 #include <Examples/ExEgamma/interface/FMCKine.h>
00077 #include <CLHEP/Vector/ThreeVector.h>
00078 
00079 #include <TROOT.h>
00080 #include <TObjectTable.h>
00081 #include <TFile.h>
00082 #include <TTree.h>
00083 #include <TBranch.h>
00084 #include <TClonesArray.h>
00085 #include <TBits.h>
00086 
00087 #include <CLHEP/Units/PhysicalConstants.h>
00088 
00089 #include <vector>
00090 #include <map>
00091 #include <string>
00092 #include <iostream> 
00093 #include <iomanip>  
00094 #include <numeric>
00095 
00096 
00097 XANAEsdBuilderForFAMOS::XANAEsdBuilderForFAMOS() : 
00098  XANAEsdBuilder()
00099 {   
00100    myEventMgr = FamosEventMgr::instance();
00101 }
00102 
00103 XANAEsdBuilder::XANAEsdBuilder() : xanaEsdTree_(0),  xanaEsdEventBranch_(0), xanaEsdFile_(0), 
00104    xanaEsdEvent_(0), xanaGenEventBranch_(0), xanaGenEvent_(0), xanaGeantEventBranch_(0),
00105    xanaGeantEvent_(0), ciso_(0), xanamemlog_(0), nbClusters_(0), nbTracks_(0), nbElectrons_(0)
00106 { 
00107 
00108    // initialize observer    
00109    xana::out(1) << "Constructing an XANAEsdBuilderForFAMOS" << endl ;
00110    Observer<G3EventProxy *>::init() ; 
00111    Observer<G3SetUp *>::init() ; 
00112 
00113    // algo selection
00114    clusterizerName_ = 
00115     SimpleConfigurable<std::string>("EGBCluster","XANADOO:Clusterizer");
00116    superClusterizerName_ =
00117     SimpleConfigurable<std::string>("EGSCluster","XANADOO:SuperClusterizer");
00118    endcapSuperClusterizerName_ =
00119     SimpleConfigurable<std::string>("EGECluster","XANADOO:EndcapSuperClusterizer");
00120    hcalClusterizerName_ =
00121     SimpleConfigurable<std::string>("TowerBuilder","XANADOO:HcalClusterizer");
00122    trackReconstructorName_=
00123     SimpleConfigurable<std::string>("TransientCombinatorialTrackFinder","XANADOO:TrackReconstructor");
00124    vertexReconstructorName_=
00125     SimpleConfigurable<std::string>("PrincipalVertexFinder","XANADOO:VertexReconstructor");
00126    muonReconstructorName_=
00127     SimpleConfigurable<std::string>("GlobalMuonReconstructor","XANADOO:MuonReconstructor");
00128    electronReconstructorName_=
00129     SimpleConfigurable<std::string>("EgammaCandidate","XANADOO:ElectronReconstructor");
00130    gsfForwardFit_=
00131     SimpleConfigurable<std::string>("KF","XANADOO:gsfForwardFit");
00132 
00133    // xanadoo steering
00134    writeTree_ = SimpleConfigurable<bool>(false,"XANADOO:writeEsd");
00135    writeCaloRecHits_ = SimpleConfigurable<bool>(false,"XANADOO:EsdWithCaloRecHits");
00136    writeTrackHits_ = SimpleConfigurable<bool>(false,"XANADOO:EsdWithTrackHits");   writeAllTrackHits_ = SimpleConfigurable<bool>(false,"XANADOO:EsdWithAllTrackHits");
00137    writeEleTrackHits_ = SimpleConfigurable<bool>(false,"XANADOO:EsdWithEleTrackHits");
00138      // jets and mets
00139    jetsFromTrueInformation_ = SimpleConfigurable<bool>(false,"XANADOO:jetsFromTrueInformation");
00140    metsFromTrueInformation_ = SimpleConfigurable<bool>(false,"XANADOO:metsFromTrueInformation");
00141 
00142      // memory log
00143    memlog_ = SimpleConfigurable<bool>(false,"XANADOO:memlog");
00144    if (memlog_) xanamemlog_=new XANAMemLog();
00145 
00146    // creates output file
00147    fileName_ = SimpleConfigurable<std::string>("single_eminus_pt100.root","XANADOO:FileName");
00148    xanaEsdFile_ = new TFile(fileName_.value().c_str(),"RECREATE");
00149    compressionLevel_ = SimpleConfigurable<int>(1,"XANADOO:CompressionLevel");        
00150    xanaEsdFile_->SetCompressionLevel(compressionLevel_);
00151 
00152    // create TTree
00153    autosave_ = SimpleConfigurable<int>(1000000000,"XANADOO:AutoSaveSize");        
00154    bufferSize_ = SimpleConfigurable<int>(256000,"XANADOO:BufferSize");        
00155    splitMode_ = SimpleConfigurable<int>(99,"XANADOO:SplitMode");        
00156    xana::out(1) << "creating xanaEsdTree  .."<<endl; 
00157    treeName_ = SimpleConfigurable<std::string>("XANAESD","XANADOO:TreeName");
00158    xanaEsdTree_ = new TTree(treeName_.value().c_str(),"test XANADOO ESD tree");  
00159    xanaEsdTree_->SetAutoSave(autosave_);
00160    int bufferSizeValue = bufferSize_.value();
00161    if (splitMode_.value()) bufferSizeValue /= 4; //from ROOT example
00162 
00163    // initialize event and create branches
00164    xanaEsdEvent_ = new XANAEsdEvent();
00165    xana::out(1) << "creating xanaEsdEventBranch  .."<<endl; 
00166    xanaEsdEventBranch_ = 
00167    xanaEsdTree_->Branch("Event.","XANAEsdEvent",&xanaEsdEvent_,bufferSizeValue,splitMode_.value());
00168    xana::out(1) << "xanaEsdEventBranch  branch created .."<<endl; 
00169 
00170    xanaGenEvent_ = new XANAGeneratorEvent();
00171    xana::out(1) << "creating xanaGeneratorEventBranch  .."<<endl; 
00172    xanaGenEventBranch_ = 
00173    xanaEsdTree_->Branch("GeneratorEvent.","XANAGeneratorEvent",&xanaGenEvent_,bufferSizeValue,splitMode_.value());
00174    xana::out(1) << "xanaGeneratorEventBranch  branch created .."<<endl; 
00175 
00176    xanaGeantEvent_ = new XANAGeantEvent();
00177    xana::out(1) << "creating xanaGeantEventBranch  .."<<endl; 
00178    xanaGeantEventBranch_ = 
00179    xanaEsdTree_->Branch("GeantEvent.","XANAGeantEvent",&xanaGeantEvent_,bufferSizeValue,splitMode_.value());
00180    xana::out(1) << "xanaGeantEventBranch  branch created .."<<endl; 
00181 
00182    // create links between generator-geant event and ESD-geant Event
00183    xanaGeantEvent_->addGeneratorEvent(xanaGenEvent_);
00184    xanaGeantEvent_->addEsdEvent(xanaEsdEvent_);
00185 
00186   //CC is this still needed?   
00187    ciso_ = Singleton<IsoCalculator>::instance();
00188   
00189 }
00190    
00191          
00192 XANAEsdBuilderForFAMOS::~XANAEsdBuilderForFAMOS() 
00193 { 
00194   delete myEventMgr;
00195 }
00196    
00197 XANAEsdBuilder::~XANAEsdBuilder() 
00198 {
00199    
00200   xana::out(1) << "Deleting XANAEsdBuilderForFAMOS" << std::endl;  
00201   xanaEsdFile_->Write();
00202   xanaEsdTree_->Print();
00203   delete xanaEsdFile_;
00204   if (memlog_) xanamemlog_->writeHistogram();
00205 
00206 }
00207 
00208 void XANAEsdBuilder::upDate( G3EventProxy * ev ) 
00209 { }
00210 
00211 void XANAEsdBuilderForFAMOS::upDate( G3EventProxy * ev ) 
00212 { 
00213     
00214   printf("Starting update for EventProxy \n");fflush(stdout);
00215   xana::out(1) << "[update] updating XANAEsdBuilderForFAMOS for event " << std::endl;  
00216   if (!ev) return ;    
00217 
00218   //    UInt_t saveNumber=TProcessID::GetObjectCount();
00219      
00220   setHeader(ev);
00221   setEmClusters(ev);
00222   setSuperClusters(ev);
00223   setEndcapClusters(ev);
00224   setElectronTracks(ev);
00225   setElectrons(ev);
00226   setMcTruth(ev);
00227 
00228   if (writeTree_) writeTree();
00229   
00230   if (memlog_) xanamemlog_->fill();
00231 
00232   xanaEsdEvent_->clear();
00233   xanaGenEvent_->clear();
00234   xanaGeantEvent_->clear();
00235 
00236 //      TProcessID::SetObjectCount(saveNumber);
00237 
00238 }      
00239 
00240 void XANAEsdBuilder::setHeader( G3EventProxy * ev ) 
00241 {
00242 /*
00243   xana::out(8) << "[update] ============================  " << std::endl;
00244   xana::out(8) << "[update] Run #" << ev->simSignal()->id().runNumber() << "; ";
00245   xana::out(8) << " Event #" << ev->simSignal()->id().eventInRun() << std::endl;     
00246   xana::out(8) << "[update] ============================  " << std::endl;
00247   xanaEsdEvent_->setHeader(ev->simSignal()->id().runNumber(),ev->simSignal()->id().eventInRun());
00248 */
00249 }
00250 
00251 void XANAEsdBuilder::setEmClusters( G3EventProxy * ev ) 
00252 {
00253   xana::out(2) << "[update] em clusters " << std::endl;
00254   
00255   RecCollection<EgammaBasicCluster> clusters (RecQuery("EGBFCluster","1.0"));
00256   RecCollection<EgammaBasicCluster> ::const_iterator cluster ;
00257   xana::out(5) << "[update] " << clusters.size() << " " << clusterizerName_.value() <<" clusters found." << std::endl;
00258   XANAEmCluster *xcluster;
00259   clusterizedCaloRecHits_.clear();
00260   xana::out(2) << "[update] before cluster loop " << std::endl;     
00261   for (cluster=clusters.begin(); cluster!=clusters.end(); cluster++) {
00262     // hybrid is algo>=100   
00263  // ORIGINAL
00264      if ( (cluster->whichAlgo()>99 && fabs(cluster->eta())>1.4791) || 
00265       (cluster->whichAlgo()<99 && fabs(cluster->eta())<1.4791) ) continue;
00266  // MODIF CARO : ISLAND Partout
00267  //   if ( (cluster->whichAlgo()>99 && fabs(cluster->eta())>1.4791) || 
00268  //      (cluster->whichAlgo()>99 && fabs(cluster->eta())<1.4791) ) continue;
00269     xana::out(2) << "[update] creating new xanacluster "<< std::endl;
00270     xana::out(2) << "[update] cluster algoName is  "<< clusterizerName_.value() << " " << getAlgoName(cluster->whichAlgo())<<std::endl;
00271     HepPoint3D uncorrectedPositionAtFrontFace(0.,0.,0.);
00272     float esum=0;
00273     std::vector<const CaloRecHit *> rhits=cluster->rHits();
00274     for (std::vector<const CaloRecHit *>::const_iterator
00275            it=rhits.begin(); it != rhits.end(); it++) {
00276       esum +=   (*it)->energy();   
00277       uncorrectedPositionAtFrontFace += (*it)->energy() * (*it)->getMyCell().position();   
00278     }
00279     if (esum != 0.) uncorrectedPositionAtFrontFace /= esum ;
00280     xcluster = new XANAEmCluster(cluster->energy(), cluster->position(),
00281                                  uncorrectedPositionAtFrontFace,
00282                                  cluster->chi2(),clusterizerName_.value());
00283     xana::out(2) << "[update] adding cluster " << std::endl;
00284     xanaEsdEvent_->addEmCluster(xcluster);
00285     for (std::vector<const CaloRecHit *>::const_iterator
00286            it=rhits.begin(); it != rhits.end(); it++) {
00287       XANACaloRecHit *xhit = new XANACaloRecHit((*it)->energy(), (*it)->time(),
00288                                                 (*it)->chi2(), XANACellID((*it)->getMyCell().position()));
00289       xanaEsdEvent_->addEmClusterHit(xhit, xcluster);
00290       clusterizedCaloRecHits_.insert(pair<CellID,XANACaloRecHit>((*it)->getMyCell(),*xhit));
00291       clusterizedCells_.push_back((*it)->getMyCell());
00292       delete xhit;
00293     } 
00294     delete xcluster;
00295   }
00296 
00297 }
00298 
00299 
00300 void XANAEsdBuilder::setSuperClusters( G3EventProxy * ev ) 
00301 {
00302            
00303   xana::out(2) << "[update] superclusters " << std::endl;
00304   
00305   // Barrel superclusters (Hybrid clustering algo chosen for the moment)
00306   RecQuery q1 = FastBasicReco::fastSCConfig();
00307   RecCollection<EgammaSuperCluster> superClusters(q1);
00308   RecCollection<EgammaSuperCluster>::const_iterator superCluster ; 
00309   xana::out(5) << "[update] " << superClusters.size() << superClusterizerName_.value() << " superclusters found." << std::endl;
00310   XANASuperCluster *xscluster;    
00311 
00312   for (superCluster=superClusters.begin(); superCluster!=superClusters.end(); superCluster++) {
00313 // ORIGINAL
00314     if ( (fabs(superCluster->eta())>1.4791) ||
00315        (superCluster->seed()->getBasicCluster()->whichAlgo()<99 &&
00316         fabs(superCluster->eta())<1.4791) ) continue;
00317 // MODIF CARO : ISLAND Partout!
00318 //    if ( (fabs(superCluster->eta())>1.4791) || 
00319 //       (superCluster->seed()->getBasicCluster()->whichAlgo()>99 && 
00320 //        fabs(superCluster->eta())<1.4791) ) continue;
00321     xana::out(2) << "[update] constructing new supercluster from algo " << superCluster->seed()->WhichAlgo()
00322               << " at eta " << superCluster->eta() << " with energy " <<superCluster->energy() << std::endl;
00323     xscluster = new XANASuperCluster;
00324     xscluster->setEnergy(superCluster->energy()); 
00325     xscluster->setPosition(superCluster->position()); 
00326     xscluster->setSum1(superCluster->seed()->getS1()); 
00327     xscluster->setSum4(superCluster->seed()->getS4()); 
00328     xscluster->setSum9(superCluster->seed()->getS9()); 
00329     xscluster->setSum25(superCluster->seed()->getS25()); 
00330     xscluster->setHadronicOverEm(superCluster->seed()->getHoe());
00331     xscluster->setCaloIsolation(superCluster->iso()); 
00332     xscluster->setDisc1(superCluster->seed()->getDisc1()); 
00333     xscluster->setDisc2(superCluster->seed()->getDisc2()); 
00334     xscluster->setDisc3(superCluster->seed()->getDisc3()); 
00335     HepSymMatrix cov(2);
00336     cov(1,1) = superCluster->seed()->getCovEE();
00337     cov(1,2) = superCluster->seed()->getCovEP();
00338     cov(2,2) = superCluster->seed()->getCovPP();
00339     xscluster->setPositionCovarianceMatrix(cov);
00340     xscluster->setEnergyScaleFactor(superCluster->energy()/superCluster->energyUncorrected());
00341     xana::out(2) << "[update] supercluster algoName is " << superClusterizerName_.value() << std::endl;
00342     xscluster->setAlgoName(superClusterizerName_.value());
00343     xana::out(2) << "[update] adding super cluster " << std::endl;
00344     xanaEsdEvent_->addSuperCluster(xscluster);
00345 
00346     // add seed basic cluster
00347     // the seed is the first in BasicClusters list
00348     EgammaVSuperCluster::const_iterator icf = superCluster->basicClusters().first;
00349     HepPoint3D uncorrectedPositionAtFrontFace(0.,0.,0.);
00350     float esum=0;
00351     std::vector<const CaloRecHit *> rhits=(*icf)->rHits();
00352     for (std::vector<const CaloRecHit *>::const_iterator
00353            it=rhits.begin(); it != rhits.end(); it++) {
00354       esum +=   (*it)->energy();   
00355       uncorrectedPositionAtFrontFace += (*it)->energy() * (*it)->getMyCell().position();   
00356     }
00357     if (esum != 0.) uncorrectedPositionAtFrontFace /= esum ;    
00358     XANAEmCluster *seed 
00359       = new XANAEmCluster((*icf)->energy(),
00360                           (*icf)->position(),
00361                           uncorrectedPositionAtFrontFace,
00362                           (*icf)->chi2(),
00363                           clusterizerName_.value()); 
00364     xana::out(2) << "[update] adding seed " << std::endl;
00365     xanaEsdEvent_->addSeedCluster(seed, xscluster);
00366     delete seed;
00367 
00368     // now add brems clusters
00369     xana::out(2) << "[update] adding brems " << std::endl;
00370     XANAEmCluster *brem = 0;
00371     // beware, the superCluster::vector<EgammaBasicCluster> list includes the seed cluster
00372     EgammaVSuperCluster::const_iterator ic = superCluster->basicClusters().first+1;
00373     for (;ic != superCluster->basicClusters().second; ic++) {
00374       HepPoint3D uncorrectedPositionAtFrontFace(0.,0.,0.);
00375       esum=0;
00376       rhits=(*ic)->rHits();
00377       for (std::vector<const CaloRecHit *>::const_iterator
00378         it=rhits.begin(); it != rhits.end(); it++) {
00379         esum +=         (*it)->energy();   
00380         uncorrectedPositionAtFrontFace += (*it)->energy() * (*it)->getMyCell().position();   
00381       }
00382       if (esum != 0.) uncorrectedPositionAtFrontFace /= esum ;                
00383               brem = new XANAEmCluster((*ic)->energy(),
00384                                       (*ic)->position(),
00385                                       uncorrectedPositionAtFrontFace,
00386                                       (*ic)->chi2(),
00387                                       clusterizerName_.value());            
00388       xana::out(2) << "[upDate] brem energy " << brem->getEnergy() << std::endl;
00389       xanaEsdEvent_->addBremCluster(brem, xscluster);
00390       delete brem;
00391     } 
00392     delete xscluster;
00393   }    
00394 
00395 }
00396 
00397 
00398 void XANAEsdBuilder::setEndcapClusters( G3EventProxy * ev ) 
00399 {
00400            
00401   xana::out(2) << "[update] endcap superclusters " << std::endl;
00402   
00403   XANASuperCluster *xscluster;    
00404   RecCollection<EgammaEndcapCluster> endcapSuperClusters(RecQuery("EGEFCluster","1.0"));
00405   RecCollection<EgammaEndcapCluster>::const_iterator endcapSuperCluster ;
00406   xana::out(5) << "[update] " << endcapSuperClusters.size() << endcapSuperClusterizerName_.value() << " endcapSuperClusters found." << std::endl;
00407   for (endcapSuperCluster=endcapSuperClusters.begin(); endcapSuperCluster!=endcapSuperClusters.end(); endcapSuperCluster++) {
00408     //    nbClusters_++;
00409     xana::out(2) << "[update] constructing new supercluster from algo " << endcapSuperCluster->seed()->WhichAlgo()
00410               << " at eta " << endcapSuperCluster->eta() << " with energy " <<
00411       endcapSuperCluster->energy() << std::endl;
00412     xscluster = new XANASuperCluster;
00413     xscluster->setEnergy(endcapSuperCluster->energy()); 
00414     xscluster->setPosition(endcapSuperCluster->position()); 
00415     xscluster->setSum1(endcapSuperCluster->seed()->getS1()); 
00416     xscluster->setSum4(endcapSuperCluster->seed()->getS4()); 
00417     xscluster->setSum9(endcapSuperCluster->seed()->getS9()); 
00418     xscluster->setSum25(endcapSuperCluster->seed()->getS25()); 
00419     xscluster->setHadronicOverEm(endcapSuperCluster->seed()->getHoe());
00420     // iso not available from endcap super clusters
00421     // compute from IsoCalculator
00422     float iso = ciso_->getIso(endcapSuperCluster->eta(),endcapSuperCluster->phi());
00423     xscluster->setCaloIsolation(iso); 
00424     xscluster->setDisc1(endcapSuperCluster->seed()->getDisc1()); 
00425     xscluster->setDisc2(endcapSuperCluster->seed()->getDisc2()); 
00426     xscluster->setDisc3(endcapSuperCluster->seed()->getDisc3()); 
00427     HepSymMatrix cov(2);
00428     cov(1,1) = endcapSuperCluster->seed()->getCovEE();
00429     cov(1,2) = endcapSuperCluster->seed()->getCovEP();
00430     cov(2,2) = endcapSuperCluster->seed()->getCovPP();
00431     xscluster->setPositionCovarianceMatrix(cov);
00432     xscluster->setEnergyScaleFactor(endcapSuperCluster->energy()/endcapSuperCluster->energyUncorrected());
00433     xana::out(2) << "[update] supercluster algoName is " << endcapSuperClusterizerName_.value() << std::endl;
00434     xscluster->setAlgoName(endcapSuperClusterizerName_.value());
00435     xana::out(2) << "[update] adding super cluster " << std::endl;
00436     xanaEsdEvent_->addSuperCluster(xscluster);
00437 
00438     // add seed basic cluster
00439     // the seed is the first in BasicClusters list
00440     EgammaVSuperCluster::const_iterator icf = endcapSuperCluster->basicClusters().first;
00441     HepPoint3D uncorrectedPositionAtFrontFace(0.,0.,0.);
00442     float esum=0;
00443     std::vector<const CaloRecHit *> rhits=(*icf)->rHits();
00444     for (std::vector<const CaloRecHit *>::const_iterator
00445            it=rhits.begin(); it != rhits.end(); it++) {
00446       esum +=   (*it)->energy();   
00447       uncorrectedPositionAtFrontFace += (*it)->energy() * (*it)->getMyCell().position();   
00448     }
00449     if (esum != 0.) uncorrectedPositionAtFrontFace /= esum ;    
00450     XANAEmCluster *seed 
00451       = new XANAEmCluster((*icf)->energy(),
00452                           (*icf)->position(),
00453                           uncorrectedPositionAtFrontFace,
00454                           (*icf)->chi2(),
00455                           clusterizerName_.value()); 
00456     xana::out(2) << "[update] adding seed " << std::endl;
00457     xanaEsdEvent_->addSeedCluster(seed, xscluster);
00458     delete seed;
00459 
00460     //  // now add brems clusters
00461     xana::out(2) << "[update] adding brems " << std::endl;
00462     XANAEmCluster *brem = 0;
00463     // beware, the endcapSuperCluster::vector<EgammaBasicCluster> list includes the seed cluster
00464     EgammaVSuperCluster::const_iterator ic = endcapSuperCluster->basicClusters().first+1;
00465     for (;ic != endcapSuperCluster->basicClusters().second; ic++) {
00466       HepPoint3D uncorrectedPositionAtFrontFace(0.,0.,0.);
00467       esum=0;
00468       rhits=(*ic)->rHits();
00469       for (std::vector<const CaloRecHit *>::const_iterator
00470         it=rhits.begin(); it != rhits.end(); it++) {
00471         esum +=         (*it)->energy();   
00472         uncorrectedPositionAtFrontFace += (*it)->energy() * (*it)->getMyCell().position();   
00473       }
00474       if (esum != 0.) uncorrectedPositionAtFrontFace /= esum ;                
00475               brem = new XANAEmCluster((*ic)->energy(),
00476                                       (*ic)->position(),
00477                                       uncorrectedPositionAtFrontFace,
00478                                       (*ic)->chi2(),
00479                                       clusterizerName_.value());            
00480       xana::out(2) << "[upDate] brem energy " << brem->getEnergy() << std::endl;
00481       xanaEsdEvent_->addBremCluster(brem, xscluster);
00482       delete brem;
00483     } 
00484     //  finally add preshower info when relevant (endcap)
00485     xana::out(2) << "[update] adding presh info " << std::endl;
00486     //  XANAPreshowerInfo *preshowerInfo = 0;
00487     // FIXME: currently not stored
00488 //      if (endcapSuperCluster->energy_P1()+endcapSuperCluster->energy_P2() != 0. ) {
00489 //           cout << "[update] constructing a  XANAPreshowerInfo" << std::endl;
00490 //        preshowerInfo = new XANAPreshowerInfo(endcapSuperCluster->energy_P1(), 
00491 //         endcapSuperCluster->energy_P2(), endcapSuperCluster->position_P1(),
00492 //         endcapSuperCluster->position_P2(), endcapSuperCluster->calibP1(), endcapSuperCluster->calibP2(), 
00493 //         endcapSuperCluster->calibSC());
00494 //        xanaEsdEvent_->addPreshCluster(preshowerInfo, xscluster);
00495 //      } 
00496     delete xscluster;
00497   }
00498 
00499 }
00500 
00501 
00502 void XANAEsdBuilder::setElectrons( G3EventProxy * ev ) 
00503 {
00504      
00505   xana::out(2) << "[update] electrons " << std::endl;
00506 
00507   RecCollection<ElectronCandidate> electrons(RecQuery("EGFElectron","1.0"));
00508   RecCollection<ElectronCandidate>::const_iterator electron;
00509   xana::out(5) << "[update] "<<electrons.size()<<" EGFElectron found " << std::endl;
00510   XANAElectronCandidate *xElectron;
00511   nbElectrons_=0;
00512   for (electron=electrons.begin(); electron!=electrons.end(); electron++) {
00513     nbElectrons_++ ;
00514     xana::out(1)<<"[OVAL] new electron "<<std::endl ;
00515     xElectron = new XANAElectronCandidate(); 
00516 
00517     // printing electron quantities from ElectronCandidate
00518     xana::out(8)<<"[OVAL] electron super cluster energy "<<electron->SCE()<<std::endl ;
00519     xana::out(8)<<"[OVAL] electron track px "<<electron->TRPx()<<std::endl ;
00520     xana::out(8)<<"[OVAL] electron track py "<<electron->TRPy()<<std::endl ;
00521     xana::out(8)<<"[OVAL] electron track pz "<<electron->TRPz()<<std::endl ;
00522     xana::out(8)<<"[OVAL] electron super cluster Et "<<electron->SCEt()<<std::endl ;
00523     xana::out(8)<<"[OVAL] electron super cluster eta "<<electron->SCEta()<<std::endl ;
00524     xana::out(8)<<"[OVAL] electron super cluster phi "<<electron->SCPhi()<<std::endl ;
00525     xana::out(8)<<"[OVAL] electron track eta "<<electron->TREta()<<std::endl ;
00526     xana::out(8)<<"[OVAL] electron track phi "<<electron->TRPhi()<<std::endl ;
00527     xana::out(8)<<"[OVAL] electron e/p "<<electron->ele_EoP()<<std::endl ;
00528     xana::out(8)<<"[OVAL] electron delta eta "<<electron->DeltaEta()<<std::endl ;
00529 
00530     const TTrack * eltrack=electron->GetTRElectron();
00531     const EgammaVSuperCluster *sc=electron->GetEGCandidate()->data();
00532     // calculate from sc and track
00533     HepPoint3D superClusterPosition;
00534     if (sc) {
00535       const EgammaCluster *seed = sc->seed();
00536 
00537       float superClusterEnergy = sc->energy();; 
00538       xana::out(8)<<"[OVAL] written electron super cluster energy "<<superClusterEnergy<<std::endl ;
00539       xElectron->setSuperClusterEnergy(superClusterEnergy);
00540       float theta = sc->theta();
00541       float x = sc->radius()*sin(theta)*cos(sc->phi());
00542       float y = sc->radius()*sin(theta)*sin(sc->phi());
00543       float z = sc->radius()*cos(theta);
00544       superClusterPosition=HepPoint3D(x,y,z);
00545       xana::out(8)<<"[OVAL] written electron super cluster position "<<superClusterPosition<<std::endl ;
00546       xElectron->setSuperClusterPosition(superClusterPosition);
00547       HepLorentzVector momentum;
00548       momentum.setPx(eltrack->momentumAtVertex().x());   
00549       momentum.setPy(eltrack->momentumAtVertex().y());   
00550       momentum.setPz(eltrack->momentumAtVertex().z());   
00551       momentum *= electron->SCE() / eltrack->momentumAtVertex().mag();
00552       momentum.setE(superClusterEnergy);        
00553       xana::out(8)<<"[OVAL] written electron momentum at vertex "<<momentum<<std::endl ;
00554       xElectron->setMomentumAtVertex(momentum); 
00555 
00556       // extrapolations from last track point
00557       TrajectoryStateOnSurface last = eltrack->outermostState(); 
00558       if (!last.isValid())
00559         throw Genexception("[writeESDForFamos::update], invalid tsos") ;
00560       // use 3D impactPointExtrapolator
00561       ImpactPointExtrapolator theExtrapolator ;
00562       // and an analytical (helix) propagator with material effect
00563       const MultipleScatteringUpdator aElectronMaterialEffectUpdator(0.000511);
00564       AnalyticalPropagator *ap = new AnalyticalPropagator(alongMomentum);
00565       PropagatorWithMaterial *thePropagator = new PropagatorWithMaterial(*(ap),aElectronMaterialEffectUpdator);
00566       delete ap;
00567 
00568       // extrapolation to super cluster
00569       GlobalPoint pos(x, y, z) ;
00570       TrajectoryStateOnSurface tssuper =
00571         theExtrapolator.extrapolate(*last.freeTrajectoryState(), pos, *thePropagator) ;  
00572       if (!tssuper.isValid()) tssuper=last;
00573       HepLorentzVector *trackMomentumAtCalo=new HepLorentzVector(tssuper.globalMomentum().x(),tssuper.globalMomentum().y(),tssuper.globalMomentum().z(),tssuper.globalMomentum().mag());
00574       xana::out(8)<<"[OVAL] written track momentum at super cluster "<<*trackMomentumAtCalo<<std::endl ;
00575       xElectron->setTrackMomentumAtCalo(*trackMomentumAtCalo);
00576       HepPoint3D *trackPositionAtCalo=new HepPoint3D(tssuper.globalPosition().x(),tssuper.globalPosition().y(),tssuper.globalPosition().z());
00577       xana::out(8)<<"[OVAL] written track position at super cluster "<<*trackPositionAtCalo<<std::endl ;
00578       xElectron->setTrackPositionAtCalo(*trackPositionAtCalo);
00579 
00580       float eSuperClusterOverP = -1.;   
00581       if (trackMomentumAtCalo->e()>0.) eSuperClusterOverP = superClusterEnergy / trackMomentumAtCalo->e();
00582       xana::out(8)<<"[OVAL] written E super cluster / p "<<eSuperClusterOverP<<std::endl ;
00583       xElectron->setESuperClusterOverP(eSuperClusterOverP);
00584       delete trackMomentumAtCalo;
00585 
00586       float deta = superClusterPosition.pseudoRapidity() - trackPositionAtCalo->pseudoRapidity();
00587       xana::out(8)<<"[OVAL] written deta (super - track) "<<deta<<std::endl ;
00588       xElectron->setDeltaEtaSuperClusterAtCalo(deta);
00589       float dphi=superClusterPosition.phi()-trackPositionAtCalo->phi();
00590       if (fabs(dphi) > pi) dphi = dphi < 0 ? 2.*pi + dphi : dphi - 2.*pi ;
00591       xana::out(8)<<"[OVAL] written dphi (super - track) "<<dphi<<std::endl ;
00592       xElectron->setDeltaPhiSuperClusterAtCalo(dphi);
00593       delete trackPositionAtCalo;
00594 
00595       // extrapolation to seed cluster
00596       Float_t thetas=2.0*atan(exp(-1*seed->Eta()));
00597       Float_t seedx=seed->Radius()*cos(seed->Phi())*sin(thetas);
00598       Float_t seedy=seed->Radius()*sin(seed->Phi())*sin(thetas);
00599       Float_t seedz=seed->Radius()*cos(thetas);
00600       GlobalPoint seedpos(seedx, seedy, seedz) ;
00601       TrajectoryStateOnSurface tsseed =
00602         theExtrapolator.extrapolate(*last.freeTrajectoryState(), seedpos, *thePropagator) ;  
00603       delete thePropagator;
00604       if (!tsseed.isValid()) tsseed=last;
00605       HepPoint3D *trackPositionAtSeed=new HepPoint3D(tsseed.globalPosition().x(),tsseed.globalPosition().y(),tsseed.globalPosition().z());
00606 
00607       float eSeedClusterOverP= -1. ;
00608       float thetaseed = 2.*atan(exp(-seed->Eta()));     
00609       float seedEnergy = 0.;
00610       if (thetaseed != 0.) seedEnergy = seed->Et() / sin(thetaseed) ;
00611       if( tsseed.globalMomentum().mag() > 0.) eSeedClusterOverP = seedEnergy / tsseed.globalMomentum().mag();
00612       xana::out(8)<<"[OVAL] written E seed cluster / p "<<eSeedClusterOverP<<std::endl ;
00613       xElectron->setESeedClusterOverP(eSeedClusterOverP);
00614 
00615       deta = seed->Eta() - trackPositionAtSeed->pseudoRapidity();
00616       xana::out(8)<<"[OVAL] written deta (seed - track) "<<deta<<std::endl ;
00617       xElectron->setDeltaEtaSeedClusterAtCalo(deta);
00618       dphi=seed->Phi()-trackPositionAtSeed->phi();
00619       if (fabs(dphi) > pi) dphi = dphi < 0 ? 2.*pi + dphi : dphi - 2.*pi ;
00620       xana::out(8)<<"[OVAL] written dphi (seed - track) "<<dphi<<std::endl ;
00621       xElectron->setDeltaPhiSeedClusterAtCalo(dphi);
00622       delete trackPositionAtSeed;
00623 
00624       // super cluster momentum, assuming perfect helix and vertex at origin
00625       GlobalPoint assumedVertex(0.,0.,0.);
00626       GlobalPoint measuredPoint(superClusterPosition.x(), superClusterPosition.y(), superClusterPosition.z()) ;                                             
00627       FTSFromVertexToPointFactory ftsFactory;
00628       FreeTrajectoryState ftsAtCalo = 
00629         ftsFactory(measuredPoint,assumedVertex,superClusterEnergy,eltrack->charge());
00630       HepLorentzVector superClusterMomentum;
00631       superClusterMomentum.setPx(ftsAtCalo.momentum().x());   
00632       superClusterMomentum.setPy(ftsAtCalo.momentum().y());   
00633       superClusterMomentum.setPz(ftsAtCalo.momentum().z());   
00634       superClusterMomentum.setE(superClusterEnergy);    
00635       xana::out(8)<<"[OVAL] written super cluster momentum  "<<superClusterMomentum<<std::endl ;
00636       xElectron->setSuperClusterMomentum(superClusterMomentum);
00637 
00638       float caloIsolation = ciso_->getIso(sc->eta(),sc->phi());
00639       xElectron->setCaloIsolation(caloIsolation);
00640     }
00641     int isTrackIsolated = electron->TRisISO();
00642     xElectron->setTrackIsolation(isTrackIsolated);
00643     float hadronicOverEm = electron->ele_HoE();
00644     xElectron->setHadronicOverEm(hadronicOverEm);
00645 
00646     xElectron->setAlgoName("EgammaCandidate");
00647 
00648     xana::out(2) << "[update] adding electron candidate No " << nbElectrons_ << std::endl;
00649     xanaEsdEvent_->addElectronCandidate(xElectron); 
00650 
00651     if (sc) {
00652       XANASuperCluster *xscluster=new XANASuperCluster();
00653       xscluster->setEnergy(sc->energy()); 
00654       xscluster->setPosition(sc->position()); 
00655       xscluster->setSum1(sc->seed()->getS1()); 
00656       xscluster->setSum4(sc->seed()->getS4()); 
00657       xscluster->setSum9(sc->seed()->getS9()); 
00658       xscluster->setSum25(sc->seed()->getS25()); 
00659       xscluster->setHadronicOverEm(sc->seed()->getHoe());
00660       // iso not available from endcap super clusters
00661       // compute from IsoCalculator
00662       float iso = ciso_->getIso(sc->eta(),sc->phi());
00663       xscluster->setCaloIsolation(iso); 
00664       xscluster->setDisc1(sc->seed()->getDisc1()); 
00665       xscluster->setDisc2(sc->seed()->getDisc2()); 
00666       xscluster->setDisc3(sc->seed()->getDisc3()); 
00667       HepSymMatrix cov(2);
00668       cov(1,1) = sc->seed()->getCovEE();
00669       cov(1,2) = sc->seed()->getCovEP();
00670       cov(2,2) = sc->seed()->getCovPP();
00671       xscluster->setPositionCovarianceMatrix(cov);
00672       xscluster->setEnergyScaleFactor(sc->energy()/sc->energyUncorrected());
00673       if ( fabs(sc->seed()->Eta())<1.4791 ) xscluster->setAlgoName(superClusterizerName_.value());
00674       else xscluster->setAlgoName(endcapSuperClusterizerName_.value());
00675       xanaEsdEvent_->addElectronSuperCluster(xscluster,xElectron);
00676       delete xscluster;
00677     }
00678 
00679 
00680       // electron seed
00681       XANAElectronSeed *xelseed= new XANAElectronSeed();
00682       const EPHLTElectronSeed * trseed = dynamic_cast<const EPHLTElectronSeed *>(&eltrack->seed().basicSeed()); 
00683       Short_t dir = 1;
00684       if (trseed->direction() == "oppositeToMomentum") dir = -1;
00685       xelseed->setDirection(dir);
00686       xelseed->setNumberOfHits(trseed->measurements().size());
00687       TrajectoryStateOnSurface tsseedfirst =
00688         trseed->measurements()[0].updatedState();
00689       HepPoint3D seedPositionAtFirst(tsseedfirst.globalPosition().x(),tsseedfirst.globalPosition().y(),
00690                                      tsseedfirst.globalPosition().z()); 
00691       HepVector3D seedMomentumAtFirst(tsseedfirst.globalMomentum().x(),tsseedfirst.globalMomentum().y(),
00692                                       tsseedfirst.globalMomentum().z());
00693       xelseed->setMomentumAtFirstPoint(seedMomentumAtFirst);
00694       xelseed->setPositionAtFirstPoint(seedPositionAtFirst);      
00695       TrajectoryStateOnSurface tsseedlast =
00696         trseed->measurements()[trseed->measurements().size()-1].updatedState();
00697       HepPoint3D seedPositionAtLast(tsseedlast.globalPosition().x(),tsseedlast.globalPosition().y(),
00698                                     tsseedlast.globalPosition().z()); 
00699       HepVector3D seedMomentumAtLast(tsseedlast.globalMomentum().x(),tsseedlast.globalMomentum().y(),
00700                                      tsseedlast.globalMomentum().z());
00701       xelseed->setMomentumAtLastPoint(seedMomentumAtLast);
00702       xelseed->setPositionAtLastPoint(seedPositionAtLast);
00703       xelseed->setAlgoName("Egammadefault");
00704 
00705       // we need an electron track just filled to make operator== work
00706       XANAElectronTrack *xTrack=new XANAElectronTrack;
00707       TrajectoryStateOnSurface tsi = eltrack->innermostState() ;
00708       TrajectoryStateOnSurface tso = eltrack->outermostState() ;
00709       HepPoint3D firstHit(tsi.globalPosition().x(),tsi.globalPosition().y(),
00710                           tsi.globalPosition().z()); 
00711       HepPoint3D lastHit(tso.globalPosition().x(),tso.globalPosition().y(),
00712                          tso.globalPosition().z()); 
00713       xTrack->setPositionAtFirstPoint(firstHit);
00714       xTrack->setPositionAtLastPoint(lastHit);
00715       
00716       float detavtx = superClusterPosition.pseudoRapidity() - XANAElectronTools::etaTransformation(xTrack->getMomentumAtFirstPoint().pseudoRapidity(),
00717                                                                                                 xTrack->getPositionAtFirstPoint().z(),xTrack->getPositionAtFirstPoint().perp()); 
00718       xElectron->setDeltaEtaSuperClusterAtVtx(detavtx);
00719       float dphivtx=superClusterPosition.phi()-XANAElectronTools::phiTransformation(xTrack->getMomentumAtFirstPoint().perp(),xTrack->getMomentumAtFirstPoint().pseudoRapidity(), 
00720                                                                          xTrack->getMomentumAtFirstPoint().phi(), xTrack->getCharge(),xTrack->getPositionAtFirstPoint().perp());
00721       if (fabs(dphivtx) > pi) dphivtx = dphivtx < 0 ? 2.*pi + dphivtx : dphivtx - 2.*pi ;
00722       xElectron->setDeltaPhiSuperClusterAtVtx(dphivtx);
00723       
00724       xanaEsdEvent_->addElectronSeed(xelseed,xTrack);
00725       xanaEsdEvent_->setElectronTrack(xTrack,xElectron);
00726       delete xelseed;
00727       delete xTrack;
00728  
00729       delete xElectron;
00730   }
00731 }
00732 
00733 void XANAEsdBuilderForFAMOS::setMcTruth( G3EventProxy * ev ) 
00734 {
00735 
00736   xana::out(2) << "[update] monte carlo truth " << std::endl;
00737   vector<XANAGeneratorParticle> particles;
00738 
00739   const RawHepEventParticle * aHepParticle;
00740   RawHepEvent& MyHepEvent = myEventMgr->event();
00741   
00742   XANAGeneratorParticle *aMcParticle=0;
00743   xana::out(5) << "[update] " << MyHepEvent.n()<< " generator particles in the event" <<std::endl;     
00744   // !! in RawHepEvent, i runs from 1 to n() !!
00745   for(unsigned int i=1; i<=MyHepEvent.n(); i++) {
00746     aHepParticle = MyHepEvent.getParticle(i);
00747     HepLorentzVector momentum(aHepParticle->px(), aHepParticle->py(),
00748      aHepParticle->pz(), aHepParticle->e());
00749      // fix case of back to back generation
00750      Short_t line = aHepParticle->line();
00751      if (MyHepEvent.n()==2 && aHepParticle->pid()<0) line=2 ;
00752      aMcParticle = new XANAGeneratorParticle(line, aHepParticle->status(),
00753      aHepParticle->pid(), aHepParticle->mother1(), aHepParticle->mother2(),
00754      aHepParticle->daughter1(), aHepParticle->daughter2(), momentum,
00755      aHepParticle->mass(), aHepParticle->vertex());
00756     // assuming signal particles are the first ones
00757     // !! in RawHepEvent, i runs from 1 to n() !!
00758     xana::out(2) << "[update] adding generator particle " << std::endl;
00759     if (i<=MyHepEvent.nSignal()) {
00760       xanaGenEvent_->addSignalParticle(aMcParticle);
00761     } else {
00762       xanaGenEvent_->addPileupParticle(aMcParticle);
00763     }
00764     delete aMcParticle;
00765   }
00766   // GEANT MC INFO
00767 
00768   vector<FEmbdSimVertex> GeantVertices; // this is an STL vector
00769   GeantVertices = *(myEventMgr->simSignal().vertices()); // get the vertex container
00770   vector<FEmbdSimTrack> GeantTracks; // this is an STL vector
00771   GeantTracks = *(myEventMgr->simSignal().tracks());
00772   xana::out(5) << "[update] " << GeantTracks.size()<< " geant tracks from true MC event." <<std::endl;     
00773   xana::out(5) << "[update] " << GeantVertices.size()<< " geant vertices from true MC event." <<std::endl;     
00774 
00775 /*
00776   vector<FEmbdGenParticle> GenParts; 
00777   GenParts = *(myEventMgr->simSignal().genparts());
00778 
00779   XANAGeantVertex *vtx=0;
00780   XANAGeantTrack *trk=0;
00781   int vertex=0;
00782 
00783   // UB, 03/01/05  moved into initialisation:
00784   //  xanaGeantEvent_->addGeneratorEvent(xanaGenEvent_);
00785 
00786   // to fix case where no GeantTracks[it] has link to genpart
00787   bool noGenParent = true;
00788   for(unsigned int it=0; it<GeantTracks.size(); it++) {
00789     if (! GeantTracks[it].noGenpart()) {
00790     noGenParent = false;
00791     break;
00792     }
00793   }
00794   xana::out(1) << "[update] noGenParent " << noGenParent << endl; 
00795 
00796 
00797   xana::out(2) << "[update] before geant track loop " << std::endl;     
00798   for(unsigned int it=0; it<GeantTracks.size(); it++) {
00799     // decode mother-daughter relationship
00800     xana::out(2) << "[update] new geant track " << GeantTracks[it].momentum() << std::endl;
00801     HepLorentzVector mom(GeantTracks[it].momentum());
00802     int mothindex = 0 ;
00803     if (mom.e() < 0.) {
00804       mothindex =  (int)(-GeantTracks[it].momentum().e()) ;
00805       xana::out(2) << "[update] track has daughter " << mothindex << std::endl;
00806       double energy2 = mom.vect().mag2();
00807       int type = GeantTracks[it].type();
00808       // the mass for the geantinos should be the mass of the corresponding particle
00809       // in the vertex
00810       // if there is a geantino  =>  there is a vertex AND
00811       //                             the geantino is the first particle
00812       //                             of that vertex by construction
00813       if (type==-103) type = GeantTracks[it+1].type();
00814       double mass = HepPDT::theTable().getParticleData(GeantTracks[it].type())->mass();
00815       energy2 += mass*mass ;
00816       double energy = sqrt(energy2) ;
00817       mom.setE(energy);
00818     }
00819     else {
00820       //IN STANDARD EVENTS WITHOUT  ANY MODIFICATION!!
00821       mothindex =  SimTrack(ev->simSignal()->track(it)).vertex().parent().id() + 1;
00822       xana::out(2) << "[update] track has daughter " << mothindex << std::endl;
00823     }
00824     xana::out(2) << "[update] recomputed geant track momentum " << mom << std::endl;
00825     // now create XANA track
00826     trk = new XANAGeantTrack(it+1,GeantTracks[it].type(), mom, mothindex);
00827     if (! GeantTracks[it].noVertex()) {
00828       vertex =  GeantTracks[it].vertIndex();
00829       xana::out(1) << "[update] vertex index is " << vertex << std::endl; 
00830       vtx = new XANAGeantVertex(GeantVertices[vertex].position());
00831       xanaGeantEvent_->addGeantVertex(vtx);
00832       xanaGeantEvent_->addGeantTrack(trk, vtx); 
00833       delete vtx;
00834     } else
00835       xanaGeantEvent_->addGeantTrack(trk);              
00836 
00837     xana::out(2) << "[update] linking geant tracks to generator particle" << std::endl;
00838 
00839     // to link XANAGeantTrack to XANAGeneratorParticle
00840 
00841     XANAGeneratorParticle *parent;
00842     xana::out(2) << "[update] track gen part index " << GeantTracks[it].genpartIndex() << " pid " << GeantTracks[it].type() << " noGenpart " << GeantTracks[it].noGenpart() << endl;
00843     if (! GeantTracks[it].noGenpart()) {
00844       int genpart = GeantTracks[it].genpartIndex();
00845       // create generator particle with available info
00846       parent = new XANAGeneratorParticle(genpart+1, GenParts[genpart].status(),
00847        GenParts[genpart].pid(), GenParts[genpart].mother1(), GenParts[genpart].mother2(),
00848        GenParts[genpart].daughter1(), GenParts[genpart].daughter2(), 
00849        GenParts[genpart].fourmomentum(), 
00850        GenParts[genpart].fourmomentum().m(), GeantVertices[vertex].position());
00851 //        parent->print();
00852       xanaGeantEvent_->setGeneratorParticle(trk, parent);
00853       delete parent;
00854     }
00855     // to fix case where no GeantTracks[it] has link to genpart:
00856     // geant tracks coming from generator particle have vertex number 0
00857     // thus for those tracks, set generator particle by hand
00858     if (noGenParent && GeantTracks[it].vertIndex()==0) {                
00859       xana::out(2) << "[update] fixing parent particle relationship to geant track "<< std::endl;
00860       int igen=0;
00861       // single particle generation
00862       if (MyHepEvent.n()==1) igen=1;
00863       // back-to-back particle/antiparticle generation
00864       if (MyHepEvent.n()==2) {
00865         if (GeantTracks[it].type() == MyHepEvent.getParticle(1)->pid()) igen=1;
00866         else igen=2;
00867       }   
00868       // create generator particle with available info
00869       if (igen!=0) {
00870         xana::out(2) << "[update] creating link to generator particle " << std::endl;
00871         aHepParticle = MyHepEvent.getParticle(igen);
00872         HepLorentzVector momentum(aHepParticle->px(), aHepParticle->py(),
00873          aHepParticle->pz(), aHepParticle->e());
00874         // in case of back to back, both particles comes from HepEvent line 1 ..
00875         parent = new XANAGeneratorParticle(igen, 1,
00876          aHepParticle->pid(), 0, 0, 0, 0, momentum, aHepParticle->mass(), aHepParticle->vertex());
00877         xana::out(2) << "[update] setting parent particle to geant track " << it << std::endl;
00878 //          parent->print();
00879         xanaGeantEvent_->setGeneratorParticle(trk, parent);
00880         delete parent;
00881       }
00882     }
00883 
00884     // to link XANAGeantTrack to XANACaloRecHit
00885 //    if (writeCaloRecHits_) SHDCaloRegion allCellsRegion(allCells_);
00886 //    SHDCaloRegion clusterizedCellsRegion(clusterizedCells_);
00887     
00888     //link only geantino geant tracks
00889     if (GeantTracks[it].type()==-103) {
00890       // first retrieve impact point if any for this geant track
00891       xana::out(3) << "[update] new geantino" << std::endl;
00892       xana::out(3) << "[update] impacting track is " << GeantTracks[it+1] << std::endl;
00893       if (! GeantTracks[it].noVertex()) {
00894         vertex =  GeantTracks[it].vertIndex();
00895         xana::out(3) << "[update] impact position is " << GeantVertices[vertex].position().vect() << std::endl; 
00896       } 
00897       // then search impacted CellID 
00898       CellID theImpactCell;
00899       // in allCells first, if writeCaloRecHits_ is true
00900       if (writeCaloRecHits_) {
00901 // version avec les bases
00902       if (fabs(GeantVertices[vertex].position().z())<300.) {
00903         theImpactCell=CaloRegistery::getBase("EBRY")->getClosestCell(GeantVertices[vertex].position().vect());
00904       } else {
00905         theImpactCell=CaloRegistery::getBase("EFRY")->getClosestCell(GeantVertices[vertex].position().vect());
00906       }
00907 // version avec SHDCaloRegion
00908 // ne trouve pas toujours la CellID
00909 //      theImpactCell = allCellsRegion.whichCellId(GeantVertices[vertex].position().vect());
00910       // then add link
00911       if ( theImpactCell.isZero() ){
00912         xana::out(3) << "[update] CellID not found !" <<  std::endl;
00913         xana::out(3) << "[update] Impact Out Of Calorimeter " << std::endl;
00914         continue ;
00915        } else {
00916         xana::out(3) << "[update] CellID found "<< theImpactCell.position() <<std::endl;
00917         map<CellID, XANACaloRecHit, less<CellID> >::iterator xhititr  = allCaloRecHits_.find(theImpactCell);
00918         if ( xhititr != allCaloRecHits_.end() ) {
00919           xana::out(3) << "[update] corresponding CaloRecHit found " <<  std::endl;
00920           xana::out(3) << "[update] calo rechit energy " << xhititr->second.getEnergy() << std::endl;
00921         }  
00922         //      XANAEsdEvent_->setCaloRecHit(trk, theImpactCell);
00923         xanaEsdEvent_->setCaloRecHit(&(xhititr->second), trk);
00924        }
00925        } 
00926        // then search in clusterizedCells only
00927     }
00928 
00929     delete trk;
00930 
00931   }
00932 
00933   xana::out(2) << "[update] printing geant event " << std::endl;
00934 //      xanaGeantEvent_->print();
00935 
00936 */
00937 
00938 }
00939 
00940 void XANAEsdBuilder::setMcTruth( G3EventProxy * ev ) 
00941 { }
00942 
00943 void XANAEsdBuilder::setElectronTracks(G3EventProxy *)
00944 {// fill electron track container for normal electron tracks
00945 //  RecCollection<TTrack> tracks( RecQuery("EPTracks","1.0"));
00946   RecCollection<TTrack> tracks( RecQuery("EPFTracks","1.0"));
00947   RecCollection<TTrack>::const_iterator eltrack;
00948   xana::out(5) << "[update] "<<tracks.size()<<" electron tracks found " << std::endl;
00949   for (eltrack=tracks.begin(); eltrack!=tracks.end(); eltrack++) {
00950     XANAElectronTrack *xTrack = new XANAElectronTrack;
00951     TrajectoryStateOnSurface tsi = eltrack->innermostState() ;
00952     TrajectoryStateOnSurface tso = eltrack->outermostState() ;
00953     HepVector3D momentumAtVertex(eltrack->momentumAtVertex().x(), 
00954                                  eltrack->momentumAtVertex().y(),eltrack->momentumAtVertex().z());
00955     HepPoint3D firstHit(tsi.globalPosition().x(),tsi.globalPosition().y(),
00956                         tsi.globalPosition().z()); 
00957     HepVector3D momentumAtFirst(tsi.globalMomentum().x(),tsi.globalMomentum().y(),
00958                                   tsi.globalMomentum().z());
00959       HepPoint3D lastHit(tso.globalPosition().x(),tso.globalPosition().y(),
00960                          tso.globalPosition().z()); 
00961       HepVector3D momentumAtLast(tso.globalMomentum().x(),tso.globalMomentum().y(),
00962                                  tso.globalMomentum().z());
00963       xTrack->setCharge(eltrack->charge()); 
00964       xTrack->setChi2OverDof(eltrack->normalisedChiSquared());
00965       xTrack->setNumberOfHits(eltrack->foundHits());
00966       xTrack->setNumberOfLostHits(eltrack->lostHits());
00967       xTrack->setImpactParameter(eltrack->impactParameter3D().value());
00968       xTrack->setLongImpactParameter(eltrack->zImpactParameter().value());
00969       xTrack->setTransImpactParameter(eltrack->transverseImpactParameter().value());
00970       xTrack->setMomentumAtVertex(momentumAtVertex);
00971       xTrack->setPositionAtFirstPoint(firstHit);
00972       xTrack->setMomentumAtFirstPoint(momentumAtFirst);
00973       xTrack->setPositionAtLastPoint(lastHit);
00974       xTrack->setMomentumAtLastPoint(momentumAtLast);
00975       xTrack->setAlgoName("EPFTrack"); 
00976       //      xanaEsdEvent_->addElectronTrack(xTrack,xElectron);
00977       xanaEsdEvent_->addElectronTrack(xTrack);
00978 
00979       // hits (if requested)
00980       // must be done after having added electron track!
00981       if (writeEleTrackHits_) {
00982         std::vector<TrajectoryMeasurement> MyEleTrMeas = eltrack->measurements();
00983         setEleTrackHits(MyEleTrMeas,xTrack);
00984       }
00985       delete xTrack;
00986   }
00987 }
00988 
00989 
00990 void XANAEsdBuilder::setEleTrackHits(std::vector<TrajectoryMeasurement> MyEleTrMeas, XANAElectronTrack *xTrack) 
00991 {
00992   // add electron track hits for this track, if asked for
00993 
00994   //     std::vector<TrajectoryMeasurement> MyEleTrMeas = eltrack->measurements();
00995      std::vector<TrajectoryMeasurement>::const_iterator IterEleTrMeas;
00996 
00997      xana::out(2) << "[update] adding electronic track hits " << std::endl; 
00998      xana::out(2) << "[update] number of measurements " <<MyEleTrMeas.size() << std::endl;
00999      //     xana::out(2) << "[update] number of valid hits "<<eltrack->foundHits()<< std::endl;
01000 
01001      XANATrackHit *elethit=0;
01002      int nbhits=0;
01003      for (IterEleTrMeas=MyEleTrMeas.begin(); IterEleTrMeas!=MyEleTrMeas.end(); IterEleTrMeas++)
01004        {
01005         // taking the rec hit 
01006         RecHit IterEleRecHit = IterEleTrMeas->recHit();                
01007         // checking if it is valid or not  
01008         if( !(IterEleRecHit.isValid()) ){xana::out(2) << "[update] non valid hit" << std::endl;}
01009 
01010         // taking the updated state 
01011         TrajectoryStateOnSurface IterEleTSN = IterEleTrMeas->updatedState();                
01012         // checking if it is valid or not  
01013         if( !(IterEleTSN.isValid()) ){xana::out(2) << "[update] non valid updated state" << std::endl;}
01014 
01015 
01016         // analysis only if both are valid 
01017         if (!(IterEleRecHit.isValid()) && (IterEleTSN.isValid())){xana::out(2)<<"[update] non valid hit but valid updated state"<< std::endl;}
01018         if ((IterEleRecHit.isValid()) && !(IterEleTSN.isValid())){xana::out(2)<<"[update] non valid updated state but valid hit"<< std::endl;}
01019         if(  IterEleRecHit.isValid() && IterEleTSN.isValid() )
01020         {
01021           HepPoint3D elehitposition( (IterEleRecHit).globalPosition().x(), 
01022                                      (IterEleRecHit).globalPosition().y(), 
01023                                      (IterEleRecHit).globalPosition().z());
01024 
01025           HepPoint3D eleupdposition( (IterEleTSN).globalPosition().x(),
01026                                      (IterEleTSN).globalPosition().y(),
01027                                      (IterEleTSN).globalPosition().z());
01028 
01029           HepVector3D eleupdmomentum( (IterEleTSN).globalMomentum().x(),
01030                                       (IterEleTSN).globalMomentum().y(),
01031                                       (IterEleTSN).globalMomentum().z());
01032 
01033           elethit = new XANATrackHit( elehitposition, eleupdposition, eleupdmomentum ); 
01034           elethit->setIsStereo((IterEleRecHit).det().detUnits().front()->type().isStereo());
01035           Bool_t isBarrel = (((IterEleRecHit).det().detUnits().front()->type().part() == barrel) ? true : false);
01036           elethit->setIsBarrel(isBarrel);
01037 
01038             xana::out(1) << "[update] ELE - NEW STEP" << std::endl; 
01039             xana::out(1) << "[update] new hit" << std::endl;    
01040             xana::out(1) << (IterEleRecHit).globalPosition().x() << std::endl;     
01041             xana::out(1) << (IterEleRecHit).globalPosition().y() << std::endl;          
01042             xana::out(1) << (IterEleRecHit).globalPosition().z() << std::endl;          
01043             xana::out(1) << "[update] new updated state: position" << std::endl;        
01044             xana::out(1) << (IterEleTSN).globalPosition().x() << std::endl;     
01045             xana::out(1) << (IterEleTSN).globalPosition().y() << std::endl;     
01046             xana::out(1) << (IterEleTSN).globalPosition().z() << std::endl;     
01047             xana::out(1) << "[update] new updated state: momentum" << std::endl;        
01048             xana::out(1) << (IterEleTSN).globalMomentum().x() << std::endl;     
01049             xana::out(1) << (IterEleTSN).globalMomentum().y() << std::endl;     
01050             xana::out(1) << (IterEleTSN).globalMomentum().z() << std::endl;     
01051 
01052           xanaEsdEvent_->addEleTrackHit(elethit, xTrack);           
01053           ++nbhits;
01054           delete elethit;
01055 
01056        } // valid hit and upd state
01057        } // loop on measurements
01058      xana::out(2) << "[update] number of hits that were added "<<nbhits<< std::endl;
01059 }
01060 
01061 void XANAEsdBuilder::writeTree() 
01062 {
01063 
01064   xana::out(2) << "[update] writing event " << std::endl;     
01065   
01066   xanaEsdTree_->Fill();  //fill the tree
01067 
01068 }
01069 
01070 void XANAEsdBuilder::upDate (G3SetUp *) 
01071 {
01072 /*
01073     xana::out(5) << "Starting initSetUp==================================" << std::endl;
01074     if (strcmp(electronTrackAlgo_.value().c_str(),"GSF"))     return;
01075     
01076     // Tracker properties
01077        LayerMediumPropertySetter MediumSetter(FullTracker::instance()->allLayers());
01078        //       const vector<DetLayer*>& aLayerContainer=FullTracker::instance()->allLayers();
01079        // 
01080        //    LayerMediumPropertySetter MediumSetter(aLayerContainer);
01081     // -----------------------------------
01082     // Track construction
01083 
01084     // Propagators - gsf   
01085     GsfMaterialEffectsUpdator* gsfMEUpdator_ = 
01086        GsfMaterialEffectsFactory::getBuilder()->constructComponent(0);
01087 
01088     theFwdGsfPropagator_ =    
01089     new GsfPropagatorWithMaterial(AnalyticalPropagator(alongMomentum),
01090                                  *gsfMEUpdator_);
01091 
01092     theBwdGsfPropagator_ = 
01093     new GsfPropagatorWithMaterial(AnalyticalPropagator(oppositeToMomentum),
01094                                     *gsfMEUpdator_);
01095     delete gsfMEUpdator_;
01096 
01097 
01098     // Propagator - kf
01099     MaterialEffectsUpdator* kfMEUpdator_ =
01100       MaterialEffectsFactory().constructComponent();
01101 
01102     theFwdKfPropagator_ = 
01103       new PropagatorWithMaterial(AnalyticalPropagator(alongMomentum),
01104                                  *kfMEUpdator_);
01105     
01106     // Updators 
01107     theKfUpdator_  = new KFUpdator();
01108     theGsfUpdator_ = new GsfMultiStateUpdator();   
01109 
01110     // Measurement estimators
01111     theKfFitterMEstimator_    = new Chi2MeasurementEstimator(150.);
01112     theGsfFitterMEstimator_   = new GsfChi2MeasurementEstimator(150.);   
01113     theGsfSmootherMEstimator_ = new GsfChi2MeasurementEstimator(1.e10);
01114     
01115     // Fitters 
01116     theKfFitter_ = new KFTrajectoryFitter(*theFwdKfPropagator_,
01117                                          *theKfUpdator_,
01118                                          *theKfFitterMEstimator_);
01119 
01120    
01121     theGsfFitter_ = new GsfTrajectoryFitter(*theFwdGsfPropagator_,
01122                                            *theGsfUpdator_,
01123                                            *theGsfFitterMEstimator_);
01124 
01125     // Smoother
01126     theGsfSmoother_ = new GsfTrajectorySmoother(*theBwdGsfPropagator_,*theGsfUpdator_,
01127                                                *theGsfSmootherMEstimator_);
01128 
01129     // Fitting 
01130     theGsfFittingSmoother_forwardKf_  = new KFFittingSmoother(theKfFitter_,theGsfSmoother_);
01131     theGsfFittingSmoother_forwardGsf_ = new KFFittingSmoother(theGsfFitter_,theGsfSmoother_);
01132 
01133 
01134     // Trajectory builder
01135     RecQuery GsfBuilderQuery("CombinatorialTrajectoryBuilder");
01136     GsfBuilderQuery.setParameter("chiSquarCut",100000.);
01137     GsfBuilderQuery.setParameter("mass",0.000511);
01138     GsfBuilderQuery.setParameter("maxCand",2);
01139     GsfBuilderQuery.setParameter("maxLostHit",1);
01140     GsfBuilderQuery.setParameter("maxConsecLostHit",1);
01141     GsfBuilderQuery.setParameter("lostHitPenalty",50.);
01142     GsfBuilderQuery.setParameter("intermediateCleaning",true);
01143     GsfBuilderQuery.setParameter("minimumNumberOfHits",5);
01144     GsfBuilderQuery.setParameter("ptCut", 3.);
01145     GsfBuilderQuery.setParameter("alwaysUseInvalidHits",true);
01146     myGsfTrajectoryBuilder_ = createAlgo<TrajectoryBuilder>(RecConfig(GsfBuilderQuery));
01147 
01148     // Track finder     
01149     if (!strcmp(gsfForwardFit_.value().c_str(),"KF") ) 
01150       {
01151        theGsfTrackFinder_ = 
01152        new ModularKfReconstructor(new CombinatorialSeedGeneratorFromPixel(1.0,0.1,15),
01153                                   myGsfTrajectoryBuilder_,
01154                                   theGsfFittingSmoother_forwardKf_,
01155                                   new TrajectoryCleanerBySharedHits);
01156       }
01157 
01158     if (!strcmp(gsfForwardFit_.value().c_str(),"GSF") ) 
01159        {
01160        theGsfTrackFinder_ = 
01161        new ModularKfReconstructor(new CombinatorialSeedGeneratorFromPixel(1.0,0.1,15),
01162                                   myGsfTrajectoryBuilder_,
01163                                   theGsfFittingSmoother_forwardGsf_,
01164                                   new TrajectoryCleanerBySharedHits);
01165       }
01166 
01167      // Extrapolator
01168      theTIPExtrapolatorGSF_ = new TransverseImpactPointExtrapolator(*theBwdGsfPropagator_);
01169 */
01170 }
01171 
01172 
01173 std::string XANAEsdBuilder::getAlgoName(int algo) {
01174 
01175   if (algo==1) return std::string("EgammaBasicIsland");
01176   else if (algo>=100) return std::string("EgammaBasicHybrid");
01177   else return std::string("");
01178 
01179 }      
01180  
01181  
01182 static PKBuilder<XANAEsdBuilderForFAMOS> eventObserver("XANAEsdBuilderForFAMOS") ;