HLT/TPCLib/tracking-ca/AliTPCtrackerCA.cxx

   1 // $Id$
   2 //***************************************************************************
   3 // This file is property of and copyright by the ALICE HLT Project          *
   4 // ALICE Experiment at CERN, All rights reserved.                           *
   5 //                                                                          *
   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.                              *
   9 //                                                                          *
  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 //***************************************************************************
  18
  19 #include "AliTPCtrackerCA.h"
  20
  21 #include <TTree.h>
  22 #include <Riostream.h>
  23 #include "AliCluster.h"
  24 #include "AliTPCClustersRow.h"
  25 #include "AliTPCParam.h"
  26 #include "AliRun.h"
  27 #include "AliRunLoader.h"
  28 #include "AliStack.h"
  29
  30 #include "AliHLTTPCCATracker.h"
  31 #include "AliHLTTPCCAGBHit.h"
  32 #include "AliHLTTPCCAGBTracker.h"
  33 #include "AliHLTTPCCAGBTrack.h"
  34 #include "AliHLTTPCCAMCTrack.h"
  35 #include "AliHLTTPCCAOutTrack.h"
  36 #include "AliHLTTPCCAPerformance.h"
  37 #include "AliHLTTPCCAParam.h"
  38 #include "AliHLTTPCCATrackConvertor.h"
  39
  40 #include "TMath.h"
  41 #include "AliTPCLoader.h"
  42 #include "AliTPC.h"
  43 #include "AliTPCclusterMI.h"
  44 #include "AliTPCTransform.h"
  45 #include "AliTPCcalibDB.h"
  46 #include "AliTPCReconstructor.h"
  47 #include "AliTPCtrack.h"
  48 #include "AliESDtrack.h"
  49 #include "AliESDEvent.h"
  50 #include "AliTrackReference.h"
  51
  52 //#include <fstream.h>
  53
  54 ClassImp(AliTPCtrackerCA)
  55
  56 AliTPCtrackerCA::AliTPCtrackerCA()
  57   :AliTracker(),fParam(0), fClusters(0), fNClusters(0), fHLTTracker(0),fHLTPerformance(0),fDoHLTPerformance(0),fDoHLTPerformanceClusters(0),fStatNEvents(0)
  58 {
  59   //* default constructor
  60 }
  61
  62 AliTPCtrackerCA::AliTPCtrackerCA(const AliTPCtrackerCA &):
  63   AliTracker(),fParam(0), fClusters(0), fNClusters(0), fHLTTracker(0),fHLTPerformance(0),fDoHLTPerformance(0),fDoHLTPerformanceClusters(0),fStatNEvents(0)
  64 {
  65   //* dummy
  66 }
  67
  68 AliTPCtrackerCA & AliTPCtrackerCA::operator=(const AliTPCtrackerCA& )
  69 {
  70   //* dummy
  71   return *this;
  72 }
  73
  74
  75 AliTPCtrackerCA::~AliTPCtrackerCA()
  76 {
  77   //* destructor
  78   if( fClusters ) delete[] fClusters;
  79   if( fHLTTracker ) delete fHLTTracker;
  80   if( fHLTPerformance ) delete fHLTPerformance;
  81 }
  82
  83 AliTPCtrackerCA::AliTPCtrackerCA(const AliTPCParam *par):
  84   AliTracker(),fParam(par), fClusters(0), fNClusters(0), fHLTTracker(0), fHLTPerformance(0),fDoHLTPerformance(0),fDoHLTPerformanceClusters(0),fStatNEvents(0)
  85 {
  86   //* constructor
  87
  88   DoHLTPerformance() = 1;
  89   DoHLTPerformanceClusters() = 1;
  90
  91   fHLTTracker = new AliHLTTPCCAGBTracker;
  92   fHLTTracker->SetNSlices( fParam->GetNSector()/2 );
  93
  94   if( fDoHLTPerformance ){
  95     fHLTPerformance = new AliHLTTPCCAPerformance;
  96     fHLTPerformance->SetTracker( fHLTTracker );
  97   }
  98
  99   for( int iSlice=0; iSlice<fHLTTracker->NSlices(); iSlice++ ){
 100
 101     Float_t bz = AliTracker::GetBz();
 102
 103     Float_t inRmin = fParam->GetInnerRadiusLow();
 104     //Float_t inRmax = fParam->GetInnerRadiusUp();
 105     //Float_t outRmin = fParam->GetOuterRadiusLow();
 106     Float_t outRmax = fParam->GetOuterRadiusUp();
 107     Float_t plusZmin = 0.0529937;
 108     Float_t plusZmax = 249.778;
 109     Float_t minusZmin = -249.645;
 110     Float_t minusZmax = -0.0799937;
 111     Float_t dalpha = 0.349066;
 112     Float_t alpha = 0.174533 + dalpha*iSlice;
 113
 114     Bool_t zPlus = (iSlice<18 );
 115     Float_t zMin =  zPlus ?plusZmin :minusZmin;
 116     Float_t zMax =  zPlus ?plusZmax :minusZmax;
 117     //TPCZmin = -249.645, ZMax = 249.778
 118     //Float_t rMin =  inRmin;
 119     //Float_t rMax =  outRmax;
 120
 121     Float_t padPitch = 0.4;
 122     Float_t sigmaZ = 0.228808;
 123
 124     Int_t NRows = fParam->GetNRowLow()+fParam->GetNRowUp();
 125
 126     Float_t rowX[200];
 127     for( Int_t irow=0; irow<fParam->GetNRowLow(); irow++){
 128       rowX[irow] = fParam->GetPadRowRadiiLow(irow);
 129     }
 130     for( Int_t irow=0; irow<fParam->GetNRowUp(); irow++){
 131       rowX[fParam->GetNRowLow()+irow] = fParam->GetPadRowRadiiUp(irow);
 132     }
 133     AliHLTTPCCAParam param;
 134     param.Initialize( iSlice, NRows, rowX, alpha, dalpha,
 135                       inRmin, outRmax, zMin, zMax, padPitch, sigmaZ, bz );
 136     param.YErrorCorrection() = 1.;//.33;
 137     param.ZErrorCorrection() = 1.;//.33;
 138     param.MaxTrackMatchDRow() = 5;
 139     param.TrackConnectionFactor() = 3.5;
 140     fHLTTracker->Slices()[iSlice].Initialize( param );
 141   }
 142 }
 143
 144
 145
 146 Int_t AliTPCtrackerCA::LoadClusters (TTree * tree)
 147 {
 148   fNClusters = 0;
 149   if( fClusters ) delete[] fClusters;
 150
 151   fHLTTracker->StartEvent();
 152   if( fDoHLTPerformance ) fHLTPerformance->StartEvent();
 153
 154   if( !fParam ) return 1;
 155
 156   // load mc tracks
 157   while( fDoHLTPerformance ){
 158     if( !gAlice ) break;
 159     AliRunLoader *rl = gAlice->GetRunLoader();
 160     if( !rl ) break;
 161     rl->LoadKinematics();
 162     AliStack *stack = rl->Stack();
 163     if( !stack ) break;
 164
 165     fHLTPerformance->SetNMCTracks( stack->GetNtrack() );
 166
 167     for( Int_t itr=0; itr<stack->GetNtrack(); itr++ ){
 168       TParticle *part = stack->Particle(itr);
 169       fHLTPerformance->ReadMCTrack( itr, part );
 170     }
 171
 172     { // check for MC tracks at the TPC entrance
 173       Bool_t *isTPC = 0;
 174       isTPC = new Bool_t [stack->GetNtrack()];
 175       for( Int_t i=0; i<stack->GetNtrack(); i++ ) isTPC[i] = 0;
 176       rl->LoadTrackRefs();
 177       TTree *TR = rl->TreeTR();
 178       if( !TR ) break;
 179       TBranch *branch=TR->GetBranch("TrackReferences");
 180       if (!branch ) break;
 181       TClonesArray tpcdummy("AliTrackReference",1000), *tpcRefs=&tpcdummy;
 182       branch->SetAddress(&tpcRefs);
 183       Int_t nr=(Int_t)TR->GetEntries();
 184       for (Int_t r=0; r<nr; r++) {
 185         TR->GetEvent(r);
 186         Int_t nref = tpcRefs->GetEntriesFast();
 187         if (!nref) continue;
 188         AliTrackReference *tpcRef= 0x0;
 189         for (Int_t iref=0; iref<nref; ++iref) {
 190           tpcRef = (AliTrackReference*)tpcRefs->UncheckedAt(iref);
 191           if (tpcRef->DetectorId() == AliTrackReference::kTPC) break;
 192           tpcRef = 0x0;
 193         }
 194         if (!tpcRef) continue;
 195
 196         if( isTPC[tpcRef->Label()] ) continue;
 197
 198         fHLTPerformance->ReadMCTPCTrack(tpcRef->Label(),
 199                                         tpcRef->X(),tpcRef->Y(),tpcRef->Z(),
 200                                         tpcRef->Px(),tpcRef->Py(),tpcRef->Pz() );
 201         isTPC[tpcRef->Label()] = 1;
 202         tpcRefs->Clear();
 203       }
 204       if( isTPC ) delete[] isTPC;
 205     }
 206
 207     while( fDoHLTPerformanceClusters ){
 208       AliTPCLoader *tpcl = (AliTPCLoader*) rl->GetDetectorLoader("TPC");
 209       if( !tpcl ) break;
 210       if( tpcl->TreeH() == 0x0 ){
 211         if( tpcl->LoadHits() ) break;
 212       }
 213       if( tpcl->TreeH() == 0x0 ) break;
 214
 215       AliTPC *tpc = (AliTPC*) gAlice->GetDetector("TPC");
 216       Int_t nEnt=(Int_t)tpcl->TreeH()->GetEntries();
 217       Int_t nPoints = 0;
 218       for (Int_t iEnt=0; iEnt<nEnt; iEnt++) {
 219         tpc->ResetHits();
 220         tpcl->TreeH()->GetEvent(iEnt);
 221         AliTPChit *phit = (AliTPChit*)tpc->FirstHit(-1);
 222         for ( ; phit; phit=(AliTPChit*)tpc->NextHit() ) nPoints++;
 223       }
 224       fHLTPerformance->SetNMCPoints( nPoints );
 225
 226       for (Int_t iEnt=0; iEnt<nEnt; iEnt++) {
 227         tpc->ResetHits();
 228         tpcl->TreeH()->GetEvent(iEnt);
 229         AliTPChit *phit = (AliTPChit*)tpc->FirstHit(-1);
 230         for ( ; phit; phit=(AliTPChit*)tpc->NextHit() ){
 231           fHLTPerformance->ReadMCPoint( phit->GetTrack(),phit->X(), phit->Y(),phit->Z(),phit->Time(), phit->fSector%36);
 232         }
 233       }
 234       break;
 235     }
 236     break;
 237   }
 238
 239   TBranch * br = tree->GetBranch("Segment");
 240   if( !br ) return 1;
 241
 242   AliTPCClustersRow *clrow = new AliTPCClustersRow;
 243   clrow->SetClass("AliTPCclusterMI");
 244   clrow->SetArray(0);
 245   clrow->GetArray()->ExpandCreateFast(10000);
 246
 247   br->SetAddress(&clrow);
 248
 249   //
 250   Int_t NEnt=Int_t(tree->GetEntries());
 251
 252   fNClusters = 0;
 253   for (Int_t i=0; i<NEnt; i++) {
 254     br->GetEntry(i);
 255     Int_t sec,row;
 256     fParam->AdjustSectorRow(clrow->GetID(),sec,row);
 257     fNClusters += clrow->GetArray()->GetEntriesFast();
 258   }
 259
 260   fClusters = new AliTPCclusterMI [fNClusters];
 261   fHLTTracker->SetNHits( fNClusters );
 262   if( fDoHLTPerformance ) fHLTPerformance->SetNHits( fNClusters );
 263   int ind=0;
 264   for (Int_t i=0; i<NEnt; i++) {
 265     br->GetEntry(i);
 266     Int_t sec,row;
 267     fParam->AdjustSectorRow(clrow->GetID(),sec,row);
 268     int NClu = clrow->GetArray()->GetEntriesFast();
 269     Float_t x = fParam->GetPadRowRadii(sec,row);
 270     for (Int_t icl=0; icl<NClu; icl++){
 271       Int_t lab0 = -1;
 272       Int_t lab1 = -1;
 273       Int_t lab2 = -1;
 274       AliTPCclusterMI* cluster = (AliTPCclusterMI*)(clrow->GetArray()->At(icl));
 275       if( !cluster ) continue;
 276       lab0 = cluster->GetLabel(0);
 277       lab1 = cluster->GetLabel(1);
 278       lab2 = cluster->GetLabel(2);
 279
 280       AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
 281       if (!transform) {
 282         AliFatal("Tranformations not in calibDB");
 283       }
 284       Double_t xx[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
 285       Int_t id[1]={cluster->GetDetector()};
 286       transform->Transform(xx,id,0,1);
 287       //if (!AliTPCReconstructor::GetRecoParam()->GetBYMirror()){
 288       if (cluster->GetDetector()%36>17){
 289         xx[1]*=-1;
 290       }
 291       //}
 292
 293       cluster->SetX(xx[0]);
 294       cluster->SetY(xx[1]);
 295       cluster->SetZ(xx[2]);
 296
 297       TGeoHMatrix  *mat = fParam->GetClusterMatrix(cluster->GetDetector());
 298       Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
 299       Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
 300       if (mat) mat->LocalToMaster(pos,posC);
 301       else{
 302         // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
 303       }
 304       cluster->SetX(posC[0]);
 305       cluster->SetY(posC[1]);
 306       cluster->SetZ(posC[2]);
 307
 308       Float_t y = cluster->GetY();
 309       Float_t z = cluster->GetZ();
 310
 311       if( sec>=36 ){
 312         sec = sec - 36;
 313         row = row + fParam->GetNRowLow();
 314       }
 315
 316       Int_t index = ind++;
 317       fClusters[index] = *cluster;
 318       fHLTTracker->ReadHit( x, y, z,
 319                             TMath::Sqrt(cluster->GetSigmaY2()), TMath::Sqrt(cluster->GetSigmaZ2()),
 320                             cluster->GetMax(), index, sec, row );
 321       if( fDoHLTPerformance ) fHLTPerformance->ReadHitLabel(index, lab0, lab1, lab2 );
 322     }
 323   }
 324   delete clrow;
 325   return 0;
 326 }
 327
 328 AliCluster * AliTPCtrackerCA::GetCluster(Int_t index) const
 329 {
 330   return &(fClusters[index]);
 331 }
 332
 333 Int_t AliTPCtrackerCA::Clusters2Tracks( AliESDEvent *event )
 334 {
 335   //cout<<"Start of AliTPCtrackerCA"<<endl;
 336
 337   fHLTTracker->FindTracks();
 338   if( fDoHLTPerformance ) fHLTPerformance->Performance();
 339
 340   if( 0 ) {// Write Event
 341     if( fStatNEvents == 0 ){
 342       fstream geo;
 343       geo.open("CAEvents/settings.dat", ios::out);
 344       if( geo.is_open() ){
 345         fHLTTracker->WriteSettings(geo);
 346       }
 347       geo.close();
 348     }
 349
 350     fstream event;
 351     char name[255];
 352     sprintf( name,"CAEvents/%i.event.dat",fStatNEvents );
 353     event.open(name, ios::out);
 354     if( event.is_open() ){
 355       fHLTTracker->WriteEvent(event);
 356       fstream tracks;
 357       sprintf( name,"CAEvents/%i.tracks.dat",fStatNEvents );
 358       tracks.open(name, ios::out);
 359       fHLTTracker->WriteTracks(tracks);
 360     }
 361     event.close();
 362     if( fDoHLTPerformance ){
 363       fstream mcevent, mcpoints;
 364       char mcname[255];
 365       sprintf( mcname,"CAEvents/%i.mcevent.dat",fStatNEvents );
 366       mcevent.open(mcname, ios::out);
 367       if( mcevent.is_open() ){
 368         fHLTPerformance->WriteMCEvent(mcevent);
 369       }
 370       if(fDoHLTPerformanceClusters ){
 371         sprintf( mcname,"CAEvents/%i.mcpoints.dat",fStatNEvents );
 372         mcpoints.open(mcname, ios::out);
 373         if( mcpoints.is_open() ){
 374           fHLTPerformance->WriteMCPoints(mcpoints);
 375         }
 376         mcpoints.close();
 377       }
 378       mcevent.close();
 379     }
 380   }
 381   fStatNEvents++;
 382
 383   if( event ){
 384
 385     Float_t bz = fHLTTracker->Slices()[0].Param().Bz();
 386
 387     for( Int_t itr=0; itr<fHLTTracker->NTracks(); itr++ ){
 388       AliTPCtrack tTPC;
 389       AliHLTTPCCAGBTrack &tCA = fHLTTracker->Tracks()[itr];
 390       AliHLTTPCCATrackParam &par = tCA.Param();
 391       AliHLTTPCCATrackConvertor::GetExtParam( par, tTPC, tCA.Alpha(), bz );
 392       tTPC.SetMass(0.13957);
 393       tTPC.SetdEdx( tCA.DeDx() );
 394       if( TMath::Abs(tTPC.GetSigned1Pt())>1./0.1 ) continue;
 395       int nhits = tCA.NHits();
 396       if( nhits>kMaxRow ) nhits = kMaxRow;
 397       tTPC.SetNumberOfClusters(nhits);
 398       for( Int_t ih=0; ih<nhits; ih++ ){
 399         Int_t index = fHLTTracker->TrackHits()[tCA.FirstHitRef()+ih];
 400         Int_t ext_index = fHLTTracker->Hits()[index].ID();
 401         tTPC.SetClusterIndex(ih, ext_index);
 402       }
 403       CookLabel(&tTPC,0.1);
 404       {
 405         Double_t xTPC=fParam->GetInnerRadiusLow();
 406         Double_t dAlpha = fParam->GetInnerAngle()/180.*TMath::Pi();
 407         if (tTPC.AliExternalTrackParam::PropagateTo(xTPC,bz)) {
 408           Double_t y=tTPC.GetY();
 409           Double_t ymax=xTPC*TMath::Tan(dAlpha/2.);
 410           if (y > ymax) {
 411             if (tTPC.Rotate(dAlpha)) tTPC.AliExternalTrackParam::PropagateTo(xTPC,bz);
 412           } else if (y <-ymax) {
 413             if (tTPC.Rotate(-dAlpha)) tTPC.AliExternalTrackParam::PropagateTo(xTPC,bz);
 414           }
 415         }
 416       }
 417
 418       AliESDtrack tESD;
 419       tESD.UpdateTrackParams( &(tTPC),AliESDtrack::kTPCin);
 420       //tESD.SetStatus( AliESDtrack::kTPCrefit );
 421       //tESD.SetTPCPoints(tTPC.GetPoints());
 422       //tESD.myTPC = tTPC;
 423       event->AddTrack(&tESD);
 424     }
 425   }
 426
 427   //cout<<"End of AliTPCtrackerCA"<<endl;
 428   return 0;
 429 }
 430
 431
 432 Int_t AliTPCtrackerCA::RefitInward (AliESDEvent *event)
 433 {
 434   //* back propagation of ESD tracks (not fully functional)
 435
 436   Float_t bz = fHLTTracker->Slices()[0].Param().Bz();
 437   Float_t xTPC = fParam->GetInnerRadiusLow();
 438   Float_t dAlpha = fParam->GetInnerAngle()/180.*TMath::Pi();
 439   Float_t yMax = xTPC*TMath::Tan(dAlpha/2.);
 440
 441   Int_t nentr=event->GetNumberOfTracks();
 442
 443   for (Int_t i=0; i<nentr; i++) {
 444     AliESDtrack *esd=event->GetTrack(i);
 445     ULong_t status=esd->GetStatus();
 446     if (!(status&AliESDtrack::kTPCin)) continue;
 447     AliHLTTPCCATrackParam t0;
 448     AliHLTTPCCATrackConvertor::SetExtParam(t0,*esd, bz );
 449     Float_t alpha = esd->GetAlpha();
 450     if( t0.TransportToXWithMaterial( xTPC, bz) ){
 451       if (t0.GetY() > yMax) {
 452         if (t0.Rotate(dAlpha)){
 453           alpha+=dAlpha;
 454           t0.TransportToXWithMaterial( xTPC, bz);
 455         }
 456       } else if (t0.GetY() <-yMax) {
 457         if (t0.Rotate(-dAlpha)){
 458           alpha+=-dAlpha;
 459           t0.TransportToXWithMaterial( xTPC, bz);
 460         }
 461       }
 462     }
 463     AliTPCtrack tt(*esd);
 464     AliHLTTPCCATrackConvertor::GetExtParam(t0,tt,alpha,bz);
 465     esd->UpdateTrackParams( &tt,AliESDtrack::kTPCrefit);
 466   }
 467   return 0;
 468 }
 469
 470 Int_t AliTPCtrackerCA::PropagateBack(AliESDEvent *)
 471 {
 472   //* not implemented yet
 473   return 0;
 474 }