MVTXTracklet/src/Hit.h

0001 #ifndef HIT_H
0002 #define HIT_H
0003
0004 #include <algorithm>
0005 #include <cmath>
0006 #include <iostream>
0007 #include <stdio.h>
0008
0009 #include <TF1.h>
0010 #include <TH1F.h>
0011 #include <TObject.h>
0012 #include <TRandom3.h>
0013 #include <TVector3.h>
0014
0015 using namespace std;
0016
0017 class Hit : public TObject
0018 {
0019   public:
0020     Hit();
0021     Hit(float, float, float, float, float, float, int);
0022     Hit(float, float, float, float, float, float, int, float);
0023     Hit(float, float, float, float, float, float, int, float, unsigned int);
0024     Hit(float, float); // simple construct
0025     ~Hit();
0026
0027     float posX();
0028     float posY();
0029     float posZ();
0030     float rho();
0031     float vtxX();
0032     float vtxY();
0033     float vtxZ();
0034     float Eta();
0035     float Phi();
0036     float R();
0037     int Layer();
0038     float PhiSize() { return _phisize; }
0039     unsigned int ClusADC() { return _clusadc; }
0040     pair<float, float> Edge();
0041
0042     void Update();
0043     void SetPos(float, float, float);
0044     void SetVtx(float, float, float);
0045     void SetEdge(float, float);
0046     void SetMatchedTkl();
0047     bool IsMatchedTkl();
0048     void Print();
0049
0050     TVector3 VecPos();
0051     TVector3 VecVtx();
0052     TVector3 VecRel();
0053
0054   private:
0055     float _x;
0056     float _y;
0057     float _z;
0058     float _vtxX;
0059     float _vtxY;
0060     float _vtxZ;
0061     float _eta;
0062     float _phi;
0063     float _R;
0064     int _layer;
0065     float _phisize;
0066     unsigned int _clusadc;
0067     pair<float, float> _edge;
0068     bool _matched_tkl;
0069     TVector3 vechit;
0070     TVector3 vecvtx;
0071     TVector3 vecrel;
0072 };
0073
0074 Hit::Hit()
0075 {
0076     vechit.SetXYZ(0, 0, 0);
0077     vecvtx.SetXYZ(0, 0, 0);
0078     _layer = -999;
0079     vecrel = vechit - vecvtx;
0080     _eta = vecrel.Eta();
0081     _phi = vecrel.Phi();
0082     _R = vecrel.Mag();
0083     _matched_tkl = false;
0084 }
0085
0086 Hit::Hit(float x, float y, float z, float vtxX, float vtxY, float vtxZ, int layer)
0087 {
0088     vechit.SetXYZ(x, y, z);
0089     vecvtx.SetXYZ(vtxX, vtxY, vtxZ);
0090     _layer = layer;
0091     vecrel = vechit - vecvtx;
0092     _eta = vecrel.Eta();
0093     _phi = vecrel.Phi();
0094     _R = vecrel.Mag();
0095     _matched_tkl = false;
0096 }
0097
0098 Hit::Hit(float x, float y, float z, float vtxX, float vtxY, float vtxZ, int layer, float phisize)
0099 {
0100     vechit.SetXYZ(x, y, z);
0101     vecvtx.SetXYZ(vtxX, vtxY, vtxZ);
0102     _layer = layer;
0103     _phisize = phisize;
0104     vecrel = vechit - vecvtx;
0105     _eta = vecrel.Eta();
0106     _phi = vecrel.Phi();
0107     _R = vecrel.Mag();
0108     _matched_tkl = false;
0109 }
0110
0111 Hit::Hit(float x, float y, float z, float vtxX, float vtxY, float vtxZ, int layer, float phisize, unsigned int clusadc)
0112 {
0113     vechit.SetXYZ(x, y, z);
0114     vecvtx.SetXYZ(vtxX, vtxY, vtxZ);
0115     _layer = layer;
0116     _phisize = phisize;
0117     _clusadc = clusadc;
0118     vecrel = vechit - vecvtx;
0119     _eta = vecrel.Eta();
0120     _phi = vecrel.Phi();
0121     _R = vecrel.Mag();
0122     _matched_tkl = false;
0123 }
0124
0125 Hit::Hit(float eta, float phi)
0126 {
0127     _eta = eta;
0128     _phi = phi;
0129     _matched_tkl = false;
0130 }
0131
0132 Hit::~Hit() {}
0133
0134 float Hit::posX() { return (vechit.X()); }
0135
0136 float Hit::posY() { return (vechit.Y()); }
0137
0138 float Hit::posZ() { return (vechit.Z()); }
0139
0140 float Hit::rho() { return (sqrt(vechit.X() * vechit.X() + vechit.Y() * vechit.Y())); }
0141
0142 float Hit::vtxX() { return (vecvtx.X()); }
0143
0144 float Hit::vtxY() { return (vecvtx.Y()); }
0145
0146 float Hit::vtxZ() { return (vecvtx.Z()); }
0147
0148 float Hit::Eta() { return (_eta); } // with respect to the vertex that it associates to
0149
0150 float Hit::Phi() { return (_phi); } // with respect to the vertex that it associates to
0151
0152 float Hit::R() { return (_R); } // with respect to the vertex that it associates to
0153
0154 int Hit::Layer() { return (_layer); }
0155
0156 pair<float, float> Hit::Edge() { return (_edge); }
0157
0158 void Hit::SetPos(float x, float y, float z) { vechit.SetXYZ(x, y, z); }
0159
0160 void Hit::SetVtx(float vtxX, float vtxY, float vtxZ) { vecvtx.SetXYZ(vtxX, vtxY, vtxZ); }
0161
0162 void Hit::SetEdge(float edge1, float edge2) { _edge = make_pair(edge1, edge2); }
0163
0164 void Hit::Update()
0165 {
0166     vecrel = vechit - vecvtx;
0167     _eta = vecrel.Eta();
0168     _phi = vecrel.Phi();
0169     _R = vecrel.Mag();
0170 }
0171
0172 void Hit::SetMatchedTkl() { _matched_tkl = true; }
0173
0174 bool Hit::IsMatchedTkl() { return _matched_tkl; }
0175
0176 TVector3 Hit::VecPos() { return (vechit); }
0177
0178 TVector3 Hit::VecVtx() { return (vecvtx); }
0179
0180 TVector3 Hit::VecRel() { return (vecrel); }
0181
0182 void Hit::Print()
0183 {
0184     printf("[Hit::Print()] (posX, posY, posZ) = (%f, %f, %f), (vtxX, vtxY, vtxZ) = (%f, %f, %f), (eta, phi) = (%f, %f) \n", vechit.X(), vechit.Y(), vechit.Z(), vecvtx.X(), vecvtx.Y(), vecvtx.Z(),
0185            vecrel.Eta(), vecrel.Phi());
0186 }
0187
0188 void UpdateHits(vector<Hit *> &Hits, vector<float> PV)
0189 {
0190     for (auto &hit : Hits)
0191     {
0192         hit->SetVtx(PV[0], PV[1], PV[2]);
0193         hit->Update();
0194     }
0195 }
0196
0197 float RandomHit_fraction(int set)
0198 {
0199     float frac = 0;
0200     switch (set)
0201     {
0202     case 0:
0203         frac = 0;
0204         break;
0205     case 1:
0206         frac = 1;
0207         break;
0208     case 2:
0209         frac = 3;
0210         break;
0211     case 3:
0212         frac = 5;
0213         break;
0214     case 4:
0215         frac = 10;
0216         break;
0217     }
0218     return frac;
0219 }
0220
0221 Hit *RandomHit(float vx, float vy, float vz, int layer)
0222 {
0223     gRandom->SetSeed(0);
0224     // The 26 unique z positions for INTT
0225     vector<float> zpos = {-22.57245, -20.57245, -18.57245, -16.57245, -14.57245, -12.57245, -10.97245, -9.372450, -7.772450, -6.172450, -4.572450, -2.972450, -1.372450,
0226                           0.4275496, 2.0275495, 3.6275494, 5.2275495, 6.8275494, 8.4275493, 10.027549, 11.627549, 13.627549, 15.627549, 17.627550, 19.627550, 21.627550};
0227     int zpos_idx = gRandom->Integer(26);
0228
0229     float layer_radius[4] = {7.453, 8.046, 9.934, 10.569};
0230     float phiMin = -TMath::Pi();
0231     float phiMax = TMath::Pi();
0232     float phi = phiMin + (phiMax - phiMin) * gRandom->Rndm();
0233     float x = layer_radius[layer] * cos(phi);
0234     float y = layer_radius[layer] * sin(phi);
0235     float z = zpos[zpos_idx];                                 // layer_radius[layer] / tan(2 * atan(exp(-eta)));
0236     Hit *randhit = new Hit(x, y, z, vx, vy, vz, layer, 1, 1); // assign the phisize and clusadc to 1 for random clusters, which should be ok
0237
0238     return randhit;
0239 }
0240
0241 float theta2pseudorapidity(float theta) { return -1. * TMath::Log(TMath::Tan(theta / 2)); }
0242
0243 TF1 *ClusADCCut(int constscale, float etascale)
0244 {
0245     TF1 *f = new TF1("f", Form("%d*TMath::CosH(%f*x)", constscale, etascale), -10, 10);
0246     return f;
0247 }
0248
0249 TH1F *ClusADCCut_StepFunc_INTTPrivate()
0250 {
0251     std::vector<float> thetastep = {0.001, 15, 20, 25, 30, 35, 45, 55, 125, 135, 145, 150, 155, 160, 165, 179.999};
0252     std::reverse(thetastep.begin(), thetastep.end());
0253     std::vector<float> adccut_theta = {225, 165, 135, 120, 105, 90, 75, 60, 75, 90, 105, 120, 135, 165, 225};
0254     float etastep_array[thetastep.size()];
0255     for (int i = 0; i < thetastep.size(); i++)
0256     {
0257         etastep_array[i] = theta2pseudorapidity(thetastep[i] * TMath::Pi() / 180);
0258     }
0259
0260     TH1F *hm_cut_inttprivate = new TH1F("hm_cut_inttprivate", "hm_cut_inttprivate", thetastep.size() - 1, etastep_array);
0261     for (int j = 0; j < hm_cut_inttprivate->GetNbinsX(); j++)
0262     {
0263         hm_cut_inttprivate->SetBinContent(j + 1, adccut_theta[j]);
0264     }
0265
0266     return hm_cut_inttprivate;
0267 }
0268
0269 TH1F *ClusADCCut_StepFunc(int constscale, float etascale)
0270 {
0271     TF1 *f_cut = ClusADCCut(constscale, etascale);
0272
0273     std::vector<float> adcstep;
0274     for (int i = 0; i < 20; i++)
0275     {
0276         adcstep.push_back(20 + i * 30);
0277     }
0278     std::vector<float> etastep;
0279     for (int i = 0; i < 20; i++)
0280     {
0281         etastep.insert(etastep.begin(), f_cut->GetX(adcstep[i], -10, 0));
0282         etastep.push_back(f_cut->GetX(adcstep[i], 0, 10));
0283     }
0284
0285     etastep.erase(std::remove_if(etastep.begin(), etastep.end(), [](float x) { return std::isnan(x); }), etastep.end());
0286
0287     float etastep_array[etastep.size()];
0288     for (int i = 0; i < etastep.size(); i++)
0289     {
0290         etastep_array[i] = etastep[i];
0291     }
0292     TH1F *hm_cut = new TH1F("hm_cut", "hm_cut", etastep.size() - 1, etastep_array);
0293     for (int j = 0; j < hm_cut->GetNbinsX(); j++)
0294     {
0295         if (hm_cut->GetBinLowEdge(j + 1) < 0)
0296         {
0297             hm_cut->SetBinContent(j + 1, f_cut->Eval(hm_cut->GetBinLowEdge(j + 1)));
0298         }
0299         else
0300         {
0301             hm_cut->SetBinContent(j + 1, f_cut->Eval(hm_cut->GetBinCenter(j + 1) + hm_cut->GetBinWidth(j + 1) / 2));
0302         }
0303     }
0304
0305     return hm_cut;
0306 }
0307
0308 int ConstADCCut(int set)
0309 {
0310     float cut = 0;
0311     switch (set)
0312     {
0313     case 0:
0314         cut = 35;
0315         break;
0316     case 1:
0317         cut = 50;
0318         break;
0319     case 2:
0320         cut = 66;
0321         break;
0322     case 3:
0323         cut = -1;
0324     }
0325     return cut;
0326 }
0327
0328 #endif