HLT/src/AliHLTReconstructor.cxx

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