BmnRoot
Loading...
Searching...
No Matches
CbmL1Performance.cxx
Go to the documentation of this file.
1/*
2 *====================================================================
3 *
4 * CBM Level 1 Reconstruction
5 *
6 * Authors: I.Kisel, S.Gorbunov
7 *
8 * e-mail : ikisel@kip.uni-heidelberg.de
9 *
10 *====================================================================
11 *
12 * L1 Fit performance
13 *
14 *====================================================================
15 */
16#include "CbmL1.h"
17
18#include "CbmL1Constants.h"
19#include "L1Algo/L1Algo.h"
20#include "L1Algo/L1StsHit.h"
21#include "L1Algo/L1Extrapolation.h" // for vertex pulls
22#include "L1Algo/L1AddMaterial.h" // for vertex pulls
23#include "FairTrackParam.h" // for vertex pulls
24#include "CbmKFTrack.h" // for vertex pulls
25
26#include "CbmKF.h"
27#include "CbmKFMath.h"
28#include "CbmMvdHit.h"
29#include "CbmMvdHitMatch.h"
30#include "CbmStsDigi.h"
31#include "CbmStsSensor.h" // for field FieldCheck.
32#include "CbmStsSector.h" // for field FieldCheck.
33#include "CbmStsStation.h" // for field FieldCheck.
34#include "CbmStsPoint.h"
35#include "CbmStsHit.h"
36#include "CbmMvdPoint.h"
37#include "CbmMCTrack.h"
38
39#include "CbmL1Counters.h"
40#include "CbmL1PartEfficiencies.h"
41#include "CbmL1ParticlesFinder.h"
42
43#include "TH1.h"
44#include "TH2.h"
45#include "TProfile.h"
46#include "TFile.h"
47
48#include <iostream>
49#include <vector>
50#include <list>
51#include <map>
52
53using std::cout;
54using std::endl;
55using std::ios;
56using std::vector;
57using std::pair;
58using std::map;
59
60void CbmL1::TrackMatch(){
61 map <int, CbmL1MCTrack*> pMCTrackMap; pMCTrackMap.clear();
62
63 // fill pMCTrackMap
64 for( vector<CbmL1MCTrack>::iterator i = vMCTracks.begin(); i != vMCTracks.end(); ++i){
65 CbmL1MCTrack &MC = *i;
66
67 if (pMCTrackMap.find(MC.ID) == pMCTrackMap.end()){
68 pMCTrackMap.insert(pair<int, CbmL1MCTrack*>(MC.ID, &MC));
69 }
70 else{
71 cout << "*** L1: Track ID " << MC.ID << " encountered twice! ***" << endl;
72 }
73 }
74 // -- prepare information about reconstructed tracks
75 const int nRTracks = vRTracks.size();
76 for (int iR = 0; iR < nRTracks; iR++){
77 CbmL1Track* prtra = &(vRTracks[iR]);
78
79 int hitsum = prtra->StsHits.size(); // number of hits in track
80
81 // count how many hits from each mcTrack belong to current recoTrack
82 map<int, int > &hitmap = prtra->hitMap; // how many hits from each mcTrack belong to current recoTrack
83 for (vector<int>::iterator ih = (prtra->StsHits).begin(); ih != (prtra->StsHits).end(); ++ih){
84
85 const int nMCPoints = vStsHits[*ih].mcPointIds.size();
86 for (int iP = 0; iP < nMCPoints; iP++){
87 int iMC = vStsHits[*ih].mcPointIds[iP];
88 int ID = -1;
89 if (iMC >= 0) ID = vMCPoints[iMC].ID;
90 if(hitmap.find(ID) == hitmap.end())
91 hitmap[ID] = 1;
92 else{
93 hitmap[ID] += 1;
94 }
95 } // for iPoint
96 } // for iHit
97
98 // RTrack <-> MCTrack identification
99 double max_percent = 0.0; // [%]. maximum persent of hits, which belong to one mcTrack.
100 for( map<int, int >::iterator posIt = hitmap.begin(); posIt != hitmap.end(); posIt++ ){ // loop over all touched MCTracks
101
102 if (posIt->first < 0) continue; // not a MC track - based on fake hits
103
104 // count max-purity
105 if (double(posIt->second) > max_percent*double(hitsum))
106 max_percent = double(posIt->second)/double(hitsum);
107
108 // set relation to the mcTrack
109 if ( double(posIt->second) > CbmL1Constants::MinPurity * double(hitsum) ){ // found correspondent MCTrack
110 if (pMCTrackMap.find(posIt->first) == pMCTrackMap.end()) continue;
111 CbmL1MCTrack* pmtra = pMCTrackMap[posIt->first];
112
113 pmtra->AddRecoTrack(prtra);
114 prtra->AddMCTrack(pmtra);
115 }
116 else{
117 if (pMCTrackMap.find(posIt->first) == pMCTrackMap.end()) continue;
118 CbmL1MCTrack* pmtra = pMCTrackMap[posIt->first];
119
120 pmtra->AddTouchTrack(prtra);
121 }
122 } // for hitmap
123 prtra->SetMaxPurity(max_percent);
124 } // for reco tracks
125} //
126
127
128
129struct TL1PerfEfficiencies: public TL1Efficiencies
130{
131 TL1PerfEfficiencies():TL1Efficiencies(),
132ratio_killed(),
133ratio_clone(),
134ratio_length(),
135ratio_fakes(),
136killed(),
137clone(),
138reco_length(),
139reco_fakes(),
140mc_length(),
141mc_length_hits()
142 {
143 // add total efficiency
144 AddCounter("long_fast_prim" ,"LongRPrim efficiency");
145 AddCounter("fast_prim" ,"RefPrim efficiency");
146 AddCounter("fast_sec" ,"RefSec efficiency");
147 AddCounter("fast" ,"Refset efficiency");
148 AddCounter("total" ,"Allset efficiency");
149 AddCounter("slow_prim" ,"ExtraPrim efficiency");
150 AddCounter("slow_sec" ,"ExtraSec efficiency");
151 AddCounter("slow" ,"Extra efficiency");
152 AddCounter("d0" ,"D0 efficiency");
153 AddCounter("short" ,"Short123s efficiency");
154 AddCounter("shortPion" ,"Short123s pion eff");
155 AddCounter("shortProton" ,"Short123s proton eff");
156 AddCounter("shortKaon" ,"Short123s kaon eff");
157 AddCounter("shortE" ,"Short123s e eff");
158 AddCounter("shortRest" ,"Short123s rest eff");
159
160 AddCounter("fast_sec_e" ,"RefSec E efficiency");
161 AddCounter("fast_e" ,"Refset E efficiency");
162 AddCounter("total_e" ,"Allset E efficiency");
163 AddCounter("slow_sec_e" ,"ExtraSecE efficiency");
164 AddCounter("slow_e" ,"Extra E efficiency");
165 }
166
167 virtual ~TL1PerfEfficiencies(){};
168
169 virtual void AddCounter(TString shortname, TString name){
170 TL1Efficiencies::AddCounter(shortname, name);
171 ratio_killed.AddCounter();
172 ratio_clone.AddCounter();
173 ratio_length.AddCounter();
174 ratio_fakes.AddCounter();
175 killed.AddCounter();
176 clone.AddCounter();
177 reco_length.AddCounter();
178 reco_fakes.AddCounter();
179 mc_length.AddCounter();
180 mc_length_hits.AddCounter();
181 };
182
192
201
202 void Inc(bool isReco, bool isKilled, double _ratio_length, double _ratio_fakes, int _nclones, int _mc_length, int _mc_length_hits, TString name){
203 TL1Efficiencies::Inc(isReco, name);
204
205 const int index = indices[name];
206
207 if (isKilled) killed.counters[index]++;
208 reco_length.counters[index] += _ratio_length;
209 reco_fakes.counters[index] += _ratio_fakes;
210 clone.counters[index] += _nclones;
211 mc_length.counters[index] += _mc_length;
212 mc_length_hits.counters[index] += _mc_length_hits;
213 };
214
215 void PrintEff(){
216 L1_assert(nEvents != 0);
217
218 cout.setf(ios::fixed);
219 cout.setf(ios::showpoint);
220 cout.precision(3);
221 cout.setf(ios::right);
222 cout << "Track category : " << " Eff " <<" / "<< "Killed" <<" / "<< "Length" <<" / "<< "Fakes " <<" / "<< "Clones" <<" / "<< "All Reco" <<" | "<< " All MC " <<" / "<< "MCl(hits)" <<" / "<< "MCl(MCps)" << endl;
223
224 int NCounters = mc.NCounters;
225 for (int iC = 0; iC < NCounters; iC++){
226 if (( names[iC] != "D0 efficiency") || (mc.counters[iC] != 0))
227 cout << names[iC] << " : "
228 << ratio_reco.counters[iC]
229 << " / " << ratio_killed.counters[iC] // tracks with aren't reco because other tracks takes their hit(-s)
230 << " / " << ratio_length.counters[iC] // nRecoMCHits/nMCHits
231 << " / " << ratio_fakes.counters[iC] // nFakeHits/nRecoAllHits
232 << " / " << ratio_clone.counters[iC] // nCloneTracks/nMCTracks
233 << " / " << setw(8) << reco.counters[iC]/double(nEvents)
234 << " | " << setw(8) << mc.counters[iC]/double(nEvents)
235 << " / " << mc_length_hits.counters[iC]/double(mc.counters[iC])
236 << " / " << mc_length.counters[iC]/double(mc.counters[iC])
237 << endl;
238 }
239 cout << "Ghost probability : " << ratio_ghosts <<" | "<< ghosts << endl;
240 };
241
242 TL1TracksCatCounters<double> ratio_killed;
243 TL1TracksCatCounters<double> ratio_clone;
244 TL1TracksCatCounters<double> ratio_length;
245 TL1TracksCatCounters<double> ratio_fakes;
246
247 TL1TracksCatCounters<int> killed;
248 TL1TracksCatCounters<int> clone;
249 TL1TracksCatCounters<double> reco_length;
250 TL1TracksCatCounters<double> reco_fakes;
251 TL1TracksCatCounters<int> mc_length;
252 TL1TracksCatCounters<int> mc_length_hits;
253};
254
255
256void CbmL1::EfficienciesPerformance()
257{
258 static TL1PerfEfficiencies L1_NTRA; // average efficiencies
259
260 static int L1_NEVENTS = 0;
261 static double L1_CATIME = 0.0;
262
263 TL1PerfEfficiencies ntra; // efficiencies for current event
264
265 for (vector<CbmL1Track>::iterator rtraIt = vRTracks.begin(); rtraIt != vRTracks.end(); ++rtraIt){
266 ntra.ghosts += rtraIt->IsGhost();
267// if(rtraIt->IsGhost()){ // Debug.
268// cout << " " << rtraIt->GetNOfHits() << " " << 1./rtraIt->T[5] << " " << rtraIt->GetMaxPurity() << " | ";
269// for( map<int, int>::iterator posIt = rtraIt->hitMap.begin(); posIt != rtraIt->hitMap.end(); posIt++ ){
270// cout << " (" << posIt->second << " " << posIt->first << ") ";
271// }
272// cout << endl;
273// }
274 }
275
276 int sta_nhits[algo->NStations], sta_nfakes[algo->NStations];
277
278 for (int i = 0; i < algo->NStations; i++){
279 sta_nhits[i] = 0;
280 sta_nfakes[i] = 0;
281 }
282
283 for ( vector<CbmL1MCTrack>::iterator mtraIt = vMCTracks.begin(); mtraIt != vMCTracks.end(); mtraIt++ ) {
284 CbmL1MCTrack &mtra = *(mtraIt);
285
286// if( !( mtra.pdg == -11 && mtra.mother_ID == -1 ) ) continue; // electrons only
287 if( ! mtra.IsReconstructable() && ! mtra.IsAdditional() ) continue;
288
289 // -- find used constans --
290 // is track reconstructed
291 const bool reco = (mtra.IsReconstructed());
292 // is track killed. At least one hit of it belong to some recoTrack
293 const bool killed = !reco && mtra.IsDisturbed();
294 // ration length for current mcTrack
295 vector< CbmL1Track* >& rTracks = mtra.GetRecoTracks(); // for length calculations
296 double ratio_length = 0;
297 double ratio_fakes = 0;
298 double mc_length_hits = mtra.NStations();
299 int mc_length = mtra.NMCStations();
300 if (reco){
301 for (unsigned int irt = 0; irt < rTracks.size(); irt++) {
302 ratio_length += static_cast<double>( rTracks[irt]->GetNOfHits() )*rTracks[irt]->GetMaxPurity() / mc_length_hits;
303 ratio_fakes += 1 - rTracks[irt]->GetMaxPurity();
304 }
305 }
306 // number of clones
307 int nclones = 0;
308 if (reco) nclones = mtra.GetNClones();
309 // if (nclones){ // Debug. Look at clones
310 // for (int irt = 0; irt < rTracks.size(); irt++){
311 // const int ista = vHitStore[rTracks[irt]->StsHits[0]].iStation;
312 // cout << rTracks[irt]->GetNOfHits() << "(" << ista << ") ";
313 // }
314 // cout << mtra.NStations() << endl;
315 // }
316
317 if ( mtra.IsAdditional() ){ // short
318 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "short");
319 switch ( mtra.pdg ) {
320 case 11:
321 case -11:
322 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortE");
323 break;
324 case 211:
325 case -211:
326 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortPion");
327 break;
328 case 321:
329 case -321:
330 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortKaon");
331 break;
332 case 2212:
333 case -2212:
334 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortProton");
335 break;
336 default:
337 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortRest");
338 }
339 }
340 else { // separate all efficiecies from short eff
341
342 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "total");
343
344 if (( mtra.IsPrimary() )&&(mtra.z > 0)){ // D0
345 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "d0");
346 }
347
348 if ( mtra.p > CbmL1Constants::MinRefMom ){ // reference tracks
349 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast");
350
351 if ( mtra.IsPrimary() ){ // reference primary
352 if ( mtra.NStations() == NStation ){ // long reference primary
353 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "long_fast_prim");
354 }
355 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_prim");
356 }
357 else{ // reference secondary
358 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_sec");
359 }
360 }
361 else{ // extra set of tracks
362 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow");
363
364 if ( mtra.IsPrimary() ){ // extra primary
365 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_prim");
366 }
367 else{
368 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_sec");
369 }
370 } // if extra
371
372 if ( mtra.pdg == 11 || mtra.pdg == -11 ) {
373 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "total_e");
374
375 if ( mtra.p > CbmL1Constants::MinRefMom ){ // reference tracks
376 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_e");
377
378 if ( mtra.IsPrimary() ){ // reference primary
379 }
380 else{ // reference secondary
381 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_sec_e");
382 }
383 }
384 else{ // extra set of tracks
385 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_e");
386
387 if ( mtra.IsPrimary() ){ // extra primary
388 }
389 else{
390 ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_sec_e");
391 }
392 } // if extra
393 }
394 }
395
396 } // for mcTracks
397
398 L1_CATIME += algo->CATime;
399 L1_NEVENTS++;
400 ntra.IncNEvents();
401 L1_NTRA += ntra;
402
403 ntra.CalcEff();
404 L1_NTRA.CalcEff();
405
406 // cout.precision(3);
407 if( fVerbose ){
408 if( fVerbose > 1 ){
409 ntra.PrintEff();
410 cout << "Number of true and fake hits in stations: " << endl;
411 for (int i = 0; i < algo->NStations; i++){
412 cout << sta_nhits[i]-sta_nfakes[i] << "+" << sta_nfakes[i] << " ";
413 }
414 cout << endl;
415 } // fVerbose > 1
416 cout << endl;
417 cout << "L1 ACCUMULATED STAT : " << L1_NEVENTS << " EVENTS " << endl << endl;
418 L1_NTRA.PrintEff();
419 cout << "MC tracks/event found : " << int(double(L1_NTRA.reco.counters[L1_NTRA.indices["total"]])/double(L1_NEVENTS)) << endl;
420 cout<<endl;
421 cout<<"CA Track Finder: " << L1_CATIME/L1_NEVENTS << " s/ev" << endl << endl;
422 }
423} // void CbmL1::Performance()
424
425
426void CbmL1::GetMCParticles()
427{
428 vMCParticles.clear();
429 // convert
430 for(int i=0; i < listMCTracks->GetEntriesFast(); i++)
431 {
432 CbmMCTrack &mtra = *(L1_DYNAMIC_CAST<CbmMCTrack*>(listMCTracks->At(i)));
433 KFMCParticle part;
434 part.SetMCTrackID( i );
435 part.SetMotherId ( mtra.GetMotherId() );
436 part.SetPDG ( mtra.GetPdgCode() );
437 vMCParticles.push_back( part );
438 }
439
440 // find relations
441 const unsigned int nParticles = vMCParticles.size();
442 for ( unsigned int iP = 0; iP < nParticles; iP++ ) {
443 KFMCParticle &part = vMCParticles[iP];
444 for(unsigned int iP2 = 0; iP2 < nParticles; iP2++) {
445 KFMCParticle &part2 = vMCParticles[iP2];
446
447 if(part.GetMotherId() == part2.GetMCTrackID()) {
448 part2.AddDaughter(iP);
449 }
450 }
451 }
452 // // correct idices according to the new array
453 // for ( unsigned int iP = 0; iP < nParticles; iP++ ) {
454 // KFMCParticle &part = vMCParticles[iP];
455 // for(unsigned int iP2 = 0; iP2 < nParticles; iP2++) {
456 // KFMCParticle &part2 = vMCParticles[iP2];
457
458 // if(part.GetMotherId() == part2.GetMCTrackID()) {
459 // part.SetMotherId (iP2);
460 // }
461 // }
462 // }
463}
464
465void CbmL1::FindReconstructableMCParticles()
466{
467 const unsigned int nParticles = vMCParticles.size();
468
469 for ( unsigned int iP = 0; iP < nParticles; iP++ ) {
470 KFMCParticle &part = vMCParticles[iP];
471 CheckMCParticleIsReconstructable(part);
472 }
473}
474
475void CbmL1::CheckMCParticleIsReconstructable(KFMCParticle &part)
476{
477
478 if ( part.IsReconstructable() ) return;
479
480 const int nParticles = 5;
481 int partPDG[nParticles] = {310,3122,3312,-3312,3334};
482 vector<int> partDaughterPdg[nParticles];
483
484 partDaughterPdg[0].push_back( 211);
485 partDaughterPdg[0].push_back(-211);
486
487 partDaughterPdg[1].push_back(2212);
488 partDaughterPdg[1].push_back(-211);
489
490 partDaughterPdg[2].push_back(3122);
491 partDaughterPdg[2].push_back(-211);
492
493 partDaughterPdg[3].push_back(3122);
494 partDaughterPdg[3].push_back( 211);
495
496 partDaughterPdg[4].push_back(3122);
497 partDaughterPdg[4].push_back(-321);
498
499 // tracks
500 if ( /*part.NDaughters() == 0*/ part.GetPDG() == -211 ||
501 part.GetPDG() == 211 ||
502 part.GetPDG() == 2212 ||
503 part.GetPDG() == 321 ||
504 part.GetPDG() == -321 ) { // TODO other particles
505
506 switch ( fFindParticlesMode ) {
507 case 1:
508 part.SetAsReconstructable();
509 break;
510 case 2:
511 unsigned int i;
512 for ( i = 0; i < vMCTracks.size(); i++ ) // TODO opt
513 if (vMCTracks[i].ID == part.GetMCTrackID()) break;
514 if ( i != vMCTracks.size() )
515 if ( vMCTracks[i].IsReconstructable() )
516 part.SetAsReconstructable();
517 break;
518 case 3:
519// int i;
520 for ( i = 0; i < vMCTracks.size(); i++ )
521 if (vMCTracks[i].ID == part.GetMCTrackID()) break;
522 if ( i != vMCTracks.size() )
523 if ( vMCTracks[i].IsReconstructed() )
524 part.SetAsReconstructable();
525 break;
526 default:
527 L1_assert(0);
528 };
529 }
530
531 // mother particles
532 else
533 {
534 for(int iPart=0; iPart<nParticles; iPart++)
535 {
536 if(part.GetPDG() == partPDG[iPart])
537 {
538 const unsigned int nDaughters = partDaughterPdg[iPart].size();
539 if( part.GetDaughterIds().size() != nDaughters ) return;
540 int pdg[nDaughters];
541
542 for(unsigned int iD=0; iD<nDaughters; iD++)
543 pdg[iD] = vMCParticles[part.GetDaughterIds()[iD]].GetPDG();
544
545 bool isDaughterFound[nDaughters];
546 for(unsigned int iDMC=0; iDMC<nDaughters; iDMC++)
547 isDaughterFound[iDMC] = 0;
548
549 bool isReco = 1;
550 for(unsigned int iDMC=0; iDMC<nDaughters; iDMC++)
551 for(unsigned int iD=0; iD<nDaughters; iD++)
552 if(pdg[iD] == partDaughterPdg[iPart][iDMC]) isDaughterFound[iDMC] = 1;
553
554 for(unsigned int iDMC=0; iDMC<nDaughters; iDMC++)
555 isReco = isReco && isDaughterFound[iDMC];
556
557 if(!isReco) return;
558 }
559 }
560
561 const vector<int>& dIds = part.GetDaughterIds();
562 const unsigned int nD = dIds.size();
563 bool reco = 1;
564 for ( unsigned int iD = 0; iD < nD && reco; iD++ ) {
565 KFMCParticle &dp = vMCParticles[dIds[iD]];
566 CheckMCParticleIsReconstructable(dp);
567 reco &= dp.IsReconstructable();
568 }
569 if (reco) part.SetAsReconstructable();
570
571 }
572
573}
574
576void CbmL1::MatchParticles()
577{
578 // get all reco particles ( temp )
579 MCtoRParticleId.clear();
580 RtoMCParticleId.clear();
581 MCtoRParticleId.resize(vMCParticles.size());
582 RtoMCParticleId.resize(vRParticles.size() );
583
584 // match tracks ( particles which are direct copy of tracks )
585 for( unsigned int iRP = 0; iRP < vRParticles.size(); iRP++ ) {
586 CbmKFParticle &rPart = vRParticles[iRP];
587 L1_assert( rPart.NDaughters() > 0 );
588 L1_ASSERT( static_cast<unsigned int>(rPart.Id()) == iRP, rPart.Id() << " != " << iRP );
589 if (rPart.NDaughters() != 1) continue;
590
591 const int rTrackId = rPart.DaughterIds()[0];
592
593 if ( vRTracks[rTrackId].IsGhost() ) continue;
594 const int mcTrackId = vRTracks[rTrackId].GetMCTrack()->ID;
595
596 for ( unsigned int iMP = 0; iMP < vMCParticles.size(); iMP++ ) {
597 KFMCParticle &mPart = vMCParticles[iMP];
598 if ( mPart.GetMCTrackID() == mcTrackId ) { // match is found
599 if( mPart.GetPDG() == rPart.GetPDG() ) {
600 MCtoRParticleId[iMP].ids.push_back(iRP);
601 RtoMCParticleId[iRP].ids.push_back(iMP);
602 }
603 else {
604 MCtoRParticleId[iMP].idsMI.push_back(iRP);
605 RtoMCParticleId[iRP].idsMI.push_back(iMP);
606 }
607 }
608 }
609 }
610
611 // match created mother particles
612 for( unsigned int iRP = 0; iRP < vRParticles.size(); iRP++ ) {
613 CbmKFParticle &rPart = vRParticles[iRP];
614 const unsigned int NRDaughters = rPart.NDaughters();
615 if (NRDaughters < 2) continue;
616
617 unsigned int iD = 0;
618 int mmId = -2; // mother MC id
619 {
620 const int rdId = rPart.DaughterIds()[iD];
621 if ( !RtoMCParticleId[rdId].IsMatched() ) continue;
622 const int mdId = RtoMCParticleId[rdId].GetBestMatch();
623 mmId = vMCParticles[mdId].GetMotherId();
624 }
625 iD++;
626 for ( ; iD < NRDaughters; iD++ ) {
627 const int rdId = rPart.DaughterIds()[iD];
628 if ( !RtoMCParticleId[rdId].IsMatched() ) break;
629 const int mdId = RtoMCParticleId[rdId].GetBestMatch();
630 if( vMCParticles[mdId].GetMotherId() != mmId ) break;
631 }
632 if ( iD == NRDaughters && mmId != -1 ) { // match is found and it is not primary vertex
633 KFMCParticle &mmPart = vMCParticles[mmId];
634 if( mmPart.GetPDG() == rPart.GetPDG() &&
635 mmPart.NDaughters() == rPart.NDaughters() ) {
636 MCtoRParticleId[mmId].ids.push_back(iRP);
637 RtoMCParticleId[iRP].ids.push_back(mmId);
638 }
639 else {
640 MCtoRParticleId[mmId].idsMI.push_back(iRP);
641 RtoMCParticleId[iRP].idsMI.push_back(mmId);
642 }
643 }
644 }
645
646}
647
648void CbmL1::PartEffPerformance()
649{
650 static int NEVENTS = 0;
651
652 static CbmL1PartEfficiencies PARTEFF; // average efficiencies
653
654 CbmKFPartEfficiencies partEff; // efficiencies for current event
655
656 const int NRP = vRParticles.size();
657 for ( int iP = 0; iP < NRP; ++iP ) {
658 const CbmKFParticle &part = vRParticles[iP];
659 const int pdg = part.GetPDG();
660
661 const bool isBG = RtoMCParticleId[iP].idsMI.size() != 0;
662 const bool isGhost = !RtoMCParticleId[iP].IsMatched();
663
664 for(int iPart=0; iPart<PARTEFF.nParticles; iPart++)
665 if ( pdg == PARTEFF.partPDG[iPart] )
666 partEff.IncReco(isGhost, isBG, PARTEFF.partName[iPart].Data());
667 }
668
669
670 const int NMP = vMCParticles.size();
671 for ( int iP = 0; iP < NMP; ++iP ) {
672 const KFMCParticle &part = vMCParticles[iP];
673 if ( !part.IsReconstructable() ) continue;
674 const int pdg = part.GetPDG();
675 const int mId = part.GetMotherId();
676
677 const bool isReco = MCtoRParticleId[iP].ids.size() != 0;
678 const int nClones = MCtoRParticleId[iP].ids.size() - 1;
679
680 for(int iPart=0; iPart<PARTEFF.nParticles; iPart++)
681 {
682 if ( pdg == PARTEFF.partPDG[iPart] ) {
683 partEff.Inc(isReco, nClones, PARTEFF.partName[iPart].Data());
684 if ( mId == -1 )
685 partEff.Inc(isReco, nClones, (PARTEFF.partName[iPart]+"_prim").Data());
686 else
687 partEff.Inc(isReco, nClones, (PARTEFF.partName[iPart]+"_sec").Data());
688 }
689 }
690 }
691
692 NEVENTS++;
693
694 PARTEFF += partEff;
695
696 partEff.CalcEff();
697 PARTEFF.CalcEff();
698
699 // cout.precision(3);
700 if( fVerbose ){
701 if(NEVENTS%500 == 0)
702 {
703 cout << " ---- KF Particle finder --- " << endl;
704 // cout << "L1 STAT : " << fNEvents << " EVENT " << endl << endl;
705 //partEff.PrintEff();
706 // cout << endl;
707 cout << "ACCUMULATED STAT : " << NEVENTS << " EVENTS " << endl << endl;
708 PARTEFF.PrintEff();
709
710 cout<<endl;
711 // cout<<"CA Track Finder: " << L1_CATIME/L1_fNEvents << " s/ev" << endl << endl;
712 }
713 }
714}
715
716void CbmL1::PartHistoPerformance()
717{
718 static const int NParticles = 2; //Ks, Lambda
719
720 static const int nFitQA = 16;
721 static TH1F *hFitDaughtersQA[NParticles][nFitQA];
722 static TH1F *hFitQA[NParticles][nFitQA];
723
724 static const int nHistoPartParam = 5;
725 static TH1F *hPartParam[NParticles][nHistoPartParam]; // mass, decay length, c*tau, chi/ndf, prob ...
726 static TH1F *hPartParamBG[NParticles][nHistoPartParam];
727 static TH1F *hPartParamSignal[NParticles][nHistoPartParam];
728 static const int nHistoPartParamQA = 3;
729 //static TH1F *hPartParamQA[NParticles][nHistoPartParamQA*2]; // residuals and pulls of these parameters
730
731 static const int nHistosPV = 6;
732 static TH1F *hPVFitQa[nHistosPV];
733
734 TString patName[NParticles] = {"K0s","Lambda"};
735
736 static bool first_call = 1;
737
738 if ( first_call )
739 {
740 first_call = 0;
741
742 TDirectory *currentDir = gDirectory;
743 gDirectory = histodir;
744 gDirectory->cd("L1");
745 gDirectory->mkdir("Particles");
746 gDirectory->cd("Particles");
747 {
748 for(int iPart=0; iPart<NParticles; ++iPart)
749 {
750 gDirectory->mkdir(patName[iPart].Data());
751 gDirectory->cd(patName[iPart].Data());
752 {
753 TString res = "res";
754 TString pull = "pull";
755
756 gDirectory->mkdir("DaughtersQA");
757 gDirectory->cd("DaughtersQA");
758 {
759 TString parName[nFitQA/2] = {"X","Y","Z","Px","Py","Pz","E","M"};
760 int nBins = 100;
761 float xMax[nFitQA/2] = {0.15,0.15,0.03,0.01,0.01,0.06,0.06,0.01};
762
763 for( int iH=0; iH<nFitQA/2; iH++ ){
764 hFitDaughtersQA[iPart][iH] = new TH1F((res+parName[iH]).Data(),
765 (res+parName[iH]).Data(),
766 nBins, -xMax[iH],xMax[iH]);
767 hFitDaughtersQA[iPart][iH+8] = new TH1F((pull+parName[iH]).Data(),
768 (pull+parName[iH]).Data(),
769 nBins, -6,6);
770 }
771 }
772 gDirectory->cd(".."); //particle directory
773
774 gDirectory->mkdir("FitQA");
775 gDirectory->cd("FitQA");
776 {
777 TString parName[nFitQA/2] = {"X","Y","Z","Px","Py","Pz","E","M"};
778 int nBins = 50;
779 float xMax[nFitQA/2] = {0.15,0.15,1.2,0.02,0.02,0.15,0.15,0.006};
780
781 for( int iH=0; iH<nFitQA/2; iH++ ){
782 hFitQA[iPart][iH] = new TH1F((res+parName[iH]).Data(),
783 (res+parName[iH]).Data(),
784 nBins, -xMax[iH],xMax[iH]);
785 hFitQA[iPart][iH+8] = new TH1F((pull+parName[iH]).Data(),
786 (pull+parName[iH]).Data(),
787 nBins, -6,6);
788 }
789 }
790 gDirectory->cd(".."); //particle directory
791
792 gDirectory->mkdir("Parameters");
793 gDirectory->cd("Parameters");
794 {
795 TString parName[nHistoPartParam] = {"M","DL","c#tau","chi2ndf","prob"};
796 int nBins[nHistoPartParam] = {1000,100,100,100,100};
797 float xMin[nHistoPartParam] = {0.3f,0. , 0., 0.,0.};
798 float xMax[nHistoPartParam] = {1.3f,30.,30.,20.,1.};
799 if(iPart == 1)
800 {
801 xMin[0] = 1.0;
802 xMax[0] = 2.0;
803 }
804
805 for(int iH=0; iH<nHistoPartParam; iH++)
806 hPartParam[iPart][iH] = new TH1F(parName[iH].Data(),parName[iH].Data(),
807 nBins[iH],xMin[iH],xMax[iH]);
808
809 gDirectory->mkdir("Signal");
810 gDirectory->cd("Signal");
811 {
812 for(int iH=0; iH<nHistoPartParam; iH++)
813 hPartParamSignal[iPart][iH] = new TH1F((parName[iH]).Data(),(parName[iH]).Data(),
814 nBins[iH],xMin[iH],xMax[iH]);
815
816 gDirectory->mkdir("QA");
817 gDirectory->cd("QA");
818 {
819 //int nBinsQA = 50;
820 //float xMaxQA[nHistoPartParamQA] = {0.01,0.001,0.001};
821 for( int iH=0; iH<nHistoPartParamQA; iH++ ){
822 //hPartParamQA[iPart][iH] =
823 // new TH1F((res+parName[iH]).Data(), (res+parName[iH]).Data(), nBinsQA, -xMaxQA[iH],xMaxQA[iH]);
824 //hPartParamQA[iPart][iH+nHistoPartParamQA] =
825 // new TH1F((pull+parName[iH]).Data(), (pull+parName[iH]).Data(), nBinsQA, -6,6);
826 }
827 }
828 gDirectory->cd(".."); // particle directory / Parameters / Signal
829
830 }
831 gDirectory->cd(".."); // particle directory / Parameters
832 gDirectory->mkdir("Background");
833 gDirectory->cd("Background");
834 {
835 for(int iH=0; iH<nHistoPartParam; iH++)
836 hPartParamBG[iPart][iH] = new TH1F((parName[iH]).Data(),(parName[iH]).Data(),
837 nBins[iH],xMin[iH],xMax[iH]);
838 }
839 gDirectory->cd(".."); // particle directory
840 }
841 gDirectory->cd(".."); //particle directory
842 }
843 gDirectory->cd(".."); //L1
844 }
845 }
846 gDirectory->cd(".."); //L1
847 gDirectory->mkdir("PrimaryVertexQA");
848 gDirectory->cd("PrimaryVertexQA");
849 {
850 struct {
851 string name;
852 string title;
853 Int_t n;
854 Double_t l,r;
855 } Table[nHistosPV]=
856 {
857 {"PVResX", "x_{rec}-x_{mc}, cm", 100, -0.006f, 0.006f},
858 {"PVResY", "y_{rec}-y_{mc}, cm", 100, -0.006f, 0.006f},
859 {"PVResZ", "z_{rec}-z_{mc}, cm", 100, -0.06f, 0.06f},
860 {"PVPullX", "Pull X", 100, -6.f, 6.f},
861 {"PVPullY", "Pull Y", 100, -6.f, 6.f},
862 {"PVPullZ", "Pull Z", 100, -6.f, 6.f}
863 };
864 for(int iHPV=0; iHPV<nHistosPV; ++iHPV){
865 hPVFitQa[iHPV] = new TH1F(Table[iHPV].name.data(),Table[iHPV].title.data(),
866 Table[iHPV].n, Table[iHPV].l, Table[iHPV].r);
867 }
868 }
869 gDirectory = currentDir;
870 }
871
872 for(unsigned int iP=0; iP < vRParticles.size(); iP++)
873 {
874 int iParticle = -1;
875 switch( vRParticles[iP].GetPDG() ) {
876 case 310:
877 {
878 iParticle = 0;
879 }
880 break;
881 case 3122:
882 {
883 iParticle = 1;
884 }
885 }
886 if(iParticle < 0) continue;
887
888 Double_t M, ErrM;
889 Double_t dL, ErrdL; // decay length
890 Double_t cT, ErrcT; // c*tau
891 CbmKFParticle TempPart = vRParticles[iP];
892 TempPart.GetMass(M,ErrM);
893 TempPart.GetDecayLength(dL,ErrdL);
894 TempPart.GetLifeTime(cT,ErrcT);
895 Double_t chi2 = TempPart.GetChi2();
896 Int_t ndf = TempPart.GetNDF();
897 Double_t prob = TMath::Prob(chi2,ndf);
898
899 hPartParam[iParticle][0]->Fill(M);
900 hPartParam[iParticle][1]->Fill(dL);
901 hPartParam[iParticle][2]->Fill(cT);
902 hPartParam[iParticle][3]->Fill(chi2/ndf);
903 hPartParam[iParticle][4]->Fill(prob);
904
905 if(!RtoMCParticleId[iP].IsMatchedWithPdg()) //background
906 {
907 hPartParamBG[iParticle][0]->Fill(M);
908 hPartParamBG[iParticle][1]->Fill(dL);
909 hPartParamBG[iParticle][2]->Fill(cT);
910 hPartParamBG[iParticle][3]->Fill(chi2/ndf);
911 hPartParamBG[iParticle][4]->Fill(prob);
912 continue;
913 }
914 hPartParamSignal[iParticle][0]->Fill(M);
915 hPartParamSignal[iParticle][1]->Fill(dL);
916 hPartParamSignal[iParticle][2]->Fill(cT);
917 hPartParamSignal[iParticle][3]->Fill(chi2/ndf);
918 hPartParamSignal[iParticle][4]->Fill(prob);
919
920 int iMCPart = RtoMCParticleId[iP].GetBestMatchWithPdg();
921 KFMCParticle &mcPart = vMCParticles[iMCPart];
922 // Fit quality of the mother particle
923 {
924 int iMCTrack = mcPart.GetMCTrackID();
925 CbmMCTrack &mcTrack = *(L1_DYNAMIC_CAST<CbmMCTrack*>(listMCTracks->At(iMCTrack)));
926 int mcDaughterId = mcPart.GetDaughterIds()[0];
927 CbmMCTrack &mcDaughter = *(L1_DYNAMIC_CAST<CbmMCTrack*>(listMCTracks->At(mcDaughterId)));
928
929 Double_t decayVtx[3] = { mcDaughter.GetStartX(), mcDaughter.GetStartY(), mcDaughter.GetStartZ() };
930 Double_t recParam[8] = { 0 };
931 Double_t errParam[8] = { 0 };
932
933 for(int iPar=0; iPar<3; iPar++)
934 {
935 recParam[iPar] = TempPart.GetParameter(iPar);
936 Double_t error = TempPart.GetCovariance(iPar,iPar);
937 if(error < 0.) { error = 1.e20;}
938 errParam[iPar] = TMath::Sqrt(error);
939 }
940 TempPart.Extrapolate(TempPart.r , TempPart.GetDStoPoint(decayVtx));
941 for(int iPar=3; iPar<7; iPar++)
942 {
943 recParam[iPar] = TempPart.GetParameter(iPar);
944 Double_t error = TempPart.GetCovariance(iPar,iPar);
945 if(error < 0.) { error = 1.e20;}
946 errParam[iPar] = TMath::Sqrt(error);
947 }
948
949 Double_t Emc = sqrt(mcTrack.GetP()*mcTrack.GetP() + mcTrack.GetMass()*mcTrack.GetMass());
950 Double_t res[8] = {0},
951 pull[8] = {0},
952 mcParam[8] = { decayVtx[0], decayVtx[1], decayVtx[2],
953 mcTrack.GetPx(), mcTrack.GetPy(), mcTrack.GetPz(), Emc, mcTrack.GetMass() };
954 for(int iPar=0; iPar < 7; iPar++ )
955 {
956 res[iPar] = recParam[iPar] - mcParam[iPar];
957 if(fabs(errParam[iPar]) > 1.e-20) pull[iPar] = res[iPar]/errParam[iPar];
958 }
959 res[7] = M - mcParam[7];
960 if(fabs(ErrM) > 1.e-20) pull[7] = res[7]/ErrM;
961
962 for(int iPar=0; iPar < 8; iPar++ )
963 {
964 hFitQA[iParticle][iPar]->Fill(res[iPar]);
965 hFitQA[iParticle][iPar+8]->Fill(pull[iPar]);
966 }
967
968// Double_t mcT = mcDaughter.GetStartT() - mcTrack.GetStartT();
969
970 }
971 // Fit quality of daughters
972 for(int iD=0; iD<mcPart.NDaughters(); ++iD)
973 {
974 int mcDaughterId = mcPart.GetDaughterIds()[iD];
975 if(!MCtoRParticleId[mcDaughterId].IsMatchedWithPdg()) continue;
976 CbmMCTrack &mcTrack = *(L1_DYNAMIC_CAST<CbmMCTrack*>(listMCTracks->At(mcDaughterId)));
977 int recDaughterId = MCtoRParticleId[mcDaughterId].GetBestMatchWithPdg();
978 CbmKFParticle Daughter = vRParticles[recDaughterId];
979 Daughter.GetMass(M,ErrM);
980
981 Double_t decayVtx[3] = {mcTrack.GetStartX(), mcTrack.GetStartY(), mcTrack.GetStartZ()};
982 Daughter.Extrapolate(Daughter.r , Daughter.GetDStoPoint(decayVtx));
983
984 Double_t Emc = sqrt(mcTrack.GetP()*mcTrack.GetP() + mcTrack.GetMass()*mcTrack.GetMass());
985 Double_t res[8] = {0},
986 pull[8] = {0},
987 mcParam[8] = { mcTrack.GetStartX(), mcTrack.GetStartY(), mcTrack.GetStartZ(),
988 mcTrack.GetPx(), mcTrack.GetPy(), mcTrack.GetPz(), Emc, mcTrack.GetMass() };
989 for(int iPar=0; iPar < 7; iPar++ )
990 {
991 Double_t error = Daughter.GetCovariance(iPar,iPar);
992 if(error < 0.) { error = 1.e20;}
993 error = TMath::Sqrt(error);
994 res[iPar] = Daughter.GetParameter(iPar) - mcParam[iPar];
995 if(fabs(error) > 1.e-20) pull[iPar] = res[iPar]/error;
996 }
997 res[7] = M - mcParam[7];
998 if(fabs(ErrM) > 1.e-20) pull[7] = res[7]/ErrM;
999
1000 for(int iPar=0; iPar < 8; iPar++ )
1001 {
1002 hFitDaughtersQA[iParticle][iPar]->Fill(res[iPar]);
1003 hFitDaughtersQA[iParticle][iPar+8]->Fill(pull[iPar]);
1004 }
1005 }
1006// for(int iD=0; iD<mcPart.NDaughters(); ++iD)
1007// {
1008// int mcDaughterId = mcPart.GetDaughterIds()[iD];
1009// if(!MCtoRParticleId[mcDaughterId].IsMatchedWithPdg()) continue;
1010// CbmMCTrack &mcTrack = *(L1_DYNAMIC_CAST<CbmMCTrack*>(listMCTracks->At(mcDaughterId)));
1011// int recDaughterId = MCtoRParticleId[mcDaughterId].GetBestMatchWithPdg();
1012//
1013// CbmKFParticleInterface Daughter1;
1014// CbmKFParticle_simd Daughter = vRParticles[recDaughterId];
1015// fvec M,ErrM;
1016// Daughter.GetMass(M,ErrM);
1017//
1018// fvec decayVtx[3] = {mcTrack.GetStartX(), mcTrack.GetStartY(), mcTrack.GetStartZ()};
1019// Daughter1.Extrapolate(&Daughter, Daughter.r , Daughter1.GetDStoPoint(&Daughter,decayVtx));
1020//
1021// fvec Emc = sqrt(mcTrack.GetP()*mcTrack.GetP() + mcTrack.GetMass()*mcTrack.GetMass());
1022// fvec res[8] = {0},
1023// pull[8] = {0},
1024// mcParam[8] = { mcTrack.GetStartX(), mcTrack.GetStartY(), mcTrack.GetStartZ(),
1025// mcTrack.GetPx(), mcTrack.GetPy(), mcTrack.GetPz(), Emc, mcTrack.GetMass() };
1026// for(int iPar=0; iPar < 7; iPar++ )
1027// {
1028// fvec error = Daughter.C[Daughter1.IJ(iPar, iPar)];
1029// if(error[0] < 0.) { error = 1.e20;}
1030// error = TMath::Sqrt(error[0]);
1031// res[iPar] = Daughter.r[iPar] - mcParam[iPar];
1032// if(fabs(error)[0] > 1.e-20) pull[iPar] = res[iPar]/error;
1033// }
1034// res[7] = M - mcParam[7];
1035// if(fabs(ErrM)[0] > 1.e-20) pull[7] = res[7]/ErrM;
1036//
1037// for(int iPar=0; iPar < 8; iPar++ )
1038// {
1039// hFitDaughtersQA[iParticle][iPar]->Fill(res[iPar][0]);
1040// hFitDaughtersQA[iParticle][iPar+8]->Fill(pull[iPar][0]);
1041// }
1042// }
1043 }
1044 //Find MC parameters of the primary vertex
1045 float mcPVx[3]={0.f};
1046 for(unsigned iMC=0; iMC<vMCTracks.size(); ++iMC)
1047 {
1048 if(vMCTracks[iMC].IsPrimary())
1049 {
1050 mcPVx[0]=vMCTracks[iMC].x;
1051 mcPVx[1]=vMCTracks[iMC].y;
1052 mcPVx[2]=vMCTracks[iMC].z;
1053 break;
1054 }
1055 }
1056 float dRPVr[3] = {(float)(PF->GetPV()->GetRefX()-mcPVx[0]),
1057 (float)(PF->GetPV()->GetRefY()-mcPVx[1]),
1058 (float)(PF->GetPV()->GetRefZ()-mcPVx[2])};
1059 float dRPVp[3] = {(float)(dRPVr[0]/sqrt(PF->GetPV()->GetCovMatrix()[0])),
1060 (float)(dRPVr[1]/sqrt(PF->GetPV()->GetCovMatrix()[2])),
1061 (float)(dRPVr[2]/sqrt(PF->GetPV()->GetCovMatrix()[5]))};
1062 for(unsigned int iHPV=0; iHPV<3; ++iHPV)
1063 hPVFitQa[iHPV]->Fill(dRPVr[iHPV]);
1064 for(unsigned int iHPV=3; iHPV<6; ++iHPV)
1065 hPVFitQa[iHPV]->Fill(dRPVp[iHPV-3]);
1066
1067} // void CbmL1::ParticlesEfficienciesPerformance()
1068
1069void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRef on match data in CbmL1**Track classes
1070{
1071
1072 //CbmKF &KF = *CbmKF::Instance();
1073
1074 static TProfile *p_eff_all_vs_mom, *p_eff_prim_vs_mom, *p_eff_sec_vs_mom, *p_eff_d0_vs_mom, *p_eff_prim_vs_theta, *p_eff_all_vs_pt, *p_eff_prim_vs_pt, *p_eff_all_vs_nhits, *p_eff_prim_vs_nhits, *p_eff_sec_vs_nhits;
1075
1076 static TH1F *h_reg_MCmom, *h_acc_MCmom, *h_reco_MCmom, *h_ghost_Rmom;
1077 static TH1F *h_reg_prim_MCmom, *h_acc_prim_MCmom, *h_reco_prim_MCmom;
1078 static TH1F *h_reg_sec_MCmom, *h_acc_sec_MCmom, *h_reco_sec_MCmom;
1079
1080 static TH1F *h_acc_mom_short123s;
1081
1082 static TH1F *h_reg_mom_prim, *h_reg_mom_sec, *h_reg_nhits_prim, *h_reg_nhits_sec;
1083 static TH1F *h_acc_mom_prim, *h_acc_mom_sec, *h_acc_nhits_prim, *h_acc_nhits_sec;
1084 static TH1F *h_reco_mom_prim, *h_reco_mom_sec, *h_reco_nhits_prim, *h_reco_nhits_sec;
1085 static TH1F *h_rest_mom_prim, *h_rest_mom_sec, *h_rest_nhits_prim, *h_rest_nhits_sec;
1086
1087 //static TH1F *h_hit_density[10];
1088
1089 static TH1F
1090 *h_ghost_mom, *h_ghost_nhits, *h_ghost_fstation,
1091 *h_ghost_chi2, *h_ghost_prob, *h_ghost_tx, *h_ghost_ty;
1092 static TH1F *h_reco_mom, *h_reco_d0_mom, *h_reco_nhits, *h_reco_station,
1093 *h_reco_chi2, *h_reco_prob, *h_rest_prob, *h_reco_clean, *h_reco_time;
1094 static TProfile *h_reco_timeNtr, *h_reco_timeNhit ;
1095 static TProfile *h_reco_fakeNtr, *h_reco_fakeNhit ;
1096 static TH1F *h_tx, *h_ty, *h_sec_r, *h_ghost_r;
1097
1098 static TH1F *h_notfound_mom, *h_notfound_nhits, *h_notfound_station, *h_notfound_r, *h_notfound_tx, *h_notfound_ty;
1099
1100 static TH1F *h_mcp, *h_nmchits;//, *h_chi2, *h_prob, *MC_vx, *MC_vy, *MC_vz, *VtxChiPrim, *VtxChiSec;
1101
1102 //static TH2F *P_vs_P ;
1103
1104 static TH2F *h2_vertex, *h2_prim_vertex, *h2_sec_vertex;
1105 //static TH3F *h3_vertex, *h3_prim_vertex, *h3_sec_vertex;
1106
1107 static TH2F *h2_reg_nhits_vs_mom_prim, *h2_reg_nhits_vs_mom_sec,
1108 *h2_reg_fstation_vs_mom_prim, *h2_reg_fstation_vs_mom_sec, *h2_reg_lstation_vs_fstation_prim, *h2_reg_lstation_vs_fstation_sec;
1109 static TH2F *h2_acc_nhits_vs_mom_prim, *h2_acc_nhits_vs_mom_sec,
1110 *h2_acc_fstation_vs_mom_prim, *h2_acc_fstation_vs_mom_sec, *h2_acc_lstation_vs_fstation_prim, *h2_acc_lstation_vs_fstation_sec;
1111 static TH2F *h2_reco_nhits_vs_mom_prim, *h2_reco_nhits_vs_mom_sec,
1112 *h2_reco_fstation_vs_mom_prim, *h2_reco_fstation_vs_mom_sec, *h2_reco_lstation_vs_fstation_prim, *h2_reco_lstation_vs_fstation_sec;
1113 static TH2F *h2_ghost_nhits_vs_mom, *h2_ghost_fstation_vs_mom, *h2_ghost_lstation_vs_fstation;
1114 static TH2F *h2_rest_nhits_vs_mom_prim, *h2_rest_nhits_vs_mom_sec,
1115 *h2_rest_fstation_vs_mom_prim, *h2_rest_fstation_vs_mom_sec, *h2_rest_lstation_vs_fstation_prim, *h2_rest_lstation_vs_fstation_sec;
1116
1117 static bool first_call = 1;
1118
1119 if ( first_call )
1120 {
1121 first_call = 0;
1122
1123 TDirectory *curdir = gDirectory;
1124 gDirectory = histodir;
1125 gDirectory->cd("L1");
1126
1127 p_eff_all_vs_mom = new TProfile("p_eff_all_vs_mom", "AllSet Efficiency vs Momentum",
1128 100, 0.0, 5.0, 0.0, 100.0);
1129 p_eff_prim_vs_mom = new TProfile("p_eff_prim_vs_mom", "Primary Set Efficiency vs Momentum",
1130 100, 0.0, 5.0, 0.0, 100.0);
1131 p_eff_sec_vs_mom = new TProfile("p_eff_sec_vs_mom", "Secondary Set Efficiency vs Momentum",
1132 100, 0.0, 5.0, 0.0, 100.0);
1133 p_eff_d0_vs_mom = new TProfile("p_eff_d0_vs_mom", "D0 Secondary Tracks Efficiency vs Momentum",
1134 150, 0.0, 15.0, 0.0, 100.0);
1135 p_eff_prim_vs_theta = new TProfile("p_eff_prim_vs_theta", "All Primary Set Efficiency vs Theta",
1136 100, 0.0, 30.0, 0.0, 100.0);
1137 p_eff_all_vs_pt = new TProfile("p_eff_all_vs_pt", "AllSet Efficiency vs Pt",
1138 100, 0.0, 5.0, 0.0, 100.0);
1139 p_eff_prim_vs_pt = new TProfile("p_eff_prim_vs_pt", "Primary Set Efficiency vs Pt",
1140 100, 0.0, 5.0, 0.0, 100.0);
1141
1142 p_eff_all_vs_nhits = new TProfile("p_eff_all_vs_nhits", "AllSet Efficiency vs NMCHits",
1143 8, 3.0, 11.0, 0.0, 100.0);
1144 p_eff_prim_vs_nhits = new TProfile("p_eff_prim_vs_nhits", "PrimSet Efficiency vs NMCHits",
1145 8, 3.0, 11.0, 0.0, 100.0);
1146 p_eff_sec_vs_nhits = new TProfile("p_eff_sec_vs_nhits", "SecSet Efficiency vs NMCHits",
1147 8, 3.0, 11.0, 0.0, 100.0);
1148
1149 h_reg_MCmom = new TH1F("h_reg_MCmom", "Momentum of registered tracks", 100, 0.0, 5.0);
1150 h_acc_MCmom = new TH1F("h_acc_MCmom", "Reconstructable tracks", 100, 0.0, 5.0);
1151 h_reco_MCmom = new TH1F("h_reco_MCmom", "Reconstructed tracks", 100, 0.0, 5.0);
1152 h_ghost_Rmom = new TH1F("h_ghost_Rmom", "Ghost tracks", 100, 0.0, 5.0);
1153 h_reg_prim_MCmom = new TH1F("h_reg_prim_MCmom", "Momentum of registered tracks", 100, 0.0, 5.0);
1154 h_acc_prim_MCmom = new TH1F("h_acc_prim_MCmom", "Reconstructable tracks", 100, 0.0, 5.0);
1155 h_reco_prim_MCmom = new TH1F("h_reco_prim_MCmom", "Reconstructed tracks", 100, 0.0, 5.0);
1156 h_reg_sec_MCmom = new TH1F("h_reg_sec_MCmom", "Momentum of registered tracks", 100, 0.0, 5.0);
1157 h_acc_sec_MCmom = new TH1F("h_acc_sec_MCmom", "Reconstructable tracks", 100, 0.0, 5.0);
1158 h_reco_sec_MCmom = new TH1F("h_reco_sec_MCmom", "Reconstructed tracks", 100, 0.0, 5.0);
1159
1160 h_acc_mom_short123s = new TH1F("h_acc_mom_short123s", "Momentum of accepted tracks with 3 hits on first stations", 500, 0.0, 5.0);
1161
1162 h_reg_mom_prim = new TH1F("h_reg_mom_prim", "Momentum of registered primary tracks", 500, 0.0, 5.0);
1163 h_reg_mom_sec = new TH1F("h_reg_mom_sec", "Momentum of registered secondary tracks", 500, 0.0, 5.0);
1164 h_acc_mom_prim = new TH1F("h_acc_mom_prim", "Momentum of accepted primary tracks", 500, 0.0, 5.0);
1165 h_acc_mom_sec = new TH1F("h_acc_mom_sec", "Momentum of accepted secondary tracks", 500, 0.0, 5.0);
1166 h_reco_mom_prim = new TH1F("h_reco_mom_prim", "Momentum of reconstructed primary tracks", 500, 0.0, 5.0);
1167 h_reco_mom_sec = new TH1F("h_reco_mom_sec", "Momentum of reconstructed secondary tracks", 500, 0.0, 5.0);
1168 h_rest_mom_prim = new TH1F("h_rest_mom_prim", "Momentum of not found primary tracks", 500, 0.0, 5.0);
1169 h_rest_mom_sec = new TH1F("h_rest_mom_sec", "Momentum of not found secondary tracks", 500, 0.0, 5.0);
1170
1171 h_reg_nhits_prim = new TH1F("h_reg_nhits_prim", "Number of hits in registered primary track", 51, -0.1, 10.1);
1172 h_reg_nhits_sec = new TH1F("h_reg_nhits_sec", "Number of hits in registered secondary track", 51, -0.1, 10.1);
1173 h_acc_nhits_prim = new TH1F("h_acc_nhits_prim", "Number of hits in accepted primary track", 51, -0.1, 10.1);
1174 h_acc_nhits_sec = new TH1F("h_acc_nhits_sec", "Number of hits in accepted secondary track", 51, -0.1, 10.1);
1175 h_reco_nhits_prim = new TH1F("h_reco_nhits_prim", "Number of hits in reconstructed primary track", 51, -0.1, 10.1);
1176 h_reco_nhits_sec = new TH1F("h_reco_nhits_sec", "Number of hits in reconstructed secondary track", 51, -0.1, 10.1);
1177 h_rest_nhits_prim = new TH1F("h_rest_nhits_prim", "Number of hits in not found primary track", 51, -0.1, 10.1);
1178 h_rest_nhits_sec = new TH1F("h_rest_nhits_sec", "Number of hits in not found secondary track", 51, -0.1, 10.1);
1179
1180 h_ghost_mom = new TH1F("h_ghost_mom", "Momentum of ghost tracks", 500, 0.0, 5.0);
1181 h_ghost_nhits = new TH1F("h_ghost_nhits", "Number of hits in ghost track", 51, -0.1, 10.1);
1182 h_ghost_fstation = new TH1F("h_ghost_fstation", "First station of ghost track", 50, -0.5, 10.0);
1183 h_ghost_chi2 = new TH1F("h_ghost_chi2", "Chi2/NDF of ghost track", 50, -0.5, 10.0);
1184 h_ghost_prob = new TH1F("h_ghost_prob", "Prob of ghost track", 505, 0., 1.01);
1185 h_ghost_r = new TH1F("h_ghost_r", "R of ghost track at the first hit", 50, 0.0, 15.0);
1186 h_ghost_tx = new TH1F("h_ghost_tx", "tx of ghost track at the first hit", 50, -5.0, 5.0);
1187 h_ghost_ty = new TH1F("h_ghost_ty", "ty of ghost track at the first hit", 50, -1.0, 1.0);
1188
1189 h_reco_mom = new TH1F("h_reco_mom", "Momentum of reco track", 50, 0.0, 5.0);
1190 h_reco_nhits = new TH1F("h_reco_nhits", "Number of hits in reco track", 50, 0.0, 10.0);
1191 h_reco_station = new TH1F("h_reco_station", "First station of reco track", 50, -0.5, 10.0);
1192 h_reco_chi2 = new TH1F("h_reco_chi2", "Chi2/NDF of reco track", 50, -0.5, 10.0);
1193 h_reco_prob = new TH1F("h_reco_prob", "Prob of reco track", 505, 0., 1.01);
1194 h_rest_prob = new TH1F("h_rest_prob", "Prob of reco rest track", 505, 0., 1.01);
1195 h_reco_clean = new TH1F("h_reco_clean", "Percentage of correct hits", 100, -0.5, 100.5);
1196 h_reco_time = new TH1F("h_reco_time", "CA Track Finder Time (s/ev)", 20, 0.0, 20.0);
1197 h_reco_timeNtr = new TProfile("h_reco_timeNtr",
1198 "CA Track Finder Time (s/ev) vs N Tracks",
1199 200, 0.0, 1000.0);
1200 h_reco_timeNhit = new TProfile("h_reco_timeNhit",
1201 "CA Track Finder Time (s/ev) vs N Hits",
1202 200, 0.0, 30000.0);
1203 h_reco_fakeNtr = new TProfile("h_reco_fakeNtr",
1204 "N Fake hits vs N Tracks",
1205 200, 0.0, 1000.0);
1206 h_reco_fakeNhit = new TProfile("h_reco_fakeNhit",
1207 "N Fake hits vs N Hits",
1208 200, 0.0, 30000.0);
1209
1210 h_reco_d0_mom = new TH1F("h_reco_d0_mom", "Momentum of reco D0 track", 150, 0.0, 15.0);
1211
1212// h_hit_density[0] = new TH1F("h_hit_density0", "Hit density station 1", 50, 0.0, 5.00);
1213// h_hit_density[1] = new TH1F("h_hit_density1", "Hit density station 2", 100, 0.0, 10.00);
1214// h_hit_density[2] = new TH1F("h_hit_density2", "Hit density station 3", 140, 0.0, 13.99);
1215// h_hit_density[3] = new TH1F("h_hit_density3", "Hit density station 4", 180, 0.0, 18.65);
1216// h_hit_density[4] = new TH1F("h_hit_density4", "Hit density station 5", 230, 0.0, 23.32);
1217// h_hit_density[5] = new TH1F("h_hit_density5", "Hit density station 6", 270, 0.0, 27.98);
1218// h_hit_density[6] = new TH1F("h_hit_density6", "Hit density station 7", 340, 0.0, 34.97);
1219// h_hit_density[7] = new TH1F("h_hit_density7", "Hit density station 8", 460, 0.0, 46.63);
1220// h_hit_density[8] = new TH1F("h_hit_density8", "Hit density station 9", 500, 0.0, 50.00);
1221// h_hit_density[9] = new TH1F("h_hit_density8", "Hit density station 9", 500, 0.0, 50.00);
1222// h_hit_density[10] = new TH1F("h_hit_density8", "Hit density station 9", 500, 0.0, 50.00);
1223
1224 h_tx = new TH1F("h_tx", "tx of MC track at the vertex", 50, -0.5, 0.5);
1225 h_ty = new TH1F("h_ty", "ty of MC track at the vertex", 50, -0.5, 0.5);
1226
1227 h_sec_r = new TH1F("h_sec_r", "R of sec MC track at the first hit", 50, 0.0, 15.0);
1228
1229 h_notfound_mom = new TH1F("h_notfound_mom", "Momentum of not found track", 50, 0.0, 5.0);
1230 h_notfound_nhits = new TH1F("h_notfound_nhits", "Num hits of not found track", 50, 0.0, 10.0);
1231 h_notfound_station = new TH1F("h_notfound_station", "First station of not found track", 50, 0.0, 10.0);
1232 h_notfound_r = new TH1F("h_notfound_r", "R at first hit of not found track", 50, 0.0, 15.0);
1233 h_notfound_tx = new TH1F("h_notfound_tx", "tx of not found track at the first hit", 50, -5.0, 5.0);
1234 h_notfound_ty = new TH1F("h_notfound_ty", "ty of not found track at the first hit", 50, -1.0, 1.0);
1235
1236 //h_chi2 = new TH1F("chi2", "Chi^2", 100, 0.0, 10.);
1237 //h_prob = new TH1F("prob", "Prob", 100, 0.0, 1.01);
1238 //VtxChiPrim = new TH1F("vtxChiPrim", "Chi to primary vtx for primary tracks", 100, 0.0, 10.);
1239 //VtxChiSec = new TH1F("vtxChiSec", "Chi to primary vtx for secondary tracks", 100, 0.0, 10.);
1240 h_mcp = new TH1F("h_mcp", "MC P ", 500, 0.0, 5.0);
1241 //MC_vx = new TH1F("MC_vx", "MC Vertex X", 100, -.05, .05);
1242 //MC_vy = new TH1F("MC_vy", "MC Vertex Y", 100, -.05, .05);
1243 //MC_vz = new TH1F("MC_vz", "MC Vertex Z", 100, -.05, .05);
1244 h_nmchits = new TH1F("h_nmchits", "N STS hits on MC track", 50, 0.0, 10.0);
1245
1246 //P_vs_P = new TH2F("P_vs_P", "Resolution P/Q vs P", 20, 0., 20.,100, -.05, .05);
1247
1248 h2_vertex = new TH2F("h2_vertex", "2D vertex distribution", 105, -5., 100., 100, -50., 50.);
1249 h2_prim_vertex = new TH2F("h2_primvertex", "2D primary vertex distribution", 105, -5., 100., 100, -50., 50.);
1250 h2_sec_vertex = new TH2F("h2_sec_vertex", "2D secondary vertex distribution", 105, -5., 100., 100, -50., 50.);
1251
1252 //h3_vertex = new TH3F("h3_vertex", "3D vertex distribution", 20, -5., 100., 100, -50., 50., 100, -50., 50.);
1253 //h3_prim_vertex = new TH3F("h3_primvertex", "3D vertex distribution", 20, -5., 100., 100, -50., 50., 100, -50., 50.);
1254 //h3_sec_vertex = new TH3F("h3_sec_vertex", "3D vertex distribution", 20, -5., 100., 100, -50., 50., 100, -50., 50.);
1255
1256 h2_reg_nhits_vs_mom_prim = new TH2F("h2_reg_nhits_vs_mom_prim", "NHits vs. Momentum (Reg. Primary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1257 h2_reg_nhits_vs_mom_sec = new TH2F("h2_reg_nhits_vs_mom_sec", "NHits vs. Momentum (Reg. Secondary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1258 h2_acc_nhits_vs_mom_prim = new TH2F("h2_acc_nhits_vs_mom_prim", "NHits vs. Momentum (Acc. Primary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1259 h2_acc_nhits_vs_mom_sec = new TH2F("h2_acc_nhits_vs_mom_sec", "NHits vs. Momentum (Acc. Secondary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1260 h2_reco_nhits_vs_mom_prim = new TH2F("h2_reco_nhits_vs_mom_prim", "NHits vs. Momentum (Reco Primary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1261 h2_reco_nhits_vs_mom_sec = new TH2F("h2_reco_nhits_vs_mom_sec", "NHits vs. Momentum (Reco Secondary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1262 h2_ghost_nhits_vs_mom = new TH2F("h2_ghost_nhits_vs_mom", "NHits vs. Momentum (Ghost Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1263 h2_rest_nhits_vs_mom_prim = new TH2F("h2_rest_nhits_vs_mom_prim", "NHits vs. Momentum (!Found Primary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1264 h2_rest_nhits_vs_mom_sec = new TH2F("h2_rest_nhits_vs_mom_sec", "NHits vs. Momentum (!Found Secondary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1265
1266 h2_reg_fstation_vs_mom_prim = new TH2F("h2_reg_fstation_vs_mom_prim", "First Station vs. Momentum (Reg. Primary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1267 h2_reg_fstation_vs_mom_sec = new TH2F("h2_reg_fstation_vs_mom_sec", "First Station vs. Momentum (Reg. Secondary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1268 h2_acc_fstation_vs_mom_prim = new TH2F("h2_acc_fstation_vs_mom_prim", "First Station vs. Momentum (Acc. Primary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1269 h2_acc_fstation_vs_mom_sec = new TH2F("h2_acc_fstation_vs_mom_sec", "First Station vs. Momentum (Acc. Secondary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1270 h2_reco_fstation_vs_mom_prim = new TH2F("h2_reco_fstation_vs_mom_prim", "First Station vs. Momentum (Reco Primary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1271 h2_reco_fstation_vs_mom_sec = new TH2F("h2_reco_fstation_vs_mom_sec", "First Station vs. Momentum (Reco Secondary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1272 h2_ghost_fstation_vs_mom = new TH2F("h2_ghost_fstation_vs_mom", "First Station vs. Momentum (Ghost Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1273 h2_rest_fstation_vs_mom_prim = new TH2F("h2_rest_fstation_vs_mom_prim", "First Station vs. Momentum (!Found Primary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1274 h2_rest_fstation_vs_mom_sec = new TH2F("h2_rest_fstation_vs_mom_sec", "First Station vs. Momentum (!Found Secondary Tracks)", 51, -0.05, 5.05, 11, -0.5, 10.5);
1275
1276 h2_reg_lstation_vs_fstation_prim = new TH2F("h2_reg_lstation_vs_fstation_prim", "Last vs. First Station (Reg. Primary Tracks)", 11, -0.5, 10.5, 11, -0.5, 10.5);
1277 h2_reg_lstation_vs_fstation_sec = new TH2F("h2_reg_lstation_vs_fstation_sec", "Last vs. First Station (Reg. Secondary Tracks)", 11, -0.5, 10.5, 11, -0.5, 10.5);
1278 h2_acc_lstation_vs_fstation_prim = new TH2F("h2_acc_lstation_vs_fstation_prim", "Last vs. First Station (Acc. Primary Tracks)", 11, -0.5, 10.5, 11, -0.5, 10.5);
1279 h2_acc_lstation_vs_fstation_sec = new TH2F("h2_acc_lstation_vs_fstation_sec", "Last vs. First Station (Acc. Secondary Tracks)", 11, -0.5, 10.5, 11, -0.5, 10.5);
1280 h2_reco_lstation_vs_fstation_prim = new TH2F("h2_reco_lstation_vs_fstation_prim", "Last vs. First Station (Reco Primary Tracks)", 11, -0.5, 10.5, 11, -0.5, 10.5);
1281 h2_reco_lstation_vs_fstation_sec = new TH2F("h2_reco_lstation_vs_fstation_sec", "Last vs. First Station (Reco Secondary Tracks)", 11, -0.5, 10.5, 11, -0.5, 10.5);
1282 h2_ghost_lstation_vs_fstation = new TH2F("h2_ghost_lstation_vs_fstation", "Last vs. First Station (Ghost Tracks)", 11, -0.5, 10.5, 11, -0.5, 10.5);
1283 h2_rest_lstation_vs_fstation_prim = new TH2F("h2_rest_lstation_vs_fstation_prim", "Last vs. First Station (!Found Primary Tracks)", 11, -0.5, 10.5, 11, -0.5, 10.5);
1284 h2_rest_lstation_vs_fstation_sec = new TH2F("h2_rest_lstation_vs_fstation_sec", "Last vs. First Station (!Found Secondary Tracks)", 11, -0.5, 10.5, 11, -0.5, 10.5);
1285
1286 //maindir->cd();
1287
1288 // ----- Create list of all histogram pointers
1289
1290 gDirectory = curdir;
1291
1292 }// first_call
1293
1294
1295 static int NEvents = 0;
1296 if ( NEvents > 0 ) {
1297 h_reg_MCmom->Scale(NEvents);
1298 h_acc_MCmom->Scale(NEvents);
1299 h_reco_MCmom->Scale(NEvents);
1300 h_ghost_Rmom->Scale(NEvents);
1301 h_reg_prim_MCmom->Scale(NEvents);
1302 h_acc_prim_MCmom->Scale(NEvents);
1303 h_reco_prim_MCmom->Scale(NEvents);
1304 h_reg_sec_MCmom->Scale(NEvents);
1305 h_acc_sec_MCmom->Scale(NEvents);
1306 h_reco_sec_MCmom->Scale(NEvents);
1307 }
1308 NEvents++;
1309
1310// //hit density calculation: h_hit_density[10]
1311 //
1312// for (vector<CbmL1HitStore>::iterator hIt = vHitStore.begin(); hIt != vHitStore.end(); ++hIt){
1313// float x = hIt->x;
1314// float y = hIt->y;
1315// float r = sqrt(x*x+y*y);
1316// h_hit_density[hIt->iStation]->Fill(r, 1.0/(2.0*3.1415*r));
1317// }
1318
1319
1320 //
1321 for (vector<CbmL1Track>::iterator rtraIt = vRTracks.begin(); rtraIt != vRTracks.end(); ++rtraIt){
1322 CbmL1Track* prtra = &(*rtraIt);
1323 if((prtra->StsHits).size() < 1) continue;
1324 { // fill histos
1325 if( fabs(prtra->T[4])>1.e-10 ) h_reco_mom->Fill(fabs(1.0/prtra->T[4]));
1326 h_reco_nhits->Fill((prtra->StsHits).size());
1327 CbmL1HitStore &mh = vHitStore[prtra->StsHits[0]];
1328 h_reco_station->Fill(mh.iStation);
1329
1330 }
1331
1332 h_reco_clean->Fill( prtra->GetMaxPurity() );
1333
1334 if (prtra->NDF > 0) {
1335 if ( prtra->IsGhost() ){
1336 h_ghost_chi2->Fill(prtra->chi2/prtra->NDF);
1337 h_ghost_prob->Fill(TMath::Prob(prtra->chi2,prtra->NDF));
1338 }
1339 else {
1340 if (prtra->GetMCTrack()[0].IsReconstructable()){
1341 h_reco_chi2->Fill(prtra->chi2/prtra->NDF);
1342 h_reco_prob->Fill(TMath::Prob(prtra->chi2,prtra->NDF));
1343 } else {
1344// h_rest_chi2->Fill(prtra->chi2/prtra->NDF);
1345 h_rest_prob->Fill(TMath::Prob(prtra->chi2,prtra->NDF));
1346 }
1347 }
1348 }
1349
1350
1351 // fill ghost histos
1352 if ( prtra->IsGhost() ){
1353 if( fabs(prtra->T[4])>1.e-10) {
1354 h_ghost_mom->Fill(fabs(1.0/prtra->T[4]));
1355 h_ghost_Rmom->Fill(fabs(1.0/prtra->T[4]));
1356 }
1357 h_ghost_nhits->Fill((prtra->StsHits).size());
1358 CbmL1HitStore &h1 = vHitStore[prtra->StsHits[0]];
1359 CbmL1HitStore &h2 = vHitStore[prtra->StsHits[1]];
1360 h_ghost_fstation->Fill(h1.iStation);
1361 h_ghost_r->Fill(sqrt(fabs(h1.x*h1.x+h1.y*h1.y)));
1362 double z1 = algo->vStations[h1.iStation].z[0];
1363 double z2 = algo->vStations[h2.iStation].z[0];
1364 if( fabs(z2-z1)>1.e-4 ){
1365 h_ghost_tx->Fill((h2.x-h1.x)/(z2-z1));
1366 h_ghost_ty->Fill((h2.y-h1.y)/(z2-z1));
1367 }
1368
1369 if( fabs(prtra->T[4])>1.e-10)
1370 h2_ghost_nhits_vs_mom->Fill(fabs(1.0/prtra->T[4]), (prtra->StsHits).size());
1371 CbmL1HitStore &hf = vHitStore[prtra->StsHits[0]];
1372 CbmL1HitStore &hl = vHitStore[prtra->StsHits[(prtra->StsHits).size()-1]];
1373 if( fabs(prtra->T[4])>1.e-10)
1374 h2_ghost_fstation_vs_mom->Fill(fabs(1.0/prtra->T[4]), hf.iStation+1);
1375 if (hl.iStation >= hf.iStation)
1376 h2_ghost_lstation_vs_fstation->Fill(hf.iStation+1, hl.iStation+1);
1377 }
1378
1379 } // for reco tracks
1380
1381 int mc_total = 0; // total amount of reconstructable mcTracks
1382 for ( vector<CbmL1MCTrack>::iterator mtraIt = vMCTracks.begin(); mtraIt != vMCTracks.end(); mtraIt++ ) {
1383 CbmL1MCTrack &mtra = *(mtraIt);
1384// if( !( mtra.pdg == -11 && mtra.mother_ID == -1 ) ) continue; // electrons only
1385
1386 // No Sts hits?
1387 int nmchits = mtra.StsHits.size();
1388 if (nmchits == 0) continue;
1389
1390 double momentum = mtra.p;
1391 double pt = sqrt(mtra.px*mtra.px+mtra.py*mtra.py);
1392 double theta = acos(mtra.pz/mtra.p)*180/3.1415;
1393
1394 h_mcp->Fill(momentum);
1395 h_nmchits->Fill(nmchits);
1396
1397 int nSta = mtra.NStations();
1398
1399 CbmL1HitStore &fh = vHitStore[*(mtra.StsHits.begin())];
1400 CbmL1HitStore &lh = vHitStore[*(mtra.StsHits.rbegin())];
1401
1402 h_reg_MCmom->Fill(momentum);
1403 if ( mtra.IsPrimary() ){
1404 h_reg_mom_prim->Fill(momentum);
1405 h_reg_prim_MCmom->Fill(momentum);
1406 h_reg_nhits_prim->Fill(nSta);
1407 h2_reg_nhits_vs_mom_prim->Fill(momentum, nSta);
1408 h2_reg_fstation_vs_mom_prim->Fill(momentum, fh.iStation+1);
1409 if (lh.iStation >= fh.iStation)
1410 h2_reg_lstation_vs_fstation_prim->Fill(fh.iStation+1, lh.iStation+1);
1411 }else{
1412 h_reg_mom_sec->Fill(momentum);
1413 h_reg_sec_MCmom->Fill(momentum);
1414 h_reg_nhits_sec->Fill(nSta);
1415 if (momentum > 0.01){
1416 h2_reg_nhits_vs_mom_sec->Fill(momentum, nSta);
1417 h2_reg_fstation_vs_mom_sec->Fill(momentum, fh.iStation+1);
1418 if (lh.iStation >= fh.iStation)
1419 h2_reg_lstation_vs_fstation_sec->Fill(fh.iStation+1, lh.iStation+1);
1420 }
1421 }
1422
1423 if ( mtra.IsAdditional() ){
1424 h_acc_mom_short123s->Fill(momentum);
1425 }
1426
1427 if( ! mtra.IsReconstructable() ) continue;
1428 mc_total++;
1429
1430 h_acc_MCmom->Fill(momentum);
1431 if ( mtra.IsPrimary() ){
1432 h_acc_mom_prim->Fill(momentum);
1433 h_acc_prim_MCmom->Fill(momentum);
1434 h_acc_nhits_prim->Fill(nSta);
1435 h2_acc_nhits_vs_mom_prim->Fill(momentum, nSta);
1436 h2_acc_fstation_vs_mom_prim->Fill(momentum, fh.iStation+1);
1437 if (lh.iStation >= fh.iStation)
1438 h2_acc_lstation_vs_fstation_prim->Fill(fh.iStation+1, lh.iStation+1);
1439 }else{
1440 h_acc_mom_sec->Fill(momentum);
1441 h_acc_sec_MCmom->Fill(momentum);
1442 h_acc_nhits_sec->Fill(nSta);
1443 if (momentum > 0.01){
1444 h2_acc_nhits_vs_mom_sec->Fill(momentum, nSta);
1445 h2_acc_fstation_vs_mom_sec->Fill(momentum, fh.iStation+1);
1446 if (lh.iStation >= fh.iStation)
1447 h2_acc_lstation_vs_fstation_sec->Fill(fh.iStation+1, lh.iStation+1);
1448 }
1449 }
1450
1451
1452 // vertex distribution of primary and secondary tracks
1453 h2_vertex->Fill(mtra.z, mtra.y);
1454 //h3_vertex->Fill(mtra.z, mtra.x, mtra.y);
1455 if (mtra.mother_ID < 0){ // primary
1456 h2_prim_vertex->Fill(mtra.z, mtra.y);
1457 //h3_prim_vertex->Fill(mtra.z, mtra.x, mtra.y);
1458 }else{
1459 h2_sec_vertex->Fill(mtra.z, mtra.y);
1460 //h3_sec_vertex->Fill(mtra.z, mtra.x, mtra.y);
1461 }
1462
1463
1464 int iph = mtra.StsHits[0];
1465 CbmL1HitStore &ph = vHitStore[iph];
1466
1467 h_sec_r->Fill(sqrt(fabs(ph.x*ph.x+ph.y*ph.y)));
1468 if( fabs( mtra.pz )>1.e-8 ){
1469 h_tx->Fill(mtra.px/mtra.pz);
1470 h_ty->Fill(mtra.py/mtra.pz);
1471 }
1472
1473 // reconstructed tracks
1474 bool reco = (mtra.IsReconstructed());
1475 int nMCHits = mtra.NStations();
1476
1477 if (reco){
1478 p_eff_all_vs_mom->Fill(momentum, 100.0);
1479 p_eff_all_vs_nhits->Fill(nMCHits, 100.0);
1480 p_eff_all_vs_pt->Fill(pt, 100.0);
1481 h_reco_MCmom->Fill(momentum);
1482 if (mtra.mother_ID < 0){ // primary
1483 p_eff_prim_vs_mom->Fill(momentum, 100.0);
1484 p_eff_prim_vs_nhits->Fill(nMCHits, 100.0);
1485 p_eff_prim_vs_pt->Fill(pt, 100.0);
1486 p_eff_prim_vs_theta->Fill(theta, 100.0);
1487 }else{
1488 p_eff_sec_vs_mom->Fill(momentum, 100.0);
1489 p_eff_sec_vs_nhits->Fill(nMCHits, 100.0);
1490 }
1491 if (mtra.mother_ID < 0){ // primary
1492 h_reco_mom_prim->Fill(momentum);
1493 h_reco_prim_MCmom->Fill(momentum);
1494 h_reco_nhits_prim->Fill(nSta);
1495 h2_reco_nhits_vs_mom_prim->Fill(momentum, nSta);
1496 h2_reco_fstation_vs_mom_prim->Fill(momentum, fh.iStation+1);
1497 if (lh.iStation >= fh.iStation)
1498 h2_reco_lstation_vs_fstation_prim->Fill(fh.iStation+1, lh.iStation+1);
1499 }else{
1500 h_reco_mom_sec->Fill(momentum);
1501 h_reco_sec_MCmom->Fill(momentum);
1502 h_reco_nhits_sec->Fill(nSta);
1503 if (momentum > 0.01){
1504 h2_reco_nhits_vs_mom_sec->Fill(momentum, nSta);
1505 h2_reco_fstation_vs_mom_sec->Fill(momentum, fh.iStation+1);
1506 if (lh.iStation >= fh.iStation)
1507 h2_reco_lstation_vs_fstation_sec->Fill(fh.iStation+1, lh.iStation+1);
1508 }
1509 }
1510 }else{
1511 h_notfound_mom->Fill(momentum);
1512 p_eff_all_vs_mom->Fill(momentum, 0.0);
1513 p_eff_all_vs_nhits->Fill(nMCHits, 0.0);
1514 p_eff_all_vs_pt->Fill(pt, 0.0);
1515 if (mtra.mother_ID < 0){ // primary
1516 p_eff_prim_vs_mom->Fill(momentum, 0.0);
1517 p_eff_prim_vs_nhits->Fill(nMCHits, 0.0);
1518 p_eff_prim_vs_pt->Fill(pt, 0.0);
1519 p_eff_prim_vs_theta->Fill(theta, 0.0);
1520 }else{
1521 p_eff_sec_vs_mom->Fill(momentum, 0.0);
1522 p_eff_sec_vs_nhits->Fill(nMCHits, 0.0);
1523 }
1524 if (mtra.mother_ID < 0){ // primary
1525 h_rest_mom_prim->Fill(momentum);
1526 h_rest_nhits_prim->Fill(nSta);
1527 h2_rest_nhits_vs_mom_prim->Fill(momentum, nSta);
1528 h2_rest_fstation_vs_mom_prim->Fill(momentum, fh.iStation+1);
1529 if (lh.iStation >= fh.iStation)
1530 h2_rest_lstation_vs_fstation_prim->Fill(fh.iStation+1, lh.iStation+1);
1531 }else{
1532 h_rest_mom_sec->Fill(momentum);
1533 h_rest_nhits_sec->Fill(nSta);
1534 if (momentum > 0.01){
1535 h2_rest_nhits_vs_mom_sec->Fill(momentum, nSta);
1536 h2_rest_fstation_vs_mom_sec->Fill(momentum, fh.iStation+1);
1537 if (lh.iStation >= fh.iStation)
1538 h2_rest_lstation_vs_fstation_sec->Fill(fh.iStation+1, lh.iStation+1);
1539 }
1540 }
1541 }
1542
1543 // killed tracks. At least one hit of they belong to some recoTrack
1544 //bool killed = 0;
1545 if(!reco){
1546 h_notfound_nhits->Fill(nmchits);
1547 h_notfound_station->Fill(ph.iStation);
1548 h_notfound_r->Fill(sqrt(fabs(ph.x*ph.x+ph.y*ph.y)));
1549
1550 if(mtra.Points.size() > 0){
1551 CbmL1MCPoint &pMC = vMCPoints[mtra.Points[0]];
1552 h_notfound_tx->Fill(pMC.px/pMC.pz);
1553 h_notfound_ty->Fill(pMC.py/pMC.pz);
1554 }
1555
1556// CbmL1HitStore &ph21 = vHitStore[mtra.StsHits[0]];
1557// CbmL1HitStore &ph22 = vHitStore[mtra.StsHits[1]];
1558
1559// double z21 = algo->vStations[ph21.iStation].z[0];
1560// double z22 = algo->vStations[ph22.iStation].z[0];
1561// if( fabs(z22-z21)>1.e-4 ){
1562// h_notfound_tx->Fill((ph22.x-ph21.x)/(z22-z21));
1563// h_notfound_ty->Fill((ph22.y-ph21.y)/(z22-z21));
1564// }
1565
1566 //if( mtra.IsDisturbed() ) killed = 1;
1567 }
1568
1569
1570
1571
1572 if (( mtra.IsPrimary() )&&(mtra.z > 0)){ // D0
1573 h_reco_d0_mom->Fill(momentum);
1574 if (reco) p_eff_d0_vs_mom->Fill(momentum, 100.0);
1575 else p_eff_d0_vs_mom->Fill(momentum, 0.0);
1576 }
1577
1578 } // for mcTracks
1579
1580 int NFakes = 0;
1581 for( unsigned int ih=0; ih<algo->vStsHits.size(); ih++){
1582 int iMC = vHitMCRef[ih]; // TODO2: adapt to linking
1583 if (iMC < 0) NFakes++;
1584 }
1585
1586 h_reco_time->Fill(algo->CATime);
1587 h_reco_timeNtr->Fill(mc_total,algo->CATime);
1588 h_reco_timeNhit->Fill(algo->vStsHits.size(),algo->CATime);
1589
1590 h_reco_fakeNtr->Fill(mc_total,NFakes);
1591 h_reco_fakeNhit->Fill(algo->vStsHits.size()-NFakes,NFakes);
1592
1593
1594 h_reg_MCmom->Scale(1.f/NEvents);
1595 h_acc_MCmom->Scale(1.f/NEvents);
1596 h_reco_MCmom->Scale(1.f/NEvents);
1597 h_ghost_Rmom->Scale(1.f/NEvents);
1598 h_reg_prim_MCmom->Scale(1.f/NEvents);
1599 h_acc_prim_MCmom->Scale(1.f/NEvents);
1600 h_reco_prim_MCmom->Scale(1.f/NEvents);
1601 h_reg_sec_MCmom->Scale(1.f/NEvents);
1602 h_acc_sec_MCmom->Scale(1.f/NEvents);
1603 h_reco_sec_MCmom->Scale(1.f/NEvents);
1604
1605} // void CbmL1::HistoPerformance()
1606
1607
1608void CbmL1::TrackFitPerformance()
1609{
1610 const int Nh_fit = 12;
1611 static TH1F *h_fit[Nh_fit], *h_fitL[Nh_fit], *h_fitSV[Nh_fit], *h_fitPV[Nh_fit], *h_fit_chi2;
1612
1613 static TH2F *PRes2D, *PRes2DPrim, *PRes2DSec;
1614
1615 static bool first_call = 1;
1616
1617 if ( first_call )
1618 {
1619 first_call = 0;
1620
1621 TDirectory *currentDir = gDirectory;
1622 gDirectory = histodir;
1623 gDirectory->cd("L1");
1624 gDirectory->mkdir("Fit");
1625 gDirectory->cd("Fit");
1626 {
1627 PRes2D = new TH2F("PRes2D", "Resolution P vs P", 100, 0., 20.,100, -15., 15.);
1628 PRes2DPrim = new TH2F("PRes2DPrim", "Resolution P vs P", 100, 0., 20.,100, -15., 15.);
1629 PRes2DSec = new TH2F("PRes2DSec", "Resolution P vs P", 100, 0., 20.,100, -15., 15.);
1630
1631 struct {
1632 const char *name;
1633 const char *title;
1634 Int_t n;
1635 Double_t l,r;
1636 } Table[Nh_fit]=
1637 {
1638 {"x", "Residual X [#mum]", 100, -45., 45.},
1639 {"y", "Residual Y [#mum]", 100, -450., 450.},
1640 {"tx", "Residual Tx [mrad]", 100, -2., 2.},
1641 {"ty", "Residual Ty [mrad]", 100, -3.5, 3.5},
1642 {"P", "Resolution P/Q [100%]", 100, -0.1, 0.1 },
1643 {"px", "Pull X [residual/estimated_error]", 100, -6., 6.},
1644 {"py", "Pull Y [residual/estimated_error]", 100, -6., 6.},
1645 {"ptx","Pull Tx [residual/estimated_error]", 100, -6., 6.},
1646 {"pty","Pull Ty [residual/estimated_error]", 100, -6., 6.},
1647 {"pQP","Pull Q/P [residual/estimated_error]", 100, -6., 6.},
1648 {"QPreco","Reco Q/P ", 100, -10., 10.},
1649 {"QPmc","MC Q/P ", 100, -10., 10.}
1650 };
1651
1652 struct Tab{
1653 const char *name;
1654 const char *title;
1655 Int_t n;
1656 Double_t l,r;
1657 };
1658 Tab TableVertex[Nh_fit]=
1659 {
1660 {"x", "Residual X [cm]", 2000, -1., 1.},
1661 {"y", "Residual Y [cm]", 2000, -1., 1.},
1662 {"tx", "Residual Tx [mrad]", 100, -2., 2.},
1663 {"ty", "Residual Ty [mrad]", 100, -2., 2.},
1664 {"P", "Resolution P/Q [100%]", 100, -0.1, 0.1 },
1665 {"px", "Pull X [residual/estimated_error]", 100, -6., 6.},
1666 {"py", "Pull Y [residual/estimated_error]", 100, -6., 6.},
1667 {"ptx","Pull Tx [residual/estimated_error]", 100, -6., 6.},
1668 {"pty","Pull Ty [residual/estimated_error]", 100, -6., 6.},
1669 {"pQP","Pull Q/P [residual/estimated_error]", 100, -6., 6.},
1670 {"QPreco","Reco Q/P ", 100, -10., 10.},
1671 {"QPmc","MC Q/P ", 100, -10., 10.}
1672 };
1673
1674 for( int i=0; i<Nh_fit; i++ ){
1675 char n[225], t[255];
1676 sprintf(n,"fst_%s",Table[i].name);
1677 sprintf(t,"First point %s",Table[i].title);
1678 h_fit[i] = new TH1F(n,t, Table[i].n, Table[i].l, Table[i].r);
1679 sprintf(n,"lst_%s",Table[i].name);
1680 sprintf(t,"Last point %s",Table[i].title);
1681 h_fitL[i] = new TH1F(n,t, Table[i].n, Table[i].l, Table[i].r);
1682 sprintf(n,"svrt_%s",TableVertex[i].name);
1683 sprintf(t,"Secondary vertex point %s",TableVertex[i].title);
1684 h_fitSV[i] = new TH1F(n,t, TableVertex[i].n, TableVertex[i].l, TableVertex[i].r);
1685 sprintf(n,"pvrt_%s",TableVertex[i].name);
1686 sprintf(t,"Primary vertex point %s",TableVertex[i].title);
1687 h_fitPV[i] = new TH1F(n,t, TableVertex[i].n, TableVertex[i].l, TableVertex[i].r);
1688 }
1689 h_fit_chi2 = new TH1F("h_fit_chi2", "Chi2/NDF", 50, -0.5, 10.0);
1690 }
1691 gDirectory = currentDir;
1692 } // if first call
1693
1694 for (vector<CbmL1Track>::iterator it = vRTracks.begin(); it != vRTracks.end(); ++it){
1695
1696 if ( it->IsGhost() ) continue;
1697 { // first hit
1698 int iMC = vHitMCRef[it->StsHits.front()]; // TODO2: adapt to linking
1699 if (iMC < 0) continue;
1700 CbmL1MCPoint &mc = vMCPoints[iMC];
1701// if( !( mc.pdg == -11 && mc.mother_ID == -1 ) ) continue; // electrons only
1702 h_fit[0]->Fill( (mc.x-it->T[0]) *1.e4);
1703 h_fit[1]->Fill( (mc.y-it->T[1]) *1.e4);
1704 h_fit[2]->Fill((mc.px/mc.pz-it->T[2])*1.e3);
1705 h_fit[3]->Fill((mc.py/mc.pz-it->T[3])*1.e3);
1706 h_fit[4]->Fill(it->T[4]/mc.q*mc.p-1);
1707
1708 PRes2D->Fill( mc.p, (1./fabs(it->T[4]) - mc.p)/mc.p*100. );
1709
1710 CbmL1MCTrack mcTrack = *(it->GetMCTracks()[0]);
1711 if(mcTrack.IsPrimary())
1712 PRes2DPrim->Fill( mc.p, (1./fabs(it->T[4]) - mc.p)/mc.p*100. );
1713 else
1714 PRes2DSec->Fill( mc.p, (1./fabs(it->T[4]) - mc.p)/mc.p*100. );
1715
1716 if( finite(it->C[0]) && it->C[0]>0 )h_fit[5]->Fill( (mc.x-it->T[0])/sqrt(it->C[0]));
1717 if( finite(it->C[2]) && it->C[2]>0 )h_fit[6]->Fill( (mc.y-it->T[1])/sqrt(it->C[2]));
1718 if( finite(it->C[5]) && it->C[5]>0 )h_fit[7]->Fill( (mc.px/mc.pz-it->T[2])/sqrt(it->C[5]));
1719 if( finite(it->C[9]) && it->C[9]>0 )h_fit[8]->Fill( (mc.py/mc.pz-it->T[3])/sqrt(it->C[9]));
1720 if( finite(it->C[14]) && it->C[14]>0 )h_fit[9]->Fill( (mc.q/mc.p-it->T[4])/sqrt(it->C[14]));
1721 h_fit[10]->Fill(it->T[4]);
1722 h_fit[11]->Fill(mc.q/mc.p);
1723 }
1724
1725 { // last hit
1726 int iMC = vHitMCRef[it->StsHits.back()]; // TODO2: adapt to linking
1727 if (iMC < 0) continue;
1728 CbmL1MCPoint &mc = vMCPoints[iMC];
1729 h_fitL[0]->Fill( (mc.x-it->TLast[0]) *1.e4);
1730 h_fitL[1]->Fill( (mc.y-it->TLast[1]) *1.e4);
1731 h_fitL[2]->Fill((mc.px/mc.pz-it->TLast[2])*1.e3);
1732 h_fitL[3]->Fill((mc.py/mc.pz-it->TLast[3])*1.e3);
1733 h_fitL[4]->Fill(it->T[4]/mc.q*mc.p-1);
1734 if( finite(it->CLast[0]) && it->CLast[0]>0 ) h_fitL[5]->Fill( (mc.x-it->TLast[0])/sqrt(it->CLast[0]));
1735 if( finite(it->CLast[2]) && it->CLast[2]>0 ) h_fitL[6]->Fill( (mc.y-it->TLast[1])/sqrt(it->CLast[2]));
1736 if( finite(it->CLast[5]) && it->CLast[5]>0 ) h_fitL[7]->Fill( (mc.px/mc.pz-it->TLast[2])/sqrt(it->CLast[5]));
1737 if( finite(it->CLast[9]) && it->CLast[9]>0 ) h_fitL[8]->Fill( (mc.py/mc.pz-it->TLast[3])/sqrt(it->CLast[9]));
1738 if( finite(it->CLast[14]) && it->CLast[14]>0 ) h_fitL[9]->Fill( (mc.q/mc.p-it->TLast[4])/sqrt(it->CLast[14]));
1739 h_fitL[10]->Fill(it->TLast[4]);
1740 h_fitL[11]->Fill(mc.q/mc.p);
1741 }
1742
1743 { // vertex
1744 CbmL1MCTrack mc = *(it->GetMCTracks()[0]);
1745 L1TrackPar trPar(it->T,it->C);
1746
1747 if (mc.mother_ID != -1){ // secondary
1748
1749 { // extrapolate to vertex
1750 L1FieldRegion fld _fvecalignment;
1751 L1FieldValue B[3] _fvecalignment;
1752 float z[3];
1753 for (int ih = 0; ih < 3; ih++){
1754 const int iMCP = mc.Points[ih];
1755 CbmL1MCPoint &mcP = vMCPoints[iMCP];
1756 L1Station &st = algo->vStations[mcP.iStation];
1757 z[ih] = st.z[0];
1758 st.fieldSlice.GetFieldValue( mcP.x, mcP.y, B[ih] );
1759 };
1760 fld.Set( B[0], z[0], B[1], z[1], B[2], z[2] );
1761
1762 L1Extrapolate(trPar, mc.z, trPar.qp, fld);
1763 // add material
1764 const int fSta = vHitStore[it->StsHits[0]].iStation;
1765 const int dir = int((mc.z - algo->vStations[fSta].z[0])/fabs(mc.z - algo->vStations[fSta].z[0]));
1766 // if (abs(mc.z - algo->vStations[fSta].z[0]) > 10.) continue; // can't extrapolate on large distance
1767 for (int iSta = fSta/*+dir*/; (iSta >= 0) && (iSta < NStation) && (dir*(mc.z - algo->vStations[iSta].z[0]) > 0); iSta += dir){
1768 // cout << iSta << " " << dir << endl;
1769 L1AddMaterial( trPar, algo->vStations[iSta].materialInfo, trPar.qp );
1770 if (iSta+dir == NMvdStations-1) L1AddPipeMaterial( trPar, trPar.qp );
1771 }
1772 }
1773 if (mc.z != trPar.z[0]) continue;
1774 // static int good = 0;
1775 // static int bad = 0;
1776 // if (mc.z != trPar.z[0]){
1777 // bad++;
1778 // continue;
1779 // }
1780 // else good++;
1781 // cout << "bad\\good" << bad << " " << good << endl;
1782
1783 // calculate pulls
1784 //h_fitSV[0]->Fill( (mc.x-trPar.x[0]) *1.e4);
1785 //h_fitSV[1]->Fill( (mc.y-trPar.y[0]) *1.e4);
1786 h_fitSV[0]->Fill( (mc.x-trPar.x[0]) );
1787 h_fitSV[1]->Fill( (mc.y-trPar.y[0]) );
1788 h_fitSV[2]->Fill((mc.px/mc.pz-trPar.tx[0])*1.e3);
1789 h_fitSV[3]->Fill((mc.py/mc.pz-trPar.ty[0])*1.e3);
1790 h_fitSV[4]->Fill(trPar.qp[0]/mc.q*mc.p-1);
1791 if( finite(trPar.C00[0]) && trPar.C00[0]>0 ) h_fitSV[5]->Fill( (mc.x-trPar.x[0])/sqrt(trPar.C00[0]));
1792 if( finite(trPar.C11[0]) && trPar.C11[0]>0 ) h_fitSV[6]->Fill( (mc.y-trPar.y[0])/sqrt(trPar.C11[0]));
1793 if( finite(trPar.C22[0]) && trPar.C22[0]>0 ) h_fitSV[7]->Fill( (mc.px/mc.pz-trPar.tx[0])/sqrt(trPar.C22[0]));
1794 if( finite(trPar.C33[0]) && trPar.C33[0]>0 ) h_fitSV[8]->Fill( (mc.py/mc.pz-trPar.ty[0])/sqrt(trPar.C33[0]));
1795 if( finite(trPar.C44[0]) && trPar.C44[0]>0 ) h_fitSV[9]->Fill( (mc.q/mc.p-trPar.qp[0])/sqrt(trPar.C44[0]));
1796 h_fitSV[10]->Fill(trPar.qp[0]);
1797 h_fitSV[11]->Fill(mc.q/mc.p);
1798 }
1799 else{ // primary
1800 if (vHitStore[it->StsHits[0]].iStation != 0) continue; // only mvd tracks
1801
1802// #define L1EXTRAPOLATE
1803#ifdef L1EXTRAPOLATE
1804 { // extrapolate to vertex
1805 L1FieldRegion fld _fvecalignment;
1806 L1FieldValue B[3], targB _fvecalignment;
1807 float z[3];
1808 for (int ih = 0; ih < 3; ih++){
1809 const int iMCP = mc.Points[ih];
1810 CbmL1MCPoint &mcP = vMCPoints[iMCP];
1811 L1Station &st = algo->vStations[mcP.iStation];
1812 z[ih] = st.z[0];
1813// cout << mcP.z << endl;
1814 st.fieldSlice.GetFieldValue( mcP.x, mcP.y, B[ih] );
1815 };
1816// cout << mc.x << " " << mc.y << " " << mc.z << endl;
1817
1818// targetFieldSlice->GetFieldValue( mc.x, mc.y, targB );
1819 targB = algo->vtxFieldValue;
1820 fld.Set( targB, 0., B[0], z[0], B[1], z[1] );
1821
1822 L1Extrapolate(trPar, mc.z, trPar.qp, fld);
1823// L1Extrapolate0(trPar, mc.z, fld);
1824 // add material
1825 const int fSta = vHitStore[it->StsHits[0]].iStation;
1826
1827 const int dir = (mc.z - algo->vStations[fSta].z[0])/abs(mc.z - algo->vStations[fSta].z[0]);
1828 // if (abs(mc.z - algo->vStations[fSta].z[0]) > 10.) continue; // can't extrapolate on large distance
1829 for (int iSta = fSta/*+dir*/; (iSta >= 0) && (iSta < NStation) && (dir*(mc.z - algo->vStations[iSta].z[0]) > 0); iSta += dir){
1830 // cout << iSta << " " << dir << endl;
1831 L1AddMaterial( trPar, algo->vStations[iSta].materialInfo, trPar.qp );
1832 if (iSta+dir == NMvdStations-1) L1AddPipeMaterial( trPar, trPar.qp );
1833 }
1834 }
1835 if (mc.z != trPar.z[0]) continue;
1836 // static int good = 0;
1837 // static int bad = 0;
1838 // if (mc.z != trPar.z[0]){
1839 // bad++;
1840 // continue;
1841 // }
1842 // else good++;
1843 // cout << "bad\\good" << bad << " " << good << endl;
1844
1845 // calculate pulls
1846 //h_fitPV[0]->Fill( (mc.x-trPar.x[0]) *1.e4);
1847 //h_fitPV[1]->Fill( (mc.y-trPar.y[0]) *1.e4);
1848 h_fitPV[0]->Fill( (mc.x-trPar.x[0]) );
1849 h_fitPV[1]->Fill( (mc.y-trPar.y[0]) );
1850 h_fitPV[2]->Fill((mc.px/mc.pz-trPar.tx[0])*1.e3);
1851 h_fitPV[3]->Fill((mc.py/mc.pz-trPar.ty[0])*1.e3);
1852 h_fitPV[4]->Fill(trPar.qp[0]/mc.q*mc.p-1);
1853 if( finite(trPar.C00[0]) && trPar.C00[0]>0 ) h_fitPV[5]->Fill( (mc.x-trPar.x[0])/sqrt(trPar.C00[0]));
1854 if( finite(trPar.C11[0]) && trPar.C11[0]>0 ) h_fitPV[6]->Fill( (mc.y-trPar.y[0])/sqrt(trPar.C11[0]));
1855 if( finite(trPar.C22[0]) && trPar.C22[0]>0 ) h_fitPV[7]->Fill( (mc.px/mc.pz-trPar.tx[0])/sqrt(trPar.C22[0]));
1856 if( finite(trPar.C33[0]) && trPar.C33[0]>0 ) h_fitPV[8]->Fill( (mc.py/mc.pz-trPar.ty[0])/sqrt(trPar.C33[0]));
1857 if( finite(trPar.C44[0]) && trPar.C44[0]>0 ) h_fitPV[9]->Fill( (mc.q/mc.p-trPar.qp[0])/sqrt(trPar.C44[0]));
1858 h_fitPV[10]->Fill(trPar.qp[0]);
1859 h_fitPV[11]->Fill(mc.q/mc.p);
1860#else
1861 FairTrackParam fTP;
1862
1863 CbmKFMath::CopyTC2TrackParam( &fTP, it->T, it->C );
1864
1865 CbmKFTrack kfTr;
1866 kfTr.SetTrackParam(fTP);
1867
1868 kfTr.Extrapolate(mc.z);
1869 CbmL1Track it2;
1870 for (int ipar = 0; ipar < 6; ipar++) it2.T[ipar] = kfTr.GetTrack()[ipar];
1871 for (int ipar = 0; ipar < 15; ipar++) it2.C[ipar] = kfTr.GetCovMatrix()[ipar];
1872
1873
1874 // calculate pulls
1875// h_fitPV[0]->Fill( (mc.x-it2.T[0]) *1.e4);
1876// h_fitPV[1]->Fill( (mc.y-it2.T[1]) *1.e4);
1877 h_fitPV[0]->Fill( (mc.x-it2.T[0]) );
1878 h_fitPV[1]->Fill( (mc.y-it2.T[1]) );
1879 h_fitPV[2]->Fill((mc.px/mc.pz-it2.T[2])*1.e3);
1880 h_fitPV[3]->Fill((mc.py/mc.pz-it2.T[3])*1.e3);
1881 h_fitPV[4]->Fill(it2.T[4]/mc.q*mc.p-1);
1882 if( finite(it2.C[0]) && it2.C[0]>0 )h_fitPV[5]->Fill( (mc.x-it2.T[0])/sqrt(it2.C[0]));
1883 if( finite(it2.C[2]) && it2.C[2]>0 )h_fitPV[6]->Fill( (mc.y-it2.T[1])/sqrt(it2.C[2]));
1884 if( finite(it2.C[5]) && it2.C[5]>0 )h_fitPV[7]->Fill( (mc.px/mc.pz-it2.T[2])/sqrt(it2.C[5]));
1885 if( finite(it2.C[9]) && it2.C[9]>0 )h_fitPV[8]->Fill( (mc.py/mc.pz-it2.T[3])/sqrt(it2.C[9]));
1886 if( finite(it2.C[14]) && it2.C[14]>0 )h_fitPV[9]->Fill( (mc.q/mc.p-it2.T[4])/sqrt(it2.C[14]));
1887 h_fitPV[10]->Fill(it2.T[4]);
1888 h_fitPV[11]->Fill(mc.q/mc.p);
1889#endif // L1EXTRAPOLATE
1890 }
1891 }
1892
1893 h_fit_chi2->Fill(it->chi2/it->NDF);
1894 }
1895
1896} // void CbmL1::TrackFitPerformance()
1897
1898
1899
1900void CbmL1::FieldApproxCheck()
1901{
1902 TDirectory *curr = gDirectory;
1903 TFile *currentFile = gFile;
1904 TFile* fout = new TFile("FieldApprox.root","RECREATE");
1905 fout->cd();
1906
1907 FairField *MF = CbmKF::Instance()->GetMagneticField();
1908 for ( int ist = 0; ist<NStation; ist++ )
1909 {
1910 double z = 0;
1911 double Xmax=-100, Ymax=-100;
1912 if( ist<NMvdStations ){
1913 CbmKFTube &t = CbmKF::Instance()->vMvdMaterial[ist];
1914 z = t.z;
1915 Xmax = Ymax = t.R;
1916 }else{
1917 //AZ CbmStsStation *st = StsDigi.GetStation(ist - NMvdStations);
1918 CbmStsStation *st = StsDigi->GetStation(ist - NMvdStations);
1919 z = st->GetZ();
1920
1921 double x,y;
1922 for(int isec = 0; isec < st->GetNSectors(); isec++)
1923 {
1924 CbmStsSector *sect = L1_DYNAMIC_CAST<CbmStsSector*>( st->GetSector(isec) );
1925 for(int isen = 0; isen < sect->GetNSensors(); isen++)
1926 {
1927 x = sect->GetSensor(isen)->GetX0() + sect->GetSensor(isen)->GetLx()/2.;
1928 y = sect->GetSensor(isen)->GetY0() + sect->GetSensor(isen)->GetLy()/2.;
1929 if(x>Xmax) Xmax = x;
1930 if(y>Ymax) Ymax = y;
1931 }
1932 }
1933 cout << "Station "<< ist << ", Xmax " << Xmax<<", Ymax" << Ymax<<endl;
1934
1935 } // if mvd
1936
1937
1938// float step = 1.;
1939
1940 int NbinsX = 100; //static_cast<int>(2*Xmax/step);
1941 int NbinsY = 100; //static_cast<int>(2*Ymax/step);
1942 float ddx = 2*Xmax/NbinsX;
1943 float ddy = 2*Ymax/NbinsY;
1944
1945 TH2F *stB = new TH2F(Form("station %i, dB", ist+1) ,Form("station %i, dB, z = %0.f cm", ist+1,z) , static_cast<int>(NbinsX+1),-(Xmax+ddx/2.),(Xmax+ddx/2.), static_cast<int>(NbinsY+1),-(Ymax+ddy/2.),(Ymax+ddy/2.));
1946 TH2F *stBx = new TH2F(Form("station %i, dBx", ist+1),Form("station %i, dBx, z = %0.f cm", ist+1,z), static_cast<int>(NbinsX+1),-(Xmax+ddx/2.),(Xmax+ddx/2.), static_cast<int>(NbinsY+1),-(Ymax+ddy/2.),(Ymax+ddy/2.));
1947 TH2F *stBy = new TH2F(Form("station %i, dBy", ist+1),Form("station %i, dBy, z = %0.f cm", ist+1,z), static_cast<int>(NbinsX+1),-(Xmax+ddx/2.),(Xmax+ddx/2.), static_cast<int>(NbinsY+1),-(Ymax+ddy/2.),(Ymax+ddy/2.));
1948 TH2F *stBz = new TH2F(Form("station %i, dBz", ist+1),Form("station %i, dBz, z = %0.f cm", ist+1,z), static_cast<int>(NbinsX+1),-(Xmax+ddx/2.),(Xmax+ddx/2.), static_cast<int>(NbinsY+1),-(Ymax+ddy/2.),(Ymax+ddy/2.));
1949
1950 Double_t r[3],B[3];
1951 L1FieldSlice FSl;
1952 L1FieldValue B_L1;
1953 Double_t bbb, bbb_L1;
1954
1955 const int M=5; // polinom order
1956 const int N=(M+1)*(M+2)/2;
1957 L1Station &st = algo->vStations[ist];
1958 for(int i=0; i<N; i++)
1959 {
1960 FSl.cx[i] = st.fieldSlice.cx[i][0];
1961 FSl.cy[i] = st.fieldSlice.cy[i][0];
1962 FSl.cz[i] = st.fieldSlice.cz[i][0];
1963 }
1964
1965 Int_t i=1,j=1;
1966
1967 double x,y;
1968 for(int ii = 1; ii <=NbinsX+1; ii++)
1969 {
1970 j=1;
1971 x = -Xmax+(ii-1)*ddx;
1972 for(int jj = 1; jj <=NbinsY+1; jj++)
1973 {
1974 y = -Ymax+(jj-1)*ddy;
1975 double rrr = sqrt(fabs(x*x/Xmax/Xmax+y/Ymax*y/Ymax));
1976 if(rrr>1. )
1977 {
1978 j++;
1979 continue;
1980 }
1981 r[2] = z; r[0] = x; r[1] = y;
1982 MF->GetFieldValue( r, B );
1983 bbb = sqrt(B[0]*B[0]+B[1]*B[1]+B[2]*B[2]);
1984
1985 FSl.GetFieldValue(x,y,B_L1);
1986 bbb_L1 = sqrt(B_L1.x[0]*B_L1.x[0] + B_L1.y[0]*B_L1.y[0] + B_L1.z[0]*B_L1.z[0]);
1987
1988 stB -> SetBinContent(ii,jj,(bbb - bbb_L1));
1989 stBx -> SetBinContent(ii,jj,(B[0] - B_L1.x[0]));
1990 stBy -> SetBinContent(ii,jj,(B[1] - B_L1.y[0]));
1991 stBz -> SetBinContent(ii,jj,(B[2] - B_L1.z[0]));
1992 j++;
1993 }
1994 i++;
1995 }
1996
1997 stB ->GetXaxis()->SetTitle("X, cm");
1998 stB ->GetYaxis()->SetTitle("Y, cm");
1999 stB ->GetXaxis()->SetTitleOffset(1);
2000 stB ->GetYaxis()->SetTitleOffset(1);
2001 stB ->GetZaxis()->SetTitle("B_map - B_L1, kGauss");
2002 stB ->GetZaxis()->SetTitleOffset(1.3);
2003
2004 stBx ->GetXaxis()->SetTitle("X, cm");
2005 stBx ->GetYaxis()->SetTitle("Y, cm");
2006 stBx ->GetXaxis()->SetTitleOffset(1);
2007 stBx ->GetYaxis()->SetTitleOffset(1);
2008 stBx ->GetZaxis()->SetTitle("Bx_map - Bx_L1, kGauss");
2009 stBx ->GetZaxis()->SetTitleOffset(1.3);
2010
2011 stBy ->GetXaxis()->SetTitle("X, cm");
2012 stBy ->GetYaxis()->SetTitle("Y, cm");
2013 stBy ->GetXaxis()->SetTitleOffset(1);
2014 stBy ->GetYaxis()->SetTitleOffset(1);
2015 stBy ->GetZaxis()->SetTitle("By_map - By_L1, kGauss");
2016 stBy ->GetZaxis()->SetTitleOffset(1.3);
2017
2018 stBz ->GetXaxis()->SetTitle("X, cm");
2019 stBz ->GetYaxis()->SetTitle("Y, cm");
2020 stBz ->GetXaxis()->SetTitleOffset(1);
2021 stBz ->GetYaxis()->SetTitleOffset(1);
2022 stBz ->GetZaxis()->SetTitle("Bz_map - Bz_L1, kGauss");
2023 stBz ->GetZaxis()->SetTitleOffset(1.3);
2024
2025 stB -> Write();
2026 stBx -> Write();
2027 stBy -> Write();
2028 stBz -> Write();
2029
2030 } // for ista
2031
2032 fout->Close();
2033 fout->Delete();
2034 gFile = currentFile;
2035 gDirectory = curr;
2036} // void CbmL1::FieldApproxCheck()
2037
2038#include "TMath.h"
2039void CbmL1::FieldIntegralCheck()
2040{
2041 TDirectory *curr = gDirectory;
2042 TFile *currentFile = gFile;
2043 TFile* fout = new TFile("FieldApprox.root","RECREATE");
2044 fout->cd();
2045
2046 FairField *MF = CbmKF::Instance()->GetMagneticField();
2047
2048 int nPointsZ = 1000;
2049 int nPointsPhi = 100;
2050 int nPointsTheta = 100;
2051 double startZ=0, endZ=100.;
2052 double startPhi=0, endPhi=2*TMath::Pi();
2053 double startTheta=-30./180.*TMath::Pi(), endTheta=30./180.*TMath::Pi();
2054
2055 double DZ=endZ-startZ;
2056 double DP=endPhi-startPhi;
2057 double DT=endTheta-startTheta;
2058
2059 float ddp = endPhi/nPointsPhi;
2060 float ddt = 2*endTheta/nPointsTheta;
2061
2062 TH2F *hSb = new TH2F("Field Integral", "Field Integral" , static_cast<int>(nPointsPhi),-(startPhi+ddp/2.),(endPhi+ddp/2.), static_cast<int>(nPointsTheta),(startTheta-ddt/2.),(endTheta+ddt/2.));
2063
2064 for(int iP=0; iP<nPointsPhi; iP++)
2065 {
2066 double phi=startPhi+iP*DP/nPointsPhi;
2067 for(int iT=0; iT<nPointsTheta; iT++)
2068 {
2069 double theta=startTheta+iT*DT/nPointsTheta;
2070
2071 double Sb=0;
2072 for(int iZ=1; iZ<nPointsZ; iZ++)
2073 {
2074 double z = startZ+ DZ*iZ/nPointsZ;
2075 double x = z*TMath::Tan(theta)*TMath::Cos(phi);
2076 double y = z*TMath::Tan(theta)*TMath::Sin(phi);
2077 double r[3] = {x,y,z};
2078 double b[3];
2079 MF->GetFieldValue( r, b );
2080 double B=sqrt(b[0]*b[0]+b[1]*b[1]+b[2]*b[2]);
2081 Sb += B*DZ/nPointsZ/100./10.;
2082 }
2083 hSb->SetBinContent(iP+1,iT+1,Sb);
2084 }
2085 }
2086
2087 hSb ->GetXaxis()->SetTitle("#phi [rad]");
2088 hSb ->GetYaxis()->SetTitle("#theta [rad]");
2089 hSb ->GetXaxis()->SetTitleOffset(1);
2090 hSb ->GetYaxis()->SetTitleOffset(1);
2091 hSb ->GetZaxis()->SetTitle("Field Integral [T#dotm]");
2092 hSb ->GetZaxis()->SetTitleOffset(1.3);
2093
2094 hSb -> Write();
2095
2096
2097 fout->Close();
2098 fout->Delete();
2099 gFile = currentFile;
2100 gDirectory = curr;
2101} // void CbmL1::FieldApproxCheck()
2102
2103void CbmL1::InputPerformance()
2104{
2105 static TH1I *nStripFHits, *nStripBHits, *nStripFMC, *nStripBMC;
2106
2107 static TH1F *resX, *resY/*, *pullZ*/;
2108 static TH1F *pullX, *pullY/*, *pullZ*/;
2109
2110 static bool first_call = 1;
2111 if ( first_call ){
2112 first_call = 0;
2113
2114 TDirectory *currentDir = gDirectory;
2115 gDirectory = histodir;
2116 gDirectory->cd("L1");
2117 gDirectory->mkdir("Input");
2118 gDirectory->cd("Input");
2119
2120 nStripFHits = new TH1I("nHits_f", "NHits On Front Strip", 10, 0, 10);
2121 nStripBHits = new TH1I("nHits_b", "NHits On Back Strip", 10, 0, 10);
2122 nStripFMC = new TH1I("nMC_f", "N MC Points On Front Strip", 10, 0, 10);
2123 nStripBMC = new TH1I("nMC_b", "N MC Points On Back Strip", 10, 0, 10);
2124
2125 pullX = new TH1F("Px", "Pull x", 50, -5, 5);
2126 pullY = new TH1F("Py", "Pull y", 50, -5, 5);
2127
2128 resX = new TH1F("x", "Residual x", 50, -50, 50);
2129 resY = new TH1F("y", "Residual y", 50, -500, 500);
2130 TH1* histo;
2131 histo = resX;
2132 histo->GetXaxis()->SetTitle("Residual x, um");
2133 histo = resY;
2134 histo->GetXaxis()->SetTitle("Residual y, um");
2135
2136 gDirectory = currentDir;
2137 } // first_call
2138
2139 std::map<unsigned int, unsigned int> stripFToNHitMap,stripBToNHitMap;
2140 std::map<unsigned int, unsigned int> stripFToNMCMap,stripBToNMCMap;
2141
2142 map<unsigned int, unsigned int>::iterator it;
2143 Int_t nMC = -1;
2144 if( listStsPts ){
2145 nMC = listStsPts->GetEntries();
2146 }
2147 if( listStsHits ){
2148 for (unsigned int iH=0; iH < vStsHits.size(); iH++ ){
2149 const CbmL1StsHit &h = vStsHits[iH];
2150
2151 if (h.extIndex < 0) continue; // mvd hit
2152 const CbmStsHit *sh = L1_DYNAMIC_CAST<CbmStsHit*>( listStsHits->At(h.extIndex) );
2153 // strip - MC correspondence
2154 int iStripF = sh->GetDigi(0);
2155 int iStripB = sh->GetDigi(1);
2156 if ( iStripF >= 0 ) stripFToNHitMap[iStripF]++;
2157 if ( iStripB >= 0 ) stripBToNHitMap[iStripB]++;
2158
2159 if ( h.mcPointIds.size() == 0 ) continue; // fake
2160 int iMC = vMCPoints[h.mcPointIds[0]].pointId; // TODO: adapt to linking
2161 if( !( iMC>=0 && iMC < nMC) )
2162 continue;
2163
2164 if ( iStripF >= 0 ) stripFToNMCMap[iStripF]++;
2165 if ( iStripB >= 0 ) stripBToNMCMap[iStripB]++;
2166
2167 // hit pulls and residuals
2168
2169 TVector3 hitPos, mcPos, hitErr;
2170 sh->Position(hitPos);
2171 sh->PositionError(hitErr);
2172 CbmStsPoint *pt = 0;
2173 pt = L1_DYNAMIC_CAST<CbmStsPoint*>( listStsPts->At(iMC) );
2174 if ( !pt ){
2175// cout << " No MC points! " << "iMC=" << iMC << endl;
2176 continue;
2177 }
2178// pt->Position(mcPos); // this is wrong!
2179 mcPos.SetX( pt->GetX( hitPos.Z() ) );
2180 mcPos.SetY( pt->GetY( hitPos.Z() ) );
2181 mcPos.SetZ( hitPos.Z() );
2182
2183#if 0 // standard errors
2184 if (hitErr.X() != 0) pullX->Fill( (hitPos.X() - mcPos.X()) / hitErr.X() ); // standard errors
2185 if (hitErr.Y() != 0) pullY->Fill( (hitPos.Y() - mcPos.Y()) / hitErr.Y() );
2186#elif 1 // qa errors
2187 if (hitErr.X() != 0) pullX->Fill( (hitPos.X() - mcPos.X()) / sh->GetDx() ); // qa errors
2188 if (hitErr.Y() != 0) pullY->Fill( (hitPos.Y() - mcPos.Y()) / sh->GetDy() );
2189#else // errors used in TF
2190 if (hitErr.X() != 0) pullX->Fill( (hitPos.X() - mcPos.X()) / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]) );
2191 if (hitErr.Y() != 0) pullY->Fill( (hitPos.Y() - mcPos.Y()) / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]) );
2192#endif
2193
2194 resX->Fill((hitPos.X() - mcPos.X())*10*1000);
2195 resY->Fill((hitPos.Y() - mcPos.Y())*10*1000);
2196
2197 }
2198 } // sts
2199 if( listMvdHits ){
2200 Int_t nEnt = listMvdHits->GetEntries();
2201// Int_t nMC = listStsPts->GetEntries();
2202 for (int j=0; j < nEnt; j++ ){
2203 CbmMvdHit *sh = L1_DYNAMIC_CAST<CbmMvdHit*>( listMvdHits->At(j) );
2204 CbmMvdHitMatch *hm = L1_DYNAMIC_CAST<CbmMvdHitMatch*>( listMvdHitMatches->At(j) );
2205
2206// int iMC = sh->GetRefIndex();
2207 int iMC = hm->GetPointId();
2208
2209 // hit pulls and residuals
2210
2211 TVector3 hitPos, mcPos, hitErr;
2212 sh->Position(hitPos);
2213 sh->PositionError(hitErr);
2214 CbmMvdPoint *pt = 0;
2215 if( iMC>=0 && iMC<nMC) pt = L1_DYNAMIC_CAST<CbmMvdPoint*>( listMvdPts->At(iMC) );
2216 if ( !pt ){
2217// cout << " No MC points! " << "iMC=" << iMC << endl;
2218 continue;
2219 }
2220// pt->Position(mcPos); // this is wrong!
2221 mcPos.SetX( ( pt->GetX() + pt->GetXOut() )/2. );
2222 mcPos.SetY( ( pt->GetY() + pt->GetYOut() )/2 );
2223 mcPos.SetZ( hitPos.Z() );
2224
2225// if (hitErr.X() != 0) pullX->Fill( (hitPos.X() - mcPos.X()) / hitErr.X() ); // standard errors
2226// if (hitErr.Y() != 0) pullY->Fill( (hitPos.Y() - mcPos.Y()) / hitErr.Y() );
2227// if (hitErr.X() != 0) pullX->Fill( (hitPos.X() - mcPos.X()) / sh->GetDx() ); // qa errors
2228// if (hitErr.Y() != 0) pullY->Fill( (hitPos.Y() - mcPos.Y()) / sh->GetDy() );
2229 if (hitErr.X() != 0) pullX->Fill( (hitPos.X() - mcPos.X()) / sqrt(algo->vStations[0].XYInfo.C00[0]) ); // errors used in TF
2230 if (hitErr.Y() != 0) pullY->Fill( (hitPos.Y() - mcPos.Y()) / sqrt(algo->vStations[0].XYInfo.C11[0]) );
2231
2232 resX->Fill((hitPos.X() - mcPos.X())*10*1000);
2233 resY->Fill((hitPos.Y() - mcPos.Y())*10*1000);
2234
2235 }
2236 } // mvd
2237
2238
2239
2240 for (it = stripFToNHitMap.begin(); it != stripFToNHitMap.end(); it++){
2241 nStripFHits->Fill(it->second);
2242 }
2243 for (it = stripBToNHitMap.begin(); it != stripBToNHitMap.end(); it++){
2244 nStripBHits->Fill(it->second);
2245 }
2246 for (it = stripFToNMCMap.begin(); it != stripFToNMCMap.end(); it++){
2247 nStripFMC->Fill(it->second);
2248 }
2249 for (it = stripBToNMCMap.begin(); it != stripBToNMCMap.end(); it++){
2250 nStripBMC->Fill(it->second);
2251 }
2252
2253 // strips Not ended
2254// if( listCbmStsDigi ){
2255// Int_t nEnt = listCbmStsDigi->GetEntries();
2256// for (int j=0; j < nEnt; j++ ){
2257// CbmStsDigi *digi = (CbmStsDigi*) listCbmStsDigi->At(j);
2258// // = sh->GetNLinks(0);
2259// // find position of mc point
2260// FairLink link = digi->GetLink(0);
2261// int iMCPoint = link.GetIndex();
2262// CbmStsPoint *mcPoint = (CbmStsPoint*) listStsPts->At(iMCPoint);
2263// TVector3 mcPos;
2264// if (digi->GetSide() == 0)
2265// mcPoint->PositionIn(mcPos);
2266// else
2267// mcPoint->PositionOut(mcPos);
2268//
2269// CbmStsStation *sta = StsDigi.GetStation(digi->GetStationNr());
2270// CbmStsSector* sector = sta->GetSector(digi->GetSectorNr());
2271// digi->GetChannelNr();
2272//
2273// }
2274// } // listCbmStsDigi
2275
2276
2277} // void CbmL1::InputPerformance()
2278
iZ
@ iZ
Definition BmnAlignDefines.h:87
L1_ASSERT
#define L1_ASSERT(v, msg)
Definition CbmL1Def.h:44
L1_assert
#define L1_assert(v)
Definition CbmL1Def.h:50
L1AddPipeMaterial
void L1AddPipeMaterial(L1TrackPar &T, fvec qp0, fvec w=1, fvec mass2=0.1395679f *0.1395679f)
Definition L1AddMaterial.h:84
L1AddMaterial
void L1AddMaterial(L1TrackPar &T, fvec radThick, fvec qp0, fvec w=1, fvec mass2=0.1395679f *0.1395679f)
Definition L1AddMaterial.h:11
L1Extrapolate
void L1Extrapolate(L1TrackPar &T, fvec z_out, fvec qp0, L1FieldRegion &F, fvec *w=0)
Definition L1Extrapolation.h:12
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition P4_F32vec4.h:34
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
_fvecalignment
#define _fvecalignment
Definition P4_F32vec4.h:236
i
int i
Definition P4_F32vec4.h:22
acos
friend F32vec4 acos(const F32vec4 &a)
Definition P4_F32vec4.h:126
CbmKFMath::CopyTC2TrackParam
static void CopyTC2TrackParam(FairTrackParam *par, Double_t T[], Double_t C[])
Definition CbmKFMath.cxx:777
CbmKFPartEfficiencies::Inc
void Inc(bool isReco, int nClones, TString name)
Definition CbmKFPartEfficiencies.h:424
CbmKFPartEfficiencies::IncReco
void IncReco(bool isGhost, bool isBg, TString name)
Definition CbmKFPartEfficiencies.h:434
CbmKFPartEfficiencies::partName
TString partName[nParticles]
Definition CbmKFPartEfficiencies.h:507
CbmKFParticle::NDaughters
int NDaughters() const
Definition CbmKFParticle.h:23
CbmKFParticle::Extrapolate
void Extrapolate(Double_t r0[], double T)
CbmKFParticle::GetMass
void GetMass(Double_t &M, Double_t &Error)
Definition CbmKFParticle.cxx:1446
CbmKFParticle::r
Double_t r[8]
Definition CbmKFParticle.h:40
CbmKFParticle::GetNDF
Int_t GetNDF() const
Definition CbmKFParticle.h:80
CbmKFParticle::GetLifeTime
void GetLifeTime(Double_t &T, Double_t &Error)
Definition CbmKFParticle.cxx:1476
CbmKFParticle::GetDecayLength
void GetDecayLength(Double_t &L, Double_t &Error)
Definition CbmKFParticle.cxx:1459
CbmKFParticle::GetParameter
Double_t GetParameter(Int_t i) const
Definition CbmKFParticle.h:87
CbmKFParticle::GetChi2
Double_t GetChi2() const
Definition CbmKFParticle.h:79
CbmKFParticle::GetPDG
int GetPDG() const
Definition CbmKFParticle.h:130
CbmKFParticle::GetCovariance
Double_t GetCovariance(Int_t i) const
Definition CbmKFParticle.h:88
CbmKFParticle::Id
int Id() const
Definition CbmKFParticle.h:22
CbmKFParticle::DaughterIds
const vector< int > & DaughterIds() const
Definition CbmKFParticle.h:24
CbmKFParticle::GetDStoPoint
Double_t GetDStoPoint(const Double_t xyz[]) const
Definition CbmKFParticle.cxx:1502
CbmKFTrackInterface::Extrapolate
Int_t Extrapolate(Double_t z, Double_t *QP0=0, Bool_t line=false)
Access to i-th hit.
Definition CbmKFTrackInterface.cxx:37
CbmKFTrack::SetTrackParam
void SetTrackParam(FairTrackParam &track)
Definition CbmKFTrack.cxx:28
CbmKFTrack::GetTrack
Double_t * GetTrack()
Is it electron.
Definition CbmKFTrack.h:51
CbmKFTrack::GetCovMatrix
Double_t * GetCovMatrix()
array[6] of track parameters(x,y,tx,ty,qp,z)
Definition CbmKFTrack.h:52
CbmKFTube::R
Double_t R
Definition CbmKFMaterial.h:86
CbmKFTube::z
Double_t z
Definition CbmKFMaterial.h:85
CbmKFVertex::GetRefZ
Double_t & GetRefZ()
Definition CbmKFVertex.h:21
CbmKFVertex::GetRefX
Double_t & GetRefX()
Definition CbmKFVertex.h:19
CbmKFVertex::GetCovMatrix
Double_t * GetCovMatrix()
Definition CbmKFVertex.h:22
CbmKFVertex::GetRefY
Double_t & GetRefY()
Definition CbmKFVertex.h:20
CbmKF::Instance
static CbmKF * Instance()
Definition CbmKF.h:35
CbmKF::GetMagneticField
BmnNewFieldMap * GetMagneticField()
Definition CbmKF.h:73
CbmKF::vMvdMaterial
std::vector< CbmKFTube > vMvdMaterial
Definition CbmKF.h:60
CbmL1HitStore::y
double y
Definition CbmL1.h:45
CbmL1HitStore::x
double x
Definition CbmL1.h:45
CbmL1HitStore::iStation
int iStation
Definition CbmL1.h:44
CbmL1MCTrack::GetNClones
int GetNClones() const
Definition CbmL1MCTrack.h:60
CbmL1MCTrack::AddRecoTrack
void AddRecoTrack(CbmL1Track *rTr)
Definition CbmL1MCTrack.h:58
CbmL1MCTrack::IsReconstructed
bool IsReconstructed() const
Definition CbmL1MCTrack.h:61
CbmL1MCTrack::IsReconstructable
bool IsReconstructable() const
Definition CbmL1MCTrack.h:50
CbmL1MCTrack::NStations
int NStations() const
Definition CbmL1MCTrack.h:52
CbmL1MCTrack::StsHits
vector< int > StsHits
Definition CbmL1MCTrack.h:35
CbmL1MCTrack::AddTouchTrack
void AddTouchTrack(CbmL1Track *tTr)
Definition CbmL1MCTrack.h:63
CbmL1MCTrack::GetRecoTracks
vector< CbmL1Track * > & GetRecoTracks()
Definition CbmL1MCTrack.h:59
CbmL1MCTrack::IsAdditional
bool IsAdditional() const
Definition CbmL1MCTrack.h:51
CbmL1MCTrack::IsDisturbed
bool IsDisturbed() const
Definition CbmL1MCTrack.h:64
CbmL1MCTrack::NMCStations
int NMCStations() const
Definition CbmL1MCTrack.h:53
CbmL1MCTrack::IsPrimary
bool IsPrimary() const
Definition CbmL1MCTrack.h:49
CbmL1MCTrack::Points
vector< int > Points
Definition CbmL1MCTrack.h:34
CbmL1StsHit::mcPointIds
vector< int > mcPointIds
Definition CbmL1StsHit.h:17
CbmL1Track::GetMaxPurity
double GetMaxPurity()
Definition CbmL1Track.h:47
CbmL1Track::GetMCTrack
CbmL1MCTrack * GetMCTrack()
Definition CbmL1Track.h:42
CbmL1Track::hitMap
map< int, int > hitMap
Definition CbmL1Track.h:57
CbmL1Track::AddMCTrack
void AddMCTrack(CbmL1MCTrack *mcTr)
Definition CbmL1Track.h:40
CbmL1Track::SetMaxPurity
void SetMaxPurity(double maxPurity_)
Definition CbmL1Track.h:46
CbmL1Track::StsHits
vector< int > StsHits
Definition CbmL1Track.h:54
CbmL1Track::IsGhost
bool IsGhost()
Definition CbmL1Track.h:44
CbmL1::algo
L1Algo * algo
Definition CbmL1.h:55
CbmL1::vHitStore
vector< CbmL1HitStore > vHitStore
Definition CbmL1.h:87
CbmL1::vRTracks
vector< CbmL1Track > vRTracks
Definition CbmL1.h:58
CbmL1::PF
CbmL1ParticlesFinder * PF
for access to L1 Algorithm from L1::Instance
Definition CbmL1.h:56
CbmMCTrack::GetMass
Double_t GetMass() const
Definition CbmMCTrack.cxx:122
CbmMCTrack::GetPy
Double_t GetPy() const
Definition CbmMCTrack.h:59
CbmMCTrack::GetPz
Double_t GetPz() const
Definition CbmMCTrack.h:60
CbmMCTrack::GetPx
Double_t GetPx() const
Definition CbmMCTrack.h:58
CbmMCTrack::GetStartZ
Double_t GetStartZ() const
Definition CbmMCTrack.h:63
CbmMCTrack::GetStartY
Double_t GetStartY() const
Definition CbmMCTrack.h:62
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition CbmMCTrack.h:57
CbmMCTrack::GetP
Double_t GetP() const
Definition CbmMCTrack.h:68
CbmMCTrack::GetStartX
Double_t GetStartX() const
Definition CbmMCTrack.h:61
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition CbmMCTrack.h:56
CbmMvdHitMatch::GetPointId
Int_t GetPointId() const
Definition CbmMvdHitMatch.h:45
CbmMvdPoint::GetYOut
Double_t GetYOut() const
Definition CbmMvdPoint.h:63
CbmMvdPoint::GetXOut
Double_t GetXOut() const
Definition CbmMvdPoint.h:62
CbmStsDigiScheme::GetStation
CbmStsStation * GetStation(Int_t iStation)
Definition CbmStsDigiScheme.cxx:184
CbmStsHit::GetDigi
Int_t GetDigi(Int_t iSide) const
Definition CbmStsHit.cxx:116
CbmStsPoint::GetY
Double_t GetY(Double_t z) const
Definition CbmStsPoint.cxx:112
CbmStsPoint::GetX
Double_t GetX(Double_t z) const
Definition CbmStsPoint.cxx:101
CbmStsSector::GetSensor
CbmStsSensor * GetSensor(Int_t iSensor)
Definition CbmStsSector.h:85
CbmStsSector::GetNSensors
Int_t GetNSensors() const
Definition CbmStsSector.h:82
CbmStsSensor::GetX0
Double_t GetX0() const
Definition CbmStsSensor.h:94
CbmStsSensor::GetLy
Double_t GetLy() const
Definition CbmStsSensor.h:99
CbmStsSensor::GetY0
Double_t GetY0() const
Definition CbmStsSensor.h:95
CbmStsSensor::GetLx
Double_t GetLx() const
Definition CbmStsSensor.h:98
CbmStsStation::GetNSectors
Int_t GetNSectors() const
Definition CbmStsStation.h:62
CbmStsStation::GetZ
Double_t GetZ(Int_t it=0)
Definition CbmStsStation.cxx:152
CbmStsStation::GetSector
CbmStsSector * GetSector(Int_t iSector)
Definition CbmStsStation.h:65
KFMCParticle::IsReconstructable
bool IsReconstructable() const
Definition KFMCParticle.h:41
KFMCParticle::GetPDG
int GetPDG() const
Definition KFMCParticle.h:40
KFMCParticle::SetMCTrackID
void SetMCTrackID(int id)
Definition KFMCParticle.h:35
KFMCParticle::GetDaughterIds
const std::vector< int > & GetDaughterIds() const
Definition KFMCParticle.h:31
KFMCParticle::NDaughters
int NDaughters() const
Definition KFMCParticle.h:30
KFMCParticle::GetMotherId
int GetMotherId() const
Definition KFMCParticle.h:39
KFMCParticle::SetAsReconstructable
void SetAsReconstructable()
Definition KFMCParticle.h:42
KFMCParticle::SetPDG
void SetPDG(int pdg)
Definition KFMCParticle.h:34
KFMCParticle::SetMotherId
void SetMotherId(int id)
Definition KFMCParticle.h:36
KFMCParticle::GetMCTrackID
int GetMCTrackID() const
Definition KFMCParticle.h:38
KFMCParticle::AddDaughter
void AddDaughter(int i)
Definition KFMCParticle.cxx:28
L1Algo::NStations
int NStations
Definition L1Algo.h:123
L1Algo::CATime
double CATime
Definition L1Algo.h:154
L1Algo::vStsHits
vector< L1StsHit > vStsHits
Definition L1Algo.h:132
L1FieldSlice::GetFieldValue
void GetFieldValue(const fvec &x, const fvec &y, L1FieldValue &B) const
Definition L1Field.h:44
L1FieldSlice::cy
fvec cy[21]
Definition L1Field.h:37
L1FieldSlice::cx
fvec cx[21]
Definition L1Field.h:37
L1FieldSlice::cz
fvec cz[21]
Definition L1Field.h:37
L1Station::fieldSlice
L1FieldSlice fieldSlice
Definition L1Station.h:21
CbmL1Constants::MinPurity
const double MinPurity
Definition CbmL1Constants.h:8
CbmL1Constants::MinRefMom
const double MinRefMom
Definition CbmL1Constants.h:7
centrality_clusterization.float
float
Definition centrality_clusterization.py:212
centrality_clusterization.int
int
Definition centrality_clusterization.py:210
nlohmann::detail::cbor_tag_handler_t::error
@ error
throw a parse_error exception in case of a tag
setup.name
name
Definition setup.py:7
CbmL1TrackPar::C
double C[15]
Definition CbmL1TrackPar.h:18
TL1Efficiencies::reco
TL1TracksCatCounters< int > reco
Definition CbmL1Counters.h:132
TL1Efficiencies::mc
TL1TracksCatCounters< int > mc
Definition CbmL1Counters.h:131
TL1Efficiencies::names
vector< TString > names
Definition CbmL1Counters.h:124
TL1Efficiencies::indices
map< TString, int > indices
Definition CbmL1Counters.h:125
TL1Efficiencies::operator+=
TL1Efficiencies & operator+=(TL1Efficiencies &a)
Definition CbmL1Counters.h:162
TL1Efficiencies::ratio_reco
TL1TracksCatCounters< double > ratio_reco
Definition CbmL1Counters.h:127
TL1Efficiencies::Inc
void Inc(bool isReco, TString name)
Definition CbmL1Counters.h:171
TL1Efficiencies::AddCounter
virtual void AddCounter(TString shortname, TString name)
Definition CbmL1Counters.h:138
TL1PerfEfficiencies::ratio_clone
TL1TracksCatCounters< double > ratio_clone
Definition CbmL1Performance.cxx:243
TL1PerfEfficiencies::mc_length_hits
TL1TracksCatCounters< int > mc_length_hits
Definition CbmL1Performance.cxx:252
TL1PerfEfficiencies::ratio_length
TL1TracksCatCounters< double > ratio_length
Definition CbmL1Performance.cxx:244
TL1PerfEfficiencies::ratio_fakes
TL1TracksCatCounters< double > ratio_fakes
Definition CbmL1Performance.cxx:245
TL1PerfEfficiencies::AddCounter
virtual void AddCounter(TString shortname, TString name)
Definition CbmL1Performance.cxx:169
TL1PerfEfficiencies::reco_length
TL1TracksCatCounters< double > reco_length
Definition CbmL1Performance.cxx:249
TL1PerfEfficiencies::mc_length
TL1TracksCatCounters< int > mc_length
Definition CbmL1Performance.cxx:251
TL1PerfEfficiencies::reco_fakes
TL1TracksCatCounters< double > reco_fakes
Definition CbmL1Performance.cxx:250
TL1PerfEfficiencies::Inc
void Inc(bool isReco, bool isKilled, double _ratio_length, double _ratio_fakes, int _nclones, int _mc_length, int _mc_length_hits, TString name)
Definition CbmL1Performance.cxx:202
TL1PerfEfficiencies::killed
TL1TracksCatCounters< int > killed
Definition CbmL1Performance.cxx:247
TL1PerfEfficiencies::ratio_killed
TL1TracksCatCounters< double > ratio_killed
Definition CbmL1Performance.cxx:242
TL1PerfEfficiencies::clone
TL1TracksCatCounters< int > clone
Definition CbmL1Performance.cxx:248
TL1PerfEfficiencies::operator+=
TL1PerfEfficiencies & operator+=(TL1PerfEfficiencies &a)
Definition CbmL1Performance.cxx:183
TL1TracksCatCounters
counters used for efficiency calculation
Definition CbmL1Counters.h:18