Fork me on GitHub

Changeset 79a7b3e in git for modules/VertexFinderDA4D.cc


Ignore:
Timestamp:
Jan 24, 2020, 3:55:17 PM (5 years ago)
Author:
GitHub <noreply@…>
Branches:
Timing
Children:
62764fb
Parents:
364dbe1 (diff), 4ac0049 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
git-author:
Michele Selvaggi <michele.selvaggi@…> (01/24/20 15:55:17)
git-committer:
GitHub <noreply@…> (01/24/20 15:55:17)
Message:

Merge pull request #5 from kaanyuxel/master

New Time Smearing Module for Neutral Particles and FCC-hh Card

File:
1 edited

Legend:

Unmodified
Added
Removed

  • modules/VertexFinderDA4D.cc

    r364dbe1 r79a7b3e  
    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);
    … …  
    425416    candidate->ClusterIndex = k;
    426417    candidate->Position.SetXYZT(0.0, 0.0, vtx.z[k] , vtx.t[k]*1E-9*c_light);
     418    candidate->InitialPosition.SetXYZT(0.0, 0.0, vtx.z[k] , vtx.t[k]*1E-9*c_light);   
    427419    candidate->PositionError.SetXYZT(0.0, 0.0, fVertexZSize , fVertexTSize*1E-9*c_light);
    428420    candidate->SumPT2 = 0;
    … …  
    483475    fTrackOutputArray->Add(tks.tt[i]);
    484476  }
    485 
    486   if(fVerbose > 10) plot_status_end(vtx, tks);
    487477
    488478}
    … …  
    877867        {
    878868          continue;
    879           // plot_split_crush(zn, tn, vtx, tks, k);
    880869          // throw std::invalid_argument( "0 division" );
    881870        }
    … …  
    10421031  }
    10431032}
    1044 
    1045 
    1046 // -----------------------------------------------------------------------------
    1047 // Plot status
    1048 void VertexFinderDA4D::plot_status(double beta, vertex_t &vtx, tracks_t &tks, int n_it, const char* flag)
    1049 {
    1050   vector<int> vtx_color = {2,4,8,1,5,6,9,14,46,3};
    1051   while(vtx.getSize() > vtx_color.size()) vtx_color.push_back(40);
    1052 
    1053   vector<double> t_PV, dt_PV, z_PV, dz_PV;
    1054   vector<double> t_PU, dt_PU, z_PU, dz_PU;
    1055 
    1056   double ETot = 0;
    1057   vector<double> pk_exp_mBetaE = Compute_pk_exp_mBetaE(beta, vtx, tks, 0);
    1058 
    1059   for(unsigned int i = 0; i < tks.getSize(); i++)
    1060   {
    1061     for(unsigned int k = 0; k < vtx.getSize(); k++)
    1062     {
    1063       unsigned int idx = k*tks.getSize() + i;
    1064       if(pk_exp_mBetaE[idx] == 0) continue;
    1065 
    1066       double p_ygx = pk_exp_mBetaE[idx] / tks.Z[i];
    1067 
    1068       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]);
    1069     }
    1070 
    1071     if(tks.tt[i]->IsPU)
    1072     {
    1073       t_PU.push_back(tks.t[i]);
    1074       dt_PU.push_back(1./tks.dt_o[i]);
    1075       z_PU.push_back(tks.z[i]);
    1076       dz_PU.push_back(1./tks.dz_o[i]);
    1077     }
    1078     else
    1079     {
    1080       t_PV.push_back(tks.t[i]);
    1081       dt_PV.push_back(1./tks.dt_o[i]);
    1082       z_PV.push_back(tks.z[i]);
    1083       dz_PV.push_back(1./tks.dz_o[i]);
    1084     }
    1085   }
    1086 
    1087 
    1088   ETot /= tks.sum_w;
    1089   fEnergy_rec.push_back(ETot);
    1090   fBeta_rec.push_back(beta);
    1091   fNvtx_rec.push_back(vtx.getSize());
    1092 
    1093   double t_min = TMath::Min(  TMath::MinElement(t_PV.size(), &t_PV[0]), TMath::MinElement(t_PU.size(), &t_PU[0])  );
    1094   t_min = TMath::Min(t_min, TMath::MinElement(vtx.getSize(), &(vtx.t[0]))  ) - fVertexTSize;
    1095   double t_max = TMath::Max(  TMath::MaxElement(t_PV.size(), &t_PV[0]), TMath::MaxElement(t_PU.size(), &t_PU[0])  );
    1096   t_max = TMath::Max(t_max, TMath::MaxElement(vtx.getSize(), &(vtx.t[0]))  ) + fVertexTSize;
    1097 
    1098   double z_min = TMath::Min(  TMath::MinElement(z_PV.size(), &z_PV[0]), TMath::MinElement(z_PU.size(), &z_PU[0])  );
    1099   z_min = TMath::Min(z_min, TMath::MinElement(vtx.getSize(), &(vtx.z[0]))  ) - 5;
    1100   double z_max = TMath::Max(  TMath::MaxElement(z_PV.size(), &z_PV[0]), TMath::MaxElement(z_PU.size(), &z_PU[0])  );
    1101   z_max = TMath::Max(z_max, TMath::MaxElement(vtx.getSize(), &(vtx.z[0]))  ) + 5;
    1102 
    1103   auto c_2Dspace = new TCanvas("c_2Dspace", "c_2Dspace", 800, 600);
    1104 
    1105   TGraphErrors* gr_PVtks = new TGraphErrors(t_PV.size(), &t_PV[0], &z_PV[0], &dt_PV[0], &dz_PV[0]);
    1106   gr_PVtks->SetTitle(Form("Clustering space - #beta = %.6f", beta));
    1107   gr_PVtks->GetXaxis()->SetTitle("t CA [ps]");
    1108   gr_PVtks->GetXaxis()->SetLimits(t_min, t_max);
    1109   gr_PVtks->GetYaxis()->SetTitle("z CA [mm]");
    1110   gr_PVtks->GetYaxis()->SetRangeUser(z_min, z_max);
    1111   gr_PVtks->SetMarkerStyle(4);
    1112   gr_PVtks->SetMarkerColor(8);
    1113   gr_PVtks->SetLineColor(8);
    1114   gr_PVtks->Draw("APE1");
    1115 
    1116   TGraphErrors* gr_PUtks = new TGraphErrors(t_PU.size(), &t_PU[0], &z_PU[0], &dt_PU[0], &dz_PU[0]);
    1117   gr_PUtks->SetMarkerStyle(3);
    1118   gr_PUtks->Draw("PE1");
    1119 
    1120   TGraph* gr_vtx = new TGraph(vtx.getSize(), &(vtx.t[0]), &(vtx.z[0]));
    1121   gr_vtx->SetMarkerStyle(28);
    1122   gr_vtx->SetMarkerColor(2);
    1123   gr_vtx->SetMarkerSize(2.);
    1124   gr_vtx->Draw("PE1");
    1125 
    1126   fItInputGenVtx->Reset();
    1127   TGraph* gr_genvtx = new TGraph(fInputGenVtx->GetEntriesFast());
    1128   Candidate *candidate;
    1129   unsigned int k = 0;
    1130   while((candidate = static_cast<Candidate*>(fItInputGenVtx->Next())))
    1131   {
    1132     gr_genvtx->SetPoint(k, candidate->Position.T()*1E9/c_light, candidate->Position.Z());
    1133     k++;
    1134   }
    1135   gr_genvtx->SetMarkerStyle(33);
    1136   gr_genvtx->SetMarkerColor(6);
    1137   gr_genvtx->SetMarkerSize(2.);
    1138   gr_genvtx->Draw("PE1");
    1139 
    1140   // auto leg = new TLegend(0.1, 0.1);
    1141   // leg->AddEntry(gr_PVtks, "PV tks", "ep");
    1142   // leg->AddEntry(gr_PUtks, "PU tks", "ep");
    1143   // leg->AddEntry(gr_vtx, "Cluster center", "p");
    1144   // leg->Draw();
    1145 
    1146   c_2Dspace->SetGrid();
    1147   c_2Dspace->SaveAs(fFigFolderPath + Form("/c_2Dspace_beta%010.0f-%s%d.png", 1E7*beta, flag, n_it));
    1148 
    1149   delete c_2Dspace;
    1150 }
    1151 
    1152 // -----------------------------------------------------------------------------
    1153 // Plot status at the end
    1154 void VertexFinderDA4D::plot_status_end(vertex_t &vtx, tracks_t &tks)
    1155 {
    1156   unsigned int nv = vtx.getSize();
    1157 
    1158   // Define colors in a meaningfull way
    1159   vector<int> MyPalette(nv);
    1160 
    1161   const int Number = 3;
    1162   double Red[Number]    = { 1.00, 0.00, 0.00};
    1163   double Green[Number]  = { 0.00, 1.00, 0.00};
    1164   double Blue[Number]   = { 1.00, 0.00, 1.00};
    1165   double Length[Number] = { 0.00, 0.50, 1.00 };
    1166   int FI = TColor::CreateGradientColorTable(Number,Length,Red,Green,Blue,nv);
    1167   for (unsigned int i=0;i<nv;i++) MyPalette[i] = FI+i;
    1168 
    1169   TCanvas * c_out = new TCanvas("c_out", "c_out", 800, 600);
    1170   double t_min = TMath::Min( TMath::MinElement(tks.getSize(), &tks.t[0]), TMath::MinElement(vtx.getSize(), &(vtx.t[0]))  ) - 2*fVertexTSize;
    1171   double t_max = TMath::Max(TMath::MaxElement(tks.getSize(), &tks.t[0]), TMath::MaxElement(vtx.getSize(), &(vtx.t[0]))  ) + 2*fVertexTSize;
    1172 
    1173   double z_min = TMath::Min( TMath::MinElement(tks.getSize(), &tks.z[0]), TMath::MinElement(vtx.getSize(), &(vtx.z[0]))  ) - 15;
    1174   double z_max = TMath::Max( TMath::MaxElement(tks.getSize(), &tks.z[0]), TMath::MaxElement(vtx.getSize(), &(vtx.z[0]))  ) + 15;
    1175 
    1176   // Draw tracks
    1177   for(unsigned int i = 0; i < tks.getSize(); i++)
    1178   {
    1179     double dt[] = {1./tks.dt_o[i]};
    1180     double dz[] = {1./tks.dz_o[i]};
    1181     TGraphErrors* gr = new TGraphErrors(1, &(tks.t[i]), &(tks.z[i]), dt, dz);
    1182 
    1183     gr->SetNameTitle(Form("gr%d",i), Form("gr%d",i));
    1184 
    1185     int marker = tks.tt[i]->IsPU? 1 : 4;
    1186     gr->SetMarkerStyle(marker);
    1187 
    1188     int idx = tks.tt[i]->ClusterIndex;
    1189     int color = idx>=0 ? MyPalette[idx] : 13;
    1190     gr->SetMarkerColor(color);
    1191     gr->SetLineColor(color);
    1192 
    1193     int line_style = idx>=0 ? 1 : 3;
    1194     gr->SetLineStyle(line_style);
    1195 
    1196     if(i==0)
    1197     {
    1198       gr->SetTitle(Form("Clustering space - Tot Vertexes = %d", nv));
    1199       gr->GetXaxis()->SetTitle("t CA [ps]");
    1200       gr->GetXaxis()->SetLimits(t_min, t_max);
    1201       gr->GetYaxis()->SetTitle("z CA [mm]");
    1202       gr->GetYaxis()->SetRangeUser(z_min, z_max);
    1203       gr->Draw("APE1");
    1204     }
    1205     else gr->Draw("PE1");
    1206   }
    1207 
    1208   // Draw vertices
    1209   for(unsigned int k = 0; k < vtx.getSize(); k++)
    1210   {
    1211     TGraph* gr = new TGraph(1, &(vtx.t[k]), &(vtx.z[k]));
    1212 
    1213     gr->SetNameTitle(Form("grv%d",k), Form("grv%d",k));
    1214 
    1215     gr->SetMarkerStyle(41);
    1216     gr->SetMarkerSize(2.);
    1217     gr->SetMarkerColor(MyPalette[k]);
    1218 
    1219     gr->Draw("P");
    1220   }
    1221 
    1222   fItInputGenVtx->Reset();
    1223   TGraph* gr_genvtx = new TGraph(fInputGenVtx->GetEntriesFast());
    1224   TGraph* gr_genPV = new TGraph(1);
    1225   Candidate *candidate;
    1226   unsigned int k = 0;
    1227   while((candidate = static_cast<Candidate*>(fItInputGenVtx->Next())))
    1228   {
    1229     if(k == 0 ) {
    1230       gr_genPV->SetPoint(k, candidate->Position.T()*1E9/c_light, candidate->Position.Z());
    1231     }
    1232     else gr_genvtx->SetPoint(k, candidate->Position.T()*1E9/c_light, candidate->Position.Z());
    1233 
    1234     k++;
    1235   }
    1236   gr_genvtx->SetMarkerStyle(20);
    1237   gr_genvtx->SetMarkerColorAlpha(kBlack, 0.8);
    1238   gr_genvtx->SetMarkerSize(.8);
    1239   gr_genvtx->Draw("PE1");
    1240   gr_genPV->SetMarkerStyle(33);
    1241   gr_genPV->SetMarkerColorAlpha(kBlack, 1);
    1242   gr_genPV->SetMarkerSize(2.5);
    1243   gr_genPV->Draw("PE1");
    1244 
    1245   // auto note =  new TLatex();
    1246   // note->DrawLatexNDC(0.5, 0.8, Form("#splitline{Vertexes Reco = %d }{Vertexes gen = %d}", vtx.getSize(), k) );
    1247 
    1248   c_out->SetGrid();
    1249   c_out->SaveAs(fFigFolderPath + Form("/c_final.root"));
    1250   delete c_out;
    1251 }
    1252 
    1253 // -----------------------------------------------------------------------------
    1254 // Plot splitting
    1255 void VertexFinderDA4D::plot_split_crush(double zn, double tn, vertex_t &vtx, tracks_t &tks, int i_vtx)
    1256 {
    1257   vector<double> t, dt, z, dz;
    1258 
    1259   for(unsigned int i = 0; i < tks.getSize(); i++)
    1260   {
    1261       t.push_back(tks.t[i]);
    1262       dt.push_back(1./tks.dt_o[i]);
    1263       z.push_back(tks.z[i]);
    1264       dz.push_back(1./tks.dz_o[i]);
    1265   }
    1266 
    1267 
    1268   double t_min = TMath::Min(TMath::MinElement(t.size(), &t[0]), TMath::MinElement(vtx.getSize(), &(vtx.t[0]))  ) - 50;
    1269   double t_max = TMath::Max(TMath::MaxElement(t.size(), &t[0]), TMath::MaxElement(vtx.getSize(), &(vtx.t[0]))  ) + 50;
    1270 
    1271   double z_min = TMath::Min(TMath::MinElement(z.size(), &z[0]), TMath::MinElement(vtx.getSize(), &(vtx.z[0]))  ) - 5;
    1272   double z_max = TMath::Max(TMath::MaxElement(z.size(), &z[0]), TMath::MaxElement(vtx.getSize(), &(vtx.z[0]))  ) + 5;
    1273 
    1274   auto c_2Dspace = new TCanvas("c_2Dspace", "c_2Dspace", 800, 600);
    1275 
    1276   TGraphErrors* gr_PVtks = new TGraphErrors(t.size(), &t[0], &z[0], &dt[0], &dz[0]);
    1277   gr_PVtks->SetTitle(Form("Clustering space"));
    1278   gr_PVtks->GetXaxis()->SetTitle("t CA [ps]");
    1279   gr_PVtks->GetXaxis()->SetLimits(t_min, t_max);
    1280   gr_PVtks->GetYaxis()->SetTitle("z CA [mm]");
    1281   gr_PVtks->GetYaxis()->SetRangeUser(z_min, z_max);
    1282   gr_PVtks->SetMarkerStyle(4);
    1283   gr_PVtks->SetMarkerColor(1);
    1284   gr_PVtks->SetLineColor(1);
    1285   gr_PVtks->Draw("APE1");
    1286 
    1287   TGraph* gr_vtx = new TGraph(1, &(vtx.t[i_vtx]), &(vtx.z[i_vtx]));
    1288   gr_vtx->SetMarkerStyle(28);
    1289   gr_vtx->SetMarkerColor(2);
    1290   gr_vtx->SetMarkerSize(2.);
    1291   gr_vtx->Draw("PE1");
    1292 
    1293   double t_pos[] = {vtx.t[i_vtx], vtx.t[i_vtx]+100};
    1294   double t_neg[] = {vtx.t[i_vtx], vtx.t[i_vtx]-100};
    1295   double z_pos[] = {vtx.z[i_vtx], vtx.z[i_vtx]+(zn/tn)*100};
    1296   double z_neg[] = {vtx.z[i_vtx], vtx.z[i_vtx]-(zn/tn)*100};
    1297 
    1298   TGraph* gr_pos = new TGraph(2, &t_pos[0], &z_pos[0]);
    1299   gr_pos->SetLineColor(8);
    1300   gr_pos->SetMarkerColor(8);
    1301   gr_pos->Draw("PL");
    1302   TGraph* gr_neg = new TGraph(2, &t_neg[0], &z_neg[0]);
    1303   gr_neg->SetLineColor(4);
    1304   gr_neg->SetMarkerColor(4);
    1305   gr_neg->Draw("PL");
    1306 
    1307 
    1308   c_2Dspace->SetGrid();
    1309   c_2Dspace->SaveAs(fFigFolderPath + Form("/crush_splitting.png"));
    1310 
    1311   delete c_2Dspace;
    1312 }
Note: See TracChangeset for help on using the changeset viewer.