2626#include " iostream"
2727#include " TMatrixD.h"
2828#include " TVectorD.h"
29+ #include " TAxis.h"
2930#include " TMatrixDSymEigen.h"
3031#include " TDatabasePDG.h"
3132#include " TLorentzVector.h"
@@ -55,17 +56,18 @@ bool DelphesO2LutWriter::fatSolve(lutEntry_t& lutEntry,
5556 const float mass,
5657 int itof,
5758 int otof,
58- int q)
59+ int q,
60+ const float nch)
5961{
6062 lutEntry.valid = false ;
6163
6264 static TLorentzVector tlv;
6365 tlv.SetPtEtaPhiM (pt, eta, 0 ., mass);
6466 o2::track::TrackParCov trkIn;
6567 o2::upgrade::convertTLorentzVectorToO2Track (q, tlv, {0 ., 0 ., 0 .}, trkIn);
66-
6768 o2::track::TrackParCov trkOut;
68- if (fat.FastTrack (trkIn, trkOut, 1 ) < 0 ) {
69+ const int status = fat.FastTrack (trkIn, trkOut, nch);
70+ if (status <= 0 ) {
6971 Printf (" --- fatSolve: FastTrack failed --- \n " );
7072 tlv.Print ();
7173 return false ;
@@ -234,6 +236,9 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, in
234236 lutEntry_t lutEntry;
235237
236238 // write entries
239+ int nCalls = 0 ;
240+ int successfullCalls = 0 ;
241+ int failedCalls = 0 ;
237242 for (int inch = 0 ; inch < nnch; ++inch) {
238243 Printf (" --- writing nch = %d/%d" , inch, nnch);
239244 auto nch = lutHeader.nchmap .eval (inch);
@@ -242,6 +247,7 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, in
242247 for (int irad = 0 ; irad < nrad; ++irad) {
243248 Printf (" --- writing irad = %d/%d" , irad, nrad);
244249 for (int ieta = 0 ; ieta < neta; ++ieta) {
250+ nCalls++;
245251 Printf (" --- writing ieta = %d/%d" , ieta, neta);
246252 auto eta = lutHeader.etamap .eval (ieta);
247253 lutEntry.eta = lutHeader.etamap .eval (ieta);
@@ -252,6 +258,7 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, in
252258 if (std::fabs (eta) <= etaMaxBarrel) { // full lever arm ends at etaMaxBarrel
253259 Printf (" Solving in the barrel" );
254260 // printf(" --- fatSolve: pt = %f, eta = %f, mass = %f, field=%f \n", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field);
261+ successfullCalls++;
255262 if (!fatSolve (lutEntry, lutEntry.pt , lutEntry.eta , lutHeader.mass , itof, otof, q)) {
256263 // printf(" --- fatSolve: error \n");
257264 lutEntry.valid = false ;
@@ -260,13 +267,16 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, in
260267 for (int i = 0 ; i < 15 ; ++i) {
261268 lutEntry.covm [i] = 0 .;
262269 }
270+ successfullCalls--;
271+ failedCalls++;
263272 }
264273 } else {
265274 Printf (" Solving outside the barrel" );
266275 // printf(" --- fwdSolve: pt = %f, eta = %f, mass = %f, field=%f \n", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field);
267276 lutEntry.eff = 1 .;
268277 lutEntry.eff2 = 1 .;
269278 bool retval = true ;
279+ successfullCalls++;
270280 if (useFlatDipole) { // Using the parametrization at the border of the barrel
271281 retval = fatSolve (lutEntry, lutEntry.pt , etaMaxBarrel, lutHeader.mass , itof, otof, q);
272282 } else if (usePara) {
@@ -288,6 +298,8 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, in
288298 for (int i = 0 ; i < 15 ; ++i) {
289299 lutEntry.covm [i] = 0 .;
290300 }
301+ successfullCalls--;
302+ failedCalls++;
291303 }
292304 }
293305 Printf (" Diagonalizing" );
@@ -298,6 +310,8 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, in
298310 }
299311 }
300312 }
313+ Printf (" --- finished writing LUT file %s" , filename);
314+ Printf (" --- successfull calls: %d/%d, failed calls: %d/%d" , successfullCalls, nCalls, failedCalls, nCalls);
301315 lutFile.close ();
302316}
303317
@@ -331,10 +345,13 @@ void DelphesO2LutWriter::diagonalise(lutEntry_t& lutEntry)
331345
332346TGraph* DelphesO2LutWriter::lutRead (const char * filename, int pdg, int what, int vs, float nch, float radius, float eta, float pt)
333347{
348+ Printf (" --- reading LUT file %s" , filename);
349+ // vs
334350 static const int kNch = 0 ;
335351 static const int kEta = 1 ;
336352 static const int kPt = 2 ;
337353
354+ // what
338355 static const int kEfficiency = 0 ;
339356 static const int kEfficiency2 = 1 ;
340357 static const int kEfficiencyInnerTOF = 2 ;
@@ -360,6 +377,58 @@ TGraph* DelphesO2LutWriter::lutRead(const char* filename, int pdg, int what, int
360377 }
361378 auto nbins = lutMap.nbins ;
362379 auto g = new TGraph ();
380+ g->SetName (Form (" lut_%s_%d_vs_%d_what_%d" , filename, pdg, vs, what));
381+ g->SetTitle (Form (" LUT for %s, pdg %d, vs %d, what %d" , filename, pdg, vs, what));
382+ switch (vs) {
383+ case kNch :
384+ Printf (" --- vs = kNch" );
385+ g->GetXaxis ()->SetTitle (" Nch" );
386+ break ;
387+ case kEta :
388+ Printf (" --- vs = kEta" );
389+ g->GetXaxis ()->SetTitle (" #eta" );
390+ break ;
391+ case kPt :
392+ Printf (" --- vs = kPt" );
393+ g->GetXaxis ()->SetTitle (" p_{T} (GeV/c)" );
394+ break ;
395+ default :
396+ Printf (" --- error: unknown vs %d" , vs);
397+ return nullptr ;
398+ }
399+ switch (what) {
400+ case kEfficiency :
401+ Printf (" --- what = kEfficiency" );
402+ g->GetYaxis ()->SetTitle (" Efficiency (%)" );
403+ break ;
404+ case kEfficiency2 :
405+ Printf (" --- what = kEfficiency2" );
406+ g->GetYaxis ()->SetTitle (" Efficiency2 (%)" );
407+ break ;
408+ case kEfficiencyInnerTOF :
409+ Printf (" --- what = kEfficiencyInnerTOF" );
410+ g->GetYaxis ()->SetTitle (" Inner TOF Efficiency (%)" );
411+ break ;
412+ case kEfficiencyOuterTOF :
413+ Printf (" --- what = kEfficiencyOuterTOF" );
414+ g->GetYaxis ()->SetTitle (" Outer TOF Efficiency (%)" );
415+ break ;
416+ case kPtResolution :
417+ Printf (" --- what = kPtResolution" );
418+ g->GetYaxis ()->SetTitle (" p_{T} Resolution (%)" );
419+ break ;
420+ case kRPhiResolution :
421+ Printf (" --- what = kRPhiResolution" );
422+ g->GetYaxis ()->SetTitle (" R#phi Resolution (#mum)" );
423+ break ;
424+ case kZResolution :
425+ Printf (" --- what = kZResolution" );
426+ g->GetYaxis ()->SetTitle (" Z Resolution (#mum)" );
427+ break ;
428+ default :
429+ Printf (" --- error: unknown what %d" , what);
430+ return nullptr ;
431+ }
363432
364433 bool canBeInvalid = true ;
365434 for (int i = 0 ; i < nbins; ++i) {
0 commit comments