SPD vertexer component created which reconstructs the XYZ primary vertex with SPD...
authorsgorbuno <sgorbuno@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 2 Nov 2009 22:56:21 +0000 (22:56 +0000)
committersgorbuno <sgorbuno@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 2 Nov 2009 22:56:21 +0000 (22:56 +0000)
The resolutions are about 400um.
It can treat an initial beam offset up to 1.5 cm in XY.
The component can do a calibration of the beam position for every N events
The speed on pp is ~5000Hz without producing histogramms and ~800 Hz with histograms

HLT/BASE/AliHLTDataTypes.cxx
HLT/BASE/AliHLTDataTypes.h
HLT/ITS/AliHLTITSAgent.cxx
HLT/ITS/AliHLTITSVertexerSPDComponent.cxx [new file with mode: 0644]
HLT/ITS/AliHLTITSVertexerSPDComponent.h [new file with mode: 0644]
HLT/ITS/macros/rec-vtxSPD.C [new file with mode: 0644]
HLT/global/AliHLTGlobalEsdConverterComponent.cxx
HLT/libAliHLTITS.pkg

index 73b90cf..ab9c60d 100644 (file)
@@ -94,6 +94,10 @@ const AliHLTComponentDataType kAliHLTDataTypeClusters = AliHLTComponentDataTypeI
 const char kAliHLTMCObjectDataTypeIDstring[8] = kAliHLTMCObjectDataTypeID;
 const AliHLTComponentDataType kAliHLTDataTypeMCObject = AliHLTComponentDataTypeInitializer(kAliHLTMCObjectDataTypeIDstring, kAliHLTDataOriginOffline);
 
+/** ESD vertex data specification */
+const char kAliHLTESDVertexDataTypeIDstring[8] = kAliHLTESDVertexDataTypeID;
+const AliHLTComponentDataType kAliHLTDataTypeESDVertex = AliHLTComponentDataTypeInitializer(kAliHLTESDVertexDataTypeIDstring, kAliHLTDataOriginAny);
+
 /** ESD data specification */
 const char kAliHLTESDObjectDataTypeIDstring[8] = kAliHLTESDObjectDataTypeID;
 const AliHLTComponentDataType kAliHLTDataTypeESDObject = AliHLTComponentDataTypeInitializer(kAliHLTESDObjectDataTypeIDstring, kAliHLTDataOriginAny);
index 943faed..ab6eb1d 100644 (file)
@@ -276,6 +276,13 @@ const int kAliHLTComponentDataTypefIDsize=8;
  */
 # define kAliHLTMCObjectDataTypeID    {'A','L','I','M','C','_','V','0'}
 
+/** ESDVertex data block
+ * an AliESDVertex object of varying origin
+ * The 'V0' at the end allows a versioning
+ * @ingroup alihlt_component_datatypes
+ */
+# define kAliHLTESDVertexDataTypeID    {'E','S','D','V','T','X','V','0'}
+
 /** ESD data block
  * an AliESD object of varying origin
  * The 'V0' at the end allows a versioning
@@ -898,6 +905,11 @@ extern "C" {
    */
   extern const AliHLTComponentDataType kAliHLTDataTypeMCObject;
 
+  /** ESD vertex object data specification, origin is 'any' 
+   * @ingroup alihlt_component_datatypes
+   */
+  extern const AliHLTComponentDataType kAliHLTDataTypeESDVertex;
+
   /** ESD object data specification, origin is 'any' 
    * @ingroup alihlt_component_datatypes
    */
index fd27fbd..902048e 100644 (file)
@@ -37,6 +37,7 @@
 #include "AliHLTITSClusterFinderComponent.h"
 #include "AliHLTITSClusterHistoComponent.h"
 #include "AliHLTITSTrackerComponent.h"
+#include "AliHLTITSVertexerSPDComponent.h"
 
 /** global instance for agent registration */
 AliHLTITSAgent gAliHLTITSAgent;
@@ -137,6 +138,7 @@ int AliHLTITSAgent::RegisterComponents(AliHLTComponentHandler* pHandler) const
   pHandler->AddComponent(new AliHLTITSClusterFinderComponent(AliHLTITSClusterFinderComponent::kClusterFinderSSD));
   pHandler->AddComponent(new AliHLTITSClusterHistoComponent);
   pHandler->AddComponent(new AliHLTITSTrackerComponent);
+  pHandler->AddComponent(new AliHLTITSVertexerSPDComponent);
 
   return 0;
 }
diff --git a/HLT/ITS/AliHLTITSVertexerSPDComponent.cxx b/HLT/ITS/AliHLTITSVertexerSPDComponent.cxx
new file mode 100644 (file)
index 0000000..ccd5ba7
--- /dev/null
@@ -0,0 +1,633 @@
+// $Id: AliHLTITSVertexerSPDComponent.cxx 32659 2009-06-02 16:08:40Z sgorbuno $
+// **************************************************************************
+// This file is property of and copyright by the ALICE HLT Project          *
+// ALICE Experiment at CERN, All rights reserved.                           *
+//                                                                          *
+// Primary Authors: Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
+//                  Ivan Kisel <kisel@kip.uni-heidelberg.de>                *
+//                  for The ALICE HLT Project.                              *
+//                                                                          *
+// Permission to use, copy, modify and distribute this software and its     *
+// documentation strictly for non-commercial purposes is hereby granted     *
+// without fee, provided that the above copyright notice appears in all     *
+// copies and that both the copyright notice and this permission notice     *
+// appear in the supporting documentation. The authors make no claims       *
+// about the suitability of this software for any purpose. It is            *
+// provided "as is" without express or implied warranty.                    *
+//                                                                          *
+//***************************************************************************
+
+///  @file   AliHLTITSVertexerSPDComponent.cxx
+///  @author Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de>
+///  @date   June 2009
+///  @brief  An ITS tracker processing component for the HLT
+
+
+/////////////////////////////////////////////////////
+//                                                 //
+// a ITS tracker processing component for the HLT  //
+//                                                 //
+/////////////////////////////////////////////////////
+
+#if __GNUC__>= 3
+using namespace std;
+#endif
+
+#include "AliHLTITSVertexerSPDComponent.h"
+#include "AliHLTArray.h"
+#include "AliExternalTrackParam.h"
+#include "TStopwatch.h"
+#include "TMath.h"
+#include "AliCDBEntry.h"
+#include "AliCDBManager.h"
+#include "AliGeomManager.h"
+#include "TObjString.h"
+#include "TObjArray.h"
+#include "AliITStrackerHLT.h"
+#include "AliHLTITSSpacePointData.h"
+#include "AliHLTITSClusterDataFormat.h"
+#include "AliHLTDataTypes.h"
+#include "AliHLTExternalTrackParam.h"
+#include "AliHLTGlobalBarrelTrack.h"
+#include "AliGeomManager.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "AliESDVertex.h"
+
+/** ROOT macro for the implementation of ROOT specific class methods */
+ClassImp( AliHLTITSVertexerSPDComponent )
+AliHLTITSVertexerSPDComponent::AliHLTITSVertexerSPDComponent()
+    :
+    fSolenoidBz( 0 ),
+    fProduceHistos(1),
+    fAutoCalibration(1000),
+    fFullTime( 0 ),
+    fRecoTime( 0 ),
+    fNEvents( 0 ),    
+    fHistoVertexXY(0),
+    fHistoVertexX(0),
+    fHistoVertexY(0),
+  fHistoVertexZ(0)
+
+{
+  // see header file for class documentation
+  // or
+  // refer to README to build package
+  // or
+  // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
+
+  fDefRunVtx[0] = 0;
+  fDefRunVtx[1] = 0;
+  fDefRunVtx[2] = 0;
+}
+
+AliHLTITSVertexerSPDComponent::AliHLTITSVertexerSPDComponent( const AliHLTITSVertexerSPDComponent& )
+    :
+    AliHLTProcessor(),
+    fSolenoidBz( 0 ),
+    fProduceHistos(1),
+    fAutoCalibration(1000),
+    fFullTime( 0 ),
+    fRecoTime( 0 ),
+    fNEvents( 0 ),    
+    fHistoVertexXY(0),
+    fHistoVertexX(0),
+    fHistoVertexY(0),
+  fHistoVertexZ(0)
+{
+  // see header file for class documentation
+  HLTFatal( "copy constructor untested" );
+}
+
+AliHLTITSVertexerSPDComponent& AliHLTITSVertexerSPDComponent::operator=( const AliHLTITSVertexerSPDComponent& )
+{
+  // see header file for class documentation
+  HLTFatal( "assignment operator untested" );
+  return *this;
+}
+
+AliHLTITSVertexerSPDComponent::~AliHLTITSVertexerSPDComponent()
+{
+  // see header file for class documentation
+  
+  delete fHistoVertexXY;
+  delete fHistoVertexX;
+  delete fHistoVertexY;
+  delete fHistoVertexZ;
+}
+
+//
+// Public functions to implement AliHLTComponent's interface.
+// These functions are required for the registration process
+//
+
+const char* AliHLTITSVertexerSPDComponent::GetComponentID()
+{
+  // see header file for class documentation
+  return "ITSVertexerSPD";
+}
+
+void AliHLTITSVertexerSPDComponent::GetInputDataTypes( vector<AliHLTComponentDataType>& list )
+{
+  // see header file for class documentation
+  list.clear();
+  list.push_back( kAliHLTDataTypeClusters|kAliHLTDataOriginITSSPD );
+}
+
+AliHLTComponentDataType AliHLTITSVertexerSPDComponent::GetOutputDataType()
+{
+  // see header file for class documentation  
+  return kAliHLTMultipleDataType;
+}
+
+int AliHLTITSVertexerSPDComponent::GetOutputDataTypes(AliHLTComponentDataTypeList& tgtList)
+
+{
+  // see header file for class documentation
+  tgtList.clear();
+  tgtList.push_back( kAliHLTDataTypeESDVertex|kAliHLTDataOriginITS);
+  tgtList.push_back(kAliHLTDataTypeHistogram);
+  return tgtList.size();
+}
+
+void AliHLTITSVertexerSPDComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
+{
+  // define guess for the output data size
+  constBase = 10000;       // minimum size
+  inputMultiplier = 0.5; // size relative to input
+}
+
+AliHLTComponent* AliHLTITSVertexerSPDComponent::Spawn()
+{
+  // see header file for class documentation
+  return new AliHLTITSVertexerSPDComponent;
+}
+
+void AliHLTITSVertexerSPDComponent::SetDefaultConfiguration()
+{
+  // Set default configuration for the CA tracker component
+  // Some parameters can be later overwritten from the OCDB
+
+  fSolenoidBz = -5.00668;
+  fProduceHistos=1;
+  fAutoCalibration = 1000;
+  fDefRunVtx[0] = 0;
+  fDefRunVtx[1] = 0;
+  fDefRunVtx[2] = 0;
+  fFullTime = 0;
+  fRecoTime = 0;
+  fNEvents = 0;
+}
+
+int AliHLTITSVertexerSPDComponent::ReadConfigurationString(  const char* arguments )
+{
+  // Set configuration parameters for the CA tracker component from the string
+
+  int iResult = 0;
+  if ( !arguments ) return iResult;
+
+  TString allArgs = arguments;
+  TString argument;
+  int bMissingParam = 0;
+
+  TObjArray* pTokens = allArgs.Tokenize( " " );
+
+  int nArgs =  pTokens ? pTokens->GetEntries() : 0;
+
+  for ( int i = 0; i < nArgs; i++ ) {
+    argument = ( ( TObjString* )pTokens->At( i ) )->GetString();
+    if ( argument.IsNull() ) continue;
+
+    if ( argument.CompareTo( "-solenoidBz" ) == 0 ) {
+      if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
+      fSolenoidBz = ( ( TObjString* )pTokens->At( i ) )->GetString().Atof();
+      HLTInfo( "Magnetic Field set to: %f", fSolenoidBz );
+      continue;
+    }
+
+    if ( argument.CompareTo( "-produceHistos" ) == 0 ) {
+      if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
+      fProduceHistos = ( ( TObjString* )pTokens->At( i ) )->GetString().Atoi();
+      HLTInfo( "fProduceHistos set to: %d", fProduceHistos );
+      continue;
+    }
+
+    if ( argument.CompareTo( "-runVertex" ) == 0 ) {
+      if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
+      fDefRunVtx[0] = ( ( TObjString* )pTokens->At( i ) )->GetString().Atof();
+      if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
+      fDefRunVtx[1] = ( ( TObjString* )pTokens->At( i ) )->GetString().Atof();
+      if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
+      fDefRunVtx[2] = ( ( TObjString* )pTokens->At( i ) )->GetString().Atof();
+      HLTInfo( "Default run vertex is set to (%f,%f,%f)",fDefRunVtx[0],
+              fDefRunVtx[1],fDefRunVtx[2] );
+      for( int i=0; i<3; i++){
+       fRunVtx[i] = fDefRunVtx[i];
+       fRunVtxNew[i] = fDefRunVtx[i];
+      }
+      fRunVtx[3] = 0.;
+      fRunVtxNew[3] = 0.;
+      continue;
+    }
+
+    if ( argument.CompareTo( "-beamDiamondCalibration" ) == 0 ) {
+      if ( ( bMissingParam = ( ++i >= pTokens->GetEntries() ) ) ) break;
+      fAutoCalibration = ( ( TObjString* )pTokens->At( i ) )->GetString().Atoi();
+      HLTInfo( "N events for recalibration of the run vertex is set to: %d", fAutoCalibration );
+      continue;
+    }
+    
+
+    HLTError( "Unknown option \"%s\"", argument.Data() );
+    iResult = -EINVAL;
+  }
+  delete pTokens;
+
+  if ( bMissingParam ) {
+    HLTError( "Specifier missed for parameter \"%s\"", argument.Data() );
+    iResult = -EINVAL;
+  }
+
+  return iResult;
+}
+
+
+int AliHLTITSVertexerSPDComponent::ReadCDBEntry( const char* cdbEntry, const char* chainId )
+{
+  // see header file for class documentation
+
+  const char* defaultNotify = "";
+
+  if ( !cdbEntry ) {
+    return 0;// need to add the HLT/ConfigITS/ITSTracker directory to cdb SG!!!
+    cdbEntry = "HLT/ConfigITS/ITSTracker";
+    defaultNotify = " (default)";
+    chainId = 0;
+  }
+
+  HLTInfo( "configure from entry \"%s\"%s, chain id %s", cdbEntry, defaultNotify, ( chainId != NULL && chainId[0] != 0 ) ? chainId : "<none>" );
+  AliCDBEntry *pEntry = AliCDBManager::Instance()->Get( cdbEntry );//,GetRunNo());
+
+  if ( !pEntry ) {
+    HLTError( "cannot fetch object \"%s\" from CDB", cdbEntry );
+    return -EINVAL;
+  }
+
+  TObjString* pString = dynamic_cast<TObjString*>( pEntry->GetObject() );
+
+  if ( !pString ) {
+    HLTError( "configuration object \"%s\" has wrong type, required TObjString", cdbEntry );
+    return -EINVAL;
+  }
+
+  HLTInfo( "received configuration object string: \"%s\"", pString->GetString().Data() );
+
+  return  ReadConfigurationString( pString->GetString().Data() );
+}
+
+
+int AliHLTITSVertexerSPDComponent::Configure( const char* cdbEntry, const char* chainId, const char *commandLine )
+{
+  // Configure the component
+  // There are few levels of configuration,
+  // parameters which are set on one step can be overwritten on the next step
+
+  //* read hard-coded values
+
+  SetDefaultConfiguration();
+
+  //* read the default CDB entry
+
+  int iResult1 = ReadCDBEntry( NULL, chainId );
+
+  //* read magnetic field
+
+  int iResult2 = ReadCDBEntry( kAliHLTCDBSolenoidBz, chainId );
+
+  //* read the actual CDB entry if required
+
+  int iResult3 = ( cdbEntry ) ? ReadCDBEntry( cdbEntry, chainId ) : 0;
+
+  //* read extra parameters from input (if they are)
+
+  int iResult4 = 0;
+
+  if ( commandLine && commandLine[0] != '\0' ) {
+    HLTInfo( "received configuration string from HLT framework: \"%s\"", commandLine );
+    iResult4 = ReadConfigurationString( commandLine );
+  }
+
+  // Initialise the tracker here
+
+  return iResult1 ? iResult1 : ( iResult2 ? iResult2 : ( iResult3 ? iResult3 : iResult4 ) );
+}
+
+
+
+int AliHLTITSVertexerSPDComponent::DoInit( int argc, const char** argv )
+{
+  // Configure the ITS tracker component
+  
+  fProduceHistos = 1;
+  fAutoCalibration = 1000;
+
+  if(AliGeomManager::GetGeometry()==NULL){
+    AliGeomManager::LoadGeometry("");
+  }
+  AliGeomManager::ApplyAlignObjsFromCDB("ITS");
+  
+  TString arguments = "";
+  for ( int i = 0; i < argc; i++ ) {
+    if ( !arguments.IsNull() ) arguments += " ";
+    arguments += argv[i];
+  }
+
+  int ret = Configure( NULL, NULL, arguments.Data() );
+
+  if( fProduceHistos ){
+    fHistoVertexXY = new TH2F("hITSvertexXY", "ITSSPD vertex in XY", 100,-2,2,100,-2,2);
+    fHistoVertexX = new TH1F("hITSvertexX", "ITSSPD vertex X", 100,-2,2);
+    fHistoVertexY = new TH1F("hITSvertexY", "ITSSPD vertex Y", 100,-2,2);
+    fHistoVertexZ = new TH1F("hITSvertexZ", "ITSSPD vertex Z", 100,-15,15);
+  }
+
+  for( int i=0; i<3; i++){
+    fRunVtx[i] = fDefRunVtx[i];
+    fRunVtxNew[i] = fDefRunVtx[i];
+  }
+  fRunVtx[3] = 0.;
+  fRunVtxNew[3] = 0.;
+
+  return ret;
+}
+
+
+int AliHLTITSVertexerSPDComponent::DoDeinit()
+{
+  // see header file for class documentation
+
+  delete fHistoVertexXY;
+  delete fHistoVertexX;
+  delete fHistoVertexY;
+  delete fHistoVertexZ;
+
+  fHistoVertexXY = 0;
+  fHistoVertexX = 0;
+  fHistoVertexY = 0;
+  fHistoVertexZ = 0;
+  fAutoCalibration = 1000;
+  fProduceHistos = 1;
+  return 0;
+}
+
+
+
+int AliHLTITSVertexerSPDComponent::Reconfigure( const char* cdbEntry, const char* chainId )
+{
+  // Reconfigure the component from OCDB .
+
+  return Configure( cdbEntry, chainId, NULL );
+}
+
+
+int AliHLTITSVertexerSPDComponent::DoEvent
+(
+  const AliHLTComponentEventData& evtData,
+  const AliHLTComponentBlockData* blocks,
+  AliHLTComponentTriggerData& /*trigData*/,
+  AliHLTUInt8_t* /*outputPtr*/,
+  AliHLTUInt32_t& size,
+  vector<AliHLTComponentBlockData>& /*outputBlocks*/ )
+{
+  //* process event
+
+  AliHLTUInt32_t maxBufferSize = size;
+  size = 0; // output size
+
+  if (!IsDataEvent()) return 0;
+
+  if ( evtData.fBlockCnt <= 0 ) {
+    HLTWarning( "no blocks in event" );
+    return 0;
+  }
+
+
+  TStopwatch timer;
+
+  // Event reconstruction in ITS
+
+  int iResult=0;
+
+  int nBlocks = evtData.fBlockCnt;
+  int nClustersTotal = 0;
+
+  const int kNPhiBins = 20;
+  const double kDPhi = TMath::TwoPi() / kNPhiBins;
+  double vtxX = fRunVtx[0], vtxY = fRunVtx[1], vtxZ = fRunVtx[2];
+
+
+  std::vector<AliHLTITSVertexerSPDComponent::AliHLTITSVZCluster> clusters[2][ kNPhiBins ];
+
+  for (int ndx=0; ndx<nBlocks && iResult>=0; ndx++) {
+
+    const AliHLTComponentBlockData* iter = blocks+ndx;
+    // Read ITS SPD clusters
+
+    if ( iter->fDataType == (kAliHLTDataTypeClusters|kAliHLTDataOriginITSSPD) ){
+
+      AliHLTITSClusterData *inPtr=reinterpret_cast<AliHLTITSClusterData*>( iter->fPtr );
+      int nClusters = inPtr->fSpacePointCnt;
+      nClustersTotal+=nClusters;
+      for( int icl=0; icl<nClusters; icl++ ){  
+       AliHLTITSSpacePointData &d = inPtr->fSpacePoints[icl];
+       if( d.fLayer>1 ) continue;// SPD only;
+       if( d.fLayer<0 ) continue;// SPD only;
+       Int_t lab[4] = { d.fTracks[0], d.fTracks[1], d.fTracks[2], d.fIndex };
+       Int_t info[3] = { d.fNy, d.fNz, d.fLayer };
+       Float_t hit[6] = { d.fY, d.fZ, d.fSigmaY2, d.fSigmaZ2, d.fQ, d.fSigmaYZ };
+       if( d.fLayer==4 ) hit[5] = -hit[5];
+
+
+       AliITSRecPoint p( lab, hit, info );
+
+       if (!p.Misalign()){
+         HLTWarning("Can not misalign an SPD cluster");
+       }
+       Float_t xyz[3];      
+       p.GetGlobalXYZ(xyz);      
+
+       // get phi bin
+       double phi = TMath::ATan2( xyz[1]-vtxY, xyz[0]-vtxX );
+       if( phi<0 ) phi+=TMath::TwoPi();
+       int iphi = (int ) phi/kDPhi;
+       if( iphi<0 ) iphi = 0;
+       if( iphi>=kNPhiBins ) iphi = kNPhiBins-1;
+       AliHLTITSVZCluster c;
+       c.fX = xyz[0];
+       c.fY = xyz[1];
+       c.fZ = xyz[2];
+       clusters[d.fLayer][iphi].push_back( c );        
+     }   
+    }    
+  }// end read input blocks
+  
+
+  // Reconstruct the vertex
+
+  TStopwatch timerReco;
+    
+
+  const double kZBinMin = -10;
+  const int kNZBins = 40;
+  const double kZStep = 1./kNZBins;
+
+  double histX[kNZBins];
+  double histY[kNZBins];
+  double histZ[kNZBins];
+  double histW[kNZBins];
+  int histN[kNZBins];
+  
+  // fit few times with decreasing of the radius cut
+
+  double vtxDR2[2] = { 2.*2., .2*.2};
+  double  maxW = 0;  
+  int bestBin = -1;
+  
+  for( int iter = 0; iter<2; iter++ ){
+    
+    for( int i=0; i<kNZBins; i++ ){
+      histX[i] = 0;
+      histY[i] = 0;
+      histZ[i] = 0;
+      histW[i] = 0;
+      histN[i] = 0;
+    }
+    for( int iPhi=0; iPhi<kNPhiBins; iPhi++ ){
+      int nCluUp = clusters[1][iPhi].size();    
+      int nCluDn = clusters[0][iPhi].size();    
+      for( int icUp=0; icUp<nCluUp; icUp++ ){
+       AliHLTITSVZCluster &cu = clusters[1][iPhi][icUp];
+       double x0 = cu.fX - vtxX;
+       double y0 = cu.fY - vtxY;
+       double z0 = cu.fZ - vtxZ;
+       double bestR2 = 1.e10;
+       int bestDn=-1, bestBin =0;
+       double bestV[3]={0,0,0};
+       for( int icDn=0; icDn<nCluDn; icDn++ ){
+         AliHLTITSVZCluster &cd = clusters[0][iPhi][icDn];
+         double x1 = cd.fX - vtxX;
+         double y1 = cd.fY - vtxY;
+         double z1 = cd.fZ - vtxZ;
+         double dx = x1 - x0;
+         double dy = y1 - y0;
+         double l2 = 1./(dx*dx + dy*dy);
+         double a = x1*y0 - x0*y1;
+         double r2 = a*a*l2; 
+         if( r2>vtxDR2[iter] ) continue;
+         double xv = -dy*a*l2;
+         double yv =  dx*a*l2;
+         double zv = ( (x1*z0-x0*z1)*dx + (y1*z0-y0*z1)*dy )*l2;
+         int zbin = (int)((zv - kZBinMin)*kZStep);
+         if( zbin<0 || zbin>=kNZBins ) continue;
+         if( r2<bestR2 ){
+           bestR2 = r2;
+           bestDn = icDn;
+           bestV[0] = xv;
+           bestV[1] = yv;
+           bestV[2] = zv;
+           bestBin = zbin;
+         }
+       }
+       if( bestDn>=0 ){
+         double w = 1./(1.+bestR2);
+         histX[bestBin]+=bestV[0]*w;
+         histY[bestBin]+=bestV[1]*w;
+         histZ[bestBin]+=bestV[2]*w;
+         histW[bestBin]+=w;
+         histN[bestBin]+=1;
+       }
+      }
+    }
+    
+    maxW = 0;  
+    bestBin = -1;
+    for( int i=0; i<kNZBins; i++ ){
+      if( histW[i]>maxW ){
+       bestBin = i;
+       maxW = histW[i];
+      }
+    }
+    if( bestBin<0 || histN[bestBin] <3 ){
+      bestBin = -1;
+      break;
+    }
+    vtxX +=histX[bestBin]/maxW;
+    vtxY +=histY[bestBin]/maxW;
+    vtxZ +=histZ[bestBin]/maxW;    
+  }
+
+  if( bestBin>=0 ){
+    double nv = 1;
+    double nNew = fRunVtxNew[3] + nv;
+    double v[3] = {vtxX, vtxY, vtxZ};
+    for( int i=0; i<3; i++){
+      fRunVtxNew[i] = ( fRunVtxNew[i]*fRunVtxNew[3] + v[i]*nv )/nNew;
+    }
+    fRunVtxNew[3] = nNew;
+  }  
+  
+  if( fAutoCalibration>0 && fRunVtxNew[3] >= fAutoCalibration ){
+    for( int i=0; i<4; i++ ){
+      fRunVtx[i] = fRunVtxNew[i];
+      fRunVtxNew[i] = 0;
+    }
+    //cout<<"ITSVertexerSPD: set run vtx to "<<fRunVtx[0]<<" "<<fRunVtx[1]<<" "<<fRunVtx[2]<<endl;  
+    if( fRunVtx[0]>3. ) fRunVtx[0] = 3;
+    if( fRunVtx[0]<-3. ) fRunVtx[0] = -3;
+    if( fRunVtx[1]>3. ) fRunVtx[1] = 3;
+    if( fRunVtx[1]<-3. ) fRunVtx[1] = -3;
+    if( fRunVtx[2]>30. ) fRunVtx[2] = 30;
+    if( fRunVtx[2]<30. ) fRunVtx[2] = -30;
+  }
+
+  timerReco.Stop();
+  
+  // Fill output 
+
+  if( bestBin>=0 ){
+    double pos[3] = {vtxX, vtxY, vtxZ};
+    double s = 400.E-4;
+    double cov[6] = {s*s,0,s*s,0,0,s*s};
+    AliESDVertex v(pos, cov, 0, histN[bestBin]);
+    PushBack( (TObject*) &v, kAliHLTDataTypeESDVertex|kAliHLTDataOriginITS,0 );
+
+    //cout<<"ITSVertexerSPD: vtx found "<<vtxX<<" "<<vtxY<<" "<<vtxZ<<endl;
+
+    if( fHistoVertexXY ) fHistoVertexXY->Fill( vtxX, vtxY );
+    if( fHistoVertexX ) fHistoVertexX->Fill( vtxX );
+    if( fHistoVertexY ) fHistoVertexY->Fill( vtxY );
+    if( fHistoVertexZ ) fHistoVertexZ->Fill( vtxZ );
+  }
+
+  if( fHistoVertexXY ) PushBack( (TObject*) fHistoVertexXY, kAliHLTDataTypeHistogram,0);
+  if( fHistoVertexX ) PushBack( (TObject*) fHistoVertexX, kAliHLTDataTypeHistogram,0);
+  if( fHistoVertexY ) PushBack( (TObject*) fHistoVertexY, kAliHLTDataTypeHistogram,0);
+  if( fHistoVertexZ ) PushBack( (TObject*) fHistoVertexZ, kAliHLTDataTypeHistogram,0);
+
+
+  timer.Stop();
+  fFullTime += timer.RealTime();
+  fRecoTime += timerReco.RealTime();
+  fNEvents++;
+
+  // Set log level to "Warning" for on-line system monitoring
+  int hz = ( int ) ( fFullTime > 1.e-10 ? fNEvents / fFullTime : 100000 );
+  int hz1 = ( int ) ( fRecoTime > 1.e-10 ? fNEvents / fRecoTime : 100000 );
+  HLTInfo( "ITS Z Vertexer: input %d clusters; time: full %d / reco %d Hz",
+          nClustersTotal, hz, hz1 );
+
+  return iResult;
+}
diff --git a/HLT/ITS/AliHLTITSVertexerSPDComponent.h b/HLT/ITS/AliHLTITSVertexerSPDComponent.h
new file mode 100644 (file)
index 0000000..29fb096
--- /dev/null
@@ -0,0 +1,159 @@
+//-*- Mode: C++ -*-
+// $Id$
+// ************************************************************************
+// This file is property of and copyright by the ALICE HLT Project        *
+// ALICE Experiment at CERN, All rights reserved.                         *
+// See cxx source for full Copyright notice                               *
+//                                                                        *
+//*************************************************************************
+
+///  @file   AliHLTITSVertexerSPDComponent.h
+///  @author Sergey Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de>
+///  @date   Nov 2009
+///  @brief  An ITS pixel vertexer component for the HLT
+
+#ifndef ALIHLTITSVERTEXERSPDCOMPONENT_H
+#define ALIHLTITSVERTEXERSPDCOMPONENT_H
+
+#include "AliHLTProcessor.h"
+#include "AliHLTDataTypes.h"
+
+class TH1F;
+class TH2F;
+
+
+/**
+ * @class AliHLTITSVertexerSPDComponent
+ * The HLT ITS SPD z-vertexer component.
+ * The vertexer uses approximate XY position of the beam
+ * It can treat initial beam offset up to 1.5 cm in XY
+ * The component does a calibration of the beam position 
+ * every n==fAutoCalibration events
+ * 
+ *
+ * <h2>General properties:</h2>
+ *
+ * Component ID: \b ITSVertexerSPD                            <br>
+ * Library: \b libAliHLTITS.so                              <br>
+ * Input Data Types:                                        <br> 
+ *    kAliHLTDataTypeClusters|kAliHLTDataOriginITSSPD       <br>
+ *      
+ * Output Data Types:                                       <br>
+ *    kAliHLTDataTypeESDVertex|kAliHLTDataOriginITS            <br>
+ *
+ * <h2>Mandatory arguments:</h2>
+ * <!-- NOTE: ignore the \li. <i> and </i>: it's just doxygen formatting -->
+ *
+ * <h2>Optional arguments:</h2>
+ * <!-- NOTE: ignore the \li. <i> and </i>: it's just doxygen formatting -->
+ *
+ * <h2>Configuration:</h2>
+ * <!-- NOTE: ignore the \li. <i> and </i>: it's just doxygen formatting -->
+ * \li -config1      <i> teststring   </i> <br>
+ *      a configuration argument with one parameter
+ * \li -config2                            <br>
+ *      a configuration argument without parameters
+ *
+ * <h2>Default CDB entries:</h2>
+ * TODO
+ *
+ * <h2>Performance:</h2>
+ * TODO
+ *
+ * <h2>Memory consumption:</h2>
+ * TODO
+ *
+ * <h2>Output size:</h2>
+ * TODO
+ * 
+ * @ingroup alihlt_its_components
+ */
+class AliHLTITSVertexerSPDComponent : public AliHLTProcessor
+{
+  public:
+    /** standard constructor */
+    AliHLTITSVertexerSPDComponent();
+
+    /** dummy copy constructor, defined according to effective C++ style */
+    AliHLTITSVertexerSPDComponent( const AliHLTITSVertexerSPDComponent& );
+
+    /** dummy assignment op, but defined according to effective C++ style */
+    AliHLTITSVertexerSPDComponent& operator=( const AliHLTITSVertexerSPDComponent& );
+
+    /** standard destructor */
+    virtual ~AliHLTITSVertexerSPDComponent();
+
+    // Public functions to implement AliHLTComponent's interface.
+    // These functions are required for the registration process
+
+    /** @see component interface @ref AliHLTComponent::GetComponentID */
+    const char* GetComponentID() ;
+
+    /** @see component interface @ref AliHLTComponent::GetInputDataTypes */
+    void GetInputDataTypes( vector<AliHLTComponentDataType>& list )  ;
+
+    /** @see component interface @ref AliHLTComponent::GetOutputDataType */
+    AliHLTComponentDataType GetOutputDataType() ;
+
+    /** @see component interface @ref AliHLTComponent::GetOutputDataTypes */
+    int GetOutputDataTypes(AliHLTComponentDataTypeList& tgtList);
+
+    /** @see component interface @ref AliHLTComponent::GetOutputDataSize */
+    virtual void GetOutputDataSize( unsigned long& constBase, double& inputMultiplier ) ;
+
+    /** @see component interface @ref AliHLTComponent::Spawn */
+    AliHLTComponent* Spawn() ;
+
+  protected:
+
+    // Protected functions to implement AliHLTComponent's interface.
+    // These functions provide initialization as well as the actual processing
+    // capabilities of the component.
+
+    /** @see component interface @ref AliHLTComponent::DoInit */
+    int DoInit( int argc, const char** argv );
+
+    /** @see component interface @ref AliHLTComponent::DoDeinit */
+    int DoDeinit();
+
+    /** reconfigure **/
+    int Reconfigure( const char* cdbEntry, const char* chainId );
+
+    /** @see component interface @ref AliHLTProcessor::DoEvent */
+    int DoEvent( const AliHLTComponentEventData& evtData, const AliHLTComponentBlockData* blocks,
+                 AliHLTComponentTriggerData& trigData, AliHLTUInt8_t* outputPtr,
+                 AliHLTUInt32_t& size, vector<AliHLTComponentBlockData>& outputBlocks );
+
+  private:
+
+  struct AliHLTITSVZCluster{
+    float fX, fY, fZ;
+  };
+    /** magnetic field */
+    double fSolenoidBz; // see above
+    Bool_t fProduceHistos;// flag to produce histogramms(def = 1)
+  Int_t   fAutoCalibration;//  n of events for recalibration of run vertex(-1=no, def = 1000);
+
+  double fDefRunVtx[3];// default vertex position
+    double fFullTime; //* total time for DoEvent() [s]
+    double fRecoTime; //* total reconstruction time [s]
+    Long_t    fNEvents;  //* number of reconstructed events
+    double fRunVtx[4]; //!
+    double fRunVtxNew[4]; //!
+
+   TH2F *fHistoVertexXY;//!
+   TH1F *fHistoVertexX;//!
+   TH1F *fHistoVertexY;//!
+   TH1F *fHistoVertexZ;//!
+    /** set configuration parameters **/
+    void SetDefaultConfiguration();
+    int ReadConfigurationString(  const char* arguments );
+    int ReadCDBEntry( const char* cdbEntry, const char* chainId );
+    int Configure( const char* cdbEntry, const char* chainId, const char *commandLine  );
+
+    ClassDef( AliHLTITSVertexerSPDComponent, 0 );
+
+};
+#endif
diff --git a/HLT/ITS/macros/rec-vtxSPD.C b/HLT/ITS/macros/rec-vtxSPD.C
new file mode 100644 (file)
index 0000000..0e3d404
--- /dev/null
@@ -0,0 +1,105 @@
+// $Id: rec-tpc-its.C 35299 2009-10-07 09:07:29Z kkanaki $
+/*
+ * Example macro to run the ITS tracker with the TPC reconstruction.
+ * The reconstruction is done from the TPC and ITS raw data.
+ *
+ * Usage:
+ * <pre>
+ *   aliroot -b -q rec-tpc-its.C | tee rec-tpc-its.log
+ *   aliroot -b -q rec-tpc-its.C'("./","spd")' | tee rec-tpc-its.log
+ * </pre>
+ *
+ * The macro asumes raw data to be available in the rawx folders, either
+ * simulated or real data. A different input can be specified as parameter
+ * <pre>
+ *   aliroot -b -q rec-tpc-its.C'("input.root")'
+ * </pre>
+ * 
+ * In the first section, an analysis chain is defined. The scale of the
+ * chain can be defined by choosing the range of sectors and partitions.
+ *
+ * The reconstruction is steered by the AliReconstruction object in the
+ * usual way.
+ *
+ * @ingroup alihlt_tpc
+ * @author Gaute.Ovrebekk@ift.uib.no
+ */
+void rec_vtxSPD(const char* input="./")
+{
+  
+  if(!gSystem->AccessPathName("galice.root")){
+    cerr << "please delete the galice.root or run at different place." << endl;
+    return;
+  }
+
+  if (!input) {
+    cerr << "please specify input or run without arguments" << endl;
+    return;
+  }
+  ///////////////////////////////////////////////////////////////////////////////////////////////////
+  //
+  // init the HLT system in order to define the analysis chain below
+  //
+  AliHLTSystem* gHLT=AliHLTPluginBase::GetInstance();
+  ///////////////////////////////////////////////////////////////////////////////////////////////////
+  //
+  // Setting up which output to give
+  //
+  TString option="libAliHLTUtil.so libAliHLTRCU.so libAliHLTTPC.so libAliHLTITS.so libAliHLTGlobal.so loglevel=0x7c chains=";
+
+
+  ///////////////////////////////////////////////////////////////////////////////////////////////////
+  //
+  // define the analysis chain to be run
+  //
+  int minddl=0;          //min ddl number for SPD
+  int maxddl=19;         //max ddl number for SPD
+  int spec=0x1;          //spec for ddl's
+  int ddlno=0;
+  TString cfout="";
+
+  for(ddlno=minddl;ddlno<=maxddl;ddlno++){  
+    TString arg, publisher, cf;
+    
+    arg.Form("-minid %d -datatype 'DDL_RAW ' 'ISPD ' -dataspec 0x%08x -verbose",ddlno, spec);
+    publisher.Form("DP_%d", ddlno);
+    new AliHLTConfiguration(publisher.Data(), "AliRawReaderPublisher", NULL , arg.Data());
+    
+    cf.Form("CF_%d",ddlno);
+    new AliHLTConfiguration(cf.Data(), "ITSClusterFinderSPD", publisher.Data(), "");
+    
+    if (cfout.Length()>0) cfout+=" ";
+    cfout+=cf;
+    
+    spec=spec<<1;
+  }
+
+  
+  AliHLTConfiguration itsvtx("itsvtx", "ITSVertexerSPD"   , cfout.Data(), "");
+
+  AliHLTConfiguration rootFileWriterClusters("historootfile", "ROOTFileWriter", "itsvtx" , "-datafile ITSHistograms -concatenate-events -overwrite");
+
+
+  AliHLTConfiguration globalConverter("globalConverter", "GlobalEsdConverter"   , "itsvtx", "");
+  AliHLTConfiguration sink("esdfile", "EsdCollector"   , "globalConverter", "");
+  option+="historootfile,esdfile ";
+
+  ///////////////////////////////////////////////////////////////////////////////////////////////////
+  //
+  // Init and run the reconstruction
+  // All but HLT reconstructio is switched off
+  //
+  AliReconstruction rec;
+  rec.SetInput(input);
+  rec.SetRunVertexFinder(kFALSE);
+  rec.SetRunReconstruction("HLT");
+  rec.SetLoadAlignFromCDB(0);
+  rec.SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+  rec.SetSpecificStorage("GRP/GRP/Data",Form("local://%s",gSystem->pwd())); 
+  rec.SetQARefDefaultStorage("local://$ALICE_ROOT/QAref") ;
+  rec.SetRunQA(":");
+  rec.SetOption("HLT", option);
+  //rec.SetEventRange(0,0);
+rec.Run();
+}
index 2c259bf..39b7a9f 100644 (file)
@@ -143,6 +143,7 @@ void AliHLTGlobalEsdConverterComponent::GetInputDataTypes(AliHLTComponentDataTyp
   list.push_back(kAliHLTDataTypeTrackMC);
   list.push_back(kAliHLTDataTypeCaloCluster);
   list.push_back(kAliHLTDataTypedEdx );
+  list.push_back(kAliHLTDataTypeESDVertex );
 }
 
 AliHLTComponentDataType AliHLTGlobalEsdConverterComponent::GetOutputDataType()
@@ -364,6 +365,14 @@ int AliHLTGlobalEsdConverterComponent::ProcessBlocks(TTree* pTree, AliESDEvent*
     }
   }
 
+
+  // Get ITS SPD vertex
+
+  for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDVertex|kAliHLTDataOriginITS); iter != NULL; iter = GetNextInputObject() ) {
+    AliESDVertex *vtx = dynamic_cast<AliESDVertex*>(const_cast<TObject*>( iter ) );
+    pESD->SetPrimaryVertexSPD( vtx );
+  }
+
   // now update ESD tracks with the ITS info
   for (const AliHLTComponentBlockData* pBlock=GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginITS);
        pBlock!=NULL; pBlock=GetNextInputBlock()) {
index 74a85d5..d9c76c0 100644 (file)
@@ -15,7 +15,8 @@ CLASS_HDRS:=          AliHLTITStrack.h \
                tracking/AliHLTITSTrackerComponent.h \
                tracking/AliHLTITSDetector.h \
                tracking/AliHLTITSLayer.h \
-               tracking/AliHLTITSTrack.h 
+               tracking/AliHLTITSTrack.h \
+               AliHLTITSVertexerSPDComponent.h 
 
 MODULE_SRCS=   $(CLASS_HDRS:.h=.cxx)