2 //***************************************************************************
3 // This file is property of and copyright by the ALICE HLT Project *
4 // ALICE Experiment at CERN, All rights reserved. *
6 // Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
7 // Ivan Kisel <kisel@kip.uni-heidelberg.de> *
8 // for The ALICE HLT Project. *
10 // Permission to use, copy, modify and distribute this software and its *
11 // documentation strictly for non-commercial purposes is hereby granted *
12 // without fee, provided that the above copyright notice appears in all *
13 // copies and that both the copyright notice and this permission notice *
14 // appear in the supporting documentation. The authors make no claims *
15 // about the suitability of this software for any purpose. It is *
16 // provided "as is" without express or implied warranty. *
17 //***************************************************************************
19 #include "AliHLTTPCCAPerformance.h"
20 #include "AliHLTTPCCAGBHit.h"
21 #include "AliHLTTPCCAMCTrack.h"
22 #include "AliHLTTPCCAMCPoint.h"
23 #include "AliHLTTPCCAOutTrack.h"
24 #include "AliHLTTPCCAGBTrack.h"
25 #include "AliHLTTPCCAGBTracker.h"
26 #include "AliHLTTPCCATracker.h"
31 #include "Riostream.h"
38 ClassImp(AliHLTTPCCAPerformance)
42 AliHLTTPCCAPerformance::AliHLTTPCCAPerformance()
103 AliHLTTPCCAPerformance::AliHLTTPCCAPerformance(const AliHLTTPCCAPerformance&)
128 fStatGBNClonesAll(0),
131 fStatGBNClonesRef(0),
163 AliHLTTPCCAPerformance &AliHLTTPCCAPerformance::operator=(const AliHLTTPCCAPerformance&)
169 AliHLTTPCCAPerformance::~AliHLTTPCCAPerformance()
175 void AliHLTTPCCAPerformance::SetTracker( AliHLTTPCCAGBTracker *Tracker )
177 //* set pointer to HLT CA Global tracker
181 void AliHLTTPCCAPerformance::StartEvent()
184 if( !fHistoDir ) CreateHistos();
185 if( fHitLabels ) delete[] fHitLabels;
188 if( fMCTracks ) delete[] fMCTracks;
191 if( fMCPoints ) delete[] fMCPoints;
196 void AliHLTTPCCAPerformance::SetNHits( Int_t NHits )
198 //* set number of hits
199 if( fHitLabels ) delete[] fHitLabels;
201 fHitLabels = new AliHLTTPCCAHitLabel[ NHits ];
205 void AliHLTTPCCAPerformance::SetNMCTracks( Int_t NMCTracks )
207 //* set number of MC tracks
208 if( fMCTracks ) delete[] fMCTracks;
210 fMCTracks = new AliHLTTPCCAMCTrack[ NMCTracks ];
211 fNMCTracks = NMCTracks;
214 void AliHLTTPCCAPerformance::SetNMCPoints( Int_t NMCPoints )
216 //* set number of MC points
217 if( fMCPoints ) delete[] fMCPoints;
219 fMCPoints = new AliHLTTPCCAMCPoint[ NMCPoints ];
223 void AliHLTTPCCAPerformance::ReadHitLabel( Int_t HitID,
224 Int_t lab0, Int_t lab1, Int_t lab2 )
226 //* read the hit labels
227 AliHLTTPCCAHitLabel hit;
231 fHitLabels[HitID] = hit;
234 void AliHLTTPCCAPerformance::ReadMCTrack( Int_t index, const TParticle *part )
236 //* read mc track to the local array
237 fMCTracks[index] = AliHLTTPCCAMCTrack(part);
240 void AliHLTTPCCAPerformance::ReadMCTPCTrack( Int_t index, Float_t X, Float_t Y, Float_t Z,
241 Float_t Px, Float_t Py, Float_t Pz )
243 //* read mc track parameters at TPC
244 fMCTracks[index].SetTPCPar(X,Y,Z,Px,Py,Pz);
247 void AliHLTTPCCAPerformance::ReadMCPoint( Int_t TrackID, Float_t X, Float_t Y, Float_t Z, Float_t Time, Int_t iSlice )
249 //* read mc point to the local array
250 AliHLTTPCCAMCPoint &p = fMCPoints[fNMCPoints];
251 p.TrackID() = TrackID;
257 fTracker->Slices()[iSlice].Param().Global2Slice( X, Y, Z,
258 &p.Sx(), &p.Sy(), &p.Sz() );
259 if( X*X + Y*Y>10.) fNMCPoints++;
262 void AliHLTTPCCAPerformance::CreateHistos()
264 //* create performance histogramms
265 TDirectory *curdir = gDirectory;
266 fHistoDir = gROOT->mkdir("HLTTPCCATrackerPerformance");
268 gDirectory->mkdir("TrackFit");
269 gDirectory->cd("TrackFit");
271 fhResY = new TH1D("resY", "track Y resoltion [cm]", 30, -.5, .5);
272 fhResZ = new TH1D("resZ", "track Z resoltion [cm]", 30, -.5, .5);
273 fhResSinPhi = new TH1D("resSinPhi", "track SinPhi resoltion ", 30, -.03, .03);
274 fhResDzDs = new TH1D("resDzDs", "track DzDs resoltion ", 30, -.01, .01);
275 fhResPt = new TH1D("resPt", "track telative Pt resoltion", 30, -.2, .2);
276 fhPullY = new TH1D("pullY", "track Y pull", 30, -10., 10.);
277 fhPullZ = new TH1D("pullZ", "track Z pull", 30, -10., 10.);
278 fhPullSinPhi = new TH1D("pullSinPhi", "track SinPhi pull", 30, -10., 10.);
279 fhPullDzDs = new TH1D("pullDzDs", "track DzDs pull", 30, -10., 10.);
280 fhPullQPt = new TH1D("pullQPt", "track Q/Pt pull", 30, -10., 10.);
282 gDirectory->cd("..");
284 fhEffVsP = new TProfile("EffVsP", "Eff vs P", 100, 0., 5.);
285 fhGBEffVsP = new TProfile("GBEffVsP", "Global tracker: Eff vs P", 100, 0., 5.);
287 gDirectory->mkdir("Clusters");
288 gDirectory->cd("Clusters");
289 fhHitResY = new TH1D("resHitY", "Y cluster resoltion [cm]", 100, -2., 2.);
290 fhHitResZ = new TH1D("resHitZ", "Z cluster resoltion [cm]", 100, -2., 2.);
291 fhHitPullY = new TH1D("pullHitY", "Y cluster pull", 50, -10., 10.);
292 fhHitPullZ = new TH1D("pullHitZ", "Z cluster pull", 50, -10., 10.);
294 fhHitResY1 = new TH1D("resHitY1", "Y cluster resoltion [cm]", 100, -2., 2.);
295 fhHitResZ1 = new TH1D("resHitZ1", "Z cluster resoltion [cm]", 100, -2., 2.);
296 fhHitPullY1 = new TH1D("pullHitY1", "Y cluster pull", 50, -10., 10.);
297 fhHitPullZ1 = new TH1D("pullHitZ1", "Z cluster pull", 50, -10., 10.);
299 fhHitErrY = new TH1D("HitErrY", "Y cluster error [cm]", 100, 0., 1.);
300 fhHitErrZ = new TH1D("HitErrZ", "Z cluster error [cm]", 100, 0., 1.);
302 gDirectory->cd("..");
304 gDirectory->mkdir("Cells");
305 gDirectory->cd("Cells");
306 fhCellPurity = new TH1D("CellPurity", "Cell Purity", 100, -0.1, 1.1);
307 fhCellNHits = new TH1D("CellNHits", "Cell NHits", 40, 0., 40.);
308 fhCellPurityVsN = new TProfile("CellPurityVsN", "Cell purity Vs N hits", 40, 2., 42.);
309 fhCellPurityVsPt = new TProfile("CellPurityVsPt", "Cell purity Vs Pt", 100, 0., 5.);
310 gDirectory->cd("..");
315 void AliHLTTPCCAPerformance::WriteDir2Current( TObject *obj )
317 //* recursive function to copy the directory 'obj' to the current one
318 if( !obj->IsFolder() ) obj->Write();
320 TDirectory *cur = gDirectory;
321 TDirectory *sub = cur->mkdir(obj->GetName());
323 TList *listSub = ((TDirectory*)obj)->GetList();
325 while( TObject *obj1=it() ) WriteDir2Current(obj1);
330 void AliHLTTPCCAPerformance::WriteHistos()
332 //* write histograms to the file
333 TDirectory *curr = gDirectory;
334 // Open output file and write histograms
335 TFile* outfile = new TFile("HLTTPCCATrackerPerformance.root","RECREATE");
337 WriteDir2Current(fHistoDir);
343 void AliHLTTPCCAPerformance::SlicePerformance( Int_t iSlice, Bool_t PrintFlag )
345 //* calculate slice tracker performance
346 if( !fTracker ) return;
348 Int_t nRecTot = 0, nGhost=0, nRecOut=0;
349 Int_t nMCAll = 0, nRecAll=0, nClonesAll=0;
350 Int_t nMCRef = 0, nRecRef=0, nClonesRef=0;
351 AliHLTTPCCATracker &slice = fTracker->Slices()[iSlice];
353 Int_t firstSliceHit = 0;
354 for( ; firstSliceHit<fTracker->NHits(); firstSliceHit++){
355 if( fTracker->Hits()[firstSliceHit].ISlice()==iSlice ) break;
357 Int_t endSliceHit = firstSliceHit;
359 for( ; endSliceHit<fTracker->NHits(); endSliceHit++){
360 if( fTracker->Hits()[endSliceHit].ISlice()!=iSlice ) break;
363 { // Cell construction performance
365 for( Int_t iRow=0; iRow<slice.Param().NRows(); iRow++ ){
366 AliHLTTPCCARow &row = slice.Rows()[iRow];
367 for (Int_t ic = 0; ic<row.NCells(); ic++){
368 AliHLTTPCCACell &c = row.Cells()[ic];
369 Int_t *lb = new Int_t[c.NHits()*3];
371 for( Int_t j=0; j<c.NHits(); j++){
372 AliHLTTPCCAHit &h = row.GetCellHit(c,j);
373 //cout<<"hit ID="<<h.ID()<<" of"<<fTracker->NHits()<<endl;
374 //cout<<"gb hit ID="<<fTracker->Hits()[h.ID()].ID()<<endl;
375 AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[h.ID()].ID()];
376 if( l.fLab[0]>=0 ) lb[nla++]= l.fLab[0];
377 if( l.fLab[1]>=0 ) lb[nla++]= l.fLab[1];
378 if( l.fLab[2]>=0 ) lb[nla++]= l.fLab[2];
381 Int_t labmax = -1, labcur=-1, lmax = 0, lcurr=0;
382 for( Int_t i=0; i<nla; i++ ){
384 if( labcur>=0 && lmax<lcurr ){
393 if( labcur>=0 && lmax<lcurr ){
398 Int_t label = labmax;
400 for( Int_t j=0; j<c.NHits(); j++){
401 AliHLTTPCCAHit &h = row.GetCellHit(c,j);
402 AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[h.ID()].ID()];
403 if( l.fLab[0]==label || l.fLab[1]==label || l.fLab[2]==label ) lmax++;
407 if( nla>0 && label>=0 ){
408 double purity = double(lmax)/double(nla);
409 fhCellPurity->Fill(purity);
410 fhCellPurityVsN->Fill(c.NHits(),purity);
411 fhCellPurityVsPt->Fill(fMCTracks[label].Pt(),purity);
413 fhCellNHits->Fill(c.NHits());
419 // Select reconstructable MC tracks
422 for (Int_t imc=0; imc<fNMCTracks; imc++) fMCTracks[imc].NHits() = 0;
424 for( Int_t ih=firstSliceHit; ih<endSliceHit; ih++){
425 AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[ih].ID()];
426 if( l.fLab[0]>=0 ) fMCTracks[l.fLab[0]].NHits()++;
427 if( l.fLab[1]>=0 ) fMCTracks[l.fLab[1]].NHits()++;
428 if( l.fLab[2]>=0 ) fMCTracks[l.fLab[2]].NHits()++;
431 for (Int_t imc=0; imc<fNMCTracks; imc++) {
432 AliHLTTPCCAMCTrack &mc = fMCTracks[imc];
434 mc.NReconstructed() = 0;
436 if( mc.NHits() >= 10 && mc.P()>=.05 ){
447 Int_t traN = slice.NOutTracks();
448 Int_t *traLabels = 0;
449 Double_t *traPurity = 0;
450 traLabels = new Int_t[traN];
451 traPurity = new Double_t[traN];
453 for (Int_t itr=0; itr<traN; itr++) {
456 AliHLTTPCCAOutTrack &tCA = slice.OutTracks()[itr];
457 Int_t nhits = tCA.NHits();
458 Int_t *lb = new Int_t[nhits*3];
460 //cout<<"\nHit labels:"<<endl;
461 for( Int_t ihit=0; ihit<nhits; ihit++){
462 Int_t index = slice.OutTrackHits()[tCA.FirstHitRef()+ihit];
463 AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[index].ID()];
464 //cout<<l.fLab[0]<<" "<<l.fLab[1]<<" "<<l.fLab[2]<<endl;
465 if(l.fLab[0]>=0 ) lb[nla++]= l.fLab[0];
466 if(l.fLab[1]>=0 ) lb[nla++]= l.fLab[1];
467 if(l.fLab[2]>=0 ) lb[nla++]= l.fLab[2];
470 Int_t labmax = -1, labcur=-1, lmax = 0, lcurr=0;
471 for( Int_t i=0; i<nla; i++ ){
473 if( labcur>=0 && lmax<lcurr ){
482 if( labcur>=0 && lmax<lcurr ){
487 for( Int_t ihit=0; ihit<nhits; ihit++){
488 Int_t index = slice.OutTrackHits()[tCA.FirstHitRef()+ihit];
489 AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[index].ID()];
490 if( l.fLab[0] == labmax || l.fLab[1] == labmax || l.fLab[2] == labmax
493 traLabels[itr] = labmax;
494 traPurity[itr] = ( (nhits>0) ?double(lmax)/double(nhits) :0 );
495 //cout<<"perf track "<<itr<<": "<<nhits<<" "<<labmax<<" "<<traPurity[itr]<<endl;
496 if( lb ) delete[] lb;
502 for(Int_t itr=0; itr<traN; itr++){
503 if( traPurity[itr]<.9 || traLabels[itr]<0 || traLabels[itr]>=fNMCTracks){
508 AliHLTTPCCAMCTrack &mc = fMCTracks[traLabels[itr]];
509 mc.NReconstructed()++;
510 if( mc.Set()== 0 ) nRecOut++;
512 if( mc.NReconstructed()==1 ) nRecAll++;
513 else if(mc.NReconstructed() > mc.NTurns() ) nClonesAll++;
515 if( mc.NReconstructed()==1 ) nRecRef++;
516 else if(mc.NReconstructed() > mc.NTurns() ) nClonesRef++;
521 for (Int_t ipart=0; ipart<fNMCTracks; ipart++) {
522 AliHLTTPCCAMCTrack &mc = fMCTracks[ipart];
523 if( mc.Set()>0 ) fhEffVsP->Fill(mc.P(), ( mc.NReconstructed()>0 ?1 :0));
526 if( traLabels ) delete[] traLabels;
527 if( traPurity ) delete[] traPurity;
529 fStatNRecTot += nRecTot;
530 fStatNRecOut += nRecOut;
531 fStatNGhost += nGhost;
532 fStatNMCAll += nMCAll;
533 fStatNRecAll += nRecAll;
534 fStatNClonesAll += nClonesAll;
535 fStatNMCRef += nMCRef;
536 fStatNRecRef += nRecRef;
537 fStatNClonesRef += nClonesRef;
539 if( nMCAll ==0 ) return;
542 cout<<"Performance for slice "<<iSlice<<" : "<<endl;
544 <<nMCAll<<" mc all, "
545 <<nMCRef<<" mc ref, "
546 <<nRecTot<<" rec total, "
547 <<nRecAll<<" rec all, "
548 <<nClonesAll<<" clones all, "
549 <<nRecRef<<" rec ref, "
550 <<nClonesRef<<" clones ref, "
552 <<nGhost<<" ghost"<<endl;
554 Int_t nRecExtr = nRecAll - nRecRef;
555 Int_t nMCExtr = nMCAll - nMCRef;
556 Int_t nClonesExtr = nClonesAll - nClonesRef;
558 Double_t dRecTot = (nRecTot>0 ) ? nRecTot :1;
559 Double_t dMCAll = (nMCAll>0 ) ? nMCAll :1;
560 Double_t dMCRef = (nMCRef>0 ) ? nMCRef :1;
561 Double_t dMCExtr = (nMCExtr>0 ) ? nMCExtr :1;
562 Double_t dRecAll = (nRecAll+nClonesAll>0 ) ? nRecAll+nClonesAll :1;
563 Double_t dRecRef = (nRecRef+nClonesRef>0 ) ? nRecRef+nClonesRef :1;
564 Double_t dRecExtr = (nRecExtr+nClonesExtr>0 ) ? nRecExtr+nClonesExtr :1;
567 if( nMCRef>0 ) cout<<nRecRef/dMCRef; else cout<<"_";
568 cout<<", CloneRef = ";
569 if( nRecRef >0 ) cout << nClonesRef/dRecRef; else cout<<"_";
571 cout<<" EffExtra = ";
572 if( nMCExtr>0 ) cout << nRecExtr/dMCExtr; else cout<<"_";
573 cout <<", CloneExtra = ";
574 if( nRecExtr>0 ) cout << nClonesExtr/dRecExtr; else cout<<"_";
577 if( nMCAll>0 ) cout<<nRecAll/dMCAll; else cout<<"_";
578 cout <<", CloneAll = ";
579 if( nRecAll>0 ) cout << nClonesAll/dRecAll; else cout<<"_";
582 if( nRecTot>0 ) cout <<nRecOut/dRecTot; else cout<<"_";
584 if( nRecTot>0 ) cout<<nGhost/dRecTot; else cout<<"_";
590 void AliHLTTPCCAPerformance::Performance()
592 // main routine for performance calculation
607 for( Int_t islice=0; islice<fTracker->NSlices(); islice++){
608 SlicePerformance(islice,0);
611 // global tracker performance
613 if( !fTracker ) return;
615 Int_t nRecTot = 0, nGhost=0, nRecOut=0;
616 Int_t nMCAll = 0, nRecAll=0, nClonesAll=0;
617 Int_t nMCRef = 0, nRecRef=0, nClonesRef=0;
619 // Select reconstructable MC tracks
622 for (Int_t imc=0; imc<fNMCTracks; imc++) fMCTracks[imc].NHits() = 0;
624 for( Int_t ih=0; ih<fNHits; ih++){
625 AliHLTTPCCAHitLabel &l = fHitLabels[ih];
626 if( l.fLab[0]>=0 ) fMCTracks[l.fLab[0]].NHits()++;
627 if( l.fLab[1]>=0 ) fMCTracks[l.fLab[1]].NHits()++;
628 if( l.fLab[2]>=0 ) fMCTracks[l.fLab[2]].NHits()++;
631 for (Int_t imc=0; imc<fNMCTracks; imc++) {
632 AliHLTTPCCAMCTrack &mc = fMCTracks[imc];
634 mc.NReconstructed() = 0;
636 if( mc.NHits() >= 50 && mc.P()>=.05 ){
647 Int_t traN = fTracker->NTracks();
648 Int_t *traLabels = 0;
649 Double_t *traPurity = 0;
650 traLabels = new Int_t[traN];
651 traPurity = new Double_t[traN];
653 for (Int_t itr=0; itr<traN; itr++) {
656 AliHLTTPCCAGBTrack &tCA = fTracker->Tracks()[itr];
657 Int_t nhits = tCA.NHits();
658 Int_t *lb = new Int_t[nhits*3];
660 for( Int_t ihit=0; ihit<nhits; ihit++){
661 Int_t index = fTracker->TrackHits()[tCA.FirstHitRef()+ihit];
662 AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[index].ID()];
663 if(l.fLab[0]>=0 ) lb[nla++]= l.fLab[0];
664 if(l.fLab[1]>=0 ) lb[nla++]= l.fLab[1];
665 if(l.fLab[2]>=0 ) lb[nla++]= l.fLab[2];
668 Int_t labmax = -1, labcur=-1, lmax = 0, lcurr=0;
669 for( Int_t i=0; i<nla; i++ ){
671 if( labcur>=0 && lmax<lcurr ){
680 if( labcur>=0 && lmax<lcurr ){
685 for( Int_t ihit=0; ihit<nhits; ihit++){
686 Int_t index = fTracker->TrackHits()[tCA.FirstHitRef()+ihit];
687 AliHLTTPCCAHitLabel &l = fHitLabels[fTracker->Hits()[index].ID()];
688 if( l.fLab[0] == labmax || l.fLab[1] == labmax || l.fLab[2] == labmax
691 traLabels[itr] = labmax;
692 traPurity[itr] = ( (nhits>0) ?double(lmax)/double(nhits) :0 );
693 if( lb ) delete[] lb;
698 for(Int_t itr=0; itr<traN; itr++){
699 if( traPurity[itr]<.9 || traLabels[itr]<0 || traLabels[itr]>=fNMCTracks){
704 AliHLTTPCCAMCTrack &mc = fMCTracks[traLabels[itr]];
706 mc.NReconstructed()++;
707 if( mc.Set()== 0 ) nRecOut++;
709 if( mc.NReconstructed()==1 ) nRecAll++;
710 else if(mc.NReconstructed() > mc.NTurns() ) nClonesAll++;
712 if( mc.NReconstructed()==1 ) nRecRef++;
713 else if(mc.NReconstructed() > mc.NTurns() ) nClonesRef++;
718 while( TMath::Abs(mc.TPCPar()[0]) + TMath::Abs(mc.TPCPar()[1])>1 ){
719 if( traPurity[itr]<.90 ) break;
720 AliHLTTPCCAGBTrack &t = fTracker->Tracks()[itr];
721 AliHLTTPCCATrackParam p = t.Param();
722 Double_t cosA = TMath::Cos( t.Alpha() );
723 Double_t sinA = TMath::Sin( t.Alpha() );
724 Double_t mcX = mc.TPCPar()[0]*cosA + mc.TPCPar()[1]*sinA;
725 Double_t mcY = -mc.TPCPar()[0]*sinA + mc.TPCPar()[1]*cosA;
726 Double_t mcZ = mc.TPCPar()[2];
727 Double_t mcEx = mc.TPCPar()[3]*cosA + mc.TPCPar()[4]*sinA;
728 Double_t mcEy = -mc.TPCPar()[3]*sinA + mc.TPCPar()[4]*cosA;
729 Double_t mcEz = mc.TPCPar()[5];
730 Double_t mcEt = TMath::Sqrt(mcEx*mcEx + mcEy*mcEy);
731 if( TMath::Abs(mcEt)<1.e-4 ) break;
732 Double_t mcSinPhi = mcEy / mcEt;
733 Double_t mcDzDs = mcEz / mcEt;
734 Double_t mcQPt = mc.TPCPar()[6]/ mcEt;
735 if( TMath::Abs(mcQPt)<1.e-4 ) break;
736 Double_t mcPt = 1./TMath::Abs(mcQPt);
738 if( t.NHits() < 50 ) break;
739 Double_t bz = fTracker->Slices()[0].Param().Bz();
740 if( !p.TransportToXWithMaterial( mcX, bz ) ) break;
741 if( p.GetCosPhi()*mcEx < 0 ){ // change direction
742 mcSinPhi = -mcSinPhi;
746 const Double_t kCLight = 0.000299792458;
747 Double_t k2QPt = 100;
748 if( TMath::Abs(bz)>1.e-4 ) k2QPt= 1./(bz*kCLight);
749 Double_t qPt = p.GetKappa()*k2QPt;
751 if( TMath::Abs(qPt) >1.e-4 ) pt = 1./TMath::Abs(qPt);
753 fhResY->Fill( p.GetY() - mcY );
754 fhResZ->Fill( p.GetZ() - mcZ );
755 fhResSinPhi->Fill( p.GetSinPhi() - mcSinPhi );
756 fhResDzDs->Fill( p.GetDzDs() - mcDzDs );
757 fhResPt->Fill( ( pt - mcPt )/mcPt );
759 if( p.GetErr2Y()>0 ) fhPullY->Fill( (p.GetY() - mcY)/TMath::Sqrt(p.GetErr2Y()) );
760 if( p.GetErr2Z()>0 ) fhPullZ->Fill( (p.GetZ() - mcZ)/TMath::Sqrt(p.GetErr2Z()) );
761 if( p.GetErr2SinPhi()>0 ) fhPullSinPhi->Fill( (p.GetSinPhi() - mcSinPhi)/TMath::Sqrt(p.GetErr2SinPhi()) );
762 if( p.GetErr2DzDs()>0 ) fhPullDzDs->Fill( (p.DzDs() - mcDzDs)/TMath::Sqrt(p.GetErr2DzDs()) );
763 if( p.GetErr2Kappa()>0 ) fhPullQPt->Fill( (qPt - mcQPt)/TMath::Sqrt(p.GetErr2Kappa()*k2QPt*k2QPt) );
768 for (Int_t ipart=0; ipart<fNMCTracks; ipart++) {
769 AliHLTTPCCAMCTrack &mc = fMCTracks[ipart];
770 if( mc.Set()>0 ) fhGBEffVsP->Fill(mc.P(), ( mc.NReconstructed()>0 ?1 :0));
773 if( traLabels ) delete[] traLabels;
774 if( traPurity ) delete[] traPurity;
776 fStatGBNRecTot += nRecTot;
777 fStatGBNRecOut += nRecOut;
778 fStatGBNGhost += nGhost;
779 fStatGBNMCAll += nMCAll;
780 fStatGBNRecAll += nRecAll;
781 fStatGBNClonesAll += nClonesAll;
782 fStatGBNMCRef += nMCRef;
783 fStatGBNRecRef += nRecRef;
784 fStatGBNClonesRef += nClonesRef;
788 // distribution of cluster errors
791 Int_t nHits = fTracker->NHits();
792 for( Int_t ih=0; ih<nHits; ih++ ){
793 AliHLTTPCCAGBHit &hit = fTracker->Hits()[ih];
794 fhHitErrY->Fill(hit.ErrY());
795 fhHitErrZ->Fill(hit.ErrZ());
801 if( fDoClusterPulls && fNMCPoints>0 ) {
804 for (Int_t ipart=0; ipart<fNMCTracks; ipart++) {
805 AliHLTTPCCAMCTrack &mc = fMCTracks[ipart];
808 sort(fMCPoints, fMCPoints+fNMCPoints, AliHLTTPCCAMCPoint::Compare );
810 for( Int_t ip=0; ip<fNMCPoints; ip++ ){
811 AliHLTTPCCAMCPoint &p = fMCPoints[ip];
812 AliHLTTPCCAMCTrack &t = fMCTracks[p.TrackID()];
813 if( t.NMCPoints()==0 ) t.FirstMCPointID() = ip;
818 for( Int_t ih=0; ih<fNHits; ih++ ){
820 AliHLTTPCCAGBHit &hit = fTracker->Hits()[ih];
821 AliHLTTPCCAHitLabel &l = fHitLabels[ih];
823 if( l.fLab[0]<0 || l.fLab[0]>=fNMCTracks
824 || l.fLab[1]>=0 || l.fLab[2]>=0 ) continue;
826 Int_t lab = l.fLab[0];
828 AliHLTTPCCAMCTrack &track = fMCTracks[lab];
829 //if( track.Pt()<1. ) continue;
830 Int_t ip1=-1, ip2=-1;
831 Double_t d1 = 1.e20, d2=1.e20;
832 for( Int_t ip=0; ip<track.NMCPoints(); ip++ ){
833 AliHLTTPCCAMCPoint &p = fMCPoints[track.FirstMCPointID() + ip];
834 if( p.ISlice() != hit.ISlice() ) continue;
835 Double_t dx = p.Sx()-hit.X();
836 Double_t dy = p.Sy()-hit.Y();
837 Double_t dz = p.Sz()-hit.Z();
838 Double_t d = dx*dx + dy*dy + dz*dz;
839 if( p.Sx()< hit.X() ){
852 if( ip1<0 || ip2<0 ) continue;
854 AliHLTTPCCAMCPoint &p1 = fMCPoints[track.FirstMCPointID() + ip1];
855 AliHLTTPCCAMCPoint &p2 = fMCPoints[track.FirstMCPointID() + ip2];
856 Double_t dx = p2.Sx() - p1.Sx();
857 Double_t dy = p2.Sy() - p1.Sy();
858 Double_t dz = p2.Sz() - p1.Sz();
859 if( TMath::Abs(dx)>1.e-8 && TMath::Abs(p1.Sx()-hit.X())<2. && TMath::Abs(p2.Sx()-hit.X())<2. ){
860 Double_t sx = hit.X();
861 Double_t sy = p1.Sy() + dy/dx*(sx-p1.Sx());
862 Double_t sz = p1.Sz() + dz/dx*(sx-p1.Sx());
866 AliHLTTPCCATrackParam t;
868 t.SinPhi() = dy/TMath::Sqrt(dx*dx+dy*dy);
869 t.CosPhi() = dx/TMath::Sqrt(dx*dx+dy*dy);
870 t.DzDs() = dz/TMath::Sqrt(dx*dx+dy*dy);
871 fTracker->GetErrors2(hit,t,errY, errZ );
872 errY = TMath::Sqrt(errY);
873 errZ = TMath::Sqrt(errZ);
876 fhHitResY->Fill((hit.Y()-sy));
877 fhHitResZ->Fill((hit.Z()-sz));
878 fhHitPullY->Fill((hit.Y()-sy)/errY);
879 fhHitPullZ->Fill((hit.Z()-sz)/errZ);
880 if( track.Pt()>=1. ){
881 fhHitResY1->Fill((hit.Y()-sy));
882 fhHitResZ1->Fill((hit.Z()-sz));
883 fhHitPullY1->Fill((hit.Y()-sy)/errY);
884 fhHitPullZ1->Fill((hit.Z()-sz)/errZ);
891 cout<<"\nSlice tracker performance: \n"<<endl;
893 <<fStatNMCAll/fStatNEvents<<" mc all, "
894 <<fStatNMCRef/fStatNEvents<<" mc ref, "
895 <<fStatNRecTot/fStatNEvents<<" rec total, "
896 <<fStatNRecAll/fStatNEvents<<" rec all, "
897 <<fStatNClonesAll/fStatNEvents<<" clones all, "
898 <<fStatNRecRef/fStatNEvents<<" rec ref, "
899 <<fStatNClonesRef/fStatNEvents<<" clones ref, "
900 <<fStatNRecOut/fStatNEvents<<" out, "
901 <<fStatNGhost/fStatNEvents<<" ghost"<<endl;
903 Int_t nRecExtr = fStatNRecAll - fStatNRecRef;
904 Int_t nMCExtr = fStatNMCAll - fStatNMCRef;
905 Int_t nClonesExtr = fStatNClonesAll - fStatNClonesRef;
907 Double_t dRecTot = (fStatNRecTot>0 ) ? fStatNRecTot :1;
908 Double_t dMCAll = (fStatNMCAll>0 ) ? fStatNMCAll :1;
909 Double_t dMCRef = (fStatNMCRef>0 ) ? fStatNMCRef :1;
910 Double_t dMCExtr = (nMCExtr>0 ) ? nMCExtr :1;
911 Double_t dRecAll = (fStatNRecAll+fStatNClonesAll>0 ) ? fStatNRecAll+fStatNClonesAll :1;
912 Double_t dRecRef = (fStatNRecRef+fStatNClonesRef>0 ) ? fStatNRecRef+fStatNClonesRef :1;
913 Double_t dRecExtr = (nRecExtr+nClonesExtr>0 ) ? nRecExtr+nClonesExtr :1;
915 cout<<" EffRef = "<< fStatNRecRef/dMCRef
916 <<", CloneRef = " << fStatNClonesRef/dRecRef <<endl;
917 cout<<" EffExtra = "<<nRecExtr/dMCExtr
918 <<", CloneExtra = " << nClonesExtr/dRecExtr<<endl;
919 cout<<" EffAll = "<<fStatNRecAll/dMCAll
920 <<", CloneAll = " << fStatNClonesAll/dRecAll<<endl;
921 cout<<" Out = "<<fStatNRecOut/dRecTot
922 <<", Ghost = "<<fStatNGhost/dRecTot<<endl;
923 cout<<" Time = "<<fTracker->StatTime(0)/fTracker->StatNEvents()*1.e3<<" msec/event "<<endl;
924 cout<<" Local timers = "
925 <<fTracker->StatTime(1)/fTracker->StatNEvents()*1.e3<<" "
926 <<fTracker->StatTime(2)/fTracker->StatNEvents()*1.e3<<" "
927 <<fTracker->StatTime(3)/fTracker->StatNEvents()*1.e3<<" "
928 <<fTracker->StatTime(4)/fTracker->StatNEvents()*1.e3<<" "
929 <<fTracker->StatTime(5)/fTracker->StatNEvents()*1.e3<<"["
930 <<fTracker->StatTime(6)/fTracker->StatNEvents()*1.e3<<"/"
931 <<fTracker->StatTime(7)/fTracker->StatNEvents()*1.e3<<"] "
932 <<fTracker->StatTime(8)/fTracker->StatNEvents()*1.e3<<" "
933 <<" msec/event "<<endl;
937 cout<<"\nGlobal tracker performance: \n"<<endl;
939 <<fStatGBNMCAll/fStatNEvents<<" mc all, "
940 <<fStatGBNMCRef/fStatNEvents<<" mc ref, "
941 <<fStatGBNRecTot/fStatNEvents<<" rec total, "
942 <<fStatGBNRecAll/fStatNEvents<<" rec all, "
943 <<fStatGBNClonesAll/fStatNEvents<<" clones all, "
944 <<fStatGBNRecRef/fStatNEvents<<" rec ref, "
945 <<fStatGBNClonesRef/fStatNEvents<<" clones ref, "
946 <<fStatGBNRecOut/fStatNEvents<<" out, "
947 <<fStatGBNGhost/fStatNEvents<<" ghost"<<endl;
949 Int_t nRecExtr = fStatGBNRecAll - fStatGBNRecRef;
950 Int_t nMCExtr = fStatGBNMCAll - fStatGBNMCRef;
951 Int_t nClonesExtr = fStatGBNClonesAll - fStatGBNClonesRef;
953 Double_t dRecTot = (fStatGBNRecTot>0 ) ? fStatGBNRecTot :1;
954 Double_t dMCAll = (fStatGBNMCAll>0 ) ? fStatGBNMCAll :1;
955 Double_t dMCRef = (fStatGBNMCRef>0 ) ? fStatGBNMCRef :1;
956 Double_t dMCExtr = (nMCExtr>0 ) ? nMCExtr :1;
957 Double_t dRecAll = (fStatGBNRecAll+fStatGBNClonesAll>0 ) ? fStatGBNRecAll+fStatGBNClonesAll :1;
958 Double_t dRecRef = (fStatGBNRecRef+fStatGBNClonesRef>0 ) ? fStatGBNRecRef+fStatGBNClonesRef :1;
959 Double_t dRecExtr = (nRecExtr+nClonesExtr>0 ) ? nRecExtr+nClonesExtr :1;
961 cout<<" EffRef = "<< fStatGBNRecRef/dMCRef
962 <<", CloneRef = " << fStatGBNClonesRef/dRecRef <<endl;
963 cout<<" EffExtra = "<<nRecExtr/dMCExtr
964 <<", CloneExtra = " << nClonesExtr/dRecExtr<<endl;
965 cout<<" EffAll = "<<fStatGBNRecAll/dMCAll
966 <<", CloneAll = " << fStatGBNClonesAll/dRecAll<<endl;
967 cout<<" Out = "<<fStatGBNRecOut/dRecTot
968 <<", Ghost = "<<fStatGBNGhost/dRecTot<<endl;
969 cout<<" Time = "<<( fTracker->StatTime(0)+fTracker->StatTime(9) )/fTracker->StatNEvents()*1.e3<<" msec/event "<<endl;
970 cout<<" Local timers: "<<endl;
971 cout<<" slice tracker "<<fTracker->StatTime(0)/fTracker->StatNEvents()*1.e3<<": "
972 <<fTracker->StatTime(1)/fTracker->StatNEvents()*1.e3<<" "
973 <<fTracker->StatTime(2)/fTracker->StatNEvents()*1.e3<<" "
974 <<fTracker->StatTime(3)/fTracker->StatNEvents()*1.e3<<" "
975 <<fTracker->StatTime(4)/fTracker->StatNEvents()*1.e3<<" "
976 <<fTracker->StatTime(5)/fTracker->StatNEvents()*1.e3<<"["
977 <<fTracker->StatTime(6)/fTracker->StatNEvents()*1.e3<<"/"
978 <<fTracker->StatTime(7)/fTracker->StatNEvents()*1.e3<<"] "
979 <<fTracker->StatTime(8)/fTracker->StatNEvents()*1.e3
980 <<" msec/event "<<endl;
981 cout<<" GB merger "<<fTracker->StatTime(9)/fTracker->StatNEvents()*1.e3<<": "
982 <<fTracker->StatTime(10)/fTracker->StatNEvents()*1.e3<<", "
983 <<fTracker->StatTime(11)/fTracker->StatNEvents()*1.e3<<", "
984 <<fTracker->StatTime(12)/fTracker->StatNEvents()*1.e3<<" "
985 <<" msec/event "<<endl;
991 void AliHLTTPCCAPerformance::WriteMCEvent( ostream &out )
993 out<<fNMCTracks<<endl;
994 for( Int_t it=0; it<fNMCTracks; it++ ){
995 AliHLTTPCCAMCTrack &t = fMCTracks[it];
998 for( Int_t i=0; i<7; i++ ) out<<t.Par()[i]<<" ";
1000 for( Int_t i=0; i<7; i++ ) out<<t.TPCPar()[i]<<" ";
1004 out<< t.NMCPoints()<<" ";
1005 out<< t.FirstMCPointID()<<" ";
1006 out<< t.NHits()<<" ";
1007 out<< t.NReconstructed()<<" ";
1009 out<< t.NTurns()<<endl;
1013 for( Int_t ih=0; ih<fNHits; ih++ ){
1014 AliHLTTPCCAHitLabel &l = fHitLabels[ih];
1015 out<<l.fLab[0]<<" "<<l.fLab[1]<<" "<<l.fLab[2]<<endl;
1019 void AliHLTTPCCAPerformance::WriteMCPoints( ostream &out )
1021 out<<fNMCPoints<<endl;
1022 for( Int_t ip=0; ip<fNMCPoints; ip++ ){
1023 AliHLTTPCCAMCPoint &p = fMCPoints[ip];
1030 out<< p.Time()<<" ";
1031 out<< p.ISlice()<<" ";
1032 out<< p.TrackID()<<endl;
1036 void AliHLTTPCCAPerformance::ReadMCEvent( istream &in )
1039 if( fMCTracks ) delete[] fMCTracks;
1042 if( fHitLabels ) delete[] fHitLabels;
1045 if( fMCPoints ) delete[] fMCPoints;
1050 fMCTracks = new AliHLTTPCCAMCTrack[fNMCTracks];
1051 for( Int_t it=0; it<fNMCTracks; it++ ){
1052 AliHLTTPCCAMCTrack &t = fMCTracks[it];
1056 for( Int_t i=0; i<7; i++ ) in>>t.Par()[i];
1057 for( Int_t i=0; i<7; i++ ) in>>t.TPCPar()[i];
1062 in>> t.FirstMCPointID();
1063 in>> t.NReconstructed();
1069 fHitLabels = new AliHLTTPCCAHitLabel[fNHits];
1070 for( Int_t ih=0; ih<fNHits; ih++ ){
1071 AliHLTTPCCAHitLabel &l = fHitLabels[ih];
1072 in>>l.fLab[0]>>l.fLab[1]>>l.fLab[2];
1076 void AliHLTTPCCAPerformance::ReadMCPoints( istream &in )
1078 if( fMCPoints ) delete[] fMCPoints;
1083 fMCPoints = new AliHLTTPCCAMCPoint[fNMCPoints];
1084 for( Int_t ip=0; ip<fNMCPoints; ip++ ){
1085 AliHLTTPCCAMCPoint &p = fMCPoints[ip];