CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/Alignment/OfflineValidation/plugins/ValidationMisalignedTracker.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    ValidationMisalignedTracker
00004 // Class:      ValidationMisalignedTracker
00005 // 
00013 //
00014 // Original Author:  Nicola De Filippis
00015 //         Created:  Thu Dec 14 13:13:32 CET 2006
00016 // $Id: ValidationMisalignedTracker.cc,v 1.7 2011/12/20 15:11:41 mussgill Exp $
00017 //
00018 //
00019 
00020 
00021 #include "Alignment/OfflineValidation/plugins/ValidationMisalignedTracker.h"
00022 
00023 
00024 // user include files
00025 
00026 #include "DataFormats/TrackReco/interface/Track.h"
00027 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
00028 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
00029 
00030 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00031 #include "FWCore/Framework/interface/MakerMacros.h"
00032 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00033 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
00034 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00035 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
00036 
00037 //
00038 // constructors and destructor
00039 //
00040 ValidationMisalignedTracker::ValidationMisalignedTracker(const edm::ParameterSet& iConfig)
00041 {
00042 
00043   //now do what ever initialization is needed
00044   mzmu=0.,recmzmu=0.,ptzmu=0.,recptzmu=0.,etazmu=0.,recetazmu=0., thetazmu=0.,recthetazmu=0.,phizmu=0.,recphizmu=0.;
00045   recenezmu=0., enezmu=0., pLzmu=0., recpLzmu=0.,yzmu=0.,recyzmu=0.,mxptmu=0.,recmxptmu=0., minptmu=0.,recminptmu=0.;
00046   // mzele=0.,recmzele=0.
00047   
00048   flag=0,flagrec=0,count=0,countrec=0;
00049   nAssoc=0;
00050 
00051   for (int i=0;i<2;i++){
00052     countpart[i]=0;
00053     countpartrec[i]=0; 
00054     for (int j=0;j<2;j++){
00055       ene[i][j]=0.;
00056       p[i][j]=0.;
00057       px[i][j]=0.;
00058       py[i][j]=0.; 
00059       pz[i][j]=0.;
00060       ptmu[i][j]=0.; 
00061       recene[i][j]=0.;
00062       recp[i][j]=0.;
00063       recpx[i][j]=0.;
00064       recpy[i][j]=0.; 
00065       recpz[i][j]=0.;       
00066       recptmu[i][j]=0.;
00067     }
00068   }
00069 
00070 
00071   eventCount_ = 0; 
00072 
00073   selection_eff    = iConfig.getUntrackedParameter<bool>("selection_eff", false);
00074   selection_fake   = iConfig.getUntrackedParameter<bool>("selection_fake", true);
00075   ZmassSelection_  = iConfig.getUntrackedParameter<bool>("ZmassSelection", false);
00076   simobject        = iConfig.getUntrackedParameter<std::string>("simobject","g4SimHits");
00077   trackassociator  = iConfig.getUntrackedParameter<std::string>("TrackAssociator","ByHits");
00078   associators      = iConfig.getParameter< std::vector<std::string> >("associators");
00079   label            = iConfig.getParameter< std::vector<edm::InputTag> >("label");
00080   label_tp_effic   = iConfig.getParameter< edm::InputTag >("label_tp_effic");
00081   label_tp_fake    = iConfig.getParameter< edm::InputTag >("label_tp_fake");
00082 
00083   rootfile_   = iConfig.getUntrackedParameter<std::string>("rootfile","myroot.root");
00084   file_ = new TFile(rootfile_.c_str(),"RECREATE");
00085   
00086   // initialize the tree
00087   tree_eff = new TTree("EffTracks","Efficiency Tracks Tree");
00088 
00089   tree_eff->Branch("Run",&irun,"irun/i");
00090   tree_eff->Branch("Event",&ievt,"ievt/i");
00091 
00092   // SimTrack
00093   tree_eff->Branch("TrackID",&trackType,"trackType/i");
00094   tree_eff->Branch("pt",&pt,"pt/F");
00095   tree_eff->Branch("eta",&eta,"eta/F");
00096   tree_eff->Branch("CotTheta",&cottheta,"cottheta/F");
00097   tree_eff->Branch("phi",&phi,"phi/F");
00098   tree_eff->Branch("d0",&d0,"d0/F");
00099   tree_eff->Branch("z0",&z0,"z0/F"); 
00100   tree_eff->Branch("nhit",&nhit,"nhit/i");
00101     
00102   // RecTrack
00103   tree_eff->Branch("recpt",&recpt,"recpt/F");
00104   tree_eff->Branch("receta",&receta,"receta/F");
00105   tree_eff->Branch("CotRecTheta",&reccottheta,"reccottheta/F");
00106   tree_eff->Branch("recphi",&recphi,"recphi/F");
00107   tree_eff->Branch("recd0",& recd0,"recd0/F");
00108   tree_eff->Branch("recz0",& recz0,"recz0/F");
00109   tree_eff->Branch("nAssoc",&nAssoc,"nAssoc/i");
00110   tree_eff->Branch("recnhit",&recnhit,"recnhit/i");
00111   tree_eff->Branch("CHISQ",&recchiq,"recchiq/F");
00112 
00113   tree_eff->Branch("reseta",&reseta,"reseta/F");
00114   tree_eff->Branch("respt",&respt,"respt/F");
00115   tree_eff->Branch("resd0",&resd0,"resd0/F");
00116   tree_eff->Branch("resz0",&resz0,"resz0/F");
00117   tree_eff->Branch("resphi",&resphi,"resphi/F");
00118   tree_eff->Branch("rescottheta",&rescottheta,"rescottheta/F");
00119   tree_eff->Branch("eff",&eff,"eff/F");
00120 
00121   // Invariant masses, pt of Z
00122   tree_eff->Branch("mzmu",&mzmu,"mzmu/F");
00123   tree_eff->Branch("ptzmu",&ptzmu,"ptzmu/F");
00124   tree_eff->Branch("pLzmu",&pLzmu,"pLzmu/F");
00125   tree_eff->Branch("enezmu",&enezmu,"enezmu/F");
00126   tree_eff->Branch("etazmu",&etazmu,"etazmu/F");
00127   tree_eff->Branch("thetazmu",&thetazmu,"thetazmu/F");
00128   tree_eff->Branch("phizmu",&phizmu,"phizmu/F");
00129   tree_eff->Branch("yzmu",&yzmu,"yzmu/F");
00130   tree_eff->Branch("mxptmu",&mxptmu,"mxptmu/F");
00131   tree_eff->Branch("minptmu",&minptmu,"minptmu/F");
00132 
00133   tree_eff->Branch("recmzmu",&recmzmu,"recmzmu/F");
00134   tree_eff->Branch("recptzmu",&recptzmu,"recptzmu/F");
00135   tree_eff->Branch("recpLzmu",&recpLzmu,"recpLzmu/F");  
00136   tree_eff->Branch("recenezmu",&recenezmu,"recenezmu/F");
00137   tree_eff->Branch("recetazmu",&recetazmu,"recetazmu/F");
00138   tree_eff->Branch("recthetazmu",&recthetazmu,"recthetazmu/F");
00139   tree_eff->Branch("recphizmu",&recphizmu,"recphizmu/F");
00140   tree_eff->Branch("recyzmu",&recyzmu,"recyzmu/F");
00141   tree_eff->Branch("recmxptmu",&recmxptmu,"recmxptmu/F");   
00142   tree_eff->Branch("recminptmu",&recminptmu,"recminptmu/F");       
00143 
00144 
00145   //tree->Branch("mzele",&ntmzele,"ntmzele/F");
00146   //tree->Branch("recmzele",&ntmzeleRec,"ntmzeleRec/F");
00147   tree_eff->Branch("chi2Associator",&recchiq,"recchiq/F");
00148 
00149   // Fake
00150 
00151   tree_fake = new TTree("FakeTracks","Fake Rate Tracks Tree");
00152 
00153   tree_fake->Branch("Run",&irun,"irun/i");
00154   tree_fake->Branch("Event",&ievt,"ievt/i");
00155 
00156   // SimTrack
00157   tree_fake->Branch("fakeTrackID",&faketrackType,"faketrackType/i");
00158   tree_fake->Branch("fakept",&fakept,"fakept/F");
00159   tree_fake->Branch("fakeeta",&fakeeta,"fakeeta/F");
00160   tree_fake->Branch("fakeCotTheta",&fakecottheta,"fakecottheta/F");
00161   tree_fake->Branch("fakephi",&fakephi,"fakephi/F");
00162   tree_fake->Branch("faked0",&faked0,"faked0/F");
00163   tree_fake->Branch("fakez0",&fakez0,"fakez0/F"); 
00164   tree_fake->Branch("fakenhit",&fakenhit,"fakenhit/i");
00165     
00166   // RecTrack
00167   tree_fake->Branch("fakerecpt",&fakerecpt,"fakerecpt/F");
00168   tree_fake->Branch("fakereceta",&fakereceta,"fakereceta/F");
00169   tree_fake->Branch("fakeCotRecTheta",&fakereccottheta,"fakereccottheta/F");
00170   tree_fake->Branch("fakerecphi",&fakerecphi,"fakerecphi/F");
00171   tree_fake->Branch("fakerecd0",& fakerecd0,"fakerecd0/F");
00172   tree_fake->Branch("fakerecz0",& fakerecz0,"fakerecz0/F");
00173   tree_fake->Branch("fakenAssoc",&fakenAssoc,"fakenAssoc/i");
00174   tree_fake->Branch("fakerecnhit",&fakerecnhit,"fakerecnhit/i");
00175   tree_fake->Branch("fakeCHISQ",&fakerecchiq,"fakerecchiq/F");
00176 
00177   tree_fake->Branch("fakereseta",&fakereseta,"fakereseta/F");
00178   tree_fake->Branch("fakerespt",&fakerespt,"fakerespt/F");
00179   tree_fake->Branch("fakeresd0",&fakeresd0,"fakeresd0/F");
00180   tree_fake->Branch("fakeresz0",&fakeresz0,"fakeresz0/F");
00181   tree_fake->Branch("fakeresphi",&fakeresphi,"fakeresphi/F");
00182   tree_fake->Branch("fakerescottheta",&fakerescottheta,"fakerescottheta/F");
00183   tree_fake->Branch("fake",&fake,"fake/F");
00184 
00185   // Invariant masses, pt of Z
00186   tree_fake->Branch("fakemzmu",&fakemzmu,"fakemzmu/F");
00187   tree_fake->Branch("fakeptzmu",&fakeptzmu,"fakeptzmu/F");
00188   tree_fake->Branch("fakepLzmu",&fakepLzmu,"fakepLzmu/F");
00189   tree_fake->Branch("fakeenezmu",&fakeenezmu,"fakeenezmu/F");
00190   tree_fake->Branch("fakeetazmu",&fakeetazmu,"fakeetazmu/F");
00191   tree_fake->Branch("fakethetazmu",&fakethetazmu,"fakethetazmu/F");
00192   tree_fake->Branch("fakephizmu",&fakephizmu,"fakephizmu/F");
00193   tree_fake->Branch("fakeyzmu",&fakeyzmu,"fakeyzmu/F");
00194   tree_fake->Branch("fakemxptmu",&fakemxptmu,"fakemxptmu/F");
00195   tree_fake->Branch("fakeminptmu",&fakeminptmu,"fakeminptmu/F");
00196 
00197   tree_fake->Branch("fakerecmzmu",&fakerecmzmu,"fakerecmzmu/F");
00198   tree_fake->Branch("fakerecptzmu",&fakerecptzmu,"fakerecptzmu/F");
00199   tree_fake->Branch("fakerecpLzmu",&fakerecpLzmu,"fakerecpLzmu/F");  
00200   tree_fake->Branch("fakerecenezmu",&fakerecenezmu,"fakerecenezmu/F");
00201   tree_fake->Branch("fakerecetazmu",&fakerecetazmu,"fakerecetazmu/F");
00202   tree_fake->Branch("fakerecthetazmu",&fakerecthetazmu,"fakerecthetazmu/F");
00203   tree_fake->Branch("fakerecphizmu",&fakerecphizmu,"fakerecphizmu/F");
00204   tree_fake->Branch("fakerecyzmu",&fakerecyzmu,"fakerecyzmu/F");
00205   tree_fake->Branch("fakerecmxptmu",&fakerecmxptmu,"fakerecmxptmu/F");   
00206   tree_fake->Branch("fakerecminptmu",&fakerecminptmu,"fakerecminptmu/F");       
00207 
00208   tree_fake->Branch("fakechi2Associator",&fakerecchiq,"fakerecchiq/F");
00209 
00210 }
00211 
00212 
00213 ValidationMisalignedTracker::~ValidationMisalignedTracker()
00214 {
00215  
00216  
00217   std::cout << "ValidationMisalignedTracker::endJob Processed " << eventCount_
00218             << " events" << std::endl;
00219                                                                                                                    
00220   // store the tree in the output file
00221   file_->Write();
00222                                                                                                                    
00223   
00224   // Closing the file deletes the tree.
00225   file_->Close();
00226   tree_eff=0;
00227   tree_fake=0;
00228 }
00229 
00230 
00231 //
00232 // member functions
00233 //
00234 
00235 // ------------ method called to for each event  ------------
00236 void
00237 ValidationMisalignedTracker::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00238 {
00239    if (watchTrackAssociatorRecord_.check(iSetup)) {
00240      associatore.clear();
00241      edm::ESHandle<TrackAssociatorBase> theAssociator;
00242      for (unsigned int w=0;w<associators.size();w++) {
00243        iSetup.get<TrackAssociatorRecord>().get(associators[w], theAssociator);
00244        associatore.push_back( theAssociator.product() );
00245      }
00246    }
00247    
00248    edm::LogInfo("Tracker Misalignment Validation") << "\n Starting!";
00249 
00250    // Monte Carlo Z selection
00251    skip=false;
00252    std::vector<int> indmu;
00253    
00254    if ( selection_eff && ZmassSelection_){ 
00255      edm::Handle<edm::HepMCProduct> evt;
00256      iEvent.getByLabel("source", evt);
00257      bool accepted = false;
00258      bool foundmuons=false;
00259      HepMC::GenEvent * myGenEvent = new  HepMC::GenEvent(*(evt->GetEvent()));
00260      
00261      for ( HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p ) { 
00262        if ( !accepted && ( (*p)->pdg_id() == 23 ) && (*p)->status() == 3 ) { 
00263          accepted=true;
00264          for( HepMC::GenVertex::particle_iterator aDaughter=(*p)->end_vertex()->particles_begin(HepMC::descendants); aDaughter !=(*p)->end_vertex()->particles_end(HepMC::descendants);aDaughter++){
00265            if ( abs((*aDaughter)->pdg_id())==13) {
00266              foundmuons=true;
00267              if ((*aDaughter)->status()!=1 ) {
00268                for( HepMC::GenVertex::particle_iterator byaDaughter=(*aDaughter)->end_vertex()->particles_begin(HepMC::descendants); byaDaughter !=(*aDaughter)->end_vertex()->particles_end(HepMC::descendants);byaDaughter++){
00269                  if ((*byaDaughter)->status()==1 && abs((*byaDaughter)->pdg_id())==13) {                   
00270                    indmu.push_back((*byaDaughter)->barcode());
00271                    std::cout << "Stable muon from Z with charge " <<  (*byaDaughter)->pdg_id()  << " and index "<< (*byaDaughter)->barcode() << std::endl;
00272                  }
00273                }
00274              }
00275              else {            
00276                indmu.push_back((*aDaughter)->barcode());
00277                std::cout << "Stable muon from Z with charge " <<  (*aDaughter)->pdg_id() << " and index "<< (*aDaughter)->barcode() << std::endl; 
00278              }       
00279            }               
00280          }
00281          if (!foundmuons){
00282            std::cout << "No muons from Z ...skip event" << std::endl;
00283            skip=true;
00284          } 
00285        }
00286      }     
00287      if ( !accepted) {
00288        std::cout << "No Z particles in the event ...skip event" << std::endl;
00289        skip=true;
00290      }   
00291    }
00292    else {
00293      skip=false;
00294    }
00295    
00296    //
00297    // Retrieve tracker geometry from event setup
00298    //
00299    edm::ESHandle<TrackerGeometry> trackerGeometry;
00300    iSetup.get<TrackerDigiGeometryRecord>().get( trackerGeometry );
00301    GeomDet* testGeomDet = trackerGeometry->detsTOB().front();
00302    std::cout << testGeomDet->position() << std::endl;
00303    
00304    
00305    //Dump Run and Event
00306    irun=iEvent.id().run();
00307    ievt=iEvent.id().event();
00308 
00309    // Reset tree variables
00310    int countpart[2]={0,0},countpartrec[2]={0,0},flag=0,flagrec=0,count=0,countrec=0;
00311    //int countsim=0;
00312    float ene[2][2],px[2][2],py[2][2],pz[2][2],ptmu[2][2];
00313    float recene[2][2],recp[2][2],recpx[2][2],recpy[2][2],recpz[2][2],recptmu[2][2];
00314    
00315    for (int i=0;i<2;i++){
00316      for (int j=0;j<2;j++){
00317        ene[i][j]=0.;
00318        px[i][j]=0.;
00319        py[i][j]=0.; 
00320        pz[i][j]=0.; 
00321        ptmu[i][j]=0.;
00322        recene[i][j]=0.;
00323        recp[i][j]=0.;
00324        recpx[i][j]=0.;
00325        recpy[i][j]=0.; 
00326        recpz[i][j]=0.;       
00327        recptmu[i][j]=0.;
00328      }
00329    }
00330    
00331    
00332    edm::Handle<TrackingParticleCollection>  TPCollectionHeff ;
00333    iEvent.getByLabel(label_tp_effic,TPCollectionHeff);
00334    const TrackingParticleCollection tPCeff = *(TPCollectionHeff.product());
00335    
00336    edm::Handle<TrackingParticleCollection>  TPCollectionHfake ;
00337    iEvent.getByLabel(label_tp_fake,TPCollectionHfake);
00338    const TrackingParticleCollection tPCfake = *(TPCollectionHfake.product());
00339    
00340                
00341    int w=0;
00342    for (unsigned int ww=0;ww<associators.size();ww++){
00343       //
00344       //get collections from the event
00345       //
00346 
00347       edm::InputTag algo = label[0];
00348      
00349       edm::Handle<edm::View<reco::Track> > trackCollection;
00350       iEvent.getByLabel(algo, trackCollection);
00351       const edm::View<reco::Track> tC = *(trackCollection.product());
00352             
00353            
00354       //associate tracks
00355       LogTrace("TrackValidator") << "Calling associateRecoToSim method" << "\n";
00356            reco::RecoToSimCollection recSimColl=associatore[ww]->associateRecoToSim(trackCollection,
00357                                                               TPCollectionHfake,
00358                                                                       &iEvent);
00359       
00360       LogTrace("TrackValidator") << "Calling associateSimToReco method" << "\n";
00361       reco::SimToRecoCollection simRecColl=associatore[ww]->associateSimToReco(trackCollection,
00362                                                                       TPCollectionHeff, 
00363                                                                       &iEvent);
00364 
00365    
00366 
00367       //
00368       //compute number of tracks per eta interval
00369       //
00370 
00371       if (selection_eff && !skip ) {
00372         std::cout << "Computing Efficiency" << std::endl;
00373 
00374         edm::LogVerbatim("TrackValidator") << "\n# of TrackingParticles (before cuts): " << tPCeff.size() << "\n";
00375         int ats = 0;
00376         int st=0;
00377         for (TrackingParticleCollection::size_type i=0; i<tPCeff.size(); i++){
00378           
00379           // Initialize variables
00380           eta = 0.,theta=0.,phi=0.,pt=0.,cottheta=0.,costheta=0.;
00381           d0=0.,z0=0.;
00382           nhit=0;
00383           receta = 0.,rectheta = 0.,recphi = 0.,recpt = 0.,reccottheta=0.,recd0=0.,recz0=0.;
00384           respt = 0.,resd0 = 0.,resz0 = 0.,reseta = 0.,resphi=0.,rescottheta=0.;
00385           recchiq = 0.;
00386           recnhit = 0;
00387           trackType = 0;
00388           eff=0;        
00389           
00390           // typedef edm::Ref<TrackingParticleCollection> TrackingParticleRef;
00391           TrackingParticleRef tp(TPCollectionHeff, i);
00392           if (tp->charge()==0) continue;
00393           st++;
00394           //pt=sqrt(tp->momentum().perp2());
00395           //eta=tp->momentum().eta();
00396           //vpos=tp->vertex().perp2()));
00397           
00398           const SimTrack * simulatedTrack = &(*tp->g4Track_begin());
00399         
00400           edm::ESHandle<MagneticField> theMF;
00401           iSetup.get<IdealMagneticFieldRecord>().get(theMF);
00402           FreeTrajectoryState 
00403             ftsAtProduction(GlobalPoint(tp->vertex().x(),tp->vertex().y(),tp->vertex().z()),
00404                             GlobalVector(simulatedTrack->momentum().x(),simulatedTrack->momentum().y(),simulatedTrack->momentum().z()),
00405                             TrackCharge(tp->charge()),
00406                             theMF.product());
00407           TSCPBuilderNoMaterial tscpBuilder;
00408           TrajectoryStateClosestToPoint tsAtClosestApproach 
00409             = tscpBuilder(ftsAtProduction,GlobalPoint(0,0,0));//as in TrackProducerAlgorithm
00410           GlobalPoint v = tsAtClosestApproach.theState().position();
00411           GlobalVector p = tsAtClosestApproach.theState().momentum();
00412 
00413           //  double qoverpSim = tsAtClosestApproach.charge()/p.mag();
00414           //  double lambdaSim = M_PI/2-p.theta();
00415           //  double phiSim    = p.phi();
00416           double dxySim    = (-v.x()*sin(p.phi())+v.y()*cos(p.phi()));
00417           double dszSim    = v.z()*p.perp()/p.mag() - (v.x()*p.x()+v.y()*p.y())/p.perp() * p.z()/p.mag();
00418           d0     = float(-dxySim);
00419           z0     = float(dszSim*p.mag()/p.perp());
00420 
00421   
00422           if (abs(simulatedTrack->type())==13 && simulatedTrack->genpartIndex() != -1 ) {
00423             std::cout << " TRACCIA SIM DI MUONI " << std::endl;  
00424             std::cout << "Gen part " << simulatedTrack->genpartIndex()<< std::endl;
00425             trackType=simulatedTrack->type();
00426             theta=simulatedTrack->momentum().theta();
00427             costheta=cos(theta);
00428             cottheta=1./tan(theta);
00429             
00430             eta=simulatedTrack->momentum().eta();
00431             phi=simulatedTrack->momentum().phi();
00432             pt=simulatedTrack->momentum().pt();
00433             nhit=tp->matchedHit();
00434             
00435 
00436             std::cout << "3) Before assoc: SimTrack of type = " << simulatedTrack->type() 
00437                       << " ,at eta = " << eta 
00438                       << " ,with pt at vertex = " << simulatedTrack->momentum().pt() << " GeV/c"
00439                       << " ,d0 =" << d0 
00440                       << " ,z0 =" << z0 
00441                       << " ,nhit=" << nhit
00442                       << std::endl;
00443 
00444             if ( ZmassSelection_ ){
00445               if (abs(trackType)==13 && (simulatedTrack->genpartIndex()==indmu[0] || simulatedTrack->genpartIndex()==indmu[1] )) { 
00446                 std::cout << " TRACK sim of muons from Z " << std::endl;  
00447                 flag=0;
00448                 count=countpart[0];
00449                 countpart[0]++;
00450               }
00451               else if (abs(trackType)==11) {
00452                 //std::cout << " TRACCIA SIM DI ELETTRONI " << std::endl;  
00453                 flag=1;
00454                 count=countpart[1];
00455                 countpart[1]++;
00456               }
00457               
00458               
00459               px[flag][count]=simulatedTrack->momentum().x();   
00460               py[flag][count]=simulatedTrack->momentum().y();
00461               pz[flag][count]=simulatedTrack->momentum().z();
00462               ptmu[flag][count]=simulatedTrack->momentum().pt(); 
00463               ene[flag][count]=simulatedTrack->momentum().e();
00464             }
00465             
00466             
00467             std::vector<std::pair<edm::RefToBase<reco::Track>, double> > rt;
00468             if(simRecColl.find(tp) != simRecColl.end()){
00469               
00470               rt = simRecColl[tp];
00471               if (rt.size()!=0) {
00472                 
00473                 edm::RefToBase<reco::Track> t = rt.begin()->first;
00474                 ats++;
00475 
00476                 // bool flagptused=false;
00477                 // for (unsigned int j=0;j<ptused.size();j++){
00478                 //   if (fabs(t->pt()-ptused[j])<0.001) {
00479                 //     flagptused=true;
00480                 //   }
00481                 // }
00482 
00483                 edm::LogVerbatim("TrackValidator") << "TrackingParticle #" << st << " with pt=" << t->pt() 
00484                                                    << " associated with quality:" << rt.begin()->second <<"\n";
00485                 std::cout << "Reconstructed Track:" << t->pt()<< std::endl;
00486                 std::cout << "\tpT: " << t->pt()<< std::endl;
00487                 std::cout << "\timpact parameter:d0: " << t->d0()<< std::endl;
00488                 std::cout << "\timpact parameter:z0: " << t->dz()<< std::endl;
00489                 std::cout << "\tAzimuthal angle of point of closest approach:" << t->phi()<< std::endl;
00490                 std::cout << "\tcharge: " << t->charge()<< std::endl;
00491                 std::cout << "\teta: " << t->eta()<< std::endl;
00492                 std::cout << "\tnormalizedChi2: " << t->normalizedChi2()<< std::endl;
00493 
00494                 recnhit=t->numberOfValidHits();
00495                 recchiq=t->normalizedChi2();
00496                 rectheta=t->theta();
00497                 reccottheta=1./tan(rectheta);
00498                 //receta=-log(tan(rectheta/2.));
00499                 receta=t->momentum().eta();
00500                 //         reccostheta=cos(matchedrectrack->momentum().theta());
00501                 recphi=t->phi(); 
00502                 recpt=t->pt();
00503                 ptused.push_back(recpt);
00504                 recd0=t->d0();
00505                 recz0=t->dz();
00506 
00507                 std::cout << "5) After call to associator: the best match has " 
00508                           << recnhit << " hits, Chi2 = " 
00509                           << recchiq << ", pt at vertex = " 
00510                           << recpt << " GeV/c, " 
00511                           << ", recd0 = " << recd0 
00512                           << ", recz0= " << recz0
00513                           << std::endl;
00514 
00515 
00516                 respt=recpt - pt;
00517                 resd0=recd0-d0;
00518                 resz0=recz0-z0;
00519                 reseta=receta-eta;
00520                 resphi=recphi-phi;
00521                 rescottheta=reccottheta-cottheta;
00522                 eff=1;
00523 
00524                 std::cout << "6) Transverse momentum residual=" << respt 
00525                           << " ,d0 residual=" << resd0 
00526                           << " ,z0 residual=" << resz0 
00527                           << " with eff=" << eff << std::endl;
00528                 
00529                 if ( ZmassSelection_ ){
00530                   
00531                   if (abs(trackType)==13) { 
00532                     std::cout << " TRACCIA RECO DI MUONI " << std::endl;  
00533                     flagrec=0;
00534                     countrec=countpartrec[0];
00535                     countpartrec[0]++;
00536                   }
00537                   else if (abs(trackType)==11) {
00538                     std::cout << " TRACCIA RECO DI ELETTRONI " << std::endl;  
00539                     flagrec=1;
00540                     countrec=countpartrec[1];
00541                     countpartrec[1]++;
00542                   }
00543                   
00544                   recp[flagrec][countrec]=sqrt(t->momentum().mag2());
00545                   recpx[flagrec][countrec]=t->momentum().x();   
00546                   recpy[flagrec][countrec]=t->momentum().y();
00547                   recpz[flagrec][countrec]=t->momentum().z();
00548                   recptmu[flagrec][countrec]=sqrt( (t->momentum().x()*t->momentum().x()) + (t->momentum().y()*t->momentum().y()) );
00549                   if (abs(trackType)==13) recene[flagrec][countrec]=sqrt(recp[flagrec][countrec]*recp[flagrec][countrec]+0.105*0.105);
00550                   if (abs(trackType)==11) recene[flagrec][countrec]=sqrt(recp[flagrec][countrec]*recp[flagrec][countrec]+0.0005*0.0005);
00551                 }
00552                 
00553                 std::cout << "7) Transverse momentum reconstructed =" << recpt 
00554                           << " at  eta= " << receta 
00555                           << " and phi= " << recphi 
00556                           << std::endl;
00557                 
00558               }
00559             }
00560             else{
00561               edm::LogVerbatim("TrackValidator") << "TrackingParticle #" << st
00562                                                  << " with pt=" << sqrt(tp->momentum().perp2())
00563                                                  << " NOT associated to any reco::Track" << "\n";
00564               receta =-100.;
00565               recphi =-100.;
00566               recpt  =-100.;
00567               recd0  =-100.;
00568               recz0  =-100;
00569               respt  =-100.;
00570               resd0  =-100.;
00571               resz0  =-100.;
00572               resphi =-100.;
00573               reseta =-100.;
00574               rescottheta=-100.;
00575               recnhit=100;
00576               recchiq=-100;
00577               eff=0;
00578               flagrec=100;            
00579             }
00580             
00581             std::cout << "Eff=" << eff << std::endl;
00582                
00583                // simulated muons
00584                
00585                std::cout <<"Flag is" << flag << std::endl;
00586                std::cout <<"RecFlag is" << flagrec << std::endl;
00587                
00588                if (countpart[0]==2 && flag==0) {
00589                  mzmu=sqrt(
00590                            (ene[0][0]+ene[0][1])*(ene[0][0]+ene[0][1])-
00591                            (px[0][0]+px[0][1])*(px[0][0]+px[0][1])-
00592                            (py[0][0]+py[0][1])*(py[0][0]+py[0][1])-
00593                            (pz[0][0]+pz[0][1])*(pz[0][0]+pz[0][1])
00594                            );
00595                  std::cout << "Mzmu " << mzmu << std::endl;
00596                  ptzmu=sqrt(
00597                            (px[0][0]+px[0][1])*(px[0][0]+px[0][1])+
00598                            (py[0][0]+py[0][1])*(py[0][0]+py[0][1])
00599                            );
00600                  
00601                  pLzmu=pz[0][0]+pz[0][1];
00602                  enezmu=ene[0][0]+ene[0][1];
00603                  phizmu=atan2((py[0][0]+py[0][1]),(px[0][0]+px[0][1]));
00604                  thetazmu=atan2(ptzmu,(pz[0][0]+pz[0][1]));
00605                  etazmu=-log(tan(thetazmu*3.14/360.));
00606                  yzmu=0.5*log((enezmu+pLzmu)/(enezmu-pLzmu));
00607                  mxptmu=std::max( ptmu[0][0], ptmu[0][1]);
00608                  minptmu=std::min( ptmu[0][0], ptmu[0][1]);
00609                }
00610                else {
00611                  mzmu=-100.;
00612                  ptzmu=-100.;
00613                  pLzmu=-100.;
00614                  enezmu=-100.;
00615                  etazmu=-100.;
00616                  phizmu=-100.;
00617                  thetazmu=-100.;
00618                  yzmu=-100.;
00619                  mxptmu=-100.;
00620                  minptmu=-100.;
00621                }      
00622 
00623                // reconstructed muons
00624                if (countpartrec[0]==2 && flagrec==0 ){
00625                  recmzmu=sqrt(
00626                               (recene[0][0]+recene[0][1])*(recene[0][0]+recene[0][1])-
00627                               (recpx[0][0]+recpx[0][1])*(recpx[0][0]+recpx[0][1])-
00628                               (recpy[0][0]+recpy[0][1])*(recpy[0][0]+recpy[0][1])-
00629                               (recpz[0][0]+recpz[0][1])*(recpz[0][0]+recpz[0][1])
00630                               );
00631                  std::cout << "RecMzmu " << recmzmu << std::endl;
00632                  recptzmu=sqrt(
00633                                (recpx[0][0]+recpx[0][1])*(recpx[0][0]+recpx[0][1])+
00634                                (recpy[0][0]+recpy[0][1])*(recpy[0][0]+recpy[0][1])
00635                                );
00636                  
00637                  recpLzmu=recpz[0][0]+recpz[0][1];      
00638                  recenezmu=recene[0][0]+recene[0][1];
00639                  recphizmu=atan2((recpy[0][0]+recpy[0][1]),(recpx[0][0]+recpx[0][1]));
00640                  recthetazmu=atan2(recptzmu,(recpz[0][0]+recpz[0][1]));
00641                  recetazmu=-log(tan(recthetazmu*3.14/360.));
00642                  recyzmu=0.5*log((recenezmu+recpLzmu)/(recenezmu-recpLzmu));
00643                  recmxptmu=std::max(recptmu[0][0], recptmu[0][1]);
00644                  recminptmu=std::min( recptmu[0][0], recptmu[0][1]);
00645                }
00646                else {
00647                  recmzmu=-100.;
00648                  recptzmu=-100.; 
00649                  recpLzmu=-100.;                 
00650                  recenezmu=-100.;       
00651                  recetazmu=-100.;
00652                  recphizmu=-100.;
00653                  recthetazmu=-100.;
00654                  recyzmu=-100.;
00655                  recmxptmu=-100;
00656                  recminptmu=-100.;
00657                }               
00658                
00659                tree_eff->Fill();
00660 
00661           } // end of loop on muons
00662         } // end of loop for tracking particle
00663       } // end of loop for efficiency
00664 
00665       //
00666       // Fake Rate
00667       // 
00668       if (selection_fake ) {
00669         std::cout << "Computing Fake Rate" << std::endl;
00670 
00671         fakeeta = 0.,faketheta=0.,fakephi=0.,fakept=0.,fakecottheta=0.,fakecostheta=0.;
00672         faked0=0.,fakez0=0.;
00673         fakenhit=0;
00674         fakereceta = 0.,fakerectheta = 0.,fakerecphi = 0.,fakerecpt = 0.,fakereccottheta=0.,fakerecd0=0.,fakerecz0=0.;
00675         fakerespt = 0.,fakeresd0 = 0.,fakeresz0 = 0.,fakereseta = 0.,fakeresphi=0.,fakerescottheta=0.;
00676         fakerecchiq = 0.;
00677         fakerecnhit = 0;
00678         faketrackType = 0;
00679         fake=0;
00680 
00681         
00682         //      int at=0;
00683         int rT=0;
00684         for(reco::TrackCollection::size_type i=0; i<tC.size(); ++i){
00685           edm::RefToBase<reco::Track> track(trackCollection, i);
00686           rT++;
00687 
00688           fakeeta = 0.,faketheta=0.,fakephi=0.,fakept=0.,fakecottheta=0.,fakecostheta=0.;
00689           faked0=0.,fakez0=0.;
00690           fakenhit=0;
00691           fakereceta = 0.,fakerectheta = 0.,fakerecphi = 0.,fakerecpt = 0.,fakereccottheta=0.,fakerecd0=0.,fakerecz0=0.;
00692           fakerespt = 0.,fakeresd0 = 0.,fakeresz0 = 0.,fakereseta = 0.,fakeresphi=0.,fakerescottheta=0.;
00693           fakerecchiq = 0.;
00694           fakerecnhit = 0;
00695           faketrackType = 0;
00696           fake=0;
00697           
00698           fakerecnhit=track->numberOfValidHits();
00699           fakerecchiq=track->normalizedChi2();
00700           fakerectheta=track->theta();
00701           fakereccottheta=1./tan(rectheta);
00702           //fakereceta=-log(tan(rectheta/2.));
00703           fakereceta=track->momentum().eta();
00704           //       fakereccostheta=cos(track->momentum().theta());
00705           fakerecphi=track->phi(); 
00706           fakerecpt=track->pt();
00707           fakerecd0=track->d0();
00708           fakerecz0=track->dz();
00709                   
00710           std::cout << "1) Before assoc: TkRecTrack at eta = " << fakereceta << std::endl;
00711           std::cout << "Track number "<< i << std::endl ;
00712           std::cout << "\tPT: " << track->pt()<< std::endl;
00713           std::cout << "\timpact parameter:d0: " << track->d0()<< std::endl;
00714           std::cout << "\timpact parameter:z0: " << track->dz()<< std::endl;
00715           std::cout << "\tAzimuthal angle of point of closest approach:" << track->phi()<< std::endl;
00716           std::cout << "\tcharge: " << track->charge()<< std::endl;
00717           std::cout << "\teta: " << track->eta()<< std::endl;
00718           std::cout << "\tnormalizedChi2: " << track->normalizedChi2()<< std::endl;
00719            
00720            
00721           std::vector<std::pair<TrackingParticleRef, double> > tp;
00722           
00723           //Compute fake rate vs eta
00724           if(recSimColl.find(track) != recSimColl.end()){
00725             tp = recSimColl[track];
00726             if (tp.size()!=0) {
00727               edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt() 
00728                                                  << " associated with quality:" << tp.begin()->second <<"\n";
00729               
00730               
00731               TrackingParticleRef tpr = tp.begin()->first;
00732               const SimTrack * fakeassocTrack = &(*tpr->g4Track_begin());
00733 
00734               edm::ESHandle<MagneticField> theMF;
00735               iSetup.get<IdealMagneticFieldRecord>().get(theMF);
00736               FreeTrajectoryState 
00737                 ftsAtProduction(GlobalPoint(tpr->vertex().x(),tpr->vertex().y(),tpr->vertex().z()),
00738                                 GlobalVector(fakeassocTrack->momentum().x(),fakeassocTrack->momentum().y(),fakeassocTrack->momentum().z()),
00739                                 TrackCharge(tpr->charge()),
00740                                 theMF.product());
00741               TSCPBuilderNoMaterial tscpBuilder;
00742               TrajectoryStateClosestToPoint tsAtClosestApproach 
00743                 = tscpBuilder(ftsAtProduction,GlobalPoint(0,0,0));//as in TrackProducerAlgorithm
00744               GlobalPoint v = tsAtClosestApproach.theState().position();
00745               GlobalVector p = tsAtClosestApproach.theState().momentum();
00746               
00747               //  double qoverpSim = tsAtClosestApproach.charge()/p.mag();
00748               //  double lambdaSim = M_PI/2-p.theta();
00749               //  double phiSim    = p.phi();
00750               double dxySim    = (-v.x()*sin(p.phi())+v.y()*cos(p.phi()));
00751               double dszSim    = v.z()*p.perp()/p.mag() - (v.x()*p.x()+v.y()*p.y())/p.perp() * p.z()/p.mag();
00752               faked0     = float(-dxySim);
00753               fakez0     = float(dszSim*p.mag()/p.perp());
00754 
00755 
00756               faketrackType=fakeassocTrack->type();
00757               faketheta=fakeassocTrack->momentum().theta();
00758               fakecottheta=1./tan(faketheta);
00759               fakeeta=fakeassocTrack->momentum().eta();
00760               fakephi=fakeassocTrack->momentum().phi();
00761               fakept=fakeassocTrack->momentum().pt();
00762               fakenhit=tpr->matchedHit();
00763 
00764               std::cout << "4) After call to associator: the best SimTrack match is of type" << fakeassocTrack->type() 
00765                         << " ,at eta = " << fakeeta 
00766                         << " and phi = " << fakephi  
00767                         << " ,with pt at vertex = " << fakept << " GeV/c" 
00768                         << " ,d0 global = " << faked0 
00769                         << " ,z0 = " << fakez0
00770                         << std::endl;
00771               fake=1;
00772 
00773               fakerespt=fakerecpt-fakept;
00774               fakeresd0=fakerecd0-faked0;
00775               fakeresz0=fakerecz0-fakez0;
00776               fakereseta=-log(tan(fakerectheta/2.))-(-log(tan(faketheta/2.)));
00777               fakeresphi=fakerecphi-fakephi;
00778               fakerescottheta=fakereccottheta-fakecottheta;
00779               
00780             }
00781           }
00782           else{
00783             edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
00784                                                << " NOT associated to any TrackingParticle" << "\n";
00785             
00786             fakeeta =-100.;
00787             faketheta=-100;
00788             fakephi =-100.;
00789             fakept  =-100.;
00790             faked0  =-100.;
00791             fakez0  =-100;
00792             fakerespt  =-100.;
00793             fakeresd0  =-100.;
00794             fakeresz0  =-100.;
00795             fakeresphi =-100.;
00796             fakereseta =-100.;
00797             fakerescottheta=-100.;
00798             fakenhit=100;
00799             fake=0;
00800           }
00801           
00802           tree_fake->Fill();
00803         }       
00804         
00805       } // End of loop on fakerate       
00806       
00807       w++;
00808       
00809    } // End of loop on associators
00810 }
00811 
00812 // ------------ method called once each job just after ending the event loop  ------------
00813 void ValidationMisalignedTracker::endJob() {
00814 
00815   std::cout << "\t Misalignment analysis completed \n" << std::endl;  
00816 
00817 }
00818 
00819 DEFINE_FWK_MODULE(ValidationMisalignedTracker);