writeESD.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/writeESD.h>
00007 
00008 #include <XANADOO/domain/test/setHeader.icc>
00009 #include <XANADOO/domain/test/setAllTrackHits.icc>
00010 #include <XANADOO/domain/test/setCaloRecHits.icc>
00011 #include <XANADOO/domain/test/setElectrons.icc>
00012 #include <XANADOO/domain/test/setElectronsGSF.icc>
00013 #include <XANADOO/domain/test/setElectronTracks.icc>
00014 #include <XANADOO/domain/test/setEmClusters.icc>
00015 #include <XANADOO/domain/test/setEndcapSuperClusters.icc>
00016 #include <XANADOO/domain/test/setHadClusters.icc>
00017 #include <XANADOO/domain/test/setJets.icc>
00018 #include <XANADOO/domain/test/setMcTruth.icc>
00019 #include <XANADOO/domain/test/setMet.icc>
00020 #include <XANADOO/domain/test/setMuons.icc>
00021 #include <XANADOO/domain/test/setSuperClusters.icc>
00022 #include <XANADOO/domain/test/setTracks.icc>
00023 #include <XANADOO/domain/test/setTriggerInfo.icc>
00024 #include <XANADOO/domain/test/setVertices.icc>
00025 
00026 #include <XANADOO/domain/interface/XANAEsdEvent.h>
00027 #include <XANADOO/XANAClusters/interface/XANASuperCluster.h>
00028 #include <XANADOO/XANAClusters/interface/XANAHadCluster.h>
00029 #include <XANADOO/XANATracks/interface/XANATrackHit.h>
00030 #include <XANADOO/XANAVertex/interface/XANAVertex.h>
00031 #include <XANADOO/XANAMcInfo/interface/XANAGeneratorEvent.h>
00032 #include <XANADOO/XANAMcInfo/interface/XANAGeneratorParticle.h>
00033 #include <XANADOO/XANAMcInfo/interface/XANAGeantEvent.h>
00034 #include <XANADOO/XANAMcInfo/interface/XANAGeantVertex.h>
00035 #include <XANADOO/XANAMcInfo/interface/XANAGeantTrack.h>
00036 #include <XANADOO/XANATriggerInfo/interface/XANATriggerInfo.h>
00037 #include <XANADOO/XANAUtilities/interface/XANAToolkit.h>
00038 #include <XANADOO/XANAUtilities/interface/XANAMemLog.h>
00039 #include <XANADOO/XANAElectronCandidate/interface/XANAElectronCandidate.h>
00040 #include <XANADOO/XANAElectronCandidate/interface/XANAElectronGSFTrack.h>
00041 #include <XANADOO/XANAElectronCandidate/interface/XANAElectronSeed.h>
00042 #include <XANADOO/XANAMuons/interface/XANAMuonCandidate.h>
00043 #include <XANADOO/XANAMuons/interface/XANAMuonTrack.h>
00044 #include <XANADOO/XANAJets/interface/XANAJet.h>
00045 #include <XANADOO/XANAMet/interface/XANAMet.h>
00046 
00047 #include <ElectronPhoton/ClusterTools/interface/IsoCalculator.h>
00048 
00049 #include <TrackerReco/TkRecoDebugTools/interface/LayerMediumPropertySetter.h>
00050 #include <TrackerReco/GsfPattern/interface/GsfMaterialEffectsFactory.h>
00051 #include <TrackerReco/GsfPattern/interface/GsfMultiStateUpdator.h>
00052 #include <TrackerReco/GsfPattern/interface/GsfPropagatorWithMaterial.h>
00053 #include <TrackerReco/GtfPattern/interface/CombinatorialTrajectoryBuilder.h>
00054 #include <TrackerReco/GsfPattern/interface/GsfChi2MeasurementEstimator.h>
00055 #include <TrackerReco/GsfPattern/interface/GsfTrajectoryFitter.h>
00056 #include <TrackerReco/GtfPattern/interface/ModularKfReconstructor.h>
00057 #include <TrackerReco/TkSeedGenerator/interface/CombinatorialSeedGeneratorFromPixel.h>
00058 #include <TrackerReco/GsfPattern/interface/GsfTrajectorySmoother.h>
00059 
00060 #include <CommonReco/TrackFitters/interface/KFFittingSmoother.h>
00061 #include <CommonReco/PatternTools/interface/TransverseImpactPointExtrapolator.h>
00062 #include <CommonReco/GeomPropagators/interface/AnalyticalPropagator.h>
00063 #include <CommonReco/MaterialEffects/interface/PropagatorWithMaterial.h>
00064 #include <CommonReco/KalmanUpdators/interface/KFUpdator.h>
00065 #include <CommonReco/KalmanUpdators/interface/Chi2MeasurementEstimator.h>
00066 #include <CommonReco/TrackFitters/interface/KFFittingSmoother.h>
00067 #include <CommonReco/TrackFitters/interface/KFTrajectoryFitter.h>
00068 #include <CommonReco/TrackFitters/interface/KFTrajectorySmoother.h>
00069 #include <CommonReco/MaterialEffects/interface/MaterialEffectsFactory.h>
00070 #include <CommonReco/PatternTools/interface/TrajectoryCleanerBySharedHits.h>
00071 
00072 #include <CommonDet/BasicDet/interface/DetUnit.h>
00073 #include <Tracker/TkLayout/interface/CmsTracker.h>
00074 #include <Tracker/TkLayout/interface/FullTracker.h>
00075 #include <CommonDet/BasicDet/interface/DetType.h>
00076 #include <CommonDet/BasicDet/interface/Readout.h>
00077 
00078 #include <CARF/G3Event/interface/G3SetUp.h>
00079 #include <CARF/Reco/interface/AutoRecCollection.h>
00080 #include <CARF/Reco/interface/RecQuery.h>
00081 #include <CARF/Reco/interface/ConfigAlgoFactory.h>
00082 #include <CARF/Reco/interface/RecBuilder.h>
00083 
00084 #include <Utilities/Notification/interface/PackageInitializer.h>
00085 #include <Utilities/Notification/interface/Singleton.h>
00086 
00087 #include <ClassReuse/GeomVector/interface/GlobalVector.h>
00088 
00089 
00090 #include <TROOT.h>
00091 #include <TObjectTable.h>
00092 #include <TFile.h>
00093 #include <TTree.h>
00094 #include <TBranch.h>
00095 #include <TClonesArray.h>
00096 #include <TBits.h>
00097 
00098 #include <CLHEP/Units/PhysicalConstants.h>
00099 
00100 #include <vector>
00101 #include <map>
00102 #include <string>
00103 #include <iostream> 
00104 #include <iomanip>  
00105 #include <numeric>
00106 
00107 
00108 XANAEsdBuilder::XANAEsdBuilder(): clusters(0), superClusters(0), endcapSuperClusters(0), towers(0), tracks(0), vertices(0),
00109    muons(0), electrons(0), eltracks(0), jets(0), mets(0), hlt(0), xanaEsdTree_(0),  xanaEsdEventBranch_(0), xanaEsdFile_(0), 
00110    xanaEsdEvent_(0), xanaGenEventBranch_(0), xanaGenEvent_(0), xanaGeantEventBranch_(0),
00111    xanaGeantEvent_(0), ciso_(0), xanamemlog_(0), nbClusters_(0), nbTracks_(0), nbElectrons_(0)
00112    
00113 { 
00114 
00115    // initialize observer    
00116    xana::out(1) << "Constructing an XANAEsdBuilder" << endl ;
00117    Observer<G3EventProxy *>::init() ; 
00118    Observer<G3SetUp *>::init() ; 
00119 
00120    // algos configuration
00121    clusterizerName_ = 
00122     SimpleConfigurable<std::string>("EGBCluster","XANADOO:Clusterizer");
00123    superClusterizerName_ =
00124     SimpleConfigurable<std::string>("EGSCluster","XANADOO:SuperClusterizer");
00125    endcapSuperClusterizerName_ =
00126     SimpleConfigurable<std::string>("EGECluster","XANADOO:EndcapSuperClusterizer");
00127    hcalClusterizerName_ =
00128     SimpleConfigurable<std::string>("TowerBuilder","XANADOO:HcalClusterizer");
00129    trackReconstructorName_=
00130     SimpleConfigurable<std::string>("TransientCombinatorialTrackFinder","XANADOO:TrackReconstructor");
00131    vertexReconstructorName_=
00132     SimpleConfigurable<std::string>("PrincipalVertexFinder","XANADOO:VertexReconstructor");
00133    muonReconstructorName_=
00134     SimpleConfigurable<std::string>("GlobalMuonReconstructor","XANADOO:MuonReconstructor");
00135    electronReconstructorName_=
00136     SimpleConfigurable<std::string>("GSF","XANADOO:ElectronReconstructor");
00137    gsfForwardFit_=
00138     SimpleConfigurable<std::string>("GSF","XANADOO:gsfForwardFit");
00139    jetsFromTrueInformation_ = 
00140     SimpleConfigurable<bool>(false,"XANADOO:jetsFromTrueInformation");
00141    metsFromTrueInformation_ = 
00142     SimpleConfigurable<bool>(false,"XANADOO:metsFromTrueInformation");
00143    // get L1 trigger config from Data/TriggerData/L1GlobalTrigger/$lumichoice/trigger_menu.dat
00144    l1Menu_ = 
00145     SimpleConfigurable<std::string>("lumi2x1033","L1Globaltrigger:TriggerMenu");
00146    // set HLT trigger config from xml file
00147    hltMenu_ = 
00148     SimpleConfigurable<std::string>("2x1033HLT.xml","HighLevelTriggerXML:XMLfile");
00149 
00150 
00151    // xanadoo steering
00152    writeTree_ = SimpleConfigurable<bool>(false,"XANADOO:writeEsd");
00153    writeCaloRecHits_ = SimpleConfigurable<bool>(false,"XANADOO:EsdWithCaloRecHits");
00154    writeTrackHits_ = SimpleConfigurable<bool>(false,"XANADOO:EsdWithTrackHits");   writeAllTrackHits_ = SimpleConfigurable<bool>(false,"XANADOO:EsdWithAllTrackHits");
00155    writeEleTrackHits_ = SimpleConfigurable<bool>(false,"XANADOO:EsdWithEleTrackHits");
00156 
00157      // memory log
00158    memlog_ = SimpleConfigurable<bool>(false,"XANADOO:memlog");
00159    if (memlog_) xanamemlog_=new XANAMemLog();
00160 
00161    // creates output file
00162    fileName_ = SimpleConfigurable<std::string>("single_eminus_pt100.root","XANADOO:FileName");
00163    xanaEsdFile_ = new TFile(fileName_.value().c_str(),"RECREATE");
00164    compressionLevel_ = SimpleConfigurable<int>(1,"XANADOO:CompressionLevel");        
00165    xanaEsdFile_->SetCompressionLevel(compressionLevel_);
00166 
00167    // create TTree
00168    autosave_ = SimpleConfigurable<int>(1000000000,"XANADOO:AutoSaveSize");        
00169    bufferSize_ = SimpleConfigurable<int>(256000,"XANADOO:BufferSize");        
00170    splitMode_ = SimpleConfigurable<int>(99,"XANADOO:SplitMode");        
00171    xana::out(1) << "creating xanaEsdTree  .."<<endl; 
00172    treeName_ = SimpleConfigurable<std::string>("XANAESD","XANADOO:TreeName");
00173    xanaEsdTree_ = new TTree(treeName_.value().c_str(),"test XANADOO ESD tree");  
00174    xanaEsdTree_->SetAutoSave(autosave_);
00175    int bufferSizeValue = bufferSize_.value();
00176    if (splitMode_.value()) bufferSizeValue /= 4; //from ROOT example
00177 
00178    // initialize event and create branches
00179    xanaEsdEvent_ = new XANAEsdEvent();
00180    xana::out(1) << "creating xanaEsdEventBranch  .."<<endl; 
00181    xanaEsdEventBranch_ = 
00182    xanaEsdTree_->Branch("Event.","XANAEsdEvent",&xanaEsdEvent_,bufferSizeValue,splitMode_.value());
00183    xana::out(1) << "xanaEsdEventBranch  branch created .."<<endl; 
00184 
00185    xanaGenEvent_ = new XANAGeneratorEvent();
00186    xana::out(1) << "creating xanaGeneratorEventBranch  .."<<endl; 
00187    xanaGenEventBranch_ = 
00188    xanaEsdTree_->Branch("GeneratorEvent.","XANAGeneratorEvent",&xanaGenEvent_,bufferSizeValue,splitMode_.value());
00189    xana::out(1) << "xanaGeneratorEventBranch  branch created .."<<endl; 
00190 
00191    xanaGeantEvent_ = new XANAGeantEvent();
00192    xana::out(1) << "creating xanaGeantEventBranch  .."<<endl; 
00193    xanaGeantEventBranch_ = 
00194    xanaEsdTree_->Branch("GeantEvent.","XANAGeantEvent",&xanaGeantEvent_,bufferSizeValue,splitMode_.value());
00195    xana::out(1) << "xanaGeantEventBranch  branch created .."<<endl; 
00196 
00197    // create links between generator-geant event and ESD-geant Event
00198    xanaGeantEvent_->addGeneratorEvent(xanaGenEvent_);
00199    xanaGeantEvent_->addEsdEvent(xanaEsdEvent_);
00200 
00201   //CC is this still needed?   
00202    ciso_ = Singleton<IsoCalculator>::instance();
00203 
00204 }
00205    
00206          
00207 XANAEsdBuilder::~XANAEsdBuilder() 
00208 {
00209    
00210   xana::out(1) << "Deleting XANAEsdBuilder" << std::endl;  
00211   xanaEsdFile_->Write();
00212   xanaEsdTree_->Print();
00213   delete xanaEsdFile_;
00214   if (memlog_) xanamemlog_->writeHistogram();
00215 
00216 }
00217 
00218 
00219 void XANAEsdBuilder::upDate( G3EventProxy * ev ) 
00220 { 
00221     
00222   xana::out(1) << "[update] updating XANAEsdBuilder for event " << std::endl;  
00223   if (!ev) return ;    
00224 
00225   //    UInt_t saveNumber=TProcessID::GetObjectCount();
00226      
00227   setHeader(ev);
00228   if (writeCaloRecHits_) setCaloRecHits(ev);
00229   setEmClusters(ev);
00230   setSuperClusters(ev);
00231   setEndcapClusters(ev);
00232   setHadClusters(ev);
00233   setTracks(ev);
00234   if (writeAllTrackHits_) setAllTrackHits(ev);
00235   setVertices(ev);
00236   setMuons(ev);
00237   if (!strcmp(electronReconstructorName_.value().c_str(),"GSF")) {
00238     setElectronsGSF(ev);
00239   }
00240   else {
00241     setElectronTracks(ev);
00242     setElectrons(ev);
00243   }
00244   setJets(ev);
00245   setMet(ev);
00246   setMcTruth(ev);
00247   setTriggerInfo(ev);
00248 
00249   if (writeTree_) writeTree();
00250   
00251   if (memlog_) xanamemlog_->fill();
00252 
00253   xanaEsdEvent_->clear();
00254   xanaGenEvent_->clear();
00255   xanaGeantEvent_->clear();
00256 
00257 //      TProcessID::SetObjectCount(saveNumber);
00258 
00259 }      
00260 
00261 
00262 void XANAEsdBuilder::upDate (G3SetUp *) 
00263 {
00264 
00265     xana::out(5) << "Starting initSetUp==================================" << std::endl;
00266     if (strcmp(electronReconstructorName_.value().c_str(),"GSF"))     return;
00267     
00268     // new from Wolfgang
00269     
00270     // RecQuery for regional tracks
00271     RecQuery regionalQuery("RegionalCombinatorialTrackFinder");
00272 
00273     // Building:
00274     RecQuery builderQuery("CombinatorialTrajectoryBuilder");
00275     builderQuery.setParameter("chiSquarCut",100000.);
00276     builderQuery.setParameter("mass",0.000511);
00277     builderQuery.setParameter("maxCand",2);
00278     builderQuery.setParameter("maxLostHit",1);
00279     builderQuery.setParameter("maxConsecLostHit",1);
00280     builderQuery.setParameter("lostHitPenalty",50.);
00281     builderQuery.setParameter("intermediateCleaning",true);
00282     builderQuery.setParameter("minimumNumberOfHits",5);
00283     builderQuery.setParameter("ptCut", 3.);
00284     builderQuery.setParameter("alwaysUseInvalidHits",true);
00285 
00286     // Fit & smooth part
00287     RecQuery fitQuery("KFFittingSmoother");
00288     fitQuery.setParameter("mass",0.000511); 
00289     
00290     // to fit forward with the gsf:
00291     if (!strcmp(gsfForwardFit_.value().c_str(),"GSF") ) 
00292       {
00293         RecQuery fwdQuery("GsfElectronFitter");
00294         fwdQuery.setParameter("mass",0.000511);
00295         fitQuery.setComponent("Fitter",fwdQuery);
00296       }
00297 
00298     // to fit forward with the kf:
00299     if (!strcmp(gsfForwardFit_.value().c_str(),"KF") ) 
00300       {
00301         RecQuery fwdQuery("KFTrajectoryFitter");   
00302         fwdQuery.setParameter("mass",0.000511); 
00303         fitQuery.setComponent("Fitter",fwdQuery);
00304       }
00305 
00306     // Smoother
00307     RecQuery bwdQuery("GsfElectronSmoother");
00308     fitQuery.setComponent("Smoother",bwdQuery);
00309 
00310     // Complete query for regional tracks 
00311     regionalQuery.setParameter("mass",0.000511);
00312     regionalQuery.setComponent("TrajectoryBuilder",builderQuery); 
00313     regionalQuery.setComponent("Smoother",fitQuery);
00314 
00315     // Create ConfigAlgorithm
00316     theRegionalTrackFinder = createAlgo<RegionalTrackFinder>(regionalQuery);
00317     // end new from Wolfgang
00318 
00319  }
00320 
00321 
00322 void XANAEsdBuilder::setEleTrackHits(std::vector<TrajectoryMeasurement> MyEleTrMeas, XANAElectronTrack *xTrack) 
00323 {
00324   // add electron track hits for this track, if asked for
00325 
00326   //     std::vector<TrajectoryMeasurement> MyEleTrMeas = eltrack->measurements();
00327      std::vector<TrajectoryMeasurement>::const_iterator IterEleTrMeas;
00328 
00329      xana::out(2) << "[update] adding electronic track hits " << std::endl; 
00330      xana::out(2) << "[update] number of measurements " <<MyEleTrMeas.size() << std::endl;
00331      //     xana::out(2) << "[update] number of valid hits "<<eltrack->foundHits()<< std::endl;
00332 
00333      XANATrackHit *elethit=0;
00334      int nbhits=0;
00335      for (IterEleTrMeas=MyEleTrMeas.begin(); IterEleTrMeas!=MyEleTrMeas.end(); IterEleTrMeas++)
00336        {
00337         // taking the rec hit 
00338         RecHit IterEleRecHit = IterEleTrMeas->recHit();                
00339         // checking if it is valid or not  
00340         if( !(IterEleRecHit.isValid()) ){xana::out(2) << "[update] non valid hit" << std::endl;}
00341 
00342         // taking the updated state 
00343         TrajectoryStateOnSurface IterEleTSN = IterEleTrMeas->updatedState();                
00344         // checking if it is valid or not  
00345         if( !(IterEleTSN.isValid()) ){xana::out(2) << "[update] non valid updated state" << std::endl;}
00346 
00347 
00348         // analysis only if both are valid 
00349         if (!(IterEleRecHit.isValid()) && (IterEleTSN.isValid())){xana::out(2)<<"[update] non valid hit but valid updated state"<< std::endl;}
00350         if ((IterEleRecHit.isValid()) && !(IterEleTSN.isValid())){xana::out(2)<<"[update] non valid updated state but valid hit"<< std::endl;}
00351         if(  IterEleRecHit.isValid() && IterEleTSN.isValid() )
00352         {
00353           HepPoint3D elehitposition( (IterEleRecHit).globalPosition().x(), 
00354                                      (IterEleRecHit).globalPosition().y(), 
00355                                      (IterEleRecHit).globalPosition().z());
00356 
00357           HepPoint3D eleupdposition( (IterEleTSN).globalPosition().x(),
00358                                      (IterEleTSN).globalPosition().y(),
00359                                      (IterEleTSN).globalPosition().z());
00360 
00361           HepVector3D eleupdmomentum( (IterEleTSN).globalMomentum().x(),
00362                                       (IterEleTSN).globalMomentum().y(),
00363                                       (IterEleTSN).globalMomentum().z());
00364 
00365           elethit = new XANATrackHit( elehitposition, eleupdposition, eleupdmomentum ); 
00366           elethit->setIsStereo((IterEleRecHit).det().detUnits().front()->type().isStereo());
00367           Bool_t isBarrel = (((IterEleRecHit).det().detUnits().front()->type().part() == barrel) ? true : false);
00368           elethit->setIsBarrel(isBarrel);
00369 
00370             xana::out(1) << "[update] ELE - NEW STEP" << std::endl; 
00371             xana::out(1) << "[update] new hit" << std::endl;    
00372             xana::out(1) << (IterEleRecHit).globalPosition().x() << std::endl;     
00373             xana::out(1) << (IterEleRecHit).globalPosition().y() << std::endl;          
00374             xana::out(1) << (IterEleRecHit).globalPosition().z() << std::endl;          
00375             xana::out(1) << "[update] new updated state: position" << std::endl;        
00376             xana::out(1) << (IterEleTSN).globalPosition().x() << std::endl;     
00377             xana::out(1) << (IterEleTSN).globalPosition().y() << std::endl;     
00378             xana::out(1) << (IterEleTSN).globalPosition().z() << std::endl;     
00379             xana::out(1) << "[update] new updated state: momentum" << std::endl;        
00380             xana::out(1) << (IterEleTSN).globalMomentum().x() << std::endl;     
00381             xana::out(1) << (IterEleTSN).globalMomentum().y() << std::endl;     
00382             xana::out(1) << (IterEleTSN).globalMomentum().z() << std::endl;     
00383 
00384           xanaEsdEvent_->addEleTrackHit(elethit, xTrack);           
00385           ++nbhits;
00386           delete elethit;
00387 
00388        } // valid hit and upd state
00389        } // loop on measurements
00390      xana::out(2) << "[update] number of hits that were added "<<nbhits<< std::endl;
00391 }
00392 
00393 
00394 void XANAEsdBuilder::writeTree() 
00395 {
00396 
00397   xana::out(2) << "[update] writing event " << std::endl;     
00398   
00399   xanaEsdTree_->Fill();  //fill the tree
00400 
00401 }
00402 
00403 
00404 std::string XANAEsdBuilder::getAlgoName(int algo) {
00405 
00406   if (algo==1) return std::string("EgammaBasicIsland");
00407   else if (algo>=100) return std::string("EgammaBasicHybrid");
00408   else return std::string("");
00409 
00410 }      
00411  
00412  
00413 static PKBuilder<XANAEsdBuilder> eventObserver("XANAEsdBuilder") ;