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