analysis_INTT/src/Tracklet.h

0001 #ifndef TRACKLET_H
0002 #define TRACKLET_H
0003
0004 #include <set>
0005
0006 #include <TTree.h>
0007
0008 #include "GenHadron.h"
0009 #include "Hit.h"
0010 #include "Utilities.h"
0011
0012 using namespace std;
0013
0014 class Tracklet : public TObject
0015 {
0016   public:
0017     Tracklet();
0018     Tracklet(Hit *hit1, Hit *hit2);
0019     ~Tracklet();
0020
0021     void SetHits(Hit *hit1, Hit *hit2);
0022     Hit *Hit1();
0023     Hit *Hit2();
0024     float dEta();
0025     float dPhi();
0026     float dR();
0027     float Eta();
0028     float Phi();
0029     float tklVtxZ();
0030     int LayerComb(); // Layer combination: can only be 12, 23, 13
0031     void SetMatchedGenHardon();
0032     bool IsMatchedGenHadron();
0033     void SetGenHadron(GenHadron *genhad);
0034     GenHadron *MatchedGenHadron();
0035
0036   private:
0037     Hit *_hit1;
0038     Hit *_hit2;
0039     float _deta;
0040     float _dphi;
0041     bool _matched_genhadron;
0042     GenHadron *_genhad;
0043 };
0044
0045 Tracklet::Tracklet()
0046 {
0047     _hit1 = nullptr;
0048     _hit2 = nullptr;
0049     _matched_genhadron = false;
0050     _genhad = nullptr;
0051 }
0052
0053 Tracklet::Tracklet(Hit *hit1, Hit *hit2)
0054 {
0055     _hit1 = hit1;
0056     _hit2 = hit2;
0057     _deta = _hit2->Eta() - _hit1->Eta();
0058     _dphi = deltaPhi(_hit2->Phi(), _hit1->Phi());
0059     _matched_genhadron = false;
0060     _genhad = nullptr;
0061 }
0062
0063 Tracklet::~Tracklet() {}
0064
0065 void Tracklet::SetHits(Hit *hit1, Hit *hit2)
0066 {
0067     _hit1 = hit1;
0068     _hit2 = hit2;
0069     _deta = _hit2->Eta() - _hit1->Eta();
0070     _dphi = deltaPhi(_hit2->Phi(), _hit1->Phi());
0071     _matched_genhadron = false;
0072     _genhad = nullptr;
0073 }
0074
0075 Hit *Tracklet::Hit1() { return (_hit1); }
0076
0077 Hit *Tracklet::Hit2() { return (_hit2); }
0078
0079 float Tracklet::dEta() { return (_deta); }
0080
0081 float Tracklet::dPhi() { return (_dphi); }
0082
0083 float Tracklet::dR() { return (sqrt(_deta * _deta + _dphi * _dphi)); }
0084
0085 float Tracklet::Eta() { return (0.5 * (_hit1->Eta() + _hit2->Eta())); }
0086
0087 float Tracklet::Phi() { return (0.5 * (_hit1->Phi() + _hit2->Phi())); }
0088
0089 float Tracklet::tklVtxZ() { return (_hit1->posZ() - _hit1->rho() * (_hit2->posZ() - _hit1->posZ()) / (_hit2->rho() - _hit1->rho())); }
0090
0091 int Tracklet::LayerComb() { return (_hit1->Layer()) * 10 + (_hit2->Layer()); }
0092
0093 void Tracklet::SetMatchedGenHardon() { _matched_genhadron = true; }
0094
0095 bool Tracklet::IsMatchedGenHadron() { return _matched_genhadron; }
0096
0097 void Tracklet::SetGenHadron(GenHadron *genhad)
0098 {
0099     _matched_genhadron = true;
0100     _genhad = genhad;
0101 }
0102
0103 GenHadron *Tracklet::MatchedGenHadron() { return _genhad; }
0104
0105 /*------------------------------------------------------------------------------------------*/
0106
0107 class TrackletData
0108 {
0109   public:
0110     bool isdata;
0111     vector<vector<Hit *>> layers = {{}, {}};
0112     vector<Tracklet *> ProtoTkls, RecoTkls, RecoTkls_GenMatched;
0113     vector<GenHadron *> GenHadrons;
0114
0115     float vtxzwei;
0116
0117     int event, NClusLayer1, NClusLayer2, NPrototkl, NRecotkl_Raw, NRecotkl_GenMatched, NGenHadron;
0118     uint64_t INTT_BCO;
0119     float PV_x, PV_y, PV_z, TruthPV_x, TruthPV_y, TruthPV_z;
0120     float centrality_mbd;
0121     float mbd_south_charge_sum, mbd_north_charge_sum, mbd_charge_sum, mbd_charge_asymm, mbd_z_vtx;
0122     bool pu0_sel, is_min_bias;
0123     int process; // single diffractive process
0124     bool firedTrig10_MBDSNgeq2, firedTrig12_vtxle10cm, firedTrig13_vtxle30cm, MbdNSge0, MbdZvtxle10cm, validMbdVtx, InttBco_IsToBeRemoved;
0125
0126     vector<int> cluslayer;
0127     vector<float> clusphi, cluseta, clusphisize, cluszsize;
0128     vector<unsigned int> clusadc;
0129
0130     vector<float> tklclus1x, tklclus1y, tklclus1z, tklclus2x, tklclus2y, tklclus2z;                     // detailed info
0131     vector<float> tklclus1phi, tklclus1eta, tklclus2phi, tklclus2eta, tklclus1phisize, tklclus2phisize; // 1=inner, 2=outer
0132     vector<unsigned int> tklclus1adc, tklclus2adc;
0133
0134     vector<float> prototkl_eta, prototkl_phi, prototkl_deta, prototkl_dphi, prototkl_dR;
0135     vector<float> recotklraw_eta, recotklraw_phi, recotklraw_deta, recotklraw_dphi, recotklraw_dR;
0136     vector<float> recotklraw_dca3dvtx;
0137
0138     vector<bool> recotkl_isMatchedGChadron;
0139     vector<int> recotkl_clus1_matchedtrackID, recotkl_clus2_matchedtrackID;
0140     vector<int> recotkl_clus1_matchedAncestorTrackID, recotkl_clus2_matchedAncestorTrackID;
0141     vector<float> recotkl_matchedPG4P_eta, recotkl_matchedPG4P_phi;
0142     vector<float> recotkl_matchedGChadron_eta, recotkl_matchedGChadron_phi;
0143     vector<int> matchedGChadron_trackID;
0144     vector<float> matchedGChadron_pt, matchedGChadron_eta, matchedGChadron_phi;
0145
0146     vector<int> G4P_trackID;
0147     vector<float> G4P_Pt, G4P_eta, G4P_phi;
0148
0149     vector<int> PrimaryG4P_trackID, PrimaryG4P_PID;
0150     vector<float> PrimaryG4P_Pt, PrimaryG4P_eta, PrimaryG4P_phi;
0151     vector<bool> PrimaryG4P_IsRecotkl;
0152
0153     vector<int> GenHadron_PID, GenHadron_trackID;
0154     vector<float> GenHadron_Pt, GenHadron_eta, GenHadron_phi;
0155     vector<bool> GenHadron_IsRecotkl;
0156 };
0157
0158 void SetMinitree(TTree *outTree, TrackletData &tkldata, bool detailed = false)
0159 {
0160     outTree->Branch("event", &tkldata.event);
0161     outTree->Branch("INTT_BCO", &tkldata.INTT_BCO);
0162     outTree->Branch("NClusLayer1", &tkldata.NClusLayer1);
0163     outTree->Branch("NClusLayer2", &tkldata.NClusLayer2);
0164     outTree->Branch("NPrototkl", &tkldata.NPrototkl);
0165     outTree->Branch("NRecotkl_Raw", &tkldata.NRecotkl_Raw);
0166     outTree->Branch("MBD_centrality", &tkldata.centrality_mbd);
0167     outTree->Branch("MBD_south_charge_sum", &tkldata.mbd_south_charge_sum);
0168     outTree->Branch("MBD_north_charge_sum", &tkldata.mbd_north_charge_sum);
0169     outTree->Branch("MBD_charge_sum", &tkldata.mbd_charge_sum);
0170     outTree->Branch("MBD_charge_asymm", &tkldata.mbd_charge_asymm);
0171     outTree->Branch("MBD_z_vtx", &tkldata.mbd_z_vtx);
0172     outTree->Branch("PV_x", &tkldata.PV_x);
0173     outTree->Branch("PV_y", &tkldata.PV_y);
0174     outTree->Branch("PV_z", &tkldata.PV_z);
0175     outTree->Branch("vtxzwei", &tkldata.vtxzwei);
0176     outTree->Branch("firedTrig10_MBDSNgeq2", &tkldata.firedTrig10_MBDSNgeq2);
0177     outTree->Branch("firedTrig12_vtxle10cm", &tkldata.firedTrig12_vtxle10cm);
0178     outTree->Branch("firedTrig13_vtxle30cm", &tkldata.firedTrig13_vtxle30cm);
0179     outTree->Branch("MbdNSge0", &tkldata.MbdNSge0);
0180     outTree->Branch("MbdZvtxle10cm", &tkldata.MbdZvtxle10cm);
0181     outTree->Branch("validMbdVtx", &tkldata.validMbdVtx);
0182     outTree->Branch("InttBco_IsToBeRemoved", &tkldata.InttBco_IsToBeRemoved);
0183
0184     outTree->Branch("is_min_bias", &tkldata.is_min_bias);
0185     outTree->Branch("clusLayer", &tkldata.cluslayer);
0186     outTree->Branch("clusPhi", &tkldata.clusphi);
0187     outTree->Branch("clusEta", &tkldata.cluseta);
0188     outTree->Branch("clusPhiSize", &tkldata.clusphisize);
0189     outTree->Branch("clusZSize", &tkldata.cluszsize);
0190     outTree->Branch("clusADC", &tkldata.clusadc);
0191     if (detailed)
0192     {
0193         outTree->Branch("tklclus1x", &tkldata.tklclus1x);
0194         outTree->Branch("tklclus1y", &tkldata.tklclus1y);
0195         outTree->Branch("tklclus1z", &tkldata.tklclus1z);
0196         outTree->Branch("tklclus2x", &tkldata.tklclus2x);
0197         outTree->Branch("tklclus2y", &tkldata.tklclus2y);
0198         outTree->Branch("tklclus2z", &tkldata.tklclus2z);
0199     }
0200     outTree->Branch("tklclus1Phi", &tkldata.tklclus1phi);
0201     outTree->Branch("tklclus1Eta", &tkldata.tklclus1eta);
0202     outTree->Branch("tklclus2Phi", &tkldata.tklclus2phi);
0203     outTree->Branch("tklclus2Eta", &tkldata.tklclus2eta);
0204     outTree->Branch("tklclus1PhiSize", &tkldata.tklclus1phisize);
0205     outTree->Branch("tklclus2PhiSize", &tkldata.tklclus2phisize);
0206     outTree->Branch("tklclus1ADC", &tkldata.tklclus1adc);
0207     outTree->Branch("tklclus2ADC", &tkldata.tklclus2adc);
0208     outTree->Branch("recotklraw_eta", &tkldata.recotklraw_eta);
0209     outTree->Branch("recotklraw_phi", &tkldata.recotklraw_phi);
0210     outTree->Branch("recotklraw_deta", &tkldata.recotklraw_deta);
0211     outTree->Branch("recotklraw_dphi", &tkldata.recotklraw_dphi);
0212     outTree->Branch("recotklraw_dR", &tkldata.recotklraw_dR);
0213     outTree->Branch("recotklraw_dca3dvtx", &tkldata.recotklraw_dca3dvtx);
0214
0215     if (!tkldata.isdata)
0216     {
0217         outTree->Branch("recotkl_isMatchedGChadron", &tkldata.recotkl_isMatchedGChadron);
0218         outTree->Branch("recotkl_clus1_matchedtrackID", &tkldata.recotkl_clus1_matchedtrackID);
0219         outTree->Branch("recotkl_clus2_matchedtrackID", &tkldata.recotkl_clus2_matchedtrackID);
0220         outTree->Branch("recotkl_clus1_matchedAncestorTrackID", &tkldata.recotkl_clus1_matchedAncestorTrackID);
0221         outTree->Branch("recotkl_clus2_matchedAncestorTrackID", &tkldata.recotkl_clus2_matchedAncestorTrackID);
0222         outTree->Branch("matchedGChadron_pt", &tkldata.matchedGChadron_pt);
0223         outTree->Branch("matchedGChadron_eta", &tkldata.matchedGChadron_eta);
0224         outTree->Branch("matchedGChadron_phi", &tkldata.matchedGChadron_phi);
0225         outTree->Branch("TruthPV_x", &tkldata.TruthPV_x);
0226         outTree->Branch("TruthPV_y", &tkldata.TruthPV_y);
0227         outTree->Branch("TruthPV_z", &tkldata.TruthPV_z);
0228         outTree->Branch("pu0_sel", &tkldata.pu0_sel);
0229         outTree->Branch("process", &tkldata.process);
0230         outTree->Branch("PrimaryG4P_PID", &tkldata.PrimaryG4P_PID);
0231         outTree->Branch("PrimaryG4P_Pt", &tkldata.PrimaryG4P_Pt);
0232         outTree->Branch("PrimaryG4P_eta", &tkldata.PrimaryG4P_eta);
0233         outTree->Branch("PrimaryG4P_phi", &tkldata.PrimaryG4P_phi);
0234         outTree->Branch("PrimaryG4P_IsRecotkl", &tkldata.PrimaryG4P_IsRecotkl);
0235         outTree->Branch("GenHadron_PID", &tkldata.GenHadron_PID);
0236         outTree->Branch("GenHadron_Pt", &tkldata.GenHadron_Pt);
0237         outTree->Branch("GenHadron_eta", &tkldata.GenHadron_eta);
0238         outTree->Branch("GenHadron_phi", &tkldata.GenHadron_phi);
0239         outTree->Branch("GenHadron_IsRecotkl", &tkldata.GenHadron_IsRecotkl);
0240         if (detailed)
0241         {
0242             outTree->Branch("recotkl_matchedPG4P_eta", &tkldata.recotkl_matchedPG4P_eta);
0243             outTree->Branch("recotkl_matchedPG4P_phi", &tkldata.recotkl_matchedPG4P_phi);
0244             outTree->Branch("recotkl_matchedGChadron_eta", &tkldata.recotkl_matchedGChadron_eta);
0245             outTree->Branch("recotkl_matchedGChadron_phi", &tkldata.recotkl_matchedGChadron_phi);
0246         }
0247     }
0248 }
0249
0250 void ResetVec(TrackletData &tkldata)
0251 {
0252     for (size_t i = 0; i < tkldata.layers.size(); i++)
0253     {
0254         for (auto &hit : tkldata.layers[i])
0255         {
0256             delete hit;
0257         }
0258         CleanVec(tkldata.layers[i]);
0259     }
0260     // delete the pointers and then clear the vector; otherwise, the memory will is not released and will cause memory leak
0261     for (auto &tkl : tkldata.ProtoTkls)
0262     {
0263         delete tkl;
0264     }
0265     CleanVec(tkldata.ProtoTkls);
0266     CleanVec(tkldata.RecoTkls); // Don't need to delete the reco tracklet pointers because they are just pointers to the proto tracklets
0267
0268     CleanVec(tkldata.prototkl_eta);
0269     CleanVec(tkldata.prototkl_phi);
0270     CleanVec(tkldata.prototkl_deta);
0271     CleanVec(tkldata.prototkl_dphi);
0272     CleanVec(tkldata.prototkl_dR);
0273
0274     CleanVec(tkldata.cluslayer);
0275     CleanVec(tkldata.clusphi);
0276     CleanVec(tkldata.cluseta);
0277     CleanVec(tkldata.clusphisize);
0278     CleanVec(tkldata.cluszsize);
0279     CleanVec(tkldata.clusadc);
0280
0281     CleanVec(tkldata.tklclus1x);
0282     CleanVec(tkldata.tklclus1y);
0283     CleanVec(tkldata.tklclus1z);
0284     CleanVec(tkldata.tklclus2x);
0285     CleanVec(tkldata.tklclus2y);
0286     CleanVec(tkldata.tklclus2z);
0287     CleanVec(tkldata.tklclus1phi);
0288     CleanVec(tkldata.tklclus1eta);
0289     CleanVec(tkldata.tklclus2phi);
0290     CleanVec(tkldata.tklclus2eta);
0291     CleanVec(tkldata.tklclus1phisize);
0292     CleanVec(tkldata.tklclus2phisize);
0293     CleanVec(tkldata.tklclus1adc);
0294     CleanVec(tkldata.tklclus2adc);
0295
0296     CleanVec(tkldata.recotklraw_eta);
0297     CleanVec(tkldata.recotklraw_phi);
0298     CleanVec(tkldata.recotklraw_deta);
0299     CleanVec(tkldata.recotklraw_dphi);
0300     CleanVec(tkldata.recotklraw_dR);
0301     CleanVec(tkldata.recotklraw_dca3dvtx);
0302
0303     CleanVec(tkldata.recotkl_isMatchedGChadron);
0304     CleanVec(tkldata.recotkl_clus1_matchedtrackID);
0305     CleanVec(tkldata.recotkl_clus2_matchedtrackID);
0306     CleanVec(tkldata.recotkl_clus1_matchedAncestorTrackID);
0307     CleanVec(tkldata.recotkl_clus2_matchedAncestorTrackID);
0308     CleanVec(tkldata.matchedGChadron_trackID);
0309     CleanVec(tkldata.matchedGChadron_pt);
0310     CleanVec(tkldata.matchedGChadron_eta);
0311
0312     CleanVec(tkldata.PrimaryG4P_PID);
0313     CleanVec(tkldata.PrimaryG4P_trackID);
0314     CleanVec(tkldata.PrimaryG4P_Pt);
0315     CleanVec(tkldata.PrimaryG4P_eta);
0316     CleanVec(tkldata.PrimaryG4P_phi);
0317     CleanVec(tkldata.PrimaryG4P_IsRecotkl);
0318
0319     CleanVec(tkldata.G4P_trackID);
0320     CleanVec(tkldata.G4P_Pt);
0321     CleanVec(tkldata.G4P_eta);
0322     CleanVec(tkldata.G4P_phi);
0323
0324     for (auto &genhad : tkldata.GenHadrons)
0325     {
0326         delete genhad;
0327     }
0328     CleanVec(tkldata.GenHadrons);
0329     CleanVec(tkldata.GenHadron_PID);
0330     CleanVec(tkldata.GenHadron_trackID);
0331     CleanVec(tkldata.GenHadron_Pt);
0332     CleanVec(tkldata.GenHadron_eta);
0333     CleanVec(tkldata.GenHadron_phi);
0334     CleanVec(tkldata.GenHadron_IsRecotkl);
0335
0336     CleanVec(tkldata.recotkl_matchedPG4P_eta);
0337     CleanVec(tkldata.recotkl_matchedPG4P_phi);
0338     CleanVec(tkldata.recotkl_matchedGChadron_eta);
0339     CleanVec(tkldata.recotkl_matchedGChadron_phi);
0340 }
0341
0342 bool compare_dR(Tracklet *a, Tracklet *b) { return a->dR() < b->dR(); }
0343
0344 void ProtoTracklets(TrackletData &tkldata, float cutdr)
0345 {
0346     float Cut_dR = cutdr;
0347
0348     for (auto ihitl1 : tkldata.layers[0])
0349     {
0350         for (auto ihitl2 : tkldata.layers[1])
0351         {
0352             // float dEta = ihitl2->Eta() - ihitl1->Eta();
0353             // float dPhi = deltaPhi(ihitl2->Phi(), ihitl1->Phi());
0354             // float dR = sqrt(dEta * dEta + dPhi * dPhi);
0355
0356             Tracklet *tmptkl = new Tracklet(ihitl1, ihitl2);
0357
0358             if (tmptkl->dR() < Cut_dR)
0359             {
0360                 tkldata.ProtoTkls.push_back(tmptkl);
0361
0362                 // tkldata.prototkl_eta.push_back(tmptkl->Hit1()->Eta());
0363                 // tkldata.prototkl_phi.push_back(tmptkl->Hit1()->Phi());
0364                 // tkldata.prototkl_deta.push_back(tmptkl->dEta());
0365                 // tkldata.prototkl_dphi.push_back(tmptkl->dPhi());
0366                 // tkldata.prototkl_dR.push_back(tmptkl->dR());
0367             }
0368             else
0369             {
0370                 delete tmptkl;
0371             }
0372         }
0373     }
0374
0375     tkldata.NPrototkl = tkldata.ProtoTkls.size();
0376 }
0377
0378 // 3D DCA calculation
0379 struct Point3D
0380 {
0381     double x, y, z;
0382 };
0383
0384 float computeDCA(const Point3D &p1, const Point3D &p2, const Point3D &p3) // p1 and p2 are the coordinates of clusters on tracklets, p3 is the event vertex
0385 {
0386     // Direction vector of the line
0387     float dx = p2.x - p1.x;
0388     float dy = p2.y - p1.y;
0389     float dz = p2.z - p1.z;
0390
0391     // Vector from p1 to p3
0392     float vx = p3.x - p1.x;
0393     float vy = p3.y - p1.y;
0394     float vz = p3.z - p1.z;
0395
0396     float denominator = dx * dx + dy * dy + dz * dz;
0397     float numerator = vx * dx + vy * dy + vz * dz;
0398
0399     float t = numerator / denominator; // Projection parameter
0400
0401     // Compute the closest point on the line
0402     float x_c = p1.x + t * dx;
0403     float y_c = p1.y + t * dy;
0404     float z_c = p1.z + t * dz;
0405
0406     // Compute the Euclidean distance between p3 and the closest point
0407     float dca = std::sqrt((p3.x - x_c) * (p3.x - x_c) + (p3.y - y_c) * (p3.y - y_c) + (p3.z - z_c) * (p3.z - z_c));
0408
0409     return dca;
0410 }
0411
0412 void RecoTracklets(TrackletData &tkldata)
0413 {
0414     sort(tkldata.ProtoTkls.begin(), tkldata.ProtoTkls.end(), compare_dR);
0415
0416     int itkl = 0;
0417     for (auto &tkl : tkldata.ProtoTkls)
0418     {
0419         if (tkl->Hit1()->IsMatchedTkl() || tkl->Hit2()->IsMatchedTkl())
0420         {
0421             continue;
0422         }
0423         else
0424         {
0425             if (!tkldata.isdata)
0426                 tkldata.RecoTkls.push_back(tkl);
0427             tkl->Hit1()->SetMatchedTkl();
0428             tkl->Hit2()->SetMatchedTkl();
0429
0430             tkldata.recotklraw_eta.push_back((tkl->Hit1()->Eta()));
0431             tkldata.recotklraw_phi.push_back(tkl->Hit1()->Phi());
0432             tkldata.recotklraw_deta.push_back(tkl->dEta());
0433             tkldata.recotklraw_dphi.push_back(tkl->dPhi());
0434             tkldata.recotklraw_dR.push_back(tkl->dR());
0435
0436             tkldata.tklclus1x.push_back(tkl->Hit1()->posX());
0437             tkldata.tklclus1y.push_back(tkl->Hit1()->posY());
0438             tkldata.tklclus1z.push_back(tkl->Hit1()->posZ());
0439             tkldata.tklclus2x.push_back(tkl->Hit2()->posX());
0440             tkldata.tklclus2y.push_back(tkl->Hit2()->posY());
0441             tkldata.tklclus2z.push_back(tkl->Hit2()->posZ());
0442             tkldata.tklclus1phi.push_back(tkl->Hit1()->Phi());
0443             tkldata.tklclus1eta.push_back(tkl->Hit1()->Eta());
0444             tkldata.tklclus2phi.push_back(tkl->Hit2()->Phi());
0445             tkldata.tklclus2eta.push_back(tkl->Hit2()->Eta());
0446             tkldata.tklclus1phisize.push_back(tkl->Hit1()->PhiSize());
0447             tkldata.tklclus2phisize.push_back(tkl->Hit2()->PhiSize());
0448             tkldata.tklclus1adc.push_back(tkl->Hit1()->ClusADC());
0449             tkldata.tklclus2adc.push_back(tkl->Hit2()->ClusADC());
0450
0451             // calculate DCA w.r.t the event vertex PV_x, PV_y, PV_z
0452             // set up Point3D objects for the 3D DCA calculation
0453             Point3D p1 = {tkl->Hit1()->posX(), tkl->Hit1()->posY(), tkl->Hit1()->posZ()};
0454             Point3D p2 = {tkl->Hit2()->posX(), tkl->Hit2()->posY(), tkl->Hit2()->posZ()};
0455             Point3D p3 = {tkldata.PV_x, tkldata.PV_y, tkldata.PV_z};
0456             tkldata.recotklraw_dca3dvtx.push_back(computeDCA(p1, p2, p3));
0457         }
0458     }
0459
0460     tkldata.NRecotkl_Raw = tkldata.recotklraw_eta.size();
0461 }
0462
0463 // function to print the size of all the vectors in the TrackletData object
0464 void PrintVecSize(TrackletData &tkldata)
0465 {
0466     cout << "Size of the vectors in the TrackletData object:" << endl;
0467     cout << "layers[0]: " << tkldata.layers[0].size() << endl;
0468     cout << "layers[1]: " << tkldata.layers[1].size() << endl;
0469     cout << "ProtoTkls: " << tkldata.ProtoTkls.size() << endl;
0470     cout << "RecoTkls: " << tkldata.RecoTkls.size() << endl;
0471     cout << "GenHadrons: " << tkldata.GenHadrons.size() << endl;
0472     cout << "cluslayer: " << tkldata.cluslayer.size() << endl;
0473     cout << "clusphi: " << tkldata.clusphi.size() << endl;
0474     cout << "cluseta: " << tkldata.cluseta.size() << endl;
0475     cout << "clusphisize: " << tkldata.clusphisize.size() << endl;
0476     cout << "cluszsize: " << tkldata.cluszsize.size() << endl;
0477     cout << "clusadc: " << tkldata.clusadc.size() << endl;
0478     cout << "tklclus1phi: " << tkldata.tklclus1phi.size() << endl;
0479     cout << "tklclus1eta: " << tkldata.tklclus1eta.size() << endl;
0480     cout << "tklclus2phi: " << tkldata.tklclus2phi.size() << endl;
0481     cout << "tklclus2eta: " << tkldata.tklclus2eta.size() << endl;
0482     cout << "tklclus1phisize: " << tkldata.tklclus1phisize.size() << endl;
0483     cout << "tklclus2phisize: " << tkldata.tklclus2phisize.size() << endl;
0484     cout << "tklclus1adc: " << tkldata.tklclus1adc.size() << endl;
0485     cout << "tklclus2adc: " << tkldata.tklclus2adc.size() << endl;
0486     cout << "prototkl_eta: " << tkldata.prototkl_eta.size() << endl;
0487     cout << "prototkl_phi: " << tkldata.prototkl_phi.size() << endl;
0488     cout << "prototkl_deta: " << tkldata.prototkl_deta.size() << endl;
0489     cout << "prototkl_dphi: " << tkldata.prototkl_dphi.size() << endl;
0490     cout << "prototkl_dR: " << tkldata.prototkl_dR.size() << endl;
0491     cout << "recotklraw_eta: " << tkldata.recotklraw_eta.size() << endl;
0492     cout << "recotklraw_phi: " << tkldata.recotklraw_phi.size() << endl;
0493     cout << "recotklraw_deta: " << tkldata.recotklraw_deta.size() << endl;
0494     cout << "recotklraw_dphi: " << tkldata.recotklraw_dphi.size() << endl;
0495     cout << "recotklraw_dR: " << tkldata.recotklraw_dR.size() << endl;
0496     cout << "GenHadron_PID: " << tkldata.GenHadron_PID.size() << endl;
0497     cout << "GenHadron_Pt: " << tkldata.GenHadron_Pt.size() << endl;
0498     cout << "GenHadron_eta: " << tkldata.GenHadron_eta.size() << endl;
0499     cout << "GenHadron_phi: " << tkldata.GenHadron_phi.size() << endl;
0500 }
0501
0502 // check if the two clusters in a tracklet are matched to the same G4 particle (primary charged particle, GenHadron_trackID)
0503 // if they are, set the matchedGChadron values in the Tracklet object
0504 void RecoTkl_isG4P(TrackletData &tkldata)
0505 {
0506     // std::set for the matched G4P track ID
0507     std::set<int> G4PTrackID_tklmatched;
0508     std::set<int> G4PAncestorTrackID_tklmatched;
0509     for (auto &tkl : tkldata.RecoTkls)
0510     {
0511         tkldata.recotkl_clus1_matchedtrackID.push_back(tkl->Hit1()->MatchedG4P_trackID());
0512         tkldata.recotkl_clus2_matchedtrackID.push_back(tkl->Hit2()->MatchedG4P_trackID());
0513         tkldata.recotkl_clus1_matchedAncestorTrackID.push_back(tkl->Hit1()->MatchedG4P_ancestor_trackID());
0514         tkldata.recotkl_clus2_matchedAncestorTrackID.push_back(tkl->Hit2()->MatchedG4P_ancestor_trackID());
0515         // std::cout << "Tracklet (delta phi, delta eta, delta R) = (" << tkl->dPhi() << ", " << tkl->dEta() << ", " << tkl->dR() << ")" << std::endl;
0516         // std::cout << "Tracklet cluster 1: matched G4P track ID = " << tkl->Hit1()->MatchedG4P_trackID() << " ; cluster 2: matched G4P track ID = " << tkl->Hit2()->MatchedG4P_trackID() << std::endl;
0517         // std::cout << "Tracklet cluster 1: matched G4P ancestor track ID = " << tkl->Hit1()->MatchedG4P_ancestor_trackID() << " ; cluster 2: matched G4P ancestor track ID = " << tkl->Hit2()->MatchedG4P_ancestor_trackID() << std::endl //
0518         //   << "--------------------------------" << std::endl;
0519
0520         // check if the two clusters in a tracklet are matched to the same G4 particle
0521         if ((tkl->Hit1()->MatchedG4P_trackID() == tkl->Hit2()->MatchedG4P_trackID()) && (tkl->Hit1()->MatchedG4P_trackID() != std::numeric_limits<int>::max()))
0522         {
0523             // see if the match G4P track ID is in the G4P_trackID vector
0524             auto it_pg4p = std::find(tkldata.G4P_trackID.begin(), tkldata.G4P_trackID.end(), tkl->Hit1()->MatchedG4P_trackID());
0525             if (it_pg4p != tkldata.G4P_trackID.end())
0526             {
0527                 int index = std::distance(tkldata.G4P_trackID.begin(), it_pg4p);
0528                 tkldata.recotkl_matchedPG4P_eta.push_back(tkldata.G4P_eta[index]);
0529                 tkldata.recotkl_matchedPG4P_phi.push_back(tkldata.G4P_phi[index]);
0530
0531                 // add the matched G4P track ID to the set
0532                 G4PTrackID_tklmatched.insert(tkl->Hit1()->MatchedG4P_trackID());
0533                 G4PAncestorTrackID_tklmatched.insert(tkl->Hit1()->MatchedG4P_ancestor_trackID());
0534             }
0535             else
0536             {
0537                 tkldata.recotkl_matchedPG4P_eta.push_back(std::numeric_limits<float>::quiet_NaN());
0538                 tkldata.recotkl_matchedPG4P_phi.push_back(std::numeric_limits<float>::quiet_NaN());
0539             }
0540
0541             // see if the matched G4P track ID is in the GenHadron_trackID vector
0542             auto it_genhad = std::find(tkldata.GenHadron_trackID.begin(), tkldata.GenHadron_trackID.end(), tkl->Hit1()->MatchedG4P_trackID());
0543             if (it_genhad != tkldata.GenHadron_trackID.end())
0544             {
0545                 tkl->SetMatchedGenHardon();
0546                 tkldata.recotkl_isMatchedGChadron.push_back(true);
0547                 tkldata.matchedGChadron_pt.push_back(tkl->Hit1()->MatchedG4P_Pt());
0548                 tkldata.matchedGChadron_eta.push_back(tkl->Hit1()->MatchedG4P_Eta());
0549                 tkldata.matchedGChadron_phi.push_back(tkl->Hit1()->MatchedG4P_Phi());
0550
0551                 tkldata.recotkl_matchedGChadron_eta.push_back(tkl->Hit1()->MatchedG4P_Eta());
0552                 tkldata.recotkl_matchedGChadron_phi.push_back(tkl->Hit1()->MatchedG4P_Phi());
0553             }
0554             else // it's not generated charged hadron
0555             {
0556                 tkldata.recotkl_isMatchedGChadron.push_back(false);
0557                 // if the matched G4P track ID is not in the GenHadron_trackID vector, set the matchedGChadron values to NAN
0558                 tkldata.matchedGChadron_pt.push_back(std::numeric_limits<float>::quiet_NaN());
0559                 tkldata.matchedGChadron_eta.push_back(std::numeric_limits<float>::quiet_NaN());
0560                 tkldata.matchedGChadron_phi.push_back(std::numeric_limits<float>::quiet_NaN());
0561
0562                 tkldata.recotkl_matchedGChadron_eta.push_back(std::numeric_limits<float>::quiet_NaN());
0563                 tkldata.recotkl_matchedGChadron_phi.push_back(std::numeric_limits<float>::quiet_NaN());
0564             }
0565         }
0566         else
0567         {
0568             tkldata.recotkl_isMatchedGChadron.push_back(false);
0569             // if the two clusters in a tracklet are not matched to the same G4 particle, set the matchedGChadron values to NAN
0570             tkldata.matchedGChadron_pt.push_back(std::numeric_limits<float>::quiet_NaN());
0571             tkldata.matchedGChadron_eta.push_back(std::numeric_limits<float>::quiet_NaN());
0572             tkldata.matchedGChadron_phi.push_back(std::numeric_limits<float>::quiet_NaN());
0573
0574             tkldata.recotkl_matchedPG4P_eta.push_back(std::numeric_limits<float>::quiet_NaN());
0575             tkldata.recotkl_matchedPG4P_phi.push_back(std::numeric_limits<float>::quiet_NaN());
0576             tkldata.recotkl_matchedGChadron_eta.push_back(std::numeric_limits<float>::quiet_NaN());
0577             tkldata.recotkl_matchedGChadron_phi.push_back(std::numeric_limits<float>::quiet_NaN());
0578         }
0579     }
0580
0581     // Check for each element in the set if it is in the PrimaryG4P_trackID vector.
0582     // If it is, set the corresponding element in PrimaryG4P_IsRecotkl to true
0583     std::cout << "G4PTrackID_tklmatched size: " << G4PTrackID_tklmatched.size() << std::endl;
0584     for (auto &trackID : G4PTrackID_tklmatched)
0585     {
0586         auto it = std::find(tkldata.PrimaryG4P_trackID.begin(), tkldata.PrimaryG4P_trackID.end(), trackID);
0587         if (it != tkldata.PrimaryG4P_trackID.end())
0588         {
0589             int index = std::distance(tkldata.PrimaryG4P_trackID.begin(), it);
0590             tkldata.PrimaryG4P_IsRecotkl[index] = true;
0591         }
0592     }
0593
0594     // Check for each element in the set if it is in the GenHadron_trackID vector.
0595     // If it is, set the corresponding element in GenHadron_IsRecotkl to true
0596     std::cout << "G4PAncestorTrackID_tklmatched size: " << G4PAncestorTrackID_tklmatched.size() << std::endl;
0597     for (auto &trackID : G4PAncestorTrackID_tklmatched)
0598     {
0599         auto it = std::find(tkldata.GenHadron_trackID.begin(), tkldata.GenHadron_trackID.end(), trackID);
0600         if (it != tkldata.GenHadron_trackID.end())
0601         {
0602             int index = std::distance(tkldata.GenHadron_trackID.begin(), it);
0603             tkldata.GenHadron_IsRecotkl[index] = true;
0604         }
0605     }
0606
0607     // clear the set
0608     G4PTrackID_tklmatched.clear();
0609     G4PAncestorTrackID_tklmatched.clear();
0610 }
0611
0612 #endif