Fork me on GitHub

Changeset 4ac0049 in git for modules


Ignore:
Timestamp:
Jan 24, 2020, 3:51:00 PM (5 years ago)
Author:
Kaan Yüksel Oyulmaz <kaanyukseloyulmaz@…>
Branches:
Timing
Children:
79a7b3e
Parents:
6049672
Message:

VertexFinderDA4D.cc module is cleaned from graphs, FCChh_PileUpVtx.tcl card, PileUpSubtractor4D.cc and TimeSmearing.cc were updated in meeting

Location:
modules
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • modules/PileUpSubtractor4D.cc

    r6049672 r4ac0049  
    170170      t_err = particle->PositionError.T();
    171171
    172       distanceCharged = pow((zvtx - z),2)/pow((zvtx_err - z_err),2) + pow((tvtx - t),2)/pow((tvtx_err - t_err),2);
    173       distanceNeutral = pow((tvtx - t),2)/pow((tvtx_err - t_err),2);
     172      distanceCharged = TMath::Sqrt(pow((zvtx - z),2)/pow((zvtx_err),2) + pow((tvtx - t),2)/pow((tvtx_err),2));
     173      distanceNeutral = TMath::Sqrt(pow((tvtx - t),2)/pow((tvtx_err),2));
    174174
    175175      if(candidate->Charge != 0 && distanceCharged < fChargedMinSignificance)
    176176      {
    177         candidate->IsRecoPU = 1;
     177        candidate->IsRecoPU = 0;
    178178      }
    179179      else if(candidate->Charge == 0 && distanceNeutral < fNeutralMinSignificance)
    180180      {
    181         candidate->IsRecoPU = 1;
     181        candidate->IsRecoPU = 0;
    182182      } 
    183183      else
    184184      {
    185         candidate->IsRecoPU = 0;
     185        candidate->IsRecoPU = 1;
    186186        if(candidate->Momentum.Pt() > fPTMin) array->Add(candidate);
    187187      }

  • modules/TimeSmearing.cc

    r6049672 r4ac0049  
    125125      candidate = static_cast<Candidate*>(candidate->Clone());  // I am not sure that we need these lines !!!
    126126      candidate->AddCandidate(mother);
    127       candidate->InitialPosition.SetT((100+ti)*1.0E3*c_light);
     127      candidate->InitialPosition.SetT((100+ti)*1.0E3*c_light);  // Dummy Value, correct value will be computed by VertexFinderDA4D
    128128      candidate->Position.SetT(tf_smeared*1.0E3*c_light);
    129129      candidate->ErrorT = timeResolution*1.0E3*c_light;

  • modules/VertexFinderDA4D.cc

    r6049672 r4ac0049  
    275275  }
    276276
    277   if( fVerbose > 10 ) plot_status(beta, vtx, tks, 0, "Ast");
    278 
    279277  if( fVerbose > 2){cout << "Cool down untill reaching the temperature to finish increasing the number of vertexes" << endl;}
    280278
    … …  
    290288      delta2 = update(beta, tks, vtx, rho0);
    291289
    292       if( fVerbose > 10 ) plot_status(beta, vtx, tks, niter, "Bup");
    293290      if (fVerbose > 3)
    294291      {
    … …  
    312309      n_it++;
    313310
    314       if( fVerbose > 10 ) plot_status(beta, vtx, tks, n_it, "Cme");
    315311    }
    316312
    … …  
    320316    {
    321317      split(beta, vtx, tks);
    322       if( fVerbose > 10 ) plot_status(beta, vtx, tks, 0, "Asp");
    323318    }
    324319    else
    … …  
    355350    }
    356351    while (delta2 > 0.3*fD2UpdateLim &&  niter < fMaxIterations);
    357     if( fVerbose > 10 ) plot_status(beta, vtx, tks, f, "Dadout");
    358352  }
    359353
    … …  
    373367        }
    374368        while (delta2 > fD2UpdateLim &&  niter < fMaxIterations);
    375         if( fVerbose > 10 ) plot_status(beta, vtx, tks, i_pu, Form("Eprg%d",min_trk));
    376369        i_pu++;
    377370      }
    … …  
    390383      n_it++;
    391384
    392       if( fVerbose > 10 ) plot_status(beta, vtx, tks, n_it, "Cme");
    393385    }
    394386  } while( beta < fBetaPurge );
    … …  
    404396      delta2 = update(beta, tks, vtx, rho0);
    405397      niter++;
    406       if( fVerbose > 10 ) plot_status(beta, vtx, tks, 0, "Bup");
    407398    }
    408399    while (delta2 > 0.3*fD2UpdateLim &&  niter < fMaxIterations);
    … …  
    484475    fTrackOutputArray->Add(tks.tt[i]);
    485476  }
    486 
    487   if(fVerbose > 10) plot_status_end(vtx, tks);
    488477
    489478}
    … …  
    878867        {
    879868          continue;
    880           // plot_split_crush(zn, tn, vtx, tks, k);
    881869          // throw std::invalid_argument( "0 division" );
    882870        }
    … …  
    10431031  }
    10441032}
    1045 
    1046 
    1047 // -----------------------------------------------------------------------------
    1048 // Plot status
    1049 void VertexFinderDA4D::plot_status(double beta, vertex_t &vtx, tracks_t &tks, int n_it, const char* flag)
    1050 {
    1051   vector<int> vtx_color = {2,4,8,1,5,6,9,14,46,3};
    1052   while(vtx.getSize() > vtx_color.size()) vtx_color.push_back(40);
    1053 
    1054   vector<double> t_PV, dt_PV, z_PV, dz_PV;
    1055   vector<double> t_PU, dt_PU, z_PU, dz_PU;
    1056 
    1057   double ETot = 0;
    1058   vector<double> pk_exp_mBetaE = Compute_pk_exp_mBetaE(beta, vtx, tks, 0);
    1059 
    1060   for(unsigned int i = 0; i < tks.getSize(); i++)
    1061   {
    1062     for(unsigned int k = 0; k < vtx.getSize(); k++)
    1063     {
    1064       unsigned int idx = k*tks.getSize() + i;
    1065       if(pk_exp_mBetaE[idx] == 0) continue;
    1066 
    1067       double p_ygx = pk_exp_mBetaE[idx] / tks.Z[i];
    1068 
    1069       ETot += tks.w[i] * p_ygx * Energy(tks.z[i], vtx.z[k], tks.dz2_o[i], tks.t[i], vtx.t[k], tks.dt2_o[i]);
    1070     }
    1071 
    1072     if(tks.tt[i]->IsPU)
    1073     {
    1074       t_PU.push_back(tks.t[i]);
    1075       dt_PU.push_back(1./tks.dt_o[i]);
    1076       z_PU.push_back(tks.z[i]);
    1077       dz_PU.push_back(1./tks.dz_o[i]);
    1078     }
    1079     else
    1080     {
    1081       t_PV.push_back(tks.t[i]);
    1082       dt_PV.push_back(1./tks.dt_o[i]);
    1083       z_PV.push_back(tks.z[i]);
    1084       dz_PV.push_back(1./tks.dz_o[i]);
    1085     }
    1086   }
    1087 
    1088 
    1089   ETot /= tks.sum_w;
    1090   fEnergy_rec.push_back(ETot);
    1091   fBeta_rec.push_back(beta);
    1092   fNvtx_rec.push_back(vtx.getSize());
    1093 
    1094   double t_min = TMath::Min(  TMath::MinElement(t_PV.size(), &t_PV[0]), TMath::MinElement(t_PU.size(), &t_PU[0])  );
    1095   t_min = TMath::Min(t_min, TMath::MinElement(vtx.getSize(), &(vtx.t[0]))  ) - fVertexTSize;
    1096   double t_max = TMath::Max(  TMath::MaxElement(t_PV.size(), &t_PV[0]), TMath::MaxElement(t_PU.size(), &t_PU[0])  );
    1097   t_max = TMath::Max(t_max, TMath::MaxElement(vtx.getSize(), &(vtx.t[0]))  ) + fVertexTSize;
    1098 
    1099   double z_min = TMath::Min(  TMath::MinElement(z_PV.size(), &z_PV[0]), TMath::MinElement(z_PU.size(), &z_PU[0])  );
    1100   z_min = TMath::Min(z_min, TMath::MinElement(vtx.getSize(), &(vtx.z[0]))  ) - 5;
    1101   double z_max = TMath::Max(  TMath::MaxElement(z_PV.size(), &z_PV[0]), TMath::MaxElement(z_PU.size(), &z_PU[0])  );
    1102   z_max = TMath::Max(z_max, TMath::MaxElement(vtx.getSize(), &(vtx.z[0]))  ) + 5;
    1103 
    1104   auto c_2Dspace = new TCanvas("c_2Dspace", "c_2Dspace", 800, 600);
    1105 
    1106   TGraphErrors* gr_PVtks = new TGraphErrors(t_PV.size(), &t_PV[0], &z_PV[0], &dt_PV[0], &dz_PV[0]);
    1107   gr_PVtks->SetTitle(Form("Clustering space - #beta = %.6f", beta));
    1108   gr_PVtks->GetXaxis()->SetTitle("t CA [ps]");
    1109   gr_PVtks->GetXaxis()->SetLimits(t_min, t_max);
    1110   gr_PVtks->GetYaxis()->SetTitle("z CA [mm]");
    1111   gr_PVtks->GetYaxis()->SetRangeUser(z_min, z_max);
    1112   gr_PVtks->SetMarkerStyle(4);
    1113   gr_PVtks->SetMarkerColor(8);
    1114   gr_PVtks->SetLineColor(8);
    1115   gr_PVtks->Draw("APE1");
    1116 
    1117   TGraphErrors* gr_PUtks = new TGraphErrors(t_PU.size(), &t_PU[0], &z_PU[0], &dt_PU[0], &dz_PU[0]);
    1118   gr_PUtks->SetMarkerStyle(3);
    1119   gr_PUtks->Draw("PE1");
    1120 
    1121   TGraph* gr_vtx = new TGraph(vtx.getSize(), &(vtx.t[0]), &(vtx.z[0]));
    1122   gr_vtx->SetMarkerStyle(28);
    1123   gr_vtx->SetMarkerColor(2);
    1124   gr_vtx->SetMarkerSize(2.);
    1125   gr_vtx->Draw("PE1");
    1126 
    1127   fItInputGenVtx->Reset();
    1128   TGraph* gr_genvtx = new TGraph(fInputGenVtx->GetEntriesFast());
    1129   Candidate *candidate;
    1130   unsigned int k = 0;
    1131   while((candidate = static_cast<Candidate*>(fItInputGenVtx->Next())))
    1132   {
    1133     gr_genvtx->SetPoint(k, candidate->Position.T()*1E9/c_light, candidate->Position.Z());
    1134     k++;
    1135   }
    1136   gr_genvtx->SetMarkerStyle(33);
    1137   gr_genvtx->SetMarkerColor(6);
    1138   gr_genvtx->SetMarkerSize(2.);
    1139   gr_genvtx->Draw("PE1");
    1140 
    1141   // auto leg = new TLegend(0.1, 0.1);
    1142   // leg->AddEntry(gr_PVtks, "PV tks", "ep");
    1143   // leg->AddEntry(gr_PUtks, "PU tks", "ep");
    1144   // leg->AddEntry(gr_vtx, "Cluster center", "p");
    1145   // leg->Draw();
    1146 
    1147   c_2Dspace->SetGrid();
    1148   c_2Dspace->SaveAs(fFigFolderPath + Form("/c_2Dspace_beta%010.0f-%s%d.png", 1E7*beta, flag, n_it));
    1149 
    1150   delete c_2Dspace;
    1151 }
    1152 
    1153 // -----------------------------------------------------------------------------
    1154 // Plot status at the end
    1155 void VertexFinderDA4D::plot_status_end(vertex_t &vtx, tracks_t &tks)
    1156 {
    1157   unsigned int nv = vtx.getSize();
    1158 
    1159   // Define colors in a meaningfull way
    1160   vector<int> MyPalette(nv);
    1161 
    1162   const int Number = 3;
    1163   double Red[Number]    = { 1.00, 0.00, 0.00};
    1164   double Green[Number]  = { 0.00, 1.00, 0.00};
    1165   double Blue[Number]   = { 1.00, 0.00, 1.00};
    1166   double Length[Number] = { 0.00, 0.50, 1.00 };
    1167   int FI = TColor::CreateGradientColorTable(Number,Length,Red,Green,Blue,nv);
    1168   for (unsigned int i=0;i<nv;i++) MyPalette[i] = FI+i;
    1169 
    1170   TCanvas * c_out = new TCanvas("c_out", "c_out", 800, 600);
    1171   double t_min = TMath::Min( TMath::MinElement(tks.getSize(), &tks.t[0]), TMath::MinElement(vtx.getSize(), &(vtx.t[0]))  ) - 2*fVertexTSize;
    1172   double t_max = TMath::Max(TMath::MaxElement(tks.getSize(), &tks.t[0]), TMath::MaxElement(vtx.getSize(), &(vtx.t[0]))  ) + 2*fVertexTSize;
    1173 
    1174   double z_min = TMath::Min( TMath::MinElement(tks.getSize(), &tks.z[0]), TMath::MinElement(vtx.getSize(), &(vtx.z[0]))  ) - 15;
    1175   double z_max = TMath::Max( TMath::MaxElement(tks.getSize(), &tks.z[0]), TMath::MaxElement(vtx.getSize(), &(vtx.z[0]))  ) + 15;
    1176 
    1177   // Draw tracks
    1178   for(unsigned int i = 0; i < tks.getSize(); i++)
    1179   {
    1180     double dt[] = {1./tks.dt_o[i]};
    1181     double dz[] = {1./tks.dz_o[i]};
    1182     TGraphErrors* gr = new TGraphErrors(1, &(tks.t[i]), &(tks.z[i]), dt, dz);
    1183 
    1184     gr->SetNameTitle(Form("gr%d",i), Form("gr%d",i));
    1185 
    1186     int marker = tks.tt[i]->IsPU? 1 : 4;
    1187     gr->SetMarkerStyle(marker);
    1188 
    1189     int idx = tks.tt[i]->ClusterIndex;
    1190     int color = idx>=0 ? MyPalette[idx] : 13;
    1191     gr->SetMarkerColor(color);
    1192     gr->SetLineColor(color);
    1193 
    1194     int line_style = idx>=0 ? 1 : 3;
    1195     gr->SetLineStyle(line_style);
    1196 
    1197     if(i==0)
    1198     {
    1199       gr->SetTitle(Form("Clustering space - Tot Vertexes = %d", nv));
    1200       gr->GetXaxis()->SetTitle("t CA [ps]");
    1201       gr->GetXaxis()->SetLimits(t_min, t_max);
    1202       gr->GetYaxis()->SetTitle("z CA [mm]");
    1203       gr->GetYaxis()->SetRangeUser(z_min, z_max);
    1204       gr->Draw("APE1");
    1205     }
    1206     else gr->Draw("PE1");
    1207   }
    1208 
    1209   // Draw vertices
    1210   for(unsigned int k = 0; k < vtx.getSize(); k++)
    1211   {
    1212     TGraph* gr = new TGraph(1, &(vtx.t[k]), &(vtx.z[k]));
    1213 
    1214     gr->SetNameTitle(Form("grv%d",k), Form("grv%d",k));
    1215 
    1216     gr->SetMarkerStyle(41);
    1217     gr->SetMarkerSize(2.);
    1218     gr->SetMarkerColor(MyPalette[k]);
    1219 
    1220     gr->Draw("P");
    1221   }
    1222 
    1223   fItInputGenVtx->Reset();
    1224   TGraph* gr_genvtx = new TGraph(fInputGenVtx->GetEntriesFast());
    1225   TGraph* gr_genPV = new TGraph(1);
    1226   Candidate *candidate;
    1227   unsigned int k = 0;
    1228   while((candidate = static_cast<Candidate*>(fItInputGenVtx->Next())))
    1229   {
    1230     if(k == 0 ) {
    1231       gr_genPV->SetPoint(k, candidate->Position.T()*1E9/c_light, candidate->Position.Z());
    1232     }
    1233     else gr_genvtx->SetPoint(k, candidate->Position.T()*1E9/c_light, candidate->Position.Z());
    1234 
    1235     k++;
    1236   }
    1237   gr_genvtx->SetMarkerStyle(20);
    1238   gr_genvtx->SetMarkerColorAlpha(kBlack, 0.8);
    1239   gr_genvtx->SetMarkerSize(.8);
    1240   gr_genvtx->Draw("PE1");
    1241   gr_genPV->SetMarkerStyle(33);
    1242   gr_genPV->SetMarkerColorAlpha(kBlack, 1);
    1243   gr_genPV->SetMarkerSize(2.5);
    1244   gr_genPV->Draw("PE1");
    1245 
    1246   // auto note =  new TLatex();
    1247   // note->DrawLatexNDC(0.5, 0.8, Form("#splitline{Vertexes Reco = %d }{Vertexes gen = %d}", vtx.getSize(), k) );
    1248 
    1249   c_out->SetGrid();
    1250   c_out->SaveAs(fFigFolderPath + Form("/c_final.root"));
    1251   delete c_out;
    1252 }
    1253 
    1254 // -----------------------------------------------------------------------------
    1255 // Plot splitting
    1256 void VertexFinderDA4D::plot_split_crush(double zn, double tn, vertex_t &vtx, tracks_t &tks, int i_vtx)
    1257 {
    1258   vector<double> t, dt, z, dz;
    1259 
    1260   for(unsigned int i = 0; i < tks.getSize(); i++)
    1261   {
    1262       t.push_back(tks.t[i]);
    1263       dt.push_back(1./tks.dt_o[i]);
    1264       z.push_back(tks.z[i]);
    1265       dz.push_back(1./tks.dz_o[i]);
    1266   }
    1267 
    1268 
    1269   double t_min = TMath::Min(TMath::MinElement(t.size(), &t[0]), TMath::MinElement(vtx.getSize(), &(vtx.t[0]))  ) - 50;
    1270   double t_max = TMath::Max(TMath::MaxElement(t.size(), &t[0]), TMath::MaxElement(vtx.getSize(), &(vtx.t[0]))  ) + 50;
    1271 
    1272   double z_min = TMath::Min(TMath::MinElement(z.size(), &z[0]), TMath::MinElement(vtx.getSize(), &(vtx.z[0]))  ) - 5;
    1273   double z_max = TMath::Max(TMath::MaxElement(z.size(), &z[0]), TMath::MaxElement(vtx.getSize(), &(vtx.z[0]))  ) + 5;
    1274 
    1275   auto c_2Dspace = new TCanvas("c_2Dspace", "c_2Dspace", 800, 600);
    1276 
    1277   TGraphErrors* gr_PVtks = new TGraphErrors(t.size(), &t[0], &z[0], &dt[0], &dz[0]);
    1278   gr_PVtks->SetTitle(Form("Clustering space"));
    1279   gr_PVtks->GetXaxis()->SetTitle("t CA [ps]");
    1280   gr_PVtks->GetXaxis()->SetLimits(t_min, t_max);
    1281   gr_PVtks->GetYaxis()->SetTitle("z CA [mm]");
    1282   gr_PVtks->GetYaxis()->SetRangeUser(z_min, z_max);
    1283   gr_PVtks->SetMarkerStyle(4);
    1284   gr_PVtks->SetMarkerColor(1);
    1285   gr_PVtks->SetLineColor(1);
    1286   gr_PVtks->Draw("APE1");
    1287 
    1288   TGraph* gr_vtx = new TGraph(1, &(vtx.t[i_vtx]), &(vtx.z[i_vtx]));
    1289   gr_vtx->SetMarkerStyle(28);
    1290   gr_vtx->SetMarkerColor(2);
    1291   gr_vtx->SetMarkerSize(2.);
    1292   gr_vtx->Draw("PE1");
    1293 
    1294   double t_pos[] = {vtx.t[i_vtx], vtx.t[i_vtx]+100};
    1295   double t_neg[] = {vtx.t[i_vtx], vtx.t[i_vtx]-100};
    1296   double z_pos[] = {vtx.z[i_vtx], vtx.z[i_vtx]+(zn/tn)*100};
    1297   double z_neg[] = {vtx.z[i_vtx], vtx.z[i_vtx]-(zn/tn)*100};
    1298 
    1299   TGraph* gr_pos = new TGraph(2, &t_pos[0], &z_pos[0]);
    1300   gr_pos->SetLineColor(8);
    1301   gr_pos->SetMarkerColor(8);
    1302   gr_pos->Draw("PL");
    1303   TGraph* gr_neg = new TGraph(2, &t_neg[0], &z_neg[0]);
    1304   gr_neg->SetLineColor(4);
    1305   gr_neg->SetMarkerColor(4);
    1306   gr_neg->Draw("PL");
    1307 
    1308 
    1309   c_2Dspace->SetGrid();
    1310   c_2Dspace->SaveAs(fFigFolderPath + Form("/crush_splitting.png"));
    1311 
    1312   delete c_2Dspace;
    1313 }

  • modules/VertexFinderDA4D.h

    r6049672 r4ac0049  
    234234    std::vector<double> Compute_pk_exp_mBetaE(double beta, vertex_t &vtx, tracks_t &tks, double Z_init);
    235235
    236     // Plot status of tracks and Vertices
    237     void plot_status(double beta, vertex_t &vtx, tracks_t &tks, int n_it = 0, const char* flag ="");
    238 
    239     // Plot status at the end of the fitting
    240     void plot_status_end(vertex_t &vtx, tracks_t &tks);
    241 
    242     // Plot Crush
    243     void plot_split_crush(double zn, double tn, vertex_t &vtx, tracks_t &tks, int i_vtx);
    244 
    245236
    246237  private:
Note: See TracChangeset for help on using the changeset viewer.