Modifications required for marking that the 'Original Data Present' bit should be...
[u/mrichter/AliRoot.git] / HLT / TPCLib / AliHLTTPCdEdxComponent.cxx
1 // **************************************************************************
2 // This file is property of and copyright by the ALICE HLT Project          *
3 // ALICE Experiment at CERN, All rights reserved.                           *
4 //                                                                          *
5 // Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
6 //                  for The ALICE HLT Project.                              *
7 //                                                                          *
8 // Permission to use, copy, modify and distribute this software and its     *
9 // documentation strictly for non-commercial purposes is hereby granted     *
10 // without fee, provided that the above copyright notice appears in all     *
11 // copies and that both the copyright notice and this permission notice     *
12 // appear in the supporting documentation. The authors make no claims       *
13 // about the suitability of this software for any purpose. It is            *
14 // provided "as is" without express or implied warranty.                    *
15 //                                                                          *
16 //***************************************************************************
17
18 ///  @file   AliHLTTPCdEdxComponent.cxx
19 ///  @author Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de>
20 ///  @date   June 2009
21 ///  @brief  An ITS tracker processing component for the HLT
22
23
24 /////////////////////////////////////////////////////
25 //                                                 //
26 // dEdx calculation component for the HLT TPC       //
27 //                                                 //
28 /////////////////////////////////////////////////////
29
30 #include "AliHLTTPCdEdxComponent.h"
31 #include "TStopwatch.h"
32 #include "TMath.h"
33 #include "AliCDBEntry.h"
34 #include "AliCDBManager.h"
35 #include "TObjString.h"
36 #include "TObjArray.h"
37 #include "AliHLTDataTypes.h"
38 #include "AliHLTExternalTrackParam.h"
39 #include "AliHLTGlobalBarrelTrack.h"
40 #include "AliHLTTPCTransform.h"
41 #include "AliHLTTPCSpacePointData.h"
42 #include "AliHLTTPCClusterDataFormat.h"
43 #include "AliHLTTPCTrackletDataFormat.h"
44 #include "AliHLTTPCDefinitions.h"
45 #include "AliTPCseed.h"
46 #include "AliTPCclusterMI.h"
47 #include "TGeoGlobalMagField.h"
48
49
50 /** ROOT macro for the implementation of ROOT specific class methods */
51 ClassImp( AliHLTTPCdEdxComponent )
52 AliHLTTPCdEdxComponent::AliHLTTPCdEdxComponent()
53     :
54     fSolenoidBz( 0 ),
55     fStatTime( 0 ),
56     fStatNEvents( 0 )
57 {
58   // see header file for class documentation
59   // or
60   // refer to README to build package
61   // or
62   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
63   for( int i=0; i<fkNPatches; i++ ){
64     fPatchClusters[i] = 0;    
65     fNPatchClusters[i] = 0;    
66   }
67 }
68
69 AliHLTTPCdEdxComponent::AliHLTTPCdEdxComponent( const AliHLTTPCdEdxComponent& )
70     :
71     AliHLTProcessor(),
72     fSolenoidBz( 0 ),
73     fStatTime( 0 ),
74     fStatNEvents( 0 )
75 {
76   // dummy
77   for( int i=0; i<fkNPatches; i++ ){
78     fPatchClusters[i] = 0;    
79     fNPatchClusters[i] = 0;    
80   }
81   HLTFatal( "copy constructor untested" );
82 }
83
84 AliHLTTPCdEdxComponent& AliHLTTPCdEdxComponent::operator=( const AliHLTTPCdEdxComponent& )
85 {
86   // see header file for class documentation
87   HLTFatal( "assignment operator untested" );
88   for( int i=0; i<fkNPatches; i++ ){
89     fPatchClusters[i] = 0;    
90     fNPatchClusters[i] = 0;    
91   }
92   return *this;
93 }
94
95 AliHLTTPCdEdxComponent::~AliHLTTPCdEdxComponent()
96 {
97   // see header file for class documentation
98   for( int i=0; i<fkNPatches; i++ ){
99     delete[] fPatchClusters[i];
100   }
101 }
102
103 //
104 // Public functions to implement AliHLTComponent's interface.
105 // These functions are required for the registration process
106 //
107
108 const char* AliHLTTPCdEdxComponent::GetComponentID()
109 {
110   // see header file for class documentation
111   return "TPCdEdx";
112 }
113
114 void AliHLTTPCdEdxComponent::GetInputDataTypes( vector<AliHLTComponentDataType>& list )
115 {
116   // see header file for class documentation
117   list.clear();
118   list.push_back( kAliHLTDataTypeTrack|kAliHLTDataOriginTPC );
119   list.push_back( AliHLTTPCDefinitions::fgkClustersDataType );
120 }
121
122 AliHLTComponentDataType AliHLTTPCdEdxComponent::GetOutputDataType()
123 {
124   // see header file for class documentation  
125   return kAliHLTDataTypedEdx|kAliHLTDataOriginTPC;
126 }
127
128 void AliHLTTPCdEdxComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
129 {
130   // define guess for the output data size
131   constBase = 200;       // minimum size
132   inputMultiplier = 0.1; // size relative to input
133 }
134
135 AliHLTComponent* AliHLTTPCdEdxComponent::Spawn()
136 {
137   // see header file for class documentation
138   return new AliHLTTPCdEdxComponent;
139 }
140
141 void AliHLTTPCdEdxComponent::SetDefaultConfiguration()
142 {
143   // Set default configuration for the CA tracker component
144   // Some parameters can be later overwritten from the OCDB
145
146   fSolenoidBz = -5.00668;
147   fStatTime = 0;
148   fStatNEvents = 0;
149 }
150
151 int AliHLTTPCdEdxComponent::ReadConfigurationString(  const char* arguments )
152 {
153   // Set configuration parameters for the CA tracker component from the string
154
155   int iResult = 0;
156   if ( !arguments ) return iResult;
157
158   TString allArgs = arguments;
159   TString argument;
160   int bMissingParam = 0;
161
162   TObjArray* pTokens = allArgs.Tokenize( " " );
163
164   int nArgs =  pTokens ? pTokens->GetEntries() : 0;
165
166   for ( int i = 0; i < nArgs; i++ ) {
167     argument = ( ( TObjString* )pTokens->At( i ) )->GetString();
168     if ( argument.IsNull() ) continue;
169
170     if ( argument.CompareTo( "-solenoidBz" ) == 0 ) {
171       if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
172       HLTWarning("argument -solenoidBz is deprecated, magnetic field set up globally (%f)", GetBz());
173       continue;
174     }
175
176     HLTError( "Unknown option \"%s\"", argument.Data() );
177     iResult = -EINVAL;
178   }
179   delete pTokens;
180
181   if ( bMissingParam ) {
182     HLTError( "Specifier missed for parameter \"%s\"", argument.Data() );
183     iResult = -EINVAL;
184   }
185
186   return iResult;
187 }
188
189
190 int AliHLTTPCdEdxComponent::ReadCDBEntry( const char* cdbEntry, const char* chainId )
191 {
192   // see header file for class documentation
193
194   const char* defaultNotify = "";
195
196   if ( !cdbEntry ) {
197     return 0;// need to add the HLT/ConfigTPC/TPCdEdx directory to cdb SG!!!
198     //cdbEntry = "HLT/ConfigTPC/TPCdEdx";
199     //defaultNotify = " (default)";
200     //chainId = 0;
201   }
202
203   HLTInfo( "configure from entry \"%s\"%s, chain id %s", cdbEntry, defaultNotify, ( chainId != NULL && chainId[0] != 0 ) ? chainId : "<none>" );
204   AliCDBEntry *pEntry = AliCDBManager::Instance()->Get( cdbEntry );//,GetRunNo());
205
206   if ( !pEntry ) {
207     HLTError( "cannot fetch object \"%s\" from CDB", cdbEntry );
208     return -EINVAL;
209   }
210   
211   TObjString* pString = dynamic_cast<TObjString*>( pEntry->GetObject() );
212   
213   if ( !pString ) {
214     HLTError( "configuration object \"%s\" has wrong type, required TObjString", cdbEntry );
215     return -EINVAL;
216   }
217
218   HLTInfo( "received configuration object string: \"%s\"", pString->GetString().Data() );
219
220   return  ReadConfigurationString( pString->GetString().Data() );
221 }
222
223
224 int AliHLTTPCdEdxComponent::Configure( const char* cdbEntry, const char* chainId, const char *commandLine )
225 {
226   // Configure the component
227   // There are few levels of configuration,
228   // parameters which are set on one step can be overwritten on the next step
229
230   //* read hard-coded values
231
232   SetDefaultConfiguration();
233
234   //* read the default CDB entry
235
236   int iResult1 = ReadCDBEntry( NULL, chainId );
237
238   //* read magnetic field
239
240   fSolenoidBz=GetBz();
241
242   //* read the actual CDB entry if required
243
244   int iResult2 = ( cdbEntry ) ? ReadCDBEntry( cdbEntry, chainId ) : 0;
245
246   //* read extra parameters from input (if they are)
247
248   int iResult3 = 0;
249
250   if ( commandLine && commandLine[0] != '\0' ) {
251     HLTInfo( "received configuration string from HLT framework: \"%s\"", commandLine );
252     iResult3 = ReadConfigurationString( commandLine );
253   }
254
255   // Initialise the tracker here
256
257   return iResult1 ? iResult1 : ( iResult2 ? iResult2 : iResult3 );
258 }
259
260
261
262 int AliHLTTPCdEdxComponent::DoInit( int argc, const char** argv )
263 {
264   // Configure the component
265
266   TString arguments = "";
267   for ( int i = 0; i < argc; i++ ) {
268     if ( !arguments.IsNull() ) arguments += " ";
269     arguments += argv[i];
270   }
271
272   int ret = Configure( NULL, NULL, arguments.Data() );
273
274   // Check field
275   if (!TGeoGlobalMagField::Instance()) {
276     HLTError("magnetic field not initialized, please set up TGeoGlobalMagField and AliMagF");
277     return -ENODEV;
278   }
279   
280   //AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform();
281   //if( transform ){
282   //AliTPCRecoParam *reco = transform->GetCurrentRecoParam();
283   //if( ! reco ){
284   //reco = new AliTPCRecoParam;
285   //reco->SetUseTotalCharge(0);
286   //transform->SetCurrentRecoParam( reco );
287   //}
288
289   return ret;
290 }
291
292
293 int AliHLTTPCdEdxComponent::DoDeinit()
294 {
295   // see header file for class documentation
296   return 0;
297 }
298
299
300
301 int AliHLTTPCdEdxComponent::Reconfigure( const char* cdbEntry, const char* chainId )
302 {
303   // Reconfigure the component from OCDB .
304
305   return Configure( cdbEntry, chainId, NULL );
306 }
307
308
309
310 int AliHLTTPCdEdxComponent::DoEvent
311 (
312   const AliHLTComponentEventData& evtData,
313   const AliHLTComponentBlockData* blocks,
314   AliHLTComponentTriggerData& /*trigData*/,
315   AliHLTUInt8_t* outputPtr,
316   AliHLTUInt32_t& size,
317   vector<AliHLTComponentBlockData>& outputBlocks )
318 {
319   //* process the event
320
321   AliHLTUInt32_t maxBufferSize = size;
322   size = 0; // output size
323
324   if (!IsDataEvent()) return 0;
325
326   if ( evtData.fBlockCnt <= 0 ) {
327     HLTWarning( "no blocks in event" );
328     return 0;
329   }
330   if ( !outputPtr ) {
331     return -ENOSPC;
332   }
333
334   int iResult=0;
335
336   TStopwatch timer;
337
338   // Initialise the data pointers
339
340   for( int i=0; i<fkNPatches; i++ ){
341     delete[] fPatchClusters[i];    
342     fPatchClusters[i] = 0;
343     fNPatchClusters[i] = 0;    
344   }
345
346   int nBlocks = (int)evtData.fBlockCnt;
347
348   int nInputClusters = 0;
349   int nInputTracks = 0;
350
351   // first read all the clusters
352
353   for (int ndx=0; ndx<nBlocks && iResult>=0; ndx++) {
354     const AliHLTComponentBlockData* iter = blocks+ndx;
355     if ( iter->fDataType != AliHLTTPCDefinitions::fgkClustersDataType ) continue;
356     Int_t slice=AliHLTTPCDefinitions::GetMinSliceNr(iter->fSpecification);
357     Int_t patch=AliHLTTPCDefinitions::GetMinPatchNr(iter->fSpecification);
358     Int_t slicepatch=slice*6+patch;
359     if( slicepatch >= fkNPatches ){
360       HLTWarning("Wrong header of TPC cluster data, slice %d, patch %d",
361                  slice, patch );
362       continue;
363     }
364     AliHLTTPCClusterData* inPtrSP = ( AliHLTTPCClusterData* )( iter->fPtr );
365     nInputClusters += inPtrSP->fSpacePointCnt;
366
367     delete[] fPatchClusters[slicepatch];
368     fPatchClusters[slicepatch] = new AliTPCclusterMI[inPtrSP->fSpacePointCnt];
369     fNPatchClusters[slicepatch] = inPtrSP->fSpacePointCnt;
370     
371     // create  off-line clusters out of the HLT clusters
372     // todo: check which cluster information is really needed for the dEdx
373
374     for ( unsigned int i = 0; i < inPtrSP->fSpacePointCnt; i++ ) {
375       AliHLTTPCSpacePointData *chlt = &( inPtrSP->fSpacePoints[i] );
376       AliTPCclusterMI *c = fPatchClusters[slicepatch]+i;
377       c->SetX(chlt->fX);
378       c->SetY(chlt->fY);
379       c->SetZ(chlt->fZ);
380       c->SetSigmaY2(chlt->fSigmaY2);
381       c->SetSigmaYZ( 0 );
382       c->SetSigmaZ2(chlt->fSigmaZ2);
383       c->SetQ( chlt->fCharge );
384       c->SetMax( chlt->fQMax );
385       Int_t sector, row;
386       Float_t padtime[3]={0,chlt->fY,chlt->fZ};
387       AliHLTTPCTransform::Slice2Sector(slice,chlt->fPadRow, sector, row);
388       AliHLTTPCTransform::Local2Raw( padtime, sector, row);
389       c->SetDetector( sector );
390       c->SetRow( row );
391       c->SetPad( (Int_t) padtime[1] );
392       c->SetTimeBin( (Int_t) padtime[2] );
393     }
394   }
395
396
397   // loop over the input tracks: calculate dEdx and write output
398
399   unsigned int outSize = 0;
400   AliHLTFloat32_t *outPtr = ( AliHLTFloat32_t * )( outputPtr );
401
402   for (const AliHLTComponentBlockData* pBlock=GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginTPC);
403        pBlock!=NULL; pBlock=GetNextInputBlock()) {
404     
405     AliHLTComponentBlockData outBlock;
406     FillBlockData( outBlock );
407     outBlock.fOffset = outSize;
408     outBlock.fSize = 0;
409     outBlock.fDataType = kAliHLTDataTypedEdx|kAliHLTDataOriginTPC;
410     outBlock.fSpecification = pBlock->fSpecification;
411   
412     AliHLTTracksData* dataPtr = ( AliHLTTracksData* ) pBlock->fPtr;
413     int nTracks = dataPtr->fCount;
414     AliHLTExternalTrackParam* currTrack = dataPtr->fTracklets;
415     nInputTracks+=nTracks;
416
417     for( int itr=0; 
418          itr<nTracks && ( (AliHLTUInt8_t *)currTrack < ((AliHLTUInt8_t *) pBlock->fPtr)+pBlock->fSize); 
419          itr++ ){
420
421       // create an off-line track
422       AliHLTGlobalBarrelTrack gb(*currTrack);
423       AliTPCseed tTPC;
424       tTPC.Set( gb.GetX(), gb.GetAlpha(),
425                 gb.GetParameter(), gb.GetCovariance() );
426
427       // set the clusters
428       
429       for( UInt_t ic=0; ic<currTrack->fNPoints; ic++){      
430         UInt_t id = currTrack->fPointIDs[ic];
431         int iSlice = AliHLTTPCSpacePointData::GetSlice(id);
432         int iPatch = AliHLTTPCSpacePointData::GetPatch(id);
433         int iCluster = AliHLTTPCSpacePointData::GetNumber(id);
434         if( iSlice<0 || iSlice>36 || iPatch<0 || iPatch>5 ){
435           HLTError("Corrupted TPC cluster Id: slice %d, patch %d, cluster %d",
436                    iSlice, iPatch,iCluster );
437           continue;
438         }
439         AliTPCclusterMI *patchClusters = fPatchClusters[iSlice*6 + iPatch];
440         if( !patchClusters ){
441           HLTError("Clusters are missed for slice %d, patch %d",
442                    iSlice, iPatch );
443           continue;
444         }
445         if( iCluster >= fNPatchClusters[iSlice*6 + iPatch] ){
446           HLTError("TPC slice %d, patch %d: ClusterID==%d >= N Cluaters==%d ",
447                    iSlice, iPatch,iCluster, fNPatchClusters[iSlice*6 + iPatch] );
448           continue;
449         }
450         AliTPCclusterMI *c = &(patchClusters[iCluster]);          
451         
452         int sec = c->GetDetector();
453         int row = c->GetRow();
454         if ( sec >= 36 ) {
455           row = row + AliHLTTPCTransform::GetNRowLow();
456         }
457         tTPC.SetClusterPointer( row, c );
458         
459         AliTPCTrackerPoint &point = *( tTPC.GetTrackPoint( row ) );
460         tTPC.Propagate( TMath::DegToRad()*(sec%18*20.+10.), c->GetX(), fSolenoidBz );
461         Double_t angle2 = tTPC.GetSnp()*tTPC.GetSnp();
462         angle2 = (angle2<1) ?TMath::Sqrt(angle2/(1-angle2)) :10.; 
463         point.SetAngleY( angle2 );
464         point.SetAngleZ( tTPC.GetTgl() );
465       }
466
467       // Cook dEdx
468
469       if( outSize+3*sizeof( AliHLTFloat32_t ) > maxBufferSize ){
470         HLTWarning( "Output buffer size %d exceed", maxBufferSize );
471         iResult = -ENOSPC;
472         break;
473       }
474
475       //Old method
476       /*
477       tTPC.CookdEdx( 0.02, 0.6 );      
478       outPtr[0] = tTPC.GetdEdx();
479       outPtr[1] = tTPC.GetSDEDX(0);
480       outPtr[2] = tTPC.GetNCDEDX(0);
481       */
482
483       //New method
484       outPtr[0] = tTPC.CookdEdxAnalytical(0.02,0.6,1,0,159,0);
485       outPtr[1] = 0.;
486       outPtr[2] = 159;
487       
488       outPtr+=3;
489       outSize+=3*sizeof( AliHLTFloat32_t );    
490       outBlock.fSize+=3*sizeof( AliHLTFloat32_t );  
491
492       unsigned int step = sizeof( AliHLTExternalTrackParam ) + currTrack->fNPoints * sizeof( unsigned int );
493       currTrack = ( AliHLTExternalTrackParam* )( (( Byte_t * )currTrack) + step );  
494     }
495     if( iResult<0 ) break;
496     outputBlocks.push_back( outBlock );    
497     size = outSize;
498   }
499
500   timer.Stop();
501   fStatTime += timer.RealTime();
502   fStatNEvents++;
503   
504   // Set log level to "Warning" for on-line system monitoring
505   int hz = ( int ) ( fStatTime > 1.e-10 ? fStatNEvents / fStatTime : 100000 );
506
507   HLTInfo( "TPCdEdx : %d clusters and %d tracks processed; average time %d Hz",
508            nInputClusters, nInputTracks, hz );
509   
510   return iResult;
511 }