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