Using common HLT track data format for TRD tracks (Theodor)
[u/mrichter/AliRoot.git] / HLT / TRD / AliHLTTRDTrackerV1Component.cxx
1 // $Id: AliHLTTRDTrackerV1Component.cxx 23618 2008-01-29 13:07:38Z hristov $
2
3 /**************************************************************************
4  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5  *                                                                        *
6  * Authors: Matthias Richter <Matthias.Richter@ift.uib.no>                *
7  *          Timm Steinbeck <timm@kip.uni-heidelberg.de>                   *
8  *          for The ALICE Off-line 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 /** @file   AliHLTTRDTrackerV1Component.cxx
20     @author Timm Steinbeck, Matthias Richter
21     @date   
22     @brief  A TRDTrackerV1 processing component for the HLT. */
23
24 #if __GNUC__ >= 3
25 using namespace std;
26 #endif
27
28 #include "AliHLTTRDTrackerV1Component.h"
29 #include "AliHLTTRDDefinitions.h"
30 #include "AliHLTTRDCluster.h"
31 #include "AliHLTTRDTrack.h"
32 #include "AliHLTTRDUtils.h"
33
34 #include "TFile.h"
35 #include "TChain.h"
36
37 #include "AliGeomManager.h"
38 #include "AliCDBManager.h"
39 #include "AliESDEvent.h"
40 #include "AliMagF.h"
41 #include "AliESDfriend.h"
42
43 #include "AliTRDcalibDB.h"
44 #include "AliTRDReconstructor.h"
45 #include "AliTRDtrackerV1.h"
46 #include "AliTRDrecoParam.h"
47
48 #include <cstdlib>
49 #include <cerrno>
50 #include <string>
51
52 #ifdef HAVE_VALGRIND_CALLGRIND_H
53 #include <valgrind/callgrind.h>
54 #else
55 #define CALLGRIND_START_INSTRUMENTATION do { } while (0)
56 #define CALLGRIND_STOP_INSTRUMENTATION do { } while (0)
57 #endif
58
59 ClassImp(AliHLTTRDTrackerV1Component)
60     
61 AliHLTTRDTrackerV1Component::AliHLTTRDTrackerV1Component():
62   AliHLTProcessor(),
63   fOutputPercentage(100), // By default we copy to the output exactly what we got as input  
64   fStrorageDBpath("local://$ALICE_ROOT/OCDB"),
65   fCDB(NULL),
66   fGeometryFileName(""),
67   fTracker(NULL),
68   fRecoParam(NULL),
69   fReconstructor(NULL)
70 {
71   // Default constructor
72
73   fGeometryFileName = getenv("ALICE_ROOT");
74   fGeometryFileName += "/HLT/TRD/geometry.root";
75 }
76
77 AliHLTTRDTrackerV1Component::~AliHLTTRDTrackerV1Component()
78 {
79   // Destructor
80 }
81
82 const char* AliHLTTRDTrackerV1Component::GetComponentID()
83 {
84   // Return the component ID const char *
85   return "TRDTrackerV1"; // The ID of this component
86 }
87
88 void AliHLTTRDTrackerV1Component::GetInputDataTypes( vector<AliHLTComponent_DataType>& list)
89 {
90   // Get the list of input data  
91   list.clear(); // We do not have any requirements for our input data type(s).
92   list.push_back( AliHLTTRDDefinitions::fgkClusterDataType );
93 }
94
95 AliHLTComponent_DataType AliHLTTRDTrackerV1Component::GetOutputDataType()
96 {
97   // Get the output data type
98   //return AliHLTTRDDefinitions::fgkClusterDataType;
99   return  kAliHLTDataTypeTrack | kAliHLTDataOriginTRD;
100 }
101
102 void AliHLTTRDTrackerV1Component::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
103 {
104   // Get the output data size
105   constBase = 0;
106   inputMultiplier = ((double)fOutputPercentage)/100.0;
107 }
108
109 // Spawn function, return new instance of this class
110 AliHLTComponent* AliHLTTRDTrackerV1Component::Spawn()
111 {
112   return new AliHLTTRDTrackerV1Component;
113 };
114
115
116 int AliHLTTRDTrackerV1Component::DoInit( int argc, const char** argv )
117 {
118   // perform initialization. We check whether our relative output size is specified in the arguments.
119   fOutputPercentage = 100;
120   int i = 0;
121   char* cpErr;
122   
123
124   Int_t iRecoParamType = -1; // default will be the low flux
125   Int_t iNtimeBins = -1;     // number of time bins for the tracker to use
126   Int_t iMagneticField = -1; // magnetic field: 0==OFF and 1==ON
127   Int_t iPIDmethod = 1;      // 0=LikelyHood(LH) 1=NeuronalNetwork(NN) 2=TruncatedMean(TM)
128   Bool_t bHLTMode = kTRUE, bWriteClusters = kFALSE;
129   
130   while ( i < argc )
131     {
132       HLTDebug("argv[%d] == %s", i, argv[i] );
133       if ( !strcmp( argv[i], "output_percentage" ) )
134         {
135           if ( i+1>=argc )
136             {
137               HLTError("Missing output_percentage parameter");
138               return ENOTSUP;
139             }
140           HLTDebug("argv[%d+1] == %s", i, argv[i+1] );
141           fOutputPercentage = strtoul( argv[i+1], &cpErr, 0 );
142           if ( *cpErr )
143             {
144               HLTError("Cannot convert output_percentage parameter '%s'", argv[i+1] );
145               return EINVAL;
146             }
147           HLTInfo("Output percentage set to %lu %%", fOutputPercentage );
148           i += 2;
149         }
150       else if ( !strcmp( argv[i], "-NTimeBins" ) )
151         {
152           if ( i+1>=argc )
153             {
154               HLTError("Missing -NTimeBins parameter");
155               return ENOTSUP;
156             }
157           HLTDebug("Arguments", "argv[%d+1] == %s", i, argv[i+1] );
158           iNtimeBins = strtoul( argv[i+1], &cpErr, 0 );
159           if ( *cpErr )
160             {
161               HLTError("Wrong Argument. Cannot convert -NTimeBins parameter '%s'", argv[i+1] );
162               return EINVAL;
163             }     
164           i += 2; 
165         }
166       else if ( strcmp( argv[i], "-cdb" ) == 0)
167         {
168           if ( i+1 >= argc )
169             {
170               HLTError( "Missing -cdb argument");
171               return ENOTSUP;         
172             }
173           fStrorageDBpath = argv[i+1];
174           HLTInfo("DB storage is %s", fStrorageDBpath.c_str() );          
175           i += 2;
176         }
177       else if ( strcmp( argv[i], "-geometry" ) == 0)
178         {
179           if ( i+1 >= argc )
180             {
181               HLTError("Missing -geometry argument");
182               return ENOTSUP;         
183             }
184           fGeometryFileName = argv[i+1];
185           HLTInfo("GeomFile storage is %s", 
186                   fGeometryFileName.c_str() );    
187           i += 2;
188         }
189       // the flux parametrizations
190       else if ( strcmp( argv[i], "-lowflux" ) == 0)
191         {
192           iRecoParamType = 0;     
193           HLTDebug("Low flux reco selected.");
194           i++;
195         }
196       else if ( strcmp( argv[i], "-highflux" ) == 0)
197         {
198           iRecoParamType = 1;     
199           HLTDebug("Low flux reco selected.");
200           i++;
201         }
202       else if ( strcmp( argv[i], "-cosmics" ) == 0)
203         {
204           iRecoParamType = 2;     
205           HLTDebug("Cosmic test reco selected.");
206           i++;
207         }
208       else if ( strcmp( argv[i], "-magnetic_field_ON" ) == 0)
209         {
210           iMagneticField = 1;
211           i++;
212         }
213       else if ( strcmp( argv[i], "-magnetic_field_OFF" ) == 0)
214         {
215           iMagneticField = 0;
216           i++;
217         }
218       else if ( strcmp( argv[i], "-offlineMode" ) == 0)
219         {
220           bHLTMode=kFALSE;
221           HLTDebug("Using standard offline tracking.");
222           i++;
223         }
224       else if ( strcmp( argv[i], "-PIDmethod" ) == 0)
225         {
226           if ( i+1 >= argc )
227             {
228               HLTError("Missing -PIDmethod argument");
229               return ENOTSUP;         
230             }
231           if( strcmp(argv[i], "LH") )
232             iPIDmethod=0;
233           else if( strcmp(argv[i], "NN") )
234             iPIDmethod=1;
235           else if( strcmp(argv[i], "TM") )
236             iPIDmethod=2;
237           else {
238             HLTError("Unknown -PIDmethod argument");
239             return ENOTSUP;           
240           }
241           i += 2;
242         }
243
244       else {
245         HLTError("Unknown option '%s'", argv[i] );
246         return EINVAL;
247       }
248       
249     }
250
251   // THE "REAL" INIT COMES HERE
252   // offline condition data base
253   fCDB = AliCDBManager::Instance();
254   if (!fCDB)
255     {
256       HLTError("Could not get CDB instance", "fCDB 0x%x", fCDB);
257       return -1;
258     }
259   else
260     {
261       fCDB->SetRun(0); // THIS HAS TO BE RETRIEVED !!!
262       fCDB->SetDefaultStorage(fStrorageDBpath.c_str());
263       HLTDebug("CDB instance", "fCDB 0x%x", fCDB);
264     }
265
266   // check if the N of time bins make sense
267   if (iNtimeBins <= 0)
268     {
269       HLTError("Sorry. Tracker needs number of time bins. At the moment you have to provide it with -NTimeBins <value>. The simulation always had 24 and the real data 30. Take your pick. Make sure the information is correct. Ask offline to implement how to propagate this information into clusters/cluster tree.");
270       return -1;
271     }
272
273   if (iNtimeBins < 24 || iNtimeBins > 30)
274     {
275       HLTWarning("The number of time bins seems to be strange = %d. But okay. Let's try it...", iNtimeBins);
276     }
277
278   HLTDebug("The number of time bins = %d.", iNtimeBins);
279   AliTRDtrackerV1::SetNTimeBins(iNtimeBins);
280
281   // !!!! THIS IS IMPORTANT
282   // init alifield map - temporarly via parameter - should come from a DB or DCS ?
283   // !!!! 
284   if (iMagneticField < 0)
285     {
286       iMagneticField = 0;
287       HLTWarning("No magnetic field switch stated. Use -magnetic_field_ON or -magnetic_field_OFF flag. Defaulting to OFF = NO MAGNETIC FIELD");
288     }
289   
290   if (!TGeoGlobalMagField::Instance()->IsLocked()) {
291     if (iMagneticField == 0)
292       {
293         // magnetic field OFF
294         AliMagF* field = new AliMagF("Maps","Maps",2,0.,0., 10.,AliMagF::k5kGUniform);
295         TGeoGlobalMagField::Instance()->SetField(field);
296         HLTDebug("Magnetic field is OFF.");
297       }
298     
299     if (iMagneticField == 1)
300       {
301         // magnetic field ON
302         AliMagF* field = new AliMagF("Maps","Maps",2,1.,1., 10.,AliMagF::k5kG);
303         TGeoGlobalMagField::Instance()->SetField(field);
304         HLTDebug("Magnetic field is ON.");
305       }
306   }
307   else {
308     HLTError("Magnetic field is already set and locked, cannot redefine it." );
309   }
310
311   // reconstruction parameters
312   if (iRecoParamType < 0 || iRecoParamType > 2)
313     {
314       HLTWarning("No reco param selected. Use -lowflux -highflux -cosmics flags. Defaulting to low flux.");
315       iRecoParamType = 0;
316     }
317
318   if (iRecoParamType == 0)
319     {
320       fRecoParam = AliTRDrecoParam::GetLowFluxParam();
321       HLTDebug("Low flux params init.");
322     }
323
324   if (iRecoParamType == 1)
325     {
326       fRecoParam = AliTRDrecoParam::GetHighFluxParam();
327       HLTDebug("High flux params init.");
328     }
329   
330   if (iRecoParamType == 2)
331     {
332       fRecoParam = AliTRDrecoParam::GetCosmicTestParam();
333       HLTDebug("Cosmic Test params init.");
334     }
335
336   if (fRecoParam == 0)
337     {
338       HLTError("No reco params initialized. Sniffing big trouble!");
339       return -1;
340     }
341
342   fReconstructor = new AliTRDReconstructor();
343   //   fRecoParam->SetChi2Y(.1);
344   //   fRecoParam->SetChi2Z(5.);
345   fReconstructor->SetRecoParam(fRecoParam);
346   // write clusters [cw] = true
347   // track seeding (stand alone tracking) [sa] = true
348   // PID method in reconstruction (NN) [nn] = true
349   // write online tracklets [tw] = false
350   // drift gas [ar] = false
351   // sl_tr_0 = StreamLevel_task_Level
352   //  fReconstructor->SetOption("sa,!cw,hlt,sl_tr_0");
353   TString recoOptions="sa,sl_tr_0";
354   
355   if (bWriteClusters)
356     {
357       recoOptions += ",cw";
358     } 
359   else
360     {
361       recoOptions += ",!cw";
362     }
363   if (bHLTMode)
364     recoOptions += ",hlt";
365
366   switch(iPIDmethod){
367   case 0: recoOptions += ",!nn"; break;
368   case 1: recoOptions += ",nn"; break;
369   case 2: recoOptions += ",!nn"; break;
370   }
371   
372   fReconstructor->SetOption(recoOptions.Data());
373   HLTDebug("Reconstructor options are: %s",recoOptions.Data());
374   
375   if((AliGeomManager::GetGeometry()) == NULL){
376     
377     if ( TFile::Open(fGeometryFileName.c_str())) {
378       AliGeomManager::LoadGeometry(fGeometryFileName.c_str());
379     }
380     else {
381       HLTError("Cannot load geometry from file %s",fGeometryFileName.c_str());
382       return EINVAL;
383     }
384   }
385   else
386     HLTInfo("Geometry Already Loaded");
387   
388   // create the tracker
389   fTracker = new AliTRDtrackerV1();
390   fTracker->SetReconstructor(fReconstructor);
391   HLTDebug("TRDTracker at 0x%x", fTracker);
392
393   if (fTracker == 0)
394     {
395       HLTError("Unable to create the tracker!");
396       return -1;
397     }
398
399   return 0;
400 }
401
402 int AliHLTTRDTrackerV1Component::DoDeinit()
403 {
404   // Deinitialization of the component
405
406   fTracker->SetClustersOwner(kFALSE);
407   delete fTracker;
408   fTracker = 0x0;
409   
410   // We need to set clusters in Reconstructor to null to prevent from 
411   // double deleting, since we delete TClonesArray by ourself in DoEvent.
412   fReconstructor->SetClusters(0x0);
413   delete fReconstructor;
414   fReconstructor = 0x0;
415   
416   AliTRDcalibDB::Terminate();
417
418   return 0;
419 }
420
421 int AliHLTTRDTrackerV1Component::DoEvent( const AliHLTComponentEventData& evtData, 
422                                           const AliHLTComponentBlockData* blocks, 
423                                           AliHLTComponent_TriggerData& /*trigData*/, 
424                                           AliHLTUInt8_t* outputPtr, 
425                                           AliHLTUInt32_t& size, 
426                                           vector<AliHLTComponent_BlockData>& outputBlocks )
427 {
428   // Process an event
429   
430   HLTDebug("NofBlocks %lu", evtData.fBlockCnt );
431   
432   AliHLTUInt32_t totalSize = 0, offset = 0;
433   AliHLTUInt32_t dBlockSpecification = 0;
434
435   vector<AliHLTComponent_DataType> expectedDataTypes;
436   GetInputDataTypes(expectedDataTypes);
437   if (evtData.fEventID > 1)
438     CALLGRIND_START_INSTRUMENTATION;
439   for ( unsigned long iBlock = 0; iBlock < evtData.fBlockCnt; iBlock++ ) 
440     {
441       const AliHLTComponentBlockData &block = blocks[iBlock];
442       AliHLTComponentDataType inputDataType = block.fDataType;
443       Bool_t correctDataType = kFALSE;
444
445       for(UInt_t i = 0; i < expectedDataTypes.size(); i++){
446         if( expectedDataTypes.at(i) == inputDataType)
447           correctDataType = kTRUE;
448       }
449       if (!correctDataType)
450         {
451           HLTDebug( "Block # %i/%i; Event 0x%08LX (%Lu) Wrong received datatype: %s - Skipping",
452                     iBlock, evtData.fBlockCnt-1,
453                     evtData.fEventID, evtData.fEventID, 
454                     DataType2Text(inputDataType).c_str());
455           continue;
456         }
457       else {
458         HLTDebug("We get the right data type: Block # %i/%i; Event 0x%08LX (%Lu) Received datatype: %s; Block Size: %i",
459                  iBlock, evtData.fBlockCnt-1,
460                  evtData.fEventID, evtData.fEventID, 
461                  DataType2Text(inputDataType).c_str(),
462                  block.fSize);
463       }
464       
465       
466       TClonesArray* clusterArray = new TClonesArray("AliTRDcluster"); // would be nice to allocate memory for all clusters here.
467       AliHLTTRDUtils::ReadClusters(clusterArray, block.fPtr, block.fSize);
468       HLTDebug("TClonesArray of clusters: nbEntries = %i", clusterArray->GetEntriesFast());
469       fTracker->LoadClusters(clusterArray);
470
471       // maybe it is not so smart to create it each event? clear is enough ?
472       AliESDEvent *esd = new AliESDEvent();
473       esd->CreateStdContent();
474       fTracker->Clusters2Tracks(esd);
475
476       Int_t nTracks = esd->GetNumberOfTracks();
477       HLTInfo("Number of tracks  == %d ==", nTracks);  
478
479       TClonesArray* trdTracks = 0x0;
480       //trdTracks = fTracker->GetListOfTracks();
481       
482       if(nTracks>0){
483         HLTDebug("We have an output ESDEvent: 0x%x with %i tracks", esd, nTracks);
484         AliHLTUInt32_t addedSize = AliHLTTRDUtils::AddESDToOutput(esd, (AliHLTUInt8_t*)(outputPtr+offset));
485         totalSize += addedSize;
486           
487         // Fill block 
488         AliHLTComponentBlockData bd;
489         FillBlockData( bd );
490         //bd.fPtr = outputPtr;
491         bd.fOffset = offset;
492         bd.fSize = addedSize;
493         bd.fSpecification = dBlockSpecification;
494         bd.fDataType = kAliHLTDataTypeTrack | kAliHLTDataOriginTRD;
495         outputBlocks.push_back( bd );
496         HLTDebug("BD fPtr 0x%x, fOffset %i, fSize %i, fSpec 0x%x", bd.fPtr, bd.fOffset, bd.fSize, bd.fSpecification);
497         offset = totalSize;
498
499         if (trdTracks){
500           //Int_t nbTracks=trdTracks->GetEntriesFast();
501           //if (nbTracks>0){
502           HLTDebug("We have an output array: pointer to trdTracks = 0x%x, nbEntries = %i", trdTracks, trdTracks->GetEntriesFast());
503           
504           AliHLTUInt32_t addedSize = AliHLTTRDUtils::AddTracksToOutput(trdTracks, (AliHLTUInt8_t*)(outputPtr+offset));
505           totalSize += addedSize;
506           
507           // Fill block 
508           AliHLTComponentBlockData bd;
509           FillBlockData( bd );
510           //bd.fPtr = outputPtr;
511           bd.fOffset = offset;
512           bd.fSize = addedSize;
513           bd.fSpecification = dBlockSpecification;
514           bd.fDataType = AliHLTTRDDefinitions::fgkTRDSATracksDataType;
515           outputBlocks.push_back( bd );
516           HLTDebug("BD fPtr 0x%x, fOffset %i, fSize %i, fSpec 0x%x", bd.fPtr, bd.fOffset, bd.fSize, bd.fSpecification);
517           offset = totalSize;
518         }
519       }
520       
521       
522       // if (trdTracks)
523       //        totalSize += TransportTracks(trdTracks, outputPtr, outputBlocks, offset, dBlockSpecification);
524       // else {
525       //          HLTDebug("Bad array trdTracks = 0x%x", trdTracks);
526       // }
527       
528       HLTDebug("totalSize: %i", totalSize);
529       
530 //       if ( totalSize > allocSize )
531 //      {
532 //        HLTError("Too much data; Data written over allowed buffer. Amount written: %lu, allowed amount: %lu.",
533 //        totalSize, size );
534 //        return EMSGSIZE;
535 //      }
536
537       //here we are deleting clusters (but not the TClonesArray itself)
538       fTracker->UnloadClusters();
539       AliTRDReconstructor::SetClusters(0x0);
540       delete esd;
541       clusterArray->Delete();
542       delete clusterArray;
543       
544       }
545       
546   size = totalSize;
547   HLTDebug("Event is done. size written to the output is %i", size);
548   return 0;
549 }