/* *==================================================================== * * CBM Level 1 Reconstruction * * Authors: I.Kisel, S.Gorbunov * * e-mail : ikisel@kip.uni-heidelberg.de * *==================================================================== * * L1 Fit performance * *==================================================================== */ #include "CbmL1.h" #include "CbmKF.h" #include "CbmKFMath.h" #include "iostream.h" #define MAGFIELD static TList *listHisto; static int L1_NMTRA_D0 = 0; static int L1_NMTRA_TOTAL = 0; static int L1_NMTRA_FAST = 0; static int L1_NMTRA_FAST_PRIM = 0; static int L1_NMTRA_FAST_SEC = 0; static int L1_NMTRA_SLOW = 0; static int L1_NMTRA_SLOW_PRIM = 0; static int L1_NMTRA_SLOW_SEC = 0; static int L1_NMTRA_D0_RECO = 0; static int L1_NMTRA_TOTAL_RECO = 0; static int L1_NMTRA_FAST_RECO = 0; static int L1_NMTRA_FAST_PRIM_RECO = 0; static int L1_NMTRA_FAST_SEC_RECO = 0; static int L1_NMTRA_SLOW_RECO = 0; static int L1_NMTRA_SLOW_PRIM_RECO = 0; static int L1_NMTRA_SLOW_SEC_RECO = 0; static int L1_NTRACKS = 0; static int L1_NCLONES = 0; static int L1_NGHOSTS = 0; static int L1_NEVENTS = 0; static map RTRAtoMTRAmap; static int RTRAtoMTRAmap_clear = 1; void FitPerformanceFillPulls( const double MC[5], const double T[5],const double C[15], TH1F* hist[10] ) { hist[0]->Fill( T[0] - MC[0] ); hist[1]->Fill( T[1] - MC[1] ); hist[2]->Fill( T[2] - MC[2] ); hist[3]->Fill( T[3] - MC[3] ); if ( fabs( T[4] )>1.e-10 ) hist[4]->Fill( T[4]/MC[4]/*/T[4]*/ - 1. ); if ( C[ 0] >1.e-20 ) hist[5]->Fill( ( T[0] - MC[0] )/sqrt(C[ 0]) ); if ( C[ 2] >1.e-20 ) hist[6]->Fill( ( T[1] - MC[1] )/sqrt(C[ 2]) ); if ( C[ 5] >1.e-20 ) hist[7]->Fill( ( T[2] - MC[2] )/sqrt(C[ 5]) ); if ( C[ 9] >1.e-20 ) hist[8]->Fill( ( T[3] - MC[3] )/sqrt(C[ 9]) ); if ( C[14] >1.e-20 ) hist[9]->Fill( ( T[4] - MC[4] )/sqrt(C[14]) ); } void CbmL1::Performance() { CbmKF &KF = *CbmKF::Instance(); // histogramms static TProfile *p_eff_all_vs_mom, *p_eff_prim_vs_mom, *p_eff_sec_vs_mom, *p_eff_d0_vs_mom; static TH1F *h_mom_track; static TH1F *h_ghost_mom, *h_ghost_nhits, *h_ghost_station, *h_ghost_chi2, *h_ghost_prob, *h_ghost_tx, *h_ghost_ty; static TH1F *h_reco_mom, *h_reco_d0_mom, *h_reco_nhits, *h_reco_station, *h_reco_chi2, *h_reco_prob, *h_reco_clean, *h_reco_time; static TH1F *h_tx, *h_ty, *h_sec_r, *h_ghost_r; static TH1F *h_notfound_mom, *h_notfound_nhits, *h_notfound_station, *h_notfound_r, *h_notfound_tx, *h_notfound_ty; static TH1F *MC_P, *nhits, *h_chi2, *h_prob, *frst[10], *last[10], *extr[10], *fvtx[10], *vtx[9], *MC_vx, *MC_vy, *MC_vz, *VtxChiPrim, *VtxChiSec; static TH2F *P_vs_P ; static bool first_call = 1; if ( first_call ) { first_call = 0; //TDirectory *maindir = gDirectory; //histodir->cd(); TDirectory *curdir = gDirectory; TDirectory *histodir = gROOT->mkdir("histodir"); histodir->cd(); p_eff_all_vs_mom = new TProfile("p_eff_all_vs_mom", "AllSet Efficiency vs Momentum", 100, 0.0, 5.0, 0.0, 100.0); p_eff_prim_vs_mom = new TProfile("p_eff_prim_vs_mom", "All Primary Set Efficiency vs Momentum", 100, 0.0, 5.0, 0.0, 100.0); p_eff_sec_vs_mom = new TProfile("p_eff_sec_vs_mom", "All Secondary Set Efficiency vs Momentum", 100, 0.0, 5.0, 0.0, 100.0); p_eff_d0_vs_mom = new TProfile("p_eff_d0_vs_mom", "D0 Secondary Tracks Efficiency vs Momentum", 150, 0.0, 15.0, 0.0, 100.0); h_mom_track = new TH1F("h_mom_track", "Momentum of accepted tracks", 100, 0.0, 20.0); h_ghost_mom = new TH1F("h_ghost_mom", "Momentum of ghost track", 50, 0.0, 5.0); h_ghost_nhits = new TH1F("h_ghost_nhits", "Number of hits in ghost track", 50, 0.0, 10.0); h_ghost_station = new TH1F("h_ghost_station", "First station of ghost track", 50, -0.5, 10.0); h_ghost_chi2 = new TH1F("h_ghost_chi2", "Chi2/NDF of ghost track", 50, -0.5, 10.0); h_ghost_prob = new TH1F("h_ghost_prob", "Prob of ghost track", 50, 0., 1.01); h_ghost_r = new TH1F("h_ghost_r", "R of ghost track at the first hit", 50, 0.0, 15.0); h_ghost_tx = new TH1F("h_ghost_tx", "tx of ghost track at the first hit", 50, -5.0, 5.0); h_ghost_ty = new TH1F("h_ghost_ty", "ty of ghost track at the first hit", 50, -1.0, 1.0); h_reco_mom = new TH1F("h_reco_mom", "Momentum of reco track", 50, 0.0, 5.0); h_reco_nhits = new TH1F("h_reco_nhits", "Number of hits in reco track", 50, 0.0, 10.0); h_reco_station = new TH1F("h_reco_station", "First station of reco track", 50, -0.5, 10.0); h_reco_chi2 = new TH1F("h_reco_chi2", "Chi2/NDF of reco track", 50, -0.5, 10.0); h_reco_prob = new TH1F("h_reco_prob", "Prob of reco track", 50, 0., 1.01); h_reco_clean = new TH1F("h_reco_clean", "Percentage of correct hits", 100, -0.5, 100.5); h_reco_time = new TH1F("h_reco_time", "CA Track Finder Time (s/ev)", 20, 0.0, 20.0); h_reco_d0_mom = new TH1F("h_reco_d0_mom", "Momentum of reco D0 track", 150, 0.0, 15.0); h_tx = new TH1F("h_tx", "tx of MC track at the vertex", 50, -0.5, 0.5); h_ty = new TH1F("h_ty", "ty of MC track at the vertex", 50, -0.5, 0.5); h_sec_r = new TH1F("h_sec_r", "R of sec MC track at the first hit", 50, 0.0, 15.0); h_notfound_mom = new TH1F("h_notfound_mom", "Momentum of not found track", 50, 0.0, 5.0); h_notfound_nhits = new TH1F("h_notfound_nhits", "Num hits of not found track", 50, 0.0, 10.0); h_notfound_station = new TH1F("h_notfound_station", "First station of not found track", 50, 0.0, 10.0); h_notfound_r = new TH1F("h_notfound_r", "R at first hit of not found track", 50, 0.0, 15.0); h_notfound_tx = new TH1F("h_notfound_tx", "tx of not found track at the first hit", 50, -5.0, 5.0); h_notfound_ty = new TH1F("h_notfound_ty", "ty of not found track at the first hit", 50, -1.0, 1.0); h_chi2 = new TH1F("chi2", "Chi^2", 100, 0.0, 10.); h_prob = new TH1F("prob", "Prob", 100, 0.0, 1.01); VtxChiPrim = new TH1F("vtxChiPrim", "Chi to primary vtx for primary tracks", 100, 0.0, 10.); VtxChiSec = new TH1F("vtxChiSec", "Chi to primary vtx for secondary tracks", 100, 0.0, 10.); MC_P = new TH1F("MC_P", "MC P/Q ", 1000, 0., 10.); MC_vx = new TH1F("MC_vx", "MC Vertex X", 100, -.05, .05); MC_vy = new TH1F("MC_vy", "MC Vertex Y", 100, -.05, .05); MC_vz = new TH1F("MC_vz", "MC Vertex Z", 100, -.05, .05); nhits = new TH1F("nhits", "N STS hits on MC track", 20, 0., 20.); P_vs_P = new TH2F("P_vs_P", "Resolution P/Q vs P", 20, 0., 20.,100, -.05, .05); frst[0]= new TH1F("frst_x", "Residual X (cm) at First STS Hit", 100, -.01, .01); frst[1]= new TH1F("frst_y", "Residual Y (cm) at First STS Hit", 100, -.01, .01); frst[2]= new TH1F("frst_tx","Residual Tx at First STS Hit", 100, -.005, .005); frst[3]= new TH1F("frst_ty","Residual Ty at First STS Hit", 100, -.005, .005); frst[4]= new TH1F("frst_P", "Resolution P/Q at First STS Hit", 100, -.1, .1); frst[5]= new TH1F("frst_pull_x", "Pull X at First STS Hit", 100, -10., 10.); frst[6]= new TH1F("frst_pull_y", "Pull Y at First STS Hit", 100, -10., 10.); frst[7]= new TH1F("frst_pull_tx", "Pull Tx at First STS Hit", 100, -10., 10.); frst[8]= new TH1F("frst_pull_ty", "Pull Ty at First STS Hit", 100, -10., 10.); frst[9]= new TH1F("frst_pull_qp", "Pull Q/P at First STS Hit", 100, -10., 10.); last[0]= new TH1F("last_x", "Residual X (cm) at Last STS Hit", 100, -.01, .01); last[1]= new TH1F("last_y", "Residual Y (cm) at Last STS Hit", 100, -.01, .01); last[2]= new TH1F("last_tx","Residual Tx at Last STS Hit", 100, -.005, .005); last[3]= new TH1F("last_ty","Residual Ty at Last STS Hit", 100, -.005, .005); last[4]= new TH1F("last_P", "Resolution P/Q at Last STS Hit", 100, -.1, .1); last[5]= new TH1F("last_pull_x", "Pull X at Last STS Hit", 100, -10., 10.); last[6]= new TH1F("last_pull_y", "Pull Y at Last STS Hit", 100, -10., 10.); last[7]= new TH1F("last_pull_tx", "Pull Tx at Last STS Hit", 100, -10., 10.); last[8]= new TH1F("last_pull_ty", "Pull Ty at Last STS Hit", 100, -10., 10.); last[9]= new TH1F("last_pull_qp", "Pull Q/P at Last STS Hit", 100, -10., 10.); extr[0] = new TH1F("extr_x", "Residual X (cm) at Vertex ", 100, -.01, .01); extr[1] = new TH1F("extr_y", "Residual Y (cm) at Vertex ", 100, -.01, .01); extr[2] = new TH1F("extr_tx", "Residual Tx at Vertex ", 100, -.005, .005); extr[3] = new TH1F("extr_ty", "Residual Ty at Vertex ", 100, -.005, .005); extr[4] = new TH1F("extr_P", "Resolution P/Q at Vertex ", 100, -.1, .1); extr[5] = new TH1F("extr_pull_x", "Pull X at Vertex ", 100, -10., 10.); extr[6] = new TH1F("extr_pull_y", "Pull Y at Vertex ", 100, -10., 10.); extr[7] = new TH1F("extr_pull_tx", "Pull Tx at Vertex ", 100, -10., 10.); extr[8] = new TH1F("extr_pull_ty", "Pull Ty at Vertex ", 100, -10., 10.); extr[9] = new TH1F("extr_pull_qp", "Pull Q/P at Vertex ", 100, -10., 10.); fvtx[0] = new TH1F("fvtx_x", "Residual X (cm) with Vertex ", 100, -.01, .01); fvtx[1] = new TH1F("fvtx_y", "Residual Y (cm) with Vertex ", 100, -.01, .01); fvtx[2] = new TH1F("fvtx_tx", "Residual Tx with Vertex ", 100, -.005, .005); fvtx[3] = new TH1F("fvtx_ty", "Residual Ty with Vertex ", 100, -.005, .005); fvtx[4] = new TH1F("fvtx_P", "Resolution P/Q with Vertex ", 100, -.1, .1); fvtx[5] = new TH1F("fvtx_pull_x", "Pull X with Vertex ", 100, -10., 10.); fvtx[6] = new TH1F("fvtx_pull_y", "Pull Y with Vertex ", 100, -10., 10.); fvtx[7] = new TH1F("fvtx_pull_tx", "Pull Tx with Vertex ", 100, -10., 10.); fvtx[8] = new TH1F("fvtx_pull_ty", "Pull Ty with Vertex ", 100, -10., 10.); fvtx[9] = new TH1F("fvtx_pull_qp", "Pull Q/P with Vertex ", 100, -10., 10.); vtx[0] = new TH1F("vtx_x", "Residual X (cm) for Primary Vertex ",100,-.0005,.0005); vtx[1] = new TH1F("vtx_y", "Residual Y (cm) for Primary Vertex ",100,-.0005,.0005); vtx[2] = new TH1F("vtx_z", "Residual Z (cm) for Primary Vertex ",100,-.004,.004); vtx[3] = new TH1F("vtx_pull_x", "Pull X for Primary Vertex ", 100, -10., 10.); vtx[4] = new TH1F("vtx_pull_y", "Pull Y for Primary Vertex ", 100, -10., 10.); vtx[5] = new TH1F("vtx_pull_z", "Pull Z for Primary Vertex ", 100, -10., 10.); vtx[6] = new TH1F("vtx_chi2", "Chi2 for Primary Vertex ", 100, 0.0, 10.); vtx[7] = new TH1F("vtx_ntracks", "N tracks at Primary Vertex ", 1000,0.0,1000.); vtx[8] = new TH1F("vtx_prob", "Prob for Primary Vertex ", 100, 0.0, 1.01); #ifdef MAGFIELD if ( KF.GetMagneticField() ) { //TDirectory *mfdir = histodir->mkdir("Field"); //histodir->cd("Field"); const double Xmin = -90., Xmax = 90., dx = 5.; const double Ymin = -50., Ymax = 50., dy = 5.; const double Zmin = -10., Zmax = 170., dz = 5., dz_fine=0.1; int BinsX = int( (Xmax-Xmin)/dx ); int BinsY = int( (Ymax-Ymin)/dy ); int BinsZ = int( (Zmax-Zmin)/dz ); int BinsZ_fine = int( (Zmax-Zmin)/dz_fine ); if(1) { TProfile *fldX = new TProfile("field_X","X-component of Magnetic Field vs Z", BinsZ_fine, Zmin, Zmax ), *fldY = new TProfile("field_Y","Y-component of Magnetic Field vs Z", BinsZ_fine, Zmin, Zmax ), *fldZ = new TProfile("field_Z","Z-component of Magnetic Field vs Z", BinsZ_fine, Zmin, Zmax ); for( double z=Zmin; z<=Zmax; z+=dz_fine ) { Double_t r[3] = { 0., 0., z}; Double_t B[3] = {10.,10.,10.}; KF.GetMagneticField()->GetFieldValue(r, B); fldX->Fill(z,B[0]); fldY->Fill(z,B[1]); fldZ->Fill(z,B[2]); } } if(0) { TProfile2D *Bx_xy = new TProfile2D("Bx_xy","X-component of Magnetic Field in XY", BinsX, Xmin, Xmax, BinsY, Ymin, Ymax ), *By_xy = new TProfile2D("By_xy","Y-component of Magnetic Field in XY", BinsX, Xmin, Xmax, BinsY, Ymin, Ymax ), *Bz_xy = new TProfile2D("Bz_xy","Z-component of Magnetic Field in XY", BinsX, Xmin, Xmax, BinsY, Ymin, Ymax ), *Bx_zr = new TProfile2D("Bx_zr","X-component of Magnetic Field in ZR", BinsZ, Zmin, Zmax, BinsX/2, 0, Xmax ), *By_zr = new TProfile2D("By_zr","Y-component of Magnetic Field in ZR", BinsZ, Zmin, Zmax, BinsX/2, 0, Xmax ), *Bz_zr = new TProfile2D("Bz_zr","Z-component of Magnetic Field in ZR", BinsZ, Zmin, Zmax, BinsX/2, 0, Xmax ); for( double x=Xmin; x<=Xmax; x+=dx ) for( double y=Ymin; y<=Ymax; y+=dy ) for( double z=Zmin; z<=Zmax; z+=dz ) { double r = sqrt(x*x+y*y); Double_t p[3] = { x, y, z}; Double_t B[3] = {10.,10.,10.}; KF.GetMagneticField()->GetFieldValue(p, B); Bx_xy->Fill(x,y,fabs(B[0])); By_xy->Fill(x,y,fabs(B[1])); Bz_xy->Fill(x,y,fabs(B[2])); Bx_zr->Fill(z,r,fabs(B[0])); By_zr->Fill(z,r,fabs(B[1])); Bz_zr->Fill(z,r,fabs(B[2])); } } if(0) { int j=0; for( vector::iterator i=KF.vStsMaterial.begin(); i!=KF.vStsMaterial.end();i++) { int Bins = 50; char n1[225], n2[225]; sprintf(n1,"B%ix",j); sprintf(n2,"X-component of Magnetic Field at z=%f",i->z); cout<R, i->R, Bins, -i->R, i->R ); sprintf(n1,"B%iy",j); sprintf(n2,"Y-component of Magnetic Field at z=%f",i->z); cout<R, i->R, Bins, -i->R, i->R ); sprintf(n1,"B%iz",j); sprintf(n2,"Z-component of Magnetic Field at z=%f",i->z); cout<R, i->R, Bins, -i->R, i->R ); double d = 2*i->R/Bins/2.; for( double x=-i->R; x<=i->R; x+=d ) for( double y=-i->R; y<=i->R; y+=d ) { double r = sqrt(x*x+y*y); if( rr || r>i->R ) continue; Double_t p[3] = { x, y, i->z}; Double_t B[3] = {10.,10.,10.}; KF.GetMagneticField()->GetFieldValue(p, B); Bx->Fill(x,y,fabs(B[0])); By->Fill(x,y,fabs(B[1])); Bz->Fill(x,y,fabs(B[2])); } j++; } } //histodir->cd(); } #endif//MAGFIELD //maindir->cd(); // ----- Create list of all histogram pointers listHisto = gDirectory->GetList(); curdir->cd(); /* { Double_t r[3] = { 0., 1., 5.0}; Double_t B[3] = {10.,10.,10.}; KF.GetMagneticField()->GetFieldValue(r, B); cout << "Magnetic filed at (" << r[0] << ", " << r[1] << ", " << r[2] << ") -> B(" << B[0] << ", " << B[1] << ", " << B[2] << ")" << endl; } { Double_t r[3] = { 0., -1., 5.0}; Double_t B[3] = {10.,10.,10.}; KF.GetMagneticField()->GetFieldValue(r, B); cout << "Magnetic filed at (" << r[0] << ", " << r[1] << ", " << r[2] << ") -> B(" << B[0] << ", " << B[1] << ", " << B[2] << ")" << endl; } */ }// first_call //track finding part typedef map Hitmap; // typedef map::iterator HitmapIt; Hitmap hitmap; map recostamap, mcstamap; map pmtramap; // reco track + clones map pmctrackmap; pmctrackmap.clear(); h_reco_time->Fill(CATime); //track finding performance pmtramap.clear(); mtrartramap.clear(); rtramtramap.clear(); if (RTRAtoMTRAmap_clear){ RTRAtoMTRAmap.clear(); RTRAtoMTRAmap_clear = 0; } int ntracks = 0; int nghosts = 0; //cout<<"reco.."<::iterator i = vMCTracks.begin(); i != vMCTracks.end(); ++i) { CbmL1MCTrack &MC = *i; if (pmctrackmap.find(MC.UID) == pmctrackmap.end()){ pmctrackmap.insert(pair(MC.UID, &MC)); } else{ cout << "*** L1: Track ID " << MC.UID << " encountered twice! ***" << endl; } }} //reconstructed tracks for (vector::iterator rtraIt = vRTracks.begin(); rtraIt != vRTracks.end(); ++rtraIt){ ntracks++; int hitsum = 0; hitmap.clear(); for (vector::iterator phitIt = (rtraIt->StsHits).begin(); phitIt != (rtraIt->StsHits).end(); ++phitIt){ int UID = -1; if ((*phitIt)->MC_Point >= 0) UID = vMCPoints[(*phitIt)->MC_Point].UID; if(hitmap.find(UID) == hitmap.end()) hitmap.insert(pair(UID, 1)); else{ hitmap[UID] += 1; } hitsum++; } { CbmL1Track* prtra = &(*rtraIt); h_reco_mom->Fill(fabs(1.0/prtra->T[4])); h_reco_nhits->Fill((prtra->StsHits).size()); CbmL1StsHit *pmh = *((prtra->StsHits).begin()); h_reco_station->Fill(pmh->iStation); if (prtra->NDF > 0) h_reco_chi2->Fill(prtra->chi2/prtra->NDF); h_reco_prob->Fill(TMath::Erf(sqrt(prtra->chi2/prtra->NDF))); } // RTRA <-> MTRA identification int ghost = 1; double momentum = 0.0; double max_percent = 0.0; for( HitmapIt posIt = hitmap.begin(); posIt != hitmap.end(); posIt++ ){ // for Hitmap if (posIt->first < 0) continue; //not a MC track - based on fake hits if (100.0*double(posIt->second)/double(hitsum) > max_percent) max_percent = 100.0*double(posIt->second)/double(hitsum); if ( double(posIt->second)/double(hitsum) > 0.7 ){ //MTRA reconstructed if (pmctrackmap.find(posIt->first) == pmctrackmap.end()) continue; CbmL1MCTrack* pmtra = pmctrackmap[posIt->first]; //pointer to MC track CbmL1Track* prtra = &(*rtraIt); mtrartramap.insert(pair(pmtra, prtra)); rtramtramap.insert(pair(prtra, pmtra)); momentum = pmtra->p; if(pmtramap.find(pmtra) == pmtramap.end()) pmtramap.insert(pair(pmtra, 1)); else{ pmtramap[pmtra]++; } RTRAtoMTRAmap.insert(pair(prtra, pmtra)); ghost = 0; } } h_reco_clean->Fill(max_percent); nghosts += ghost; if (ghost == 1){ //ghost CbmL1Track* prtra = &(*rtraIt); h_ghost_mom->Fill(fabs(1.0/prtra->T[4])); h_ghost_nhits->Fill((prtra->StsHits).size()); CbmL1StsHit *pmh1 = *((prtra->StsHits).begin()); h_ghost_station->Fill(pmh1->iStation); if (prtra->NDF > 0){ h_ghost_chi2->Fill(prtra->chi2/prtra->NDF); h_ghost_prob->Fill(TMath::Erf(sqrt(prtra->chi2/prtra->NDF))); } h_ghost_r->Fill(sqrt(pmh1->x()*pmh1->x()+pmh1->y()*pmh1->y())); CbmL1StsHit *pmh2 = *(++((prtra->StsHits).begin())); h_ghost_tx->Fill((pmh2->x()-pmh1->x())/(pmh2->z()-pmh1->z())); h_ghost_ty->Fill((pmh2->y()-pmh1->y())/(pmh2->z()-pmh1->z())); } }//reco tracks int nmtra_d0 = 0; int nmtra_total = 0; int nmtra_fast = 0; int nmtra_fast_prim = 0; int nmtra_fast_sec = 0; int nmtra_slow = 0; int nmtra_slow_prim = 0; int nmtra_slow_sec = 0; int nmtra_d0_reco = 0; int nmtra_total_reco = 0; int nmtra_fast_reco = 0; int nmtra_fast_prim_reco = 0; int nmtra_fast_sec_reco = 0; int nmtra_slow_reco = 0; int nmtra_slow_prim_reco = 0; int nmtra_slow_sec_reco = 0; int nclones = 0; for ( map::iterator mtraIt = pmctrackmap.begin(); mtraIt != pmctrackmap.end(); mtraIt++ ) { // No Sts hits? int nmchits = (mtraIt->second)->StsHits.size(); if (nmchits == 0) continue; double momentum = (mtraIt->second)->p; if (nmchits < 4) continue;// allset 3 or 4 or 5 stations ? int sta_count = 1; vector::iterator previous_hit = (mtraIt->second)->StsHits.begin(); vector::iterator next_hit = previous_hit; next_hit++; /* {for (vector::reverse_iterator imh = previous_hit; imh != (mtraIt->second)->StsHits.rend(); ++imh){ cout << " " << (*imh)->iStation; }} */ for (vector::iterator imh = next_hit; imh != (mtraIt->second)->StsHits.end(); ++imh){ if ((*imh)->iStation - (*previous_hit)->iStation == 1){ sta_count++; } else{ sta_count = 1; } if (sta_count == 4) break; previous_hit = imh; } //cout << " => " << sta_count << endl; if (sta_count < 4) continue;//less than 4 continuous hits //continuous stations CbmL1StsHit* ph = *((mtraIt->second)->StsHits.begin()); h_sec_r->Fill(sqrt(ph->x()*ph->x()+ph->y()*ph->y())); //if (momentum < 0.2) continue;// reject very slow tracks from analysis h_tx->Fill((mtraIt->second)->px/(mtraIt->second)->pz); h_ty->Fill((mtraIt->second)->py/(mtraIt->second)->pz); h_mom_track->Fill(momentum); int reco = 0; if (pmtramap.find(mtraIt->second) != pmtramap.end()){ reco = 1; nclones += pmtramap[mtraIt->second]-1; } if(!reco){ h_notfound_mom->Fill(momentum); h_notfound_nhits->Fill(nmchits); h_notfound_station->Fill(ph->iStation); h_notfound_r->Fill(sqrt(ph->x()*ph->x()+ph->y()*ph->y())); CbmL1StsHit* ph21 = *((mtraIt->second)->StsHits.begin()); CbmL1StsHit* ph22 = *(++((mtraIt->second)->StsHits.begin())); h_notfound_tx->Fill((ph22->x()-ph21->x())/(ph22->z()-ph21->z())); h_notfound_ty->Fill((ph22->y()-ph21->y())/(ph22->z()-ph21->z())); } if (((mtraIt->second)->mother_ID < 0)&&((mtraIt->second)->z > 0)){ // D0 nmtra_d0++; if (reco) nmtra_d0_reco++; if (reco) p_eff_d0_vs_mom->Fill(momentum, 100.0); else p_eff_d0_vs_mom->Fill(momentum, 0.0); h_reco_d0_mom->Fill(momentum); } nmtra_total++; if (reco) nmtra_total_reco++; if (reco){ p_eff_all_vs_mom->Fill(momentum, 100.0); if ((mtraIt->second)->mother_ID < 0){ // primary p_eff_prim_vs_mom->Fill(momentum, 100.0); }else{ p_eff_sec_vs_mom->Fill(momentum, 100.0); } }else{ p_eff_all_vs_mom->Fill(momentum, 0.0); if ((mtraIt->second)->mother_ID < 0){ // primary p_eff_prim_vs_mom->Fill(momentum, 0.0); }else{ p_eff_sec_vs_mom->Fill(momentum, 0.0); } } if (momentum>1.0){ // reference tracks nmtra_fast++; if (reco) nmtra_fast_reco++; if ((mtraIt->second)->mother_ID < 0){ // primary nmtra_fast_prim++; if (reco) nmtra_fast_prim_reco++; } else{ nmtra_fast_sec++; if (reco) nmtra_fast_sec_reco++; } } else{ // extra set of tracks nmtra_slow++; if (reco){ nmtra_slow_reco++; } if ((mtraIt->second)->mother_ID < 0){ // primary nmtra_slow_prim++; if (reco) nmtra_slow_prim_reco++; } else{ nmtra_slow_sec++; if (reco) nmtra_slow_sec_reco++; } } } L1_NMTRA_D0 += nmtra_d0; L1_NMTRA_TOTAL += nmtra_total; L1_NMTRA_FAST += nmtra_fast; L1_NMTRA_FAST_PRIM += nmtra_fast_prim; L1_NMTRA_FAST_SEC += nmtra_fast_sec; L1_NMTRA_SLOW += nmtra_slow; L1_NMTRA_SLOW_PRIM += nmtra_slow_prim; L1_NMTRA_SLOW_SEC += nmtra_slow_sec; L1_NMTRA_D0_RECO += nmtra_d0_reco; L1_NMTRA_TOTAL_RECO += nmtra_total_reco; L1_NMTRA_FAST_RECO += nmtra_fast_reco; L1_NMTRA_FAST_PRIM_RECO += nmtra_fast_prim_reco; L1_NMTRA_FAST_SEC_RECO += nmtra_fast_sec_reco; L1_NMTRA_SLOW_RECO += nmtra_slow_reco; L1_NMTRA_SLOW_PRIM_RECO += nmtra_slow_prim_reco; L1_NMTRA_SLOW_SEC_RECO += nmtra_slow_sec_reco; L1_NCLONES += nclones; L1_NGHOSTS += nghosts; L1_NTRACKS += ntracks; L1_NEVENTS++; double p_total=0.0,p_fast=0.0,p_fast_prim=0.0,p_fast_sec=0.0,p_slow=0.0,p_slow_prim=0.0,p_slow_sec=0.0; double p_clone=0.0,p_ghost=0.0; if (nmtra_total!=0) p_total = double(nmtra_total_reco)/double(nmtra_total); if (nmtra_fast!=0) p_fast = double(nmtra_fast_reco)/double(nmtra_fast); if (nmtra_fast_prim!=0) p_fast_prim = double(nmtra_fast_prim_reco)/double(nmtra_fast_prim); if (nmtra_fast_sec!=0) p_fast_sec = double(nmtra_fast_sec_reco)/double(nmtra_fast_sec); if (nmtra_slow!=0) p_slow = double(nmtra_slow_reco)/double(nmtra_slow); if (nmtra_slow_prim!=0) p_slow_prim = double(nmtra_slow_prim_reco)/double(nmtra_slow_prim); if (nmtra_slow_sec!=0) p_slow_sec = double(nmtra_slow_sec_reco)/double(nmtra_slow_sec); if (ntracks!=0) p_clone = double(nclones)/double(ntracks); if (ntracks!=0) p_ghost = double(nghosts)/double(ntracks); double P_D0=0.0,P_TOTAL=0.0,P_FAST=0.0,P_FAST_PRIM=0.0,P_FAST_SEC=0.0,P_SLOW=0.0,P_SLOW_PRIM=0.0,P_SLOW_SEC=0.0; double P_CLONE=0.0,P_GHOST=0.0; if (L1_NMTRA_D0!=0) P_D0 = double(L1_NMTRA_D0_RECO)/double(L1_NMTRA_D0); if (L1_NMTRA_TOTAL!=0) P_TOTAL = double(L1_NMTRA_TOTAL_RECO)/double(L1_NMTRA_TOTAL); if (L1_NMTRA_FAST!=0) P_FAST = double(L1_NMTRA_FAST_RECO)/double(L1_NMTRA_FAST); if (L1_NMTRA_FAST_PRIM!=0) P_FAST_PRIM = double(L1_NMTRA_FAST_PRIM_RECO)/double(L1_NMTRA_FAST_PRIM); if (L1_NMTRA_FAST_SEC!=0) P_FAST_SEC = double(L1_NMTRA_FAST_SEC_RECO)/double(L1_NMTRA_FAST_SEC); if (L1_NMTRA_SLOW!=0) P_SLOW = double(L1_NMTRA_SLOW_RECO)/double(L1_NMTRA_SLOW); if (L1_NMTRA_SLOW_PRIM!=0) P_SLOW_PRIM = double(L1_NMTRA_SLOW_PRIM_RECO)/double(L1_NMTRA_SLOW_PRIM); if (L1_NMTRA_SLOW_SEC!=0) P_SLOW_SEC = double(L1_NMTRA_SLOW_SEC_RECO)/double(L1_NMTRA_SLOW_SEC); if (L1_NMTRA_TOTAL_RECO!=0) //P_CLONE = double(L1_NCLONES)/double(L1_NTRACKS); P_CLONE = double(L1_NCLONES)/double(L1_NMTRA_TOTAL); if (L1_NMTRA_TOTAL_RECO!=0) //P_GHOST = double(L1_NGHOSTS)/double(L1_NTRACKS); P_GHOST = double(L1_NGHOSTS)/double(L1_NMTRA_TOTAL_RECO); cout.setf(ios::fixed); cout.setf(ios::showpoint); cout.precision(4); #ifdef XXX cout << endl; cout << "L1 PER EVENT STAT : " << endl; cout << "MC Refset : " << nmtra_fast << endl; cout << "MC Extras : " << nmtra_slow << endl; cout << "ALL SIMULATED : " << nmtra_total << endl; cout << endl; cout << "RC Refset : " << nmtra_fast_reco << endl; cout << "RC Extras : " << nmtra_slow_reco << endl; cout << "ghosts : " << nghosts << endl; cout << "clones : " << nclones << endl; cout << "ALL RECONSTRUCTED : " << ntracks << endl; cout << endl; cout << "Refset efficiency : " << p_fast << endl; cout << "Allset efficiency : " << p_total << endl; cout << "Extra efficiency : " << p_slow << endl; cout << "clone probability : " << p_clone << endl; cout << "ghost probability : " << p_ghost << endl; #endif//XXX cout << endl; cout << "L1 ACCUMULATED STAT : " << L1_NEVENTS << " EVENTS " << endl << endl; if (L1_NMTRA_D0!=0) cout << "D0 efficiency : " << P_D0 <<" | "<first || !imap->second ) continue; CbmL1Track &T = *(imap->first); CbmL1MCTrack &MC = *(imap->second); if( MC.iPileUp ) continue; nhits->Fill( T.StsHits.size() ); if ( !T.is_long ) continue; if ( fabs(MC.z) >=1. ) continue; if ( MC.p < 0.050 ) continue; MC_P->Fill(fabs(MC.p/MC.q)); MC_vx->Fill( MC.x ); MC_vy->Fill( MC.y ); MC_vz->Fill( MC.z ); { // P resolution vs P if ( fabs( T.T[4] )>1.e-10 ) P_vs_P->Fill( MC.p, MC.q/MC.p/T.T[4] - 1. ); } h_chi2->Fill( T.chi2/T.NDF ); h_prob->Fill( TMath::Erf(sqrt(T.chi2/T.NDF)) ); { // Perf at first hit CbmL1MCPoint *p=0, *plast = 0; for(vector::iterator ip = MC.Points.begin(); ip!= MC.Points.end(); ip++) { if( !p || fabs( (*ip)->z - T.T[5])z - T.T[5]) ) p = *ip; if( !plast || fabs( (*ip)->z - T.TLast[5])z - T.TLast[5]) ) plast = *ip; } if( p ) { KF.Propagate( p->T, 0, T.T[5], p->T[4]); FitPerformanceFillPulls( p->T, T.T, T.C, frst ); } if( plast ) { KF.Propagate( plast->TOut, 0, T.TLast[5], plast->TOut[4]); FitPerformanceFillPulls( plast->TOut, T.TLast, T.CLast, last ); } } { // perf at vertex CbmL1TrackPar t = T; t.Extrapolate( MC.z ); FitPerformanceFillPulls( MC.T, t.T, t.C, extr ); t.Extrapolate( PrimVtx.z ); Double_t ChiPrimary = CbmKFMath::getDeviation( t.T[0], t.T[1], t.C, PrimVtx.x, PrimVtx.y, PrimVtx.C ); if( MC.mother_ID==-1) VtxChiPrim->Fill( ChiPrimary ); else VtxChiSec->Fill( ChiPrimary ); if ( MC.mother_ID==-1){ t.Fit2Vertex( PrimVtx ); FitPerformanceFillPulls( MC.T, t.T, t.C, fvtx ); } } } }// Tracks // primary vertex fit performance { CbmL1Vtx &v = PrimVtx; vtx[0]->Fill( v.x - v.MC_x ); vtx[1]->Fill( v.y - v.MC_y ); vtx[2]->Fill( v.z - v.MC_z ); if ( v.C[ 0] >1.e-20 ) vtx[3]->Fill( ( v.x - v.MC_x )/sqrt(v.C[ 0]) ); if ( v.C[ 2] >1.e-20 ) vtx[4]->Fill( ( v.y - v.MC_y )/sqrt(v.C[ 2]) ); if ( v.C[ 5] >1.e-20 ) vtx[5]->Fill( ( v.z - v.MC_z )/sqrt(v.C[ 5]) ); vtx[6]->Fill( v.chi2/v.NDF ); vtx[7]->Fill( stat_vtx_ntracks ); vtx[8]->Fill( TMath::Erf(sqrt(v.chi2/v.NDF)) ); } cout<<" Fit time per track = "<= 100){//write at every 100th event if ( L1_NEVENTS % 100 == 0) write = true; } if (write){ TDirectory *curr = gDirectory; // Open output file and write histograms TFile* outfile = new TFile("L1_histo.root","RECREATE"); TIter hiter(listHisto); while (TObject* obj = hiter()) obj->Write(); outfile->Close(); curr->cd(); } }