PhotonConversion/truthconversion/Conversion.h

0001 /**
0002  * @Conversion.h
0003  * @author Francesco Vassalli <Francesco.Vassalli@colorado.edu>
0004  * @version 1.0
0005  *
0006  * @brief Simple class to hold conversion information
0007  * This could be optimized by initializing the object with the necessary node pointers
0008  **/
0009 #ifndef CONVERSION_H__
0010 #define CONVERSION_H__
0011
0012 #include <g4main/PHG4Particle.h>
0013 #include <g4main/PHG4VtxPoint.h>
0014 #include <trackbase_historic/SvtxTrack.h>
0015 #include <trackbase_historic/SvtxClusterMap.h>
0016 #include <trackbase_historic/SvtxVertex.h>
0017 #include <g4eval/SvtxTrackEval.h>
0018 #include <GenFit/GFRaveVertex.h>
0019 #include <TLorentzVector.h>
0020 #include <TVector3.h>
0021 #include <utility>
0022
0023 namespace PHGenFit {
0024   class Track;
0025 } /* namespace PHGenFit */
0026
0027
0028 class SvtxTrackEval;
0029 class PHCompositeNode;
0030 class PHG4TruthInfoContainer;
0031 class TLorentzVector;
0032 class SvtxHitMap;
0033 class SVReco;
0034 class VtxRegressor;
0035
0036 //TODO check function to see ifthey can be made const
0037 class Conversion
0038 {
0039   public:
0040     Conversion(){verbosity=0;}
0041     Conversion(SvtxTrackEval* trackeval,int verbosity=0);
0042     Conversion(PHG4VtxPoint* vtx,int verbosity=0);
0043     Conversion(PHG4VtxPoint* vtx,SvtxTrackEval *trackeval,int verbosity=0);
0044     ~Conversion();
0045     /** sets the daughters of the conversion
0046      * use this to set the electron and positron
0047      * initializes both points but does not determine charge**/
0048     void setElectron(PHG4Particle* e);
0049     inline void setTrackEval(SvtxTrackEval *trackeval){this->trackeval=trackeval;}
0050     inline void setVtx(PHG4VtxPoint* vtx){this->vtx=vtx;}
0051     inline PHG4VtxPoint* getVtx()const{return vtx;}
0052     ///true if both electrons and photon are set
0053     inline bool isComplete()const{
0054       if (e1&&e2&&photon)
0055       {
0056         return true;
0057       }
0058       else{
0059         return false;
0060       }
0061     }
0062     ///check if both electrons are initialized
0063     inline bool hasPair()const{
0064       if (e1&&e2)
0065       {
0066         return true;
0067       }
0068       else{
0069         return false;
0070       }
0071     }
0072     ///sets the photon @return true if there are no errors
0073     bool setParent(PHG4Particle* parent);
0074
0075     ///calls @{link setParent} {@link photon} is set with @param parent and the {@link primaryPhoton} is set with the primary particle associated with @param parent
0076     void setPrimaryPhoton(PHG4Particle* parent,PHG4TruthInfoContainer* truthinfo);
0077
0078     /// @return the number of matched reco tracks
0079     inline int recoCount()const{
0080       int r=0;
0081       if (reco1)
0082       {
0083         r++;
0084       }
0085       if (reco2)
0086       {
0087         r++;
0088       }
0089       return r;
0090     }
0091     ///calls {@link setElectron()} then returns {@link e1}
0092     PHG4Particle* getElectron();
0093     ///@return {@link e1}
0094     inline PHG4Particle* getElectron()const{return e1;}
0095     ///calls {@link setElectron()} if true returns {@link e2} else returns {@link e1}
0096     PHG4Particle* getPositron();
0097     ///@return {@link e2}
0098     inline PHG4Particle* getPositron()const{return e2;}
0099     inline PHG4Particle* getPhoton()const{return photon;}
0100     inline PHG4Particle* getPrimaryPhoton()const{return primaryPhoton;}
0101     inline int getEmbed() const {return embedID;}
0102
0103     inline void setEmbed(int embedID) {this->embedID=embedID;}
0104     inline void setSourceId(int source){sourceId=source;}
0105     inline int getSourceId()const{return sourceId;}
0106     inline genfit::GFRaveVertex* getRecoVtx()const{return recoVertex;}
0107     /** Finds the delta eta of the reco tracks.
0108      * @return -1 if reco tracks are not set **/
0109     float trackDEta()const;
0110     static float trackDEta(SvtxTrack* reco1, SvtxTrack* reco2);
0111     /** Finds the delta phi of the reco tracks.
0112      * @return -1 if reco tracks are not set **/
0113     inline float trackDPhi()const{
0114       if (recoCount()==2)
0115       {
0116         return DeltaPhi(reco1->get_phi(),reco2->get_phi());
0117       }
0118       else return -1.;
0119     }
0120     inline static float trackDPhi(SvtxTrack* reco1, SvtxTrack* reco2){
0121       return DeltaPhi(reco1->get_phi(),reco2->get_phi());
0122     }
0123     ///@return the minimun reco track pT
0124     float minTrackpT();
0125     std::pair<float,float> getTrackpTs();
0126     std::pair<float,float> getTrackEtas();
0127     std::pair<float,float> getTrackPhis();
0128     inline std::pair<SvtxTrack*,SvtxTrack*> getRecoTracks()const{
0129       return std::pair<SvtxTrack*,SvtxTrack*>(reco1,reco2);
0130     }
0131     ///@return the reco track associated with the @param truthID
0132     SvtxTrack* getRecoTrack(unsigned truthID) const;
0133     /** set the reco maps used for {@link trackDEta}, {@link trackDLayer},{@link hasSilicon}
0134      *Note that all use of SvtxClusterMap is now deprecated**/
0135     inline void setClusterMap(SvtxClusterMap* cmap){
0136       _svtxClusterMap=cmap;
0137     }
0138     inline int setVerbosity(int v){
0139       verbosity=v;
0140       return verbosity;
0141     }
0142     /**Finds the cluster associated with {@link reco1}
0143      * if the trackeval was not previously set it needs to be set now
0144      * if {@link reco1} not set will attempt to set
0145      * @return -1 if no cluster is found**/
0146     int get_cluster_id();
0147     /**Finds the cluster associated with {@link reco1}
0148      * @return -1 if no cluster is found**/
0149     int get_cluster_id()const;
0150     /**Finds the cluster associated with {@link reco1}
0151      * if {@link reco1} not set will attempt to set
0152      * @return -1 if no cluster is found**/
0153     int get_cluster_id(SvtxTrackEval *trackeval);
0154     /**@return the cluster ID for both tracks if both are set.
0155      *If one is set pair.first will be the id.
0156      *the unset tracks get -1 as the cluster id**/
0157     std::pair<int,int> get_cluster_ids();
0158     /** match the truth particles to reco tracks
0159      * if the trackeval was not previously set it needs to be set now
0160      * @return the number of matched tracks [0,2]**/
0161     int setRecoTracks(SvtxTrackEval* trackeval);
0162     int setRecoTracks();
0163     /** call {@link setElectron()} then {@link setRecoTrakcs} then set the
0164      * reco track for the truth particle that corresponds to @param truthID
0165      * if it does not already have a reco track**/
0166     void setRecoTrack(int truthID, SvtxTrack* recoTrack);
0167     /**Sets the negative charge to {@link e1}
0168      * @return false otherwise
0169      * will print warnings if something is wrong with the Conversion
0170      **/
0171     bool setElectron();
0172     /** Return the difference in layers of the first hits of the reco track
0173      * @return -1 if reco tracks are not set**/
0174     int trackDLayer(SvtxClusterMap* cmap,SvtxHitMap *hitmap);
0175     int trackDLayer(TrkrClusterContainer* );
0176     static int trackDLayer(TrkrClusterContainer*,SvtxTrack*,SvtxTrack* );
0177     ///@return the first layer the associated reco clusters hit
0178     int firstLayer(SvtxClusterMap* cmap,SvtxHitMap *hitmap);
0179     int firstLayer(TrkrClusterContainer* );
0180     ///@return true if there are any silicon hits for the conversion. Is deprecated.
0181     bool hasSilicon(SvtxClusterMap* );
0182     /** distance between two closest points on the reco tracks
0183      * @return -1 if tracks are not set**/
0184     double approachDistance()const;
0185     static double approachDistance(SvtxTrack*,SvtxTrack*);
0186     ///@return distance in xyz space between the vertex and the closest track hit
0187     double dist(PHG4VtxPoint* vtx, SvtxClusterMap* cmap);
0188     double dist(TVector3* vtx, SvtxClusterMap* cmap);
0189     double dist(TVector3* vtx, TrkrClusterContainer* cmap);
0190     ///@return dca of the reco track with the lower dca
0191     float minDca();
0192     ///@return the RPhi distance between the reco vtx and further track
0193     float vtxTrackRPhi(TVector3 vertpos);
0194     static float vtxTrackRPhi(TVector3 vertpos,SvtxTrack*,SvtxTrack*);
0195     ///@return the RZ distance between the reco vtx and further track
0196     float vtxTrackRZ(TVector3 vertpos);
0197     static float vtxTrackRZ(TVector3 vertpos,SvtxTrack*,SvtxTrack*);
0198
0199     //Deprecated
0200     //float setRecoVtx(SvtxVertex* recovtx,SvtxClusterMap* cmap);
0201     TLorentzVector* setRecoPhoton();///< constructs a tlv for a photon by adding the tlvs for the reco tracks
0202     TLorentzVector* getRecoPhoton();///<@return the constructed tlv does not return ownserhip
0203     static TLorentzVector* getRecoPhoton(SvtxTrack* reco1, SvtxTrack* reco2);///<@return the addition of respective tlvs, @return ownership
0204     ///Uses the {@link _refit_phgf_tracks} to set {@link recoPhoton}. Will delete any existing {@link recoPhoton}. @return NULL if either {@link _refit_phgf_tracks} are NULL.
0205     TLorentzVector* getRefitRecoPhoton();
0206     PHG4Particle* getTruthPhoton(PHG4TruthInfoContainer*);///<@return NULL if not valid conversion else return photon
0207     ///Uses {@link recoVertex} and \param vertexer to make {@link _refit_phgf_tracks}. Sets {@link recoVertex} to \param recoVtx
0208     std::pair<PHGenFit::Track*,PHGenFit::Track*> refitTracks(SVReco* vertexer, SvtxVertex* recoVtx);
0209     ///Uses {@link recoVertex} and \param vertexer to make {@link _refit_phgf_tracks}. If {@link recoVertex} is not set call {@link getSecondaryVertex(SVReco* vertexer)}.
0210     std::pair<PHGenFit::Track*,PHGenFit::Track*> refitTracks(SVReco* vertexer);
0211
0212     ///Uses the truth vertex and {@link SVReco} to improve the fit of {@link reco1} and {@link reco2}
0213     std::pair<PHGenFit::Track*,PHGenFit::Track*> refitTracksTruthVtx(SVReco* vertexer);
0214     ///Uses the truth vertex and {@link SVReco} to improve the fit of {@link reco1} and {@link reco2}. Uses @param seedVtx to set the uncertainies on the truth vertex.
0215     std::pair<PHGenFit::Track*,PHGenFit::Track*> refitTracksTruthVtx(SVReco* vertexer,SvtxVertex* seedVtx);
0216     //TODO Set {@link recoVertex} to the vertex reconstructed from {@link SVReco} using the reco tracks and @return {@link recoVertex}. @return NULL if {@link recoCount()}!=2. Delete any existing {@link recoVertex}
0217     ///@return the secondary vertex fit with {@link reco1} and {@link reco2}
0218     genfit::GFRaveVertex* getSecondaryVertex(SVReco* vertexer);
0219     ///get the PHGF version of {@link reco1} and {@link reco2}. If not possible @return NULL for that track.
0220     std::pair<PHGenFit::Track*,PHGenFit::Track*> getPHGFTracks(SVReco* vertexer);
0221     /*Convertes {@link reco1} and {@link reco2} into a pair of TLorentzVector* using the electron mass.
0222      * references {@link recoCount()} to determine which NULL values to @return
0223      * Ownership is returned.**/
0224     std::pair<TLorentzVector*,TLorentzVector*> getRecoTlvs();
0225     /*Convertes {@link _refit_phgf_tracks} into a pair of TLorentzVector* using the electron mass.
0226      * If either refit track is not defined @return will be a pair of NULL vectors
0227      * Ownership is returned.**/
0228     std::pair<TLorentzVector*,TLorentzVector*> getRefitRecoTlvs();
0229     /** @return the truth particle(s) that are missing reco tracks
0230      * fill NULL if the particles have a reco track
0231      * in the case that only one is missing a track it will be the first value of the pair**/
0232     std::pair<PHG4Particle*,PHG4Particle*> getParticlesMissingTrack();
0233     ///Uses the specified TMVA method to update and @return {@link recoVertex}
0234     genfit::GFRaveVertex* correctSecondaryVertex(VtxRegressor* regressor);
0235     static genfit::GFRaveVertex* correctSecondaryVertex(VtxRegressor* regressor,genfit::GFRaveVertex*, SvtxTrack* reco1, SvtxTrack* reco2);
0236     /**
0237      * Returns the equivalent angle in the range 0 to 2pi.
0238      **/
0239     inline static double InTwoPi (double phi) {
0240       while (phi < 0 || 2*TMath::Pi() <= phi) {
0241         if (phi < 0) phi += 2*TMath::Pi();
0242         else phi -= 2*TMath::Pi();
0243       }
0244       return phi;
0245     }
0246
0247     /**
0248      * Returns the difference between two angles in 0 to pi.
0249      **/
0250     inline static double DeltaPhi (double phi1, double phi2) {
0251       phi1 = InTwoPi(phi1);
0252       phi2 = InTwoPi(phi2);
0253       double dphi = abs(phi1 - phi2);
0254       while (dphi > TMath::Pi()) dphi = abs (dphi - 2*TMath::Pi());
0255       return dphi;
0256     }
0257
0258     ///print the truth info calls {@link e1->identify()}, {@link e2->identify()}, {@link vtx->identify()}, {@link photon->identify()}
0259     void printTruth();
0260     ///print the reco info calls {@link reco1->identify()}, {@link reco2->identify()}, {@link recoVertex->identify()}, {@link recoPhoton->print()}
0261     void printReco();
0262     ///print the info relating to photon reconstruction
0263     void PrintPhotonRecoInfo();
0264     void PrintPhotonRecoInfo(TLorentzVector *tlv_photon,TLorentzVector *tlv_electron, TLorentzVector *tlv_positron,float mass);
0265   private:
0266     PHG4Particle* e1=NULL;
0267     PHG4Particle* e2=NULL;
0268     PHG4Particle* photon=NULL;
0269     PHG4Particle* primaryPhoton=NULL;
0270     PHG4VtxPoint* vtx=NULL;
0271     SvtxVertex* truthSvtxVtx=NULL;
0272     SvtxTrack* reco1=NULL;
0273     SvtxTrack* reco2=NULL;
0274     SvtxTrackEval* trackeval=NULL;
0275     SvtxClusterMap* _svtxClusterMap=NULL;
0276     genfit::GFRaveVertex *recoVertex=NULL;
0277     TLorentzVector *recoPhoton=NULL;
0278     std::pair<PHGenFit::Track*,PHGenFit::Track*> _refit_phgf_tracks;
0279     std::pair<unsigned,unsigned> pairTruthReco1;
0280     std::pair<unsigned,unsigned> pairTruthReco2;
0281
0282     static const int _kNSiliconLayer =7; ///<hardcoded
0283     int embedID=0;
0284     int verbosity;
0285     int sourceId;
0286     static const float _kElectronRestM;
0287     ///helper function
0288     static float vtxTrackRZ(TVector3 recoVertPos,SvtxTrack *track);
0289     ///helper function
0290     static float vtxTrackRPhi(TVector3 recoVertPos,SvtxTrack *track);
0291     /**Initializes {@link truthSvtxVtx} with x,y,x,t from the PHG4VtxPoint.
0292      *Chisq and ndof are set to 1.
0293      *The tracks are set with {@link reco1} and {@link reco2}.
0294      *The errors are set to a randomish positive-definite matrix**/
0295     void PHG4VtxPointToSvtxVertex();
0296     /**Initializes {@link truthSvtxVtx} with x,y,x,t from the PHG4VtxPoint.
0297      *Chisq and ndof are set to 1.
0298      *The tracks are set with {@link reco1} and {@link reco2}.
0299      *The errors are set to the errors of the seedVtx. The primary vtx is suggested.**/
0300     void PHG4VtxPointToSvtxVertex(SvtxVertex* seedVtx);
0301
0302 };
0303 #endif //CONVERSION_H__