StRoot: Sti/StiTPCCATrackerInterface.cxx Source File

00001 #ifdef DO_TPCCATRACKER
00002 #include "StiTPCCATrackerInterface.h"
00003 
00004 #include "TPCCATracker/code/AliHLTTPCCAGBHit.h"
00005 #include "TPCCATracker/code/AliHLTTPCCAGBTrack.h"
00006 #include "TPCCATracker/code/AliHLTTPCCAParam.h"
00007 
00008   // need for hits data
00009 #include "StTpcHit.h"                
00010 #include "StTpcDb/StTpcDb.h"
00011 #include "StDbUtilities/StTpcCoordinateTransform.hh"
00012 #include "StDbUtilities/StTpcLocalSectorCoordinate.hh"
00013   // for MCdata
00014 #include "tables/St_g2t_track_Table.h" 
00015 #include "tables/St_g2t_tpc_hit_Table.h"
00016 #include "TDatabasePDG.h"
00017 #include "StBFChain.h"
00018 #include "StiDetectorBuilder.h"
00019 #include "StiDetectorGroups.h"
00020 #include "StiGenericDetectorGroup.h"
00021 #include "StiToolkit.h"
00022   //to obtain error coefficients
00023 #include "StDetectorDbMaker/StiTpcInnerHitErrorCalculator.h"
00024 #include "StDetectorDbMaker/StiTpcOuterHitErrorCalculator.h"
00025   //to get Magnetic Field
00026 #include "StarMagField/StarMagField.h"
00027 
00028   // for sti perfo
00029 #include "TPCCATracker/code/AliHLTTPCCAStiPerformance.h"
00030 #include "TPCCATracker/code/AliHLTTPCCAMergerPerformance.h"
00031 #include "StDetectorDbMaker/St_tpcPadPlanesC.h"
00032 #include "StiKalmanTrack.h"
00033 
00034 #include <vector>
00035 #include <algorithm>
00036 using std::vector;
00037 
00038 #include <string>
00039 using std::string;
00040 
00041 
00042 //#define STORE_STANDALONE_DATA // write data in files for Standalone
00043 
00044 StiTPCCATrackerInterface &StiTPCCATrackerInterface::Instance()
00045 {
00046     // reference to static object
00047   static StiTPCCATrackerInterface g;
00048   return g;
00049 }
00050 
00051 StiTPCCATrackerInterface::StiTPCCATrackerInterface()
00052 {
00053   fOutFile = StMaker::GetChain()->GetTFile();
00054   if(!fOutFile) cout << "W StiTPCCATrackerInterface: Warning - There isn't any tag file, so histograms won't be saved!" << endl;
00055 
00056   fPerformance = &(AliHLTTPCCAPerformance::Instance());
00057   fPerformance->SetOutputFile(fOutFile);
00058 
00059   SetNewEvent();
00060 } // StiTPCCATrackerInterface::StiTPCCATrackerInterface()
00061 
00062 StiTPCCATrackerInterface::~StiTPCCATrackerInterface(  )
00063 { // never called for static object
00064 } // StiTPCCATrackerInterface::StiTPCCATrackerInterface()
00065 
00066 void StiTPCCATrackerInterface::SetNewEvent()
00067 {
00068   assert(fPerformance != 0);
00069   
00070   fHitsMap = 0;
00071   fSeeds.clear();
00072 
00073   fSeedFinder = 0;
00074   fStiTracks = 0;
00075 
00076   fIdTruth.clear(); // id of the Track, which has created CaHit
00077   fMCTracks.clear();
00078   fMCPoints.clear();
00079   fHitLabels.clear();
00080 
00081   fCaParam.clear();// settings for all sectors to give CATracker
00082   fCaHits.clear(); // hits to give CATracker
00083   fSeedHits.clear();          // hits to make seeds
00084 
00085   if (!fSeedFinder) fSeedFinder = new StiTpcSeedFinder;
00086   
00087   if (fTracker)    delete fTracker;
00088   fTracker    = new AliHLTTPCCAGBTracker;
00089   if (fStiTracker) delete fStiTracker;
00090   fStiTracker = new AliHLTTPCCAGBTracker; 
00091   
00092   fPerformance->SetTracker(fTracker);
00093 }
00094 
00095 
00097 void StiTPCCATrackerInterface::Run()
00098 {
00099   assert(fHitsMap != 0);
00100 
00101   TStopwatch timer;
00102   timer.Start();
00103 
00104   MakeSettings();
00105   MakeHits();
00106 
00107 #ifdef DO_TPCCATRACKER_EFF_PERFORMANCE
00108     // put data in performance
00109   FillPerformance(fCaHits,fIdTruth, fMCTracks,fMCPoints,fHitLabels);
00110   fPerformance->SetMCTracks(fMCTracks);
00111   fPerformance->SetMCPoints(fMCPoints);
00112   fPerformance->SetHitLabels(fHitLabels);
00113 #endif //  DO_TPCCATRACKER_EFF_PERFORMANCE
00114 
00115   
00116     // run tracker
00117   fTracker->SetSettings(fCaParam);
00118   fTracker->SetHits(fCaHits);
00119   
00120 #ifdef STORE_STANDALONE_DATA // write data in files for Standalone
00121   static int iEvent = -1;
00122   iEvent++;
00123   TString name = "./data/";
00124   if (iEvent == 0) fTracker->SaveSettingsInFile(string(name));
00125   name += "event";
00126   name += iEvent;
00127   name += "_";
00128   fTracker->SaveHitsInFile(string(name));
00129   fPerformance->SaveDataInFiles(string(name));
00130 
00131 // check
00132   if(1){
00133   if (fTracker)    delete fTracker;
00134   fTracker    = new AliHLTTPCCAGBTracker;
00135   fPerformance->SetTracker(fTracker);
00136   TString name = "./data/";
00137   fTracker->ReadSettingsFromFile(string(name));
00138   name += "event";
00139   name += iEvent;
00140   name += "_";
00141   fTracker->ReadHitsFromFile(string(name));
00142   fPerformance->ReadDataFromFiles(string(name));
00143   fTracker->SetSettings(fCaParam);
00144   fTracker->SetHits(fCaHits);
00145 
00146   }
00147 #endif // STORE_STANDALONE_DATA
00148 
00149 
00150 
00151   timer.Stop();
00152   fPreparationTime_real = timer.RealTime();
00153   fPreparationTime_cpu = timer.CpuTime();  
00154   
00155   fTracker->FindTracks();
00156 
00157     // copy hits
00158   timer.Start();
00159   
00160   MakeSeeds();
00161   
00162   timer.Stop();
00163   fPreparationTime_real += timer.RealTime();
00164   fPreparationTime_cpu += timer.CpuTime();    
00165 } // void StiTPCCATrackerInterface::Run()
00166 
00167 void StiTPCCATrackerInterface::RunPerformance()
00168 {
00169   cout << " ---- CA TPC Tracker ---- " << endl;
00170 #ifdef DO_TPCCATRACKER_EFF_PERFORMANCE
00171   TStopwatch timer;
00172   timer.Start();
00173     
00174   assert(fStiTracks  != 0);
00175   assert(fStiTracker != 0);
00176 
00177   // ----- Efficiency ---------
00178   
00179     // init all performances
00180   fPerformance->InitSubPerformances();
00181 
00182     // -- correct init of stiPerformance --
00183 
00184   FillStiPerformance(); // get fStiCaHits fStiIdTruth fStiCaTracks
00185   
00186     // prepare tracker
00187   fStiTracker->SetHits(fStiCaHits);
00188   fStiTracker->SetSettings(fCaParam);
00189 
00190   AliHLTTPCCAStiPerformance* stiPerf = (AliHLTTPCCAStiPerformance*)fPerformance->GetSubPerformance("Sti Performance");
00191   stiPerf->SetNewEvent(fStiTracker,&fStiHitLabels,&fStiMCTracks,&fStiMCPoints);
00192   stiPerf->SetRecoTracks(fStiCaTracks);  
00193   fPerformance->ExecPerformance();
00194 
00195   timer.Stop();
00196   fPreparationTime_real += timer.RealTime();
00197   fPreparationTime_cpu += timer.CpuTime();   
00198 #endif //  DO_TPCCATRACKER_EFF_PERFORMANCE
00199 
00200   // ----- Timing ---------
00201   //  TODO saparete in procedure in TPCCATracker/ ...
00202 #if 1
00203   const int printTime = 3; // 0 - show nothing; 1 - show avarege; 2 - show cur event; 3 - show both
00204   const int fullTiming = 1; // show info about each stage of tracking; same as printTime
00205   if ((printTime == 2) || (printTime == 3)){
00206         
00207     std::cout << "Sector reconstruction Time"
00208               << " Real = " << std::setw( 10 ) << fTracker->SliceTrackerTime() * 1.e3 << " ms,"
00209               << " CPU = " << std::setw( 10 ) << fTracker->SliceTrackerCpuTime() * 1.e3 << " ms,"
00210               << " parallelization speedup: " << fTracker->SliceTrackerCpuTime() / fTracker->SliceTrackerTime()
00211               << std::endl;
00212     if ((fullTiming == 2) || (fullTiming == 3)) {
00213       std::cout
00214         << " |  sum slice trackers: " << std::setw( 10 ) << fTracker->StatTime( 0 ) * 1000. << " ms\n"
00215         << " |    NeighboursFinder: " << std::setw( 10 ) << fTracker->StatTime( 1 ) * 1000. << " ms, " << std::setw( 12 ) << fTracker->StatTime( 5 ) << " cycles\n"
00216         << " |     StartHitsFinder: " << std::setw( 10 ) << fTracker->StatTime( 4 ) * 1000. << " ms\n"
00217         << " | TrackletConstructor: " << std::setw( 10 ) << fTracker->StatTime( 2 ) * 1000. << " ms, " << std::setw( 12 ) << fTracker->StatTime( 7 ) << " cycles\n"
00218         << " |    TrackletSelector: " << std::setw( 10 ) << fTracker->StatTime( 3 ) * 1000. << " ms, " << std::setw( 12 ) << fTracker->StatTime( 8 ) << " cycles\n"
00219         << " |         WriteOutput: " << std::setw( 10 ) << fTracker->StatTime( 6 ) * 1000. << " ms\n"
00220         << " |               merge: " << std::setw( 10 ) << fTracker->StatTime( 9 ) * 1000. << " ms" << std::endl;
00221     }
00222     std::cout << "Total (sector+merge) reconstuction time"
00223               << " Real = " << std::setw( 10 ) << (fTracker->SliceTrackerTime()+fTracker->StatTime( 9 )) * 1.e3 << " ms,"
00224               << " CPU = " << std::setw( 10 ) << (fTracker->SliceTrackerCpuTime()+fTracker->StatTime( 10 )) * 1.e3 << " ms"
00225               << std::endl;
00226     std::cout << "Preparation time"
00227               << " Real = " << std::setw( 10 ) << fPreparationTime_real * 1.e3 << " ms,"
00228               << " CPU = " << std::setw( 10 ) << fPreparationTime_cpu * 1.e3 << " ms"
00229               << std::endl;
00230   };
00231   if ((printTime == 1) || (printTime == 3)){
00232     const int NTimers = 20;
00233     static int statIEvent = 0;
00234     static double statTime[NTimers];
00235     static double statTime_SliceTrackerTime = 0;
00236     static double statTime_SliceTrackerCpuTime = 0;
00237     static double statPreparationTime_real = 0;
00238     static double statPreparationTime_cpu = 0;
00239     
00240     if (!statIEvent){
00241       for (int i = 0; i < NTimers; i++){
00242         statTime[i] = 0;
00243       }
00244     }
00245 
00246     statIEvent++;
00247     for (int i = 0; i < NTimers; i++){
00248       statTime[i] += fTracker->StatTime( i );
00249     }
00250     statTime_SliceTrackerTime += fTracker->SliceTrackerTime();
00251     statTime_SliceTrackerCpuTime += fTracker->SliceTrackerCpuTime();
00252     statPreparationTime_real += fPreparationTime_real;
00253     statPreparationTime_cpu  += fPreparationTime_cpu;
00254         
00255     std::cout << "Avarage sector reconstruction Time"
00256               << " Real = " << std::setw( 10 ) << 1./statIEvent*statTime_SliceTrackerTime * 1.e3 << " ms,"
00257               << " CPU = " << std::setw( 10 ) << 1./statIEvent*statTime_SliceTrackerCpuTime * 1.e3 << " ms,"
00258               << " parallelization speedup: " << statTime_SliceTrackerCpuTime / statTime_SliceTrackerTime
00259               << std::endl;
00260     if ((fullTiming == 1) || (fullTiming == 3)) {
00261       std::cout
00262         << " |  sum slice trackers: " << std::setw( 10 ) << 1./statIEvent*statTime[ 0 ] * 1000. << " ms\n"
00263         << " |    NeighboursFinder: " << std::setw( 10 ) << 1./statIEvent*statTime[ 1 ] * 1000. << " ms, " << std::setw( 12 ) << 1./statIEvent*statTime[ 5 ] << " cycles\n"
00264         << " |     StartHitsFinder: " << std::setw( 10 ) << 1./statIEvent*statTime[ 4 ] * 1000. << " ms\n"
00265         << " | TrackletConstructor: " << std::setw( 10 ) << 1./statIEvent*statTime[ 2 ] * 1000. << " ms, " << std::setw( 12 ) << 1./statIEvent*statTime[ 7 ] << " cycles\n"
00266         << " |    TrackletSelector: " << std::setw( 10 ) << 1./statIEvent*statTime[ 3 ] * 1000. << " ms, " << std::setw( 12 ) << 1./statIEvent*statTime[ 8 ] << " cycles\n"
00267         << " |         WriteOutput: " << std::setw( 10 ) << 1./statIEvent*statTime[ 6 ] * 1000. << " ms\n"
00268         << " |               merge: " << std::setw( 10 ) << 1./statIEvent*statTime[ 9 ] * 1000. << " ms" << std::endl;
00269     }
00270     std::cout << "Total (sector+merge) avarage reconstuction time" 
00271               << " Real = " << std::setw( 10 ) << (statTime_SliceTrackerTime+statTime[ 9 ])/statIEvent * 1.e3 << " ms,"
00272               << " CPU = " << std::setw( 10 ) << (statTime_SliceTrackerCpuTime+statTime[ 10 ])/statIEvent * 1.e3 << " ms"
00273               << std::endl;
00274     std::cout << "Avarage preparation time"
00275               << " Real = " << std::setw( 10 ) << statPreparationTime_real/statIEvent * 1.e3 << " ms,"
00276               << " CPU = " << std::setw( 10 ) << statPreparationTime_cpu/statIEvent * 1.e3 << " ms."
00277               << std::endl;
00278   }
00279 #endif // 0 timing
00280   
00281 #ifdef DO_TPCCATRACKER_EFF_PERFORMANCE  
00282   ((AliHLTTPCCAMergerPerformance*)(AliHLTTPCCAPerformance::Instance().GetSubPerformance("Merger")))->FillTree();
00283 #endif //DO_TPCCATRACKER_EFF_PERFORMANCE
00284 } // void StiTPCCATrackerInterface::Run()
00285 
00286 
00287 void StiTPCCATrackerInterface::MakeSettings()
00288 {
00289 
00290   const int NSlices = 24; //TODO initialize from StRoot
00291   const int NoOfInnerRows = St_tpcPadPlanesC::instance()->innerPadRows();
00292   const int NRows = St_tpcPadPlanesC::instance()->padRows();
00293   for ( int iSlice = 0; iSlice < NSlices; iSlice++ ) {
00294     AliHLTTPCCAParam SlicePar;
00295     memset(&SlicePar, 0, sizeof(AliHLTTPCCAParam));
00296 
00297     Int_t sector = iSlice+1;
00298       // Int_t sector = iSlice;
00299     SlicePar.SetISlice( iSlice );
00300     SlicePar.SetNRows ( NRows ); 
00301     Double_t beta = 0;
00302     if (sector > 12) beta = (24-sector)*2.*TMath::Pi()/12.;
00303     else             beta =     sector *2.*TMath::Pi()/12.;
00304     SlicePar.SetAlpha  ( beta );
00305     SlicePar.SetDAlpha  ( 30*TMath::DegToRad() );                        //TODO initialize from StRoot
00306     SlicePar.SetCosAlpha ( TMath::Cos(SlicePar.Alpha()) );
00307     SlicePar.SetSinAlpha ( TMath::Sin(SlicePar.Alpha()) );
00308     SlicePar.SetAngleMin ( SlicePar.Alpha() - 0.5*SlicePar.DAlpha() );
00309     SlicePar.SetAngleMax ( SlicePar.Alpha() + 0.5*SlicePar.DAlpha() );
00310     SlicePar.SetRMin     (  51. );                                        //TODO initialize from StRoot
00311     SlicePar.SetRMax     ( 194. );                                        //TODO initialize from StRoot
00312     SlicePar.SetErrX     (   0. );                                        //TODO initialize from StRoot
00313     SlicePar.SetErrY     (   0.12 ); // 0.06  for Inner                        //TODO initialize from StRoot
00314     SlicePar.SetErrZ     (   0.16 ); // 0.12  for Inner                NodePar->fitPars()        //TODO initialize from StRoot
00315       //   SlicePar.SetPadPitch (   0.675 );// 0.335 -"-
00316     float x[3]={0,0,0},b[3];
00317     StarMagField::Instance()->BField(x,b);
00318     SlicePar.SetBz       ( - b[2] );   // change sign because change z
00319     if (sector <= 12) {
00320       SlicePar.SetZMin     (   0. );                                        //TODO initialize from StRoot
00321       SlicePar.SetZMax     ( 210. );                                        //TODO initialize from StRoot
00322     } else {
00323       SlicePar.SetZMin     (-210. );                                        //TODO initialize from StRoot
00324       SlicePar.SetZMax     (   0. );                                        //TODO initialize from StRoot
00325     }
00326     for( int iR = 0; iR < NRows; iR++){
00327       if (iR < NoOfInnerRows) {
00328         SlicePar.SetRowX(iR, St_tpcPadPlanesC::instance()->innerRowRadii()[iR]);
00329       } else {
00330         SlicePar.SetRowX(iR, St_tpcPadPlanesC::instance()->outerRowRadii()[iR-NoOfInnerRows]);
00331       }
00332     }
00333 
00334     Double_t *coeffInner = StiTpcInnerHitErrorCalculator::instance()->coeff();
00335     for(int iCoef=0; iCoef<6; iCoef++)
00336     {
00337       SlicePar.SetParamS0Par(0, 0, iCoef, (float)coeffInner[iCoef] );
00338     }  
00339   
00340     SlicePar.SetParamS0Par(0, 0, 6, 0.0f );
00341     
00342     Double_t *coeffOuter =StiTpcOuterHitErrorCalculator::instance()->coeff();
00343     for(int iCoef=0; iCoef<6; iCoef++)
00344     {
00345       SlicePar.SetParamS0Par(0, 1, iCoef, (float)coeffOuter[iCoef] );
00346     }
00347   
00348     SlicePar.SetParamS0Par(0, 1, 6, 0.0f );
00349     SlicePar.SetParamS0Par(0, 2, 0, 0.0f );
00350     SlicePar.SetParamS0Par(0, 2, 1, 0.0f );
00351     SlicePar.SetParamS0Par(0, 2, 2, 0.0f );
00352     SlicePar.SetParamS0Par(0, 2, 3, 0.0f );
00353     SlicePar.SetParamS0Par(0, 2, 4, 0.0f );
00354     SlicePar.SetParamS0Par(0, 2, 5, 0.0f );
00355     SlicePar.SetParamS0Par(0, 2, 6, 0.0f );
00356     SlicePar.SetParamS0Par(1, 0, 0, 0.0f );
00357     SlicePar.SetParamS0Par(1, 0, 1, 0.0f );
00358     SlicePar.SetParamS0Par(1, 0, 2, 0.0f );
00359     SlicePar.SetParamS0Par(1, 0, 3, 0.0f );
00360     SlicePar.SetParamS0Par(1, 0, 4, 0.0f );
00361     SlicePar.SetParamS0Par(1, 0, 5, 0.0f );
00362     SlicePar.SetParamS0Par(1, 0, 6, 0.0f );
00363     SlicePar.SetParamS0Par(1, 1, 0, 0.0f );
00364     SlicePar.SetParamS0Par(1, 1, 1, 0.0f );
00365     SlicePar.SetParamS0Par(1, 1, 2, 0.0f );
00366     SlicePar.SetParamS0Par(1, 1, 3, 0.0f );
00367     SlicePar.SetParamS0Par(1, 1, 4, 0.0f );
00368     SlicePar.SetParamS0Par(1, 1, 5, 0.0f );
00369     SlicePar.SetParamS0Par(1, 1, 6, 0.0f );
00370     SlicePar.SetParamS0Par(1, 2, 0, 0.0f );
00371     SlicePar.SetParamS0Par(1, 2, 1, 0.0f );
00372     SlicePar.SetParamS0Par(1, 2, 2, 0.0f );
00373     SlicePar.SetParamS0Par(1, 2, 3, 0.0f );
00374     SlicePar.SetParamS0Par(1, 2, 4, 0.0f );
00375     SlicePar.SetParamS0Par(1, 2, 5, 0.0f );
00376     SlicePar.SetParamS0Par(1, 2, 6, 0.0f );
00377     
00378     fCaParam.push_back(SlicePar);
00379   } // for iSlice
00380 } // void StiTPCCATrackerInterface::MakeSettings()
00381 
00382 
00383 void StiTPCCATrackerInterface::MakeHits()
00384 {
00385   StTpcCoordinateTransform tran(gStTpcDb);
00386   StTpcLocalSectorCoordinate loc;
00387   for (HitMapToVectorAndEndType::iterator it= fHitsMap->begin(); it!= fHitsMap->end(); ++it){
00388     vector<StiHit*>& tempvec = (*it).second.hits();
00389     vector<StiHit*>::iterator  start = tempvec.begin();
00390     vector<StiHit*>::iterator  stop  = tempvec.end();
00391     for (vector<StiHit*>::iterator cit = start; cit != stop; cit++) {
00392 
00393         // get local coordinates. take into account distortion
00394       StiHit *hit = *cit;
00395       if (! hit->stHit()) continue;
00396       const StTpcHit *tpcHit = dynamic_cast<const StTpcHit*>(hit->stHit());
00397       if ( ! tpcHit) continue;
00398       StGlobalCoordinate glob(tpcHit->position());
00399       tran(glob,loc,tpcHit->sector(),tpcHit->padrow());
00400 
00401         // obtain seed Hit
00402       SeedHit_t hitc;
00403       hitc.mMinPad  = tpcHit->minPad();
00404       hitc.mMaxPad  = tpcHit->maxPad();
00405       hitc.mMinTmbk = tpcHit->minTmbk();
00406       hitc.mMaxTmbk = tpcHit->maxTmbk();
00407       hitc.padrow = tpcHit->padrow()-1;
00408       hitc.x = loc.position().x();
00409       hitc.y = loc.position().y();
00410       hitc.z = loc.position().z();
00411       hitc.status=0;
00412       hitc.taken=0;
00413       hitc.track_key=tpcHit->idTruth();
00414       hitc.hit  = hit;
00415       fSeedHits.push_back(hitc);
00416 
00417         // convert to CA Hit
00418       AliHLTTPCCAGBHit caHit;
00419       caHit.SetIRow( hitc.padrow );
00420 //      caHit.SetX( hit->x() );
00421       caHit.SetX( hit->position() ); // take position of the row
00422       caHit.SetY( - hit->y() );
00423       caHit.SetZ( - hit->z() );
00424         // caHit.SetErrX(   );
00425       caHit.SetErrY( 0.12 );// TODO: read parameters from somewhere 
00426       caHit.SetErrZ( 0.16 );
00427       caHit.SetISlice( tpcHit->sector() - 1 );
00428       caHit.SetIRow( hitc.padrow );
00429       caHit.SetID( fCaHits.size() );
00430       fIdTruth.push_back( hitc.track_key );
00431 
00432       fCaHits.push_back(caHit);
00433     }
00434   }
00435 
00436 } // void StiTPCCATrackerInterface::MakeHits()
00437 
00438 void StiTPCCATrackerInterface::ConvertPars(const AliHLTTPCCATrackParam& caPar, double _alpha, StiNodePars& nodePars, StiNodeErrs& nodeErrs)
00439 {
00440     // set jacobian integral coef
00441   double JI[5]; 
00442   JI[0] = -1.;                    // y
00443   JI[1] = -1.;                    // z
00444   JI[2] = -1.;         // eta
00445   JI[3] =  1.;         // ptin
00446   JI[4] = -1.;         // tanl
00447     // get parameters
00448   nodePars.x()    = caPar.GetX();
00449   nodePars.y()    = JI[0] * caPar.GetY();
00450   nodePars.z()    = JI[1] * caPar.GetZ();
00451   nodePars.eta()  = JI[2] * asin(caPar.GetSinPhi()); // (signed curvature)*(local Xc of helix axis - X current point on track)
00452   nodePars.ptin() = JI[3] * caPar.GetQPt();        // signed invert pt [sign = sign(-qB)]
00453   nodePars.tanl() = JI[4] * caPar.GetDzDs();         // tangent of the track momentum dip angle
00454   
00455     // get h & signB
00456   static const double EC = 2.99792458e-4, ZEROHZ = 2e-6;
00457 // #define USE_CA_FIELD
00458 #ifndef USE_CA_FIELD // use field in the point. Need this because of check in StiTrackNodeHelper::set()
00459   const double ca = cos(_alpha), sa = sin(_alpha); // code has been taken from StiKalmanTrackNode::getHz()
00460   double globalXYZ[3] = {
00461     ca * nodePars.x() - sa * nodePars.y(),
00462     sa * nodePars.x() + ca * nodePars.y(),
00463     nodePars.z()
00464   };
00465   double h[3];
00466   StarMagField::Instance()->BField( globalXYZ, h );
00467   double &h2 = h[2];
00468 #else  // these parameters have been obtained with that MF, so let's use it.
00469   double h2 = - fTracker->Slice(0).Param().Bz(); // change sign because change z
00470 #endif // 1
00471   h2 *= EC;
00472   if (fabs(h2) < ZEROHZ) h2 = 0;
00473 
00474     // get parameters. continue
00475   nodePars.hz() = h2;  // Z component magnetic field in units Pt(Gev) = Hz * RCurv(cm)
00476   nodePars.ready(); // set cosCA, sinCA & curv
00477 //std::cout << nodePars._ptin << std::endl;
00478     // set jacobian integral coef
00479   double J[5]; 
00480   J[0] = JI[0];                    // y
00481   J[1] = JI[1];                    // z
00482   J[2] = JI[2]/cos(nodePars.eta()); // eta
00483   J[3] = JI[3];                    // ptin
00484   J[4] = JI[4];                    // tanl
00485   
00486     // get cov matrises
00487   const float *caCov = caPar.GetCov();
00488   double nodeCov[15];
00489   for (int i1 = 0, i = 0; i1 < 5; i1++){
00490     for (int i2 = 0; i2 <= i1; i2++, i++){
00491       nodeCov[i] = J[i1]*J[i2]*caCov[i];
00492     }
00493   }
00494   // if ( (caCov[0] <= 0) || (caCov[2] <= 0) || (caCov[5] <= 0) || (caCov[9] <= 0) || (caCov[14] <= 0))
00495   //   cout << "Warrning: Bad CA Cov Matrix." << endl;
00496   // if ( (nodeCov[0] <= 0) || (nodeCov[2] <= 0) || (nodeCov[5] <= 0) || (nodeCov[9] <= 0) || (nodeCov[14] <= 0))
00497   //   cout << "Warrning: Bad Node Cov Matrix." << endl;
00498 
00499   double *A = nodeErrs.A;
00500 /*  for (int i1 = 0, i = 0; i1 < 5; i1++){
00501     for (int i2 = 0; i2 <= i1; i2++, i++){
00502       A[i+i1+2] = caCov[i];
00503     }
00504   }*/
00505   
00506   nodeErrs._cYY = caCov[ 0];
00507   nodeErrs._cZY = caCov[ 1]*J[0]*J[1];
00508   nodeErrs._cZZ = caCov[ 2]*J[0]*J[1];
00509   nodeErrs._cEY = caCov[ 3]*J[0]*J[2];
00510   nodeErrs._cEZ = caCov[ 4]*J[1]*J[2];
00511   nodeErrs._cEE = caCov[ 5]*J[2]*J[2];
00512   nodeErrs._cTY = caCov[ 6]*J[0]*J[4];
00513   nodeErrs._cTZ = caCov[ 7]*J[1]*J[4];
00514   nodeErrs._cTE = caCov[ 8]*J[2]*J[4];    
00515   nodeErrs._cTT = caCov[ 9]*J[4]*J[4];
00516   nodeErrs._cPY = caCov[10]*J[0]*J[3];
00517   nodeErrs._cPZ = caCov[11]*J[1]*J[3];
00518   nodeErrs._cPE = caCov[12]*J[2]*J[3];
00519   nodeErrs._cTP = caCov[13]*J[4]*J[3];
00520   nodeErrs._cPP = caCov[14]*J[3]*J[3];
00521   
00522   A[0] = 1; // don't use parameter X
00523   A[2] = 0;
00524   A[5] = 0;
00525   A[9] = 0;
00526   A[14] = 0;
00527   A[20] = 0;
00528 }
00529 
00530 void StiTPCCATrackerInterface::MakeSeeds()
00531 {
00532   const int NRecoTracks = fTracker->NTracks();
00533   for ( int iTr = 0; iTr < NRecoTracks; iTr++ ) {
00534       // get seed
00535     const AliHLTTPCCAGBTrack tr = fTracker->Track( iTr );
00536     Seed_t seed;
00537 
00538     const int NHits = tr.NHits();
00539     float last_x = 1e10; // for check
00540     for ( int iHit = NHits-1; iHit >= 0; iHit-- ){ 
00541       const int index = fTracker->TrackHit( tr.FirstHitRef() + iHit );
00542       const int hId   = fTracker->Hit( index ).ID();
00543 //      if ( last_x == fSeedHits[hId].hit->position() ) continue; // track can have 2 hits on 1 row because of track segments merger.
00544       seed.vhit.push_back(&(fSeedHits[hId]));
00545       assert( last_x >= fSeedHits[hId].hit->position() ); // should be back order - from outer to inner
00546       last_x = fSeedHits[hId].hit->position();
00547     }
00548 
00549     seed.total_hits = seed.vhit.size();
00550 
00551     ConvertPars( tr.OuterParam(), tr.Alpha(), seed.firstNodePars, seed.firstNodeErrs );
00552     ConvertPars( tr.InnerParam(), tr.Alpha(), seed.lastNodePars,  seed.lastNodeErrs );
00553 
00554     fSeeds.push_back(seed);
00555   }
00556 } // void StiTPCCATrackerInterface::MakeSeeds()
00557 
00558 bool myfunction (AliHLTTPCCALocalMCPoint i,AliHLTTPCCALocalMCPoint j) { 
00559   return (i.TrackI() < j.TrackI()) || ( i.TrackI()==j.TrackI() && i.IRow() < j.IRow() ); 
00560 }
00562 void StiTPCCATrackerInterface::FillPerformance(const vector<AliHLTTPCCAGBHit>& caHits, const vector<int>& idsTruth, vector<AliHLTTPCCAMCTrack>& mcTracks, vector<AliHLTTPCCALocalMCPoint>& mcPoints, vector<AliHLTTPCCAHitLabel>& hitLabels)
00563 {
00564   mcTracks.clear();
00565   mcPoints.clear();
00566   hitLabels.clear();
00567   hitLabels.resize(caHits.size());
00568 
00569   StBFChain *chain = (StBFChain *) StMaker::GetChain();
00570   St_g2t_track *Track = (St_g2t_track *) chain->FindObject("g2t_track");
00571   St_g2t_tpc_hit *mcTpcHits = (St_g2t_tpc_hit *) chain->FindObject("g2t_tpc_hit");
00572   if (Track && mcTpcHits) {
00573     static  StTpcLocalSectorCoordinate        localSect;
00574     static  StTpcLocalSectorDirection         dirSect;
00575     StTpcCoordinateTransform transform(gStTpcDb);
00576     int nMCTracks = Track->GetNRows();
00577     const g2t_track_st *track = Track->GetTable();
00578     
00579     vector<int> truthToIndex; // need for finding MCPoints to MCTracks correspondence
00580     truthToIndex.resize(nMCTracks+1,-1);
00581     
00582     const g2t_tpc_hit_st *mcTpcHit = mcTpcHits->GetTable();
00583     const Int_t nMcTpcHits = mcTpcHits->GetNRows();
00584      
00585     AliHLTTPCCAMCTrack mcTrack;
00586     for (Int_t l = 0; l < nMCTracks; l++, track++) {
00587       if (! track->n_tpc_hit) continue;
00588       if (track->ptot < 1e-3) continue;
00589       memset(&mcTrack, 0, sizeof(AliHLTTPCCAMCTrack));
00590       mcTrack.SetPDG( TDatabasePDG::Instance()->ConvertGeant3ToPdg(track->ge_pid) );
00591       TParticlePDG *part = TDatabasePDG::Instance()->GetParticle( mcTrack.PDG() );
00592       Float_t      q     = part->Charge();
00593       mcTrack.SetP ( track->ptot );
00594       mcTrack.SetPt( track->pt );
00595         //      mcTrack.SetNHits( track->n_tpc_hit ); // it is not really what we need. See below
00596       for (Int_t i = 0; i < 3; i++) { // TODO
00597         mcTrack.SetPar( i, ((i>=1)?-1:1) * 0 );
00598       }
00599       for (Int_t i = 0; i < 3; i++) {
00600         mcTrack.SetPar( 3+i, ((i>=1)?-1:1) * track->p[i]/mcTrack.P() );
00601       }
00602       mcTrack.SetPar(6, q/mcTrack.P()/3 ); // q=3q      
00603       mcTrack.SetNTurns( 1 ); // TODO: read from somewhere
00604       const int idTruth = l+1;
00605 
00606       truthToIndex[idTruth] = mcTracks.size();
00607       mcTracks.push_back(mcTrack);
00608     } // for iMCTrack
00609     
00610       // get hits
00611     const int nHits = caHits.size();
00612     for (int iHit = 0; iHit < nHits; iHit++){
00613       const int iTrack = truthToIndex[idsTruth[iHit]];
00614       if (iTrack == -1) continue;
00615       AliHLTTPCCAHitLabel &l = hitLabels[iHit];
00616       l.fLab[0] = iTrack;
00617       l.fLab[1] = -1;
00618       l.fLab[2] = -1;
00619       mcTracks[iTrack].SetNHits( 1 + mcTracks[iTrack].NHits() );
00620     } // j
00621       
00622       // get MCPoints
00623     for (Int_t iP = 0; iP < nMcTpcHits; iP++, mcTpcHit++) {
00624       AliHLTTPCCALocalMCPoint mcPoint;
00625       Int_t is     =  mcTpcHit->volume_id/100000;
00626       if (is) continue;
00627       const int iTrack = truthToIndex[mcTpcHit->track_p];
00628       if (iTrack == -1) continue;
00629       AliHLTTPCCAMCTrack& mcTrack = mcTracks[iTrack];
00630       
00631       mcPoint.SetIRow( mcTpcHit->volume_id%100 - 1 );
00632       mcPoint.SetISlice( (mcTpcHit->volume_id/100)%100 - 1 );
00633 
00634       Int_t sector = (mcTpcHit->volume_id%100000)/100;
00635       Int_t row    =  mcTpcHit->volume_id%100;
00636 
00637         // get detector
00638       StiDetectorBuilder* detectorBuilder = 0;
00639       StiDetectorGroups *groups=StiToolkit::instance()->getDetectorGroups();
00640       vector<StiGenericDetectorGroup *>::iterator it = groups->begin();
00641       for (; it != groups->end(); ++it) {
00642         StiGenericDetectorGroup *group = *it;
00643         detectorBuilder = group->getDetectorBuilder();
00644         TString builderName = (const char*)(detectorBuilder->getName().c_str());
00645         if (builderName == "TpcBuilder")
00646           break;
00647       }
00648       int sectorTpm = sector - 1;
00649       if (sectorTpm >= 12)    sectorTpm = 11 - (sectorTpm - 11)%12;
00650       StiDetector* detector = detectorBuilder->getDetector(row-1,sectorTpm);
00651 
00652       static StiHit tmpStiHit; // used for convert only
00653 
00654         // Get XYZ
00655       StTpcHit* mpTmp = new StTpcHit;
00656       double energyTmp = 0.;
00657       tmpStiHit.setGlobal(detector, mpTmp, mcTpcHit->x[0], mcTpcHit->x[1], mcTpcHit->x[2], energyTmp);
00658       mcPoint.SetX( tmpStiHit.x() );
00659       mcPoint.SetY( - tmpStiHit.y() );
00660       mcPoint.SetZ( - tmpStiHit.z() );
00661 
00662         // Get P
00663       StGlobalDirection  gDir(mcTpcHit->p[0],mcTpcHit->p[1],mcTpcHit->p[2]);
00664       transform(gDir,dirSect,sector,row);
00665       double px = mcTpcHit->p[0], py = mcTpcHit->p[1], pz_new = mcTpcHit->p[2], px_new, py_new;
00666       double angle = detector->getPlacement()->getNormalRefAngle();
00667       px_new =  cos(angle)*px + sin(angle)*py;
00668       py_new = -sin(angle)*px + cos(angle)*py;
00669 
00670       mcPoint.SetPx( px_new );
00671       mcPoint.SetPy( - py_new );
00672       mcPoint.SetPz( - pz_new );
00673       const int qpSign = fabs(mcTrack.Par(6))/mcTrack.Par(6);
00674       mcPoint.SetQP( double( qpSign ) / sqrt(px_new*px_new + py_new*py_new + pz_new*pz_new) );
00675       mcPoint.SetTrackI( iTrack );
00676       mcPoint.SetTrackID( mcTpcHit->track_p );
00677 
00678       mcPoints.push_back(mcPoint);
00679     } // for mcPoints
00680     
00681     std::sort( mcPoints.begin(), mcPoints.end(), myfunction );
00682     mcTracks[0].SetFirstMCPointID(0);
00683     int iPLast = 0;
00684     int iTr = 0;
00685     for( UInt_t iP = 0; iP < mcPoints.size(); ++iP ) {
00686       if ( mcPoints[iP].TrackI() != iTr ) {
00687         mcTracks[iTr].SetNMCPoints(iP - iPLast);
00688         iPLast = iP;
00689         iTr = mcPoints[iP].TrackI();
00690         mcTracks[iTr].SetFirstMCPointID(iP);
00691       }
00692     }
00693     mcTracks[iTr].SetNMCPoints(mcPoints.size() - iPLast);
00694   }
00695 } // void StiTPCCATrackerInterface::FillPerformance()
00696 
00697 void StiTPCCATrackerInterface::FillStiPerformance()
00698 {
00699   fStiCaHits.clear();
00700   fStiIdTruth.clear();
00701   fStiCaTracks.clear();
00702   
00703   AliHLTTPCCAStiPerformance* stiPerf = (AliHLTTPCCAStiPerformance*)fPerformance->GetSubPerformance("Sti Performance");
00704   int NGBHits = 0;
00705 
00706   
00707   StTpcCoordinateTransform tran(gStTpcDb);
00708   StTpcLocalSectorCoordinate loc;
00709 
00710   for(int iTr=0; iTr<fStiTracks->getTrackCount(0); iTr++)
00711   {
00712     StiKalmanTrack * track = (StiKalmanTrack*) fStiTracks->at(iTr);
00713     vector<StiHit*> hits_v = track-> getHits();
00714 
00715     AliHLTTPCCAGBTrack GBTrack;
00716     int nTrackHits = 0;
00717 
00718         // get local coordinates. take into account distortion
00719 
00720     for(unsigned iH=0; iH<hits_v.size(); iH++)
00721     {
00722       StiHit *hit = hits_v[iH];
00723       if (! hit->stHit()) continue;
00724       const StTpcHit *tpcHit = dynamic_cast<const StTpcHit*>(hit->stHit());
00725       if ( ! tpcHit) continue;
00726       StGlobalCoordinate glob(tpcHit->position());
00727       tran(glob,loc,tpcHit->sector(),tpcHit->padrow());
00728 
00729         // convert to CA Hit
00730       AliHLTTPCCAGBHit caHit;
00731       caHit.SetIRow( tpcHit->padrow()-1 );
00732       caHit.SetX( hit->position() ); // take position of the row
00733       caHit.SetY( - hit->y() );
00734       caHit.SetZ( - hit->z() );
00735       caHit.SetErrY( 0.12 );// TODO: read parameters from somewhere 
00736       caHit.SetErrZ( 0.16 );
00737       caHit.SetISlice( tpcHit->sector() - 1 );
00738       caHit.SetID( fStiCaHits.size() );
00739       fStiIdTruth.push_back( tpcHit->idTruth() );
00740 
00741       int iHit = NGBHits+iH;
00742       stiPerf->AddHit(iHit);
00743       nTrackHits++;
00744       
00745       fStiCaHits.push_back(caHit);
00746     }
00747             
00748     if(!(track->getInnerMostTPCHitNode(0))) continue;
00749     StiKalmanTrackNode *NodePar = track->getInnerMostTPCHitNode(0);
00750 //    std::cout <<"node x  "<< track->getInnerMostTPCHitNode(0)->getX() << std::endl << std::endl;
00751 
00752     double JI[5]; 
00753     JI[0] = -1.;                    // y
00754     JI[1] = -1.;                    // z
00755     JI[2] = -1.;                    // eta
00756     JI[3] =  1.;                    // ptin
00757     JI[4] = -1.;                    // tanl
00758 
00759     double J[5]; 
00760     J[0] = JI[0];                              // y
00761     J[1] = JI[1];                              // z
00762     J[2] = JI[2]*cos(NodePar->fitPars().eta()); // eta
00763     J[3] = JI[3];                              // ptin
00764     J[4] = JI[4];                              // tanl
00765     
00766     double Cos = cos(NodePar->fitPars().eta());
00767     float sgnCos = fabs(Cos)/Cos;
00768 
00769     AliHLTTPCCATrackParam    GBTrackParam;
00770     GBTrackParam.SetX     ( NodePar->fitPars().x() );
00771     GBTrackParam.SetY     ( JI[0] * NodePar->fitPars().y() );
00772     GBTrackParam.SetZ     ( JI[1] * NodePar->fitPars().z() );
00773     GBTrackParam.SetSinPhi( JI[2] * sin(NodePar->fitPars().eta()));
00774     GBTrackParam.SetQPt   ( JI[3] * NodePar->fitPars().ptin() );
00775     GBTrackParam.SetDzDs  ( JI[4] * NodePar->fitPars().tanl() );
00776 
00777     GBTrackParam.SetCov( 0, NodePar->fitErrs()._cYY);
00778     GBTrackParam.SetCov( 1, NodePar->fitErrs()._cZY*J[0]*J[1]);
00779     GBTrackParam.SetCov( 2, NodePar->fitErrs()._cZZ*J[0]*J[1]);
00780     GBTrackParam.SetCov( 3, NodePar->fitErrs()._cEY*J[0]*J[2]);
00781     GBTrackParam.SetCov( 4, NodePar->fitErrs()._cEZ*J[1]*J[2]);
00782     GBTrackParam.SetCov( 5, NodePar->fitErrs()._cEE*J[2]*J[2]);
00783     GBTrackParam.SetCov( 6, NodePar->fitErrs()._cTY*J[0]*J[4]);
00784     GBTrackParam.SetCov( 7, NodePar->fitErrs()._cTZ*J[1]*J[4]);
00785     GBTrackParam.SetCov( 8, NodePar->fitErrs()._cTE*J[2]*J[4]);    
00786     GBTrackParam.SetCov( 9, NodePar->fitErrs()._cTT*J[4]*J[4]);
00787     GBTrackParam.SetCov(10, NodePar->fitErrs()._cPY*J[0]*J[3]);
00788     GBTrackParam.SetCov(11, NodePar->fitErrs()._cPZ*J[1]*J[3]);
00789     GBTrackParam.SetCov(12, NodePar->fitErrs()._cPE*J[2]*J[3]);
00790     GBTrackParam.SetCov(13, NodePar->fitErrs()._cTP*J[4]*J[3]);
00791     GBTrackParam.SetCov(14, NodePar->fitErrs()._cPP*J[3]*J[3]);
00792     
00793     GBTrackParam.SetSignCosPhi(sgnCos);
00794     
00795     GBTrack.SetInnerParam(GBTrackParam);
00796     
00797     GBTrack.SetNHits( nTrackHits );
00798     GBTrack.SetFirstHitRef( NGBHits );
00799     
00800     NGBHits += nTrackHits;
00801     
00802     fStiCaTracks.push_back(GBTrack);
00803  }
00804 
00805     // prepare mc
00806   vector<AliHLTTPCCAMCTrack> stiMCTracks;
00807   vector<AliHLTTPCCALocalMCPoint> stiMCPoints;
00808   vector<AliHLTTPCCAHitLabel> stiHitLabels;
00809   FillPerformance(fStiCaHits,fStiIdTruth, stiMCTracks,stiMCPoints,stiHitLabels);
00810 
00811   fStiMCTracks.Resize(stiMCTracks.size());
00812   for (unsigned i = 0; i < stiMCTracks.size(); i++) {
00813     fStiMCTracks[i] = stiMCTracks[i];
00814   }
00815   fStiMCPoints.Resize(stiMCPoints.size());
00816   for (unsigned i = 0; i < stiMCPoints.size(); i++) {
00817     fStiMCPoints[i] = stiMCPoints[i];
00818   }
00819   fStiHitLabels.Resize(stiHitLabels.size());
00820   for (unsigned i = 0; i < stiHitLabels.size(); i++) {
00821     fStiHitLabels[i] = stiHitLabels[i];
00822   }
00823   
00824 
00825 } // void StiTPCCATrackerInterface::FillStiPerformance()
00826 #endif /* DO_TPCCATRACKER */