HLT/src/AliHLTReconstructor.cxx

   1 // $Id$
   2
   3 ///////////////////////////////////////////////////////////////////////////////
   4 //                                                                           //
   5 // class for HLT reconstruction                                              //
   6 // <Cvetan.Cheshkov@cern.ch>                                                 //
   7 // <loizides@ikf.uni-frankfurt.de>                                           //
   8 ///////////////////////////////////////////////////////////////////////////////
   9
  10 // very ugly but it has to work fast
  11 #ifdef use_reconstruction
  12
  13 #include <Riostream.h>
  14 #include <TSystem.h>
  15 #include <TArrayF.h>
  16
  17 #include "AliL3StandardIncludes.h"
  18 #include "AliL3Logging.h"
  19 #include "AliLevel3.h"
  20 #include "AliL3Evaluate.h"
  21 #include "AliHLTReconstructor.h"
  22 #include "AliL3Transform.h"
  23 #include "AliL3Hough.h"
  24 #include "AliL3FileHandler.h"
  25 #include "AliL3Track.h"
  26 #include "AliL3HoughTrack.h"
  27 #include "AliL3TrackArray.h"
  28 #include "AliRunLoader.h"
  29 #include "AliHeader.h"
  30 #include "AliGenEventHeader.h"
  31 #include "AliESD.h"
  32 #include "AliESDHLTtrack.h"
  33
  34 #if __GNUC__== 3
  35 using namespace std;
  36 #endif
  37
  38
  39 ClassImp(AliHLTReconstructor)
  40
  41 void AliHLTReconstructor::Reconstruct(AliRunLoader* runLoader) const
  42 {
  43   if(!runLoader) {
  44     LOG(AliL3Log::kFatal,"AliHLTReconstructor::Reconstruct","RunLoader")
  45       <<" Missing RunLoader! 0x0"<<ENDLOG;
  46     return;
  47   }
  48   gSystem->Exec("rm -rf hlt");
  49   gSystem->MakeDirectory("hlt");
  50   gSystem->Exec("rm -rf hough");
  51   gSystem->MakeDirectory("hough");
  52
  53   Bool_t isinit=AliL3Transform::Init(runLoader);
  54   if(!isinit){
  55     LOG(AliL3Log::kError,"AliHLTReconstructor::Reconstruct","Transformer")
  56      << "Could not create transform settings, please check log for error messages!" << ENDLOG;
  57     return;
  58   }
  59
  60   Int_t nEvents = runLoader->GetNumberOfEvents();
  61
  62   for(Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
  63     runLoader->GetEvent(iEvent);
  64
  65     if(fDoTracker) ReconstructWithConformalMapping(runLoader,iEvent);
  66     if(fDoHough) ReconstructWithHoughTransform(runLoader,iEvent);
  67   }
  68 }
  69
  70 void AliHLTReconstructor::ReconstructWithConformalMapping(AliRunLoader* runLoader,Int_t iEvent) const
  71 {
  72   AliLevel3 *fHLT = new AliLevel3(runLoader);
  73   fHLT->Init("./", AliLevel3::kRunLoader, 1);
  74
  75   Int_t phiSegments = 50;
  76   Int_t etaSegments = 100;
  77   Int_t trackletlength = 3;
  78   Int_t tracklength = 10;
  79   Int_t rowscopetracklet = 2;
  80   Int_t rowscopetrack = 10;
  81   Double_t minPtFit = 0;
  82   Double_t maxangle = 0.1745;
  83   Double_t goodDist = 5;
  84   Double_t maxphi = 0.1;
  85   Double_t maxeta = 0.1;
  86   Double_t hitChi2Cut = 20;
  87   Double_t goodHitChi2 = 5;
  88   Double_t trackChi2Cut = 10;
  89   Double_t xyerror = -1;
  90   Double_t zerror =  -1;
  91
  92   fHLT->SetClusterFinderParam(xyerror,zerror,kTRUE);
  93   fHLT->SetTrackerParam(phiSegments, etaSegments,
  94                         trackletlength, tracklength,
  95                         rowscopetracklet, rowscopetrack,
  96                         minPtFit, maxangle, goodDist, hitChi2Cut,
  97                         goodHitChi2, trackChi2Cut, 50, maxphi, maxeta, kTRUE);
  98   fHLT->SetTrackerParam(phiSegments, etaSegments,
  99                         trackletlength, tracklength,
 100                         rowscopetracklet, rowscopetrack,
 101                         minPtFit, maxangle, goodDist, hitChi2Cut,
 102                         goodHitChi2, trackChi2Cut, 50, maxphi, maxeta, kFALSE);
 103   fHLT->SetMergerParameters(2,3,0.003,0.1,0.05);
 104   fHLT->DoMc();
 105   fHLT->DoNonVertexTracking(); /*2 tracking passes, last without vertex contraint.*/
 106   fHLT->WriteFiles("./hlt/");
 107   fHLT->ProcessEvent(0, 35, iEvent);
 108   if(fDoBench){
 109     char filename[256];
 110     sprintf(filename, "confmap_%d",iEvent);
 111     fHLT->DoBench(filename);
 112   }
 113
 114   delete fHLT;
 115 }
 116
 117 void AliHLTReconstructor::ReconstructWithHoughTransform(AliRunLoader* runLoader,Int_t iEvent) const
 118 {
 119   Float_t ptmin = 0.1*AliL3Transform::GetSolenoidField();
 120
 121   Float_t zvertex = 0;
 122   TArrayF mcVertex(3);
 123   runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(mcVertex);
 124   zvertex = mcVertex[2];
 125
 126   LOG(AliL3Log::kInformational,"AliHLTReconstructor::Reconstruct","HoughTransform")
 127     <<" Hough Transform will run with ptmin="<<ptmin<<" and zvertex="<<zvertex<<ENDLOG;
 128
 129   AliL3Hough *hough = new AliL3Hough();
 130
 131   hough->SetThreshold(4);
 132   hough->SetTransformerParams(76,140,ptmin,-1);
 133   hough->SetPeakThreshold(70,-1);
 134   hough->SetRunLoader(runLoader);
 135   hough->Init("./", kFALSE, 100, kFALSE,4,0,0,zvertex);
 136   hough->SetAddHistograms();
 137
 138   for(Int_t slice=0; slice<=35; slice++)
 139     {
 140       //cout<<"Processing slice "<<slice<<endl;
 141       hough->ReadData(slice,iEvent);
 142       hough->Transform();
 143       hough->AddAllHistogramsRows();
 144       hough->FindTrackCandidatesRow();
 145       //hough->WriteTracks(slice,"./hough");
 146       hough->AddTracks();
 147     }
 148   hough->WriteTracks("./hough");
 149
 150   if(fDoBench){
 151     char filename[256];
 152     sprintf(filename, "hough_%d",iEvent);
 153     hough->DoBench(filename);
 154   }
 155   delete hough;
 156 }
 157
 158 void AliHLTReconstructor::FillESD(AliRunLoader* runLoader,
 159                                   AliESD* esd) const
 160 {
 161   Int_t iEvent = runLoader->GetEventNumber();
 162
 163   if(fDoTracker) FillESDforConformalMapping(esd,iEvent);
 164   if(fDoHough) FillESDforHoughTransform(esd,iEvent);
 165
 166   if(fDoCleanUp){
 167     char name[256];
 168     gSystem->Exec("rm -rf hlt");
 169     sprintf(name, "rm -f confmap_%d.root confmap_%d.dat",iEvent,iEvent);
 170     gSystem->Exec(name);
 171     gSystem->Exec("rm -rf hough");
 172     sprintf(name, "rm -f hough_%d.root hough_%d.dat",iEvent,iEvent);
 173     gSystem->Exec(name);
 174   }
 175 }
 176
 177 void AliHLTReconstructor::FillESDforConformalMapping(AliESD* esd,Int_t iEvent) const
 178 {
 179   //Assign MC labels for found tracks
 180   Int_t slicerange[2]={0,35};
 181   Int_t good = (int)(0.4*AliL3Transform::GetNRows());
 182   Int_t nclusters = (int)(0.4*AliL3Transform::GetNRows());
 183   Float_t ptmin = 0.;
 184   Float_t ptmax = 0.;
 185   Float_t maxfalseratio = 0.1;
 186
 187   AliL3Evaluate *fHLTEval = new AliL3Evaluate("./hlt",nclusters,good,ptmin,ptmax,slicerange);
 188   fHLTEval->SetMaxFalseClusters(maxfalseratio);
 189   fHLTEval->LoadData(iEvent,kTRUE);
 190   fHLTEval->AssignPIDs();
 191   fHLTEval->AssignIDs();
 192   AliL3TrackArray *fTracks = fHLTEval->GetTracks();
 193   if(!fTracks){
 194     delete fHLTEval;
 195     return;
 196   }
 197   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 198     {
 199       AliL3Track *tpt = (AliL3Track *)fTracks->GetCheckedTrack(i);
 200       if(!tpt) continue;
 201
 202       AliESDHLTtrack *esdtrack = new AliESDHLTtrack() ;
 203
 204       esdtrack->SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
 205       esdtrack->SetNHits(tpt->GetNHits());
 206       esdtrack->SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
 207       esdtrack->SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
 208       esdtrack->SetPt(tpt->GetPt());
 209       esdtrack->SetPsi(tpt->GetPsi());
 210       esdtrack->SetTgl(tpt->GetTgl());
 211       esdtrack->SetCharge(tpt->GetCharge());
 212       esdtrack->SetMCid(tpt->GetMCid());
 213       esdtrack->SetSector(tpt->GetSector());
 214       esdtrack->SetPID(tpt->GetPID());
 215       esdtrack->ComesFromMainVertex(tpt->ComesFromMainVertex());
 216
 217       esd->AddHLTConfMapTrack(esdtrack);
 218       delete esdtrack;
 219     }
 220   delete fHLTEval;
 221 }
 222
 223 void AliHLTReconstructor::FillESDforHoughTransform(AliESD* esd,Int_t iEvent) const
 224 {
 225   char filename[256];
 226   sprintf(filename,"./hough/tracks_%d.raw",iEvent);
 227
 228   AliL3FileHandler *tfile = new AliL3FileHandler();
 229   if(!tfile->SetBinaryInput(filename)){
 230     LOG(AliL3Log::kError,"AliHLTReconstructor::FillESDforHoughTransform","Input file")
 231       <<" Missing file "<<filename<<ENDLOG;
 232     return;
 233   }
 234
 235   AliL3TrackArray *fTracks = new AliL3TrackArray("AliL3HoughTrack");
 236   tfile->Binary2TrackArray(fTracks);
 237   tfile->CloseBinaryInput();
 238   delete tfile;
 239   if(!fTracks) return;
 240   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 241     {
 242       AliL3HoughTrack *tpt = (AliL3HoughTrack *)fTracks->GetCheckedTrack(i);
 243       if(!tpt) continue;
 244
 245       AliESDHLTtrack *esdtrack = new AliESDHLTtrack() ;
 246
 247       esdtrack->SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
 248       esdtrack->SetNHits(tpt->GetNHits());
 249       esdtrack->SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
 250       esdtrack->SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
 251       esdtrack->SetPt(tpt->GetPt());
 252       esdtrack->SetPsi(tpt->GetPsi());
 253       esdtrack->SetTgl(tpt->GetTgl());
 254       esdtrack->SetCharge(tpt->GetCharge());
 255       esdtrack->SetMCid(tpt->GetMCid());
 256       esdtrack->SetWeight(tpt->GetWeight());
 257       esdtrack->SetSector(tpt->GetSector());
 258       esdtrack->SetBinXY(tpt->GetBinX(),tpt->GetBinY(),tpt->GetSizeX(),tpt->GetSizeY());
 259       esdtrack->SetPID(tpt->GetPID());
 260       esdtrack->ComesFromMainVertex(tpt->ComesFromMainVertex());
 261
 262       esd->AddHLTHoughTrack(esdtrack);
 263       delete esdtrack;
 264     }
 265
 266   delete fTracks;
 267 }
 268 #endif