[u/mrichter/AliRoot.git] / HLT / global / physics / AliHLTV0HistoComponent.cxx

//**************************************************************************
//* This file is property of and copyright by the ALICE HLT Project        * 
//* ALICE Experiment at CERN, All rights reserved.                         *
//*                                                                        *
//* Primary Authors: S.Gorbunov <sergey.gorbunov@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   AliHLTV0HistoComponent.cxx
    @author Sergey Gorbunov
    @brief  Component for ploting charge in clusters
*/

#if __GNUC__>= 3
using namespace std;
#endif

#include "AliHLTV0HistoComponent.h"
#include "AliCDBEntry.h"
#include "AliCDBManager.h"
#include <TFile.h>
#include <TString.h>
#include "TObjString.h"
#include "TObjArray.h"
#include "AliKFParticle.h"
#include "AliKFVertex.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TSystem.h"
#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliESDv0.h"
#include "AliHLTMessage.h"
#include "TTimeStamp.h"

//#include "AliHLTTPC.h"
//#include <stdlib.h>
//#include <cerrno>

/** ROOT macro for the implementation of ROOT specific class methods */
ClassImp(AliHLTV0HistoComponent)

AliHLTV0HistoComponent::AliHLTV0HistoComponent() :
  fUID(0),
  fGamma(0),
  fKShort(0),
  fLambda(0),
  fPi0(0),
  fAP(0),
  fGammaXY(0),
  fNEvents(0),
  fNGammas(0),
  fNKShorts(0),
  fNLambdas(0),
  fNPi0s(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

}

AliHLTV0HistoComponent::~AliHLTV0HistoComponent()
{
  // see header file for class documentation
}

// Public functions to implement AliHLTComponent's interface.
// These functions are required for the registration process

const char* AliHLTV0HistoComponent::GetComponentID()
{
  // see header file for class documentation
  
  return "V0Histo";
}

void AliHLTV0HistoComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list)
{
  // see header file for class documentation
  list.clear();
  list.push_back( kAliHLTDataTypeESDObject|kAliHLTDataOriginOut );
}

AliHLTComponentDataType AliHLTV0HistoComponent::GetOutputDataType()
{
  // see header file for class documentation
  return kAliHLTDataTypeHistogram  | kAliHLTDataOriginOut;
}

void AliHLTV0HistoComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
{
  // see header file for class documentation
  // XXX TODO: Find more realistic values.
  constBase = 80000;
  inputMultiplier = 0;
}

AliHLTComponent* AliHLTV0HistoComponent::Spawn()
{
  // see header file for class documentation
  return new AliHLTV0HistoComponent;
}

int AliHLTV0HistoComponent::DoInit( int argc, const char** argv ) {
  // init

  fUID = 0;

  fGamma = new TH1F("hGamma","HLT:  #gamma inv mass",50,-.06,.2); 
  fGamma->SetFillColor(kGreen);
  fGamma->SetStats(0);
  fKShort = new TH1F("hKShort","HLT:  K_{s}^{0} inv mass",80,0.4,.6); 
  fKShort->SetFillColor(kGreen);
  fKShort->SetStats(0);
  fLambda = new TH1F("hLambda","HLT:  #Lambda^{0} inv mass",50,1.0,1.36); 
  fLambda->SetFillColor(kGreen);
  fLambda->SetStats(0);

  fPi0 = new TH1F("hPi0","HLT:  #Pi^{0} inv mass",50,0.0,0.38); 
  fPi0->SetFillColor(kGreen);
  fPi0->SetStats(0);
  
  fAP = new TH2F("hAP","HLT:  Armenteros-Podolanski plot",60,-1.,1.,60,-0.02,0.3);
  fAP->SetMarkerStyle(8);
  fAP->SetMarkerSize(0.4);
  fAP->SetYTitle("p_{t}[GeV]");
  fAP->SetXTitle("(p^{+}_{L}-p^{-}_{L})/(p^{+}_{L}+p^{-}_{L})");
  fAP->SetStats(0);
  fAP->SetOption("CONT4 Z");

  fGammaXY = new TH2F("hGammaXY","HLT:  #gamma conversions",100,-100.,100.,100,-100.,100.);
  fGammaXY->SetMarkerStyle(8);
  fGammaXY->SetMarkerSize(0.4);
  fGammaXY->SetYTitle("X[cm]");
  fGammaXY->SetXTitle("Y[cm]");
  fGammaXY->SetStats(0);
  fGammaXY->SetOption("COLZ");

  fNEvents =0;
  fNGammas = 0;
  fNKShorts = 0;
  fNLambdas = 0;
  fNPi0s = 0;

  // cuts: 
  // [0] == 0   --- N clusters on each daughter track
  // [1] == 2.5 --- (daughter-primVtx)/sigma >= cut
  // [2] == 3.5 --- (v0 - primVtx)/sigma <= cut
  // [3] == 3.0 --- (decay_length)/sigma >= cut
  // [4] == 0.0 --- (decay_length)[cm]   >= cut
  // [5] == 300.0 --- (v0 radius)[cm]    <= cut
  // [6] == 3.5  --- (v0 mass - true value)/sigma <= cut (for identification)
  // [7] == 0.05 --- (v0 mass - true value)       <= cut (for identification)

  fGammaCuts[0] = 0;
  fGammaCuts[1] = 2.5;
  fGammaCuts[2] = 3.5;
  fGammaCuts[3] = 3.0;
  fGammaCuts[4] = 0.0;
  fGammaCuts[5] = 300.0;
  fGammaCuts[6] = 3.5;
  fGammaCuts[7] = 0.05;

  fAPCuts[0] = 60;
  fAPCuts[1] = 2.5;
  fAPCuts[2] = 3.5;
  fAPCuts[3] = 3.0;
  fAPCuts[4] = 0.0;
  fAPCuts[5] = 50.0;
  fAPCuts[6] = 4.0;
  fAPCuts[7] = 0.05;

  fKsCuts[0] = 60;
  fKsCuts[1] = 2.5;
  fKsCuts[2] = 3.5;
  fKsCuts[3] = 3.0;
  fKsCuts[4] = 1.5;
  fKsCuts[5] = 50.0;
  fKsCuts[6] = 4.0;
  fKsCuts[7] = 0.015;

  fLambdaCuts[0] = 60;   // [0] == 60  --- N clusters on each daughter track
  fLambdaCuts[1] = 3.0;  // [1] == 3.0 --- (daughter-primVtx)/sigma >= cut
  fLambdaCuts[2] = 3.0;  // [2] == 3.0 --- (v0 - primVtx)/sigma <= cut
  fLambdaCuts[3] = 3.5;  // [3] == 3.5 --- (decay_length)/sigma >= cut
  fLambdaCuts[4] = 4.0;  // [4] == 0.0 --- (decay_length)[cm]   >= cut
  fLambdaCuts[5] = 50.0; // [5] == 300.0 --- (v0 radius)[cm]    <= cut
  fLambdaCuts[6] = 4.0;  // [6] == 3.5  --- (v0 mass - true value)/sigma <= cut (for identification)
  fLambdaCuts[7] = 0.03; // [7] == 0.05 --- (v0 mass - true value)       <= cut (for identification)

  int iResult=0;
  TString configuration="";
  TString argument="";
  for (int i=0; i<argc && iResult>=0; i++) {
    argument=argv[i];
    if (!configuration.IsNull()) configuration+=" ";
    configuration+=argument;
  }
  
  if (!configuration.IsNull()) {
    iResult=Configure(configuration.Data());
  }  

  return iResult; 
}
  
int AliHLTV0HistoComponent::DoDeinit()
{
  // see header file for class documentation
  delete fGamma;
  delete fKShort;
  delete fLambda;
  delete fAP;
  delete fGammaXY;
  fUID = 0;
  return 0;
}

int AliHLTV0HistoComponent::DoEvent(const AliHLTComponentEventData& evtData, AliHLTComponentTriggerData& /*trigData*/)
{
  
  if ( GetFirstInputBlock( kAliHLTDataTypeSOR ) || GetFirstInputBlock( kAliHLTDataTypeEOR ) )
    return 0;

  if( fUID == 0 ){
    TTimeStamp t;
    fUID = ( gSystem->GetPid() + t.GetNanoSec())*10 + evtData.fEventID;
  }

  fNEvents++;

  for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {

    AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
    event->GetStdContent();
    Int_t nV0 = event->GetNumberOfV0s();
    AliKFParticle::SetField( event->GetMagneticField() );

    const double kKsMass = 0.49767;
    const double kLambdaMass = 1.11568;
    const double kPi0Mass = 0.13498;

    std::vector<AliKFParticle> vGammas;

    for (Int_t iv=0; iv<nV0; iv++) {
	
      AliESDtrack *t1=event->GetTrack( event->GetV0(iv)->GetNindex());
      AliESDtrack *t2=event->GetTrack( event->GetV0(iv)->GetPindex());      

      AliKFParticle kf1( *t1, 11 );
      AliKFParticle kf2( *t2, 11 );

      AliKFVertex primVtx( *event->GetPrimaryVertexTracks() );
      double dev1 = kf1.GetDeviationFromVertex( primVtx );
      double dev2 = kf2.GetDeviationFromVertex( primVtx );
      
      AliKFParticle v0( kf1, kf2 );
      double devPrim = v0.GetDeviationFromVertex( primVtx );
      primVtx+=v0;
      v0.SetProductionVertex( primVtx );

      Double_t length, sigmaLength;
      if( v0.GetDecayLength( length, sigmaLength ) ) continue;

      double dx = v0.GetX()-primVtx.GetX();
      double dy = v0.GetY()-primVtx.GetY();
      double r = sqrt(dx*dx + dy*dy);

      // AP plot

      double pt=0, ap=0;
      {
	AliKFParticle kf01 = kf1, kf02 = kf2;
	kf01.SetProductionVertex(v0);
	kf02.SetProductionVertex(v0);
	kf01.TransportToProductionVertex();
	kf02.TransportToProductionVertex();      
	double px1=kf01.GetPx(), py1=kf01.GetPy(), pz1=kf01.GetPz();
	double px2=kf02.GetPx(), py2=kf02.GetPy(), pz2=kf02.GetPz();
	double px = px1+px2, py = py1+py2, pz = pz1+pz2;
	double p = sqrt(px*px+py*py+pz*pz);
	double l1 = (px*px1 + py*py1 + pz*pz1)/p;
	double l2 = (px*px2 + py*py2 + pz*pz2)/p;
	pt = sqrt(px1*px1+py1*py1+pz1*pz1 - l1*l1);
	ap = (l2-l1)/(l1+l2);
      }
      
      if( 
	 t1->GetTPCNcls()>=fAPCuts[0]
	 && t2->GetTPCNcls()>=fAPCuts[0]
	 && dev1>=fAPCuts[1]
	 && dev2>=fAPCuts[1]
	 && devPrim <= fAPCuts[2]
	 && length >= fAPCuts[3]*sigmaLength
	 && length >= fAPCuts[4]
	 && r <= fAPCuts[5]
	 ){	
	if( fAP ) fAP->Fill( ap, pt );
      } 

      // Gamma finder

      bool isGamma = 0;
      
      if( 
	 t1->GetTPCNcls()>=fGammaCuts[0]
	 && t2->GetTPCNcls()>=fGammaCuts[0]
	 && dev1>=fGammaCuts[1]
	 && dev2>=fGammaCuts[1]
	 && devPrim <= fGammaCuts[2]
	 && length >= fGammaCuts[3]*sigmaLength
	 && length >= fGammaCuts[4]
	 && r <= fGammaCuts[5]
	 ){
	double mass, error;       
	v0.GetMass(mass,error);	
	if( fGamma ) fGamma->Fill( mass );

	if( TMath::Abs(mass)<=fGammaCuts[6]*error || TMath::Abs(mass)<=fGammaCuts[7] ){	  
	  AliKFParticle gamma = v0;
	  gamma.SetMassConstraint(0);
	  if( fGammaXY
	      &&  t1->GetTPCNcls()>=60
	      && t2->GetTPCNcls()>=60
	      ) fGammaXY->Fill(gamma.GetX(), gamma.GetY());
	  isGamma = 1;
	  fNGammas++;
	  vGammas.push_back( gamma );
	}	     
      }
      
      if( isGamma ) continue;


      // KShort finder
      
      bool isKs = 0;
      
      if( 
	 t1->GetTPCNcls()>=fKsCuts[0]
	 && t2->GetTPCNcls()>=fKsCuts[0]
	 && dev1>=fKsCuts[1]
	 && dev2>=fKsCuts[1]
	 && devPrim <= fKsCuts[2]
	 && length >= fKsCuts[3]*sigmaLength
	 && length >= fKsCuts[4]
	 && r <= fKsCuts[5]
	 ){	
    
	AliKFParticle piN( *t1, 211 );	
	AliKFParticle piP( *t2, 211 );	
	
	AliKFParticle Ks( piN, piP );
	Ks.SetProductionVertex( primVtx );
	
	double mass, error;
	Ks.GetMass( mass, error);
	if( fKShort ) fKShort->Fill( mass );	
	if( TMath::Abs( mass - kKsMass )<=fKsCuts[6]*error || TMath::Abs( mass - kKsMass )<=fKsCuts[7] ){  
	  isKs = 1;
	  fNKShorts++;
	}
      }
      
      if( isKs ) continue;
      
      // Lambda finder 
     
      if( 
	 t1->GetTPCNcls()>=fLambdaCuts[0]
	 && t2->GetTPCNcls()>=fLambdaCuts[0]
	 && dev1>=fLambdaCuts[1]
	 && dev2>=fLambdaCuts[1]
	 && devPrim <= fLambdaCuts[2]
	 && length >= fLambdaCuts[3]*sigmaLength
	 && length >= fLambdaCuts[4]
	 && r <= fLambdaCuts[5]
	 && TMath::Abs( ap )>.4
	 ){
	
	AliKFParticle kP, kpi;
	if( ap<0 ){ 
	  kP = AliKFParticle( *t2, 2212 );
	  kpi = AliKFParticle( *t1, 211 );
	} else {
	  kP = AliKFParticle( *t1, 2212 );
	  kpi = AliKFParticle( *t2, 211 );
	}

	AliKFParticle lambda = AliKFParticle( kP, kpi );
	lambda.SetProductionVertex( primVtx );	
	double mass, error;
	lambda.GetMass( mass, error);
	if( fLambda ) fLambda->Fill( mass );
	if( TMath::Abs( mass - kLambdaMass )<=fLambdaCuts[6]*error || TMath::Abs( mass - kLambdaMass )<=fLambdaCuts[7] ){
	  fNLambdas++;
	}
      }

    }// V0's


    // Pi0 finder 

    for(UInt_t g1=0;g1<vGammas.size();g1++){
      for(UInt_t g2=g1+1;g2<vGammas.size();g2++){
	AliKFParticle pi0(vGammas.at(g1),vGammas.at(g2));
	double mass, error;
	pi0.GetMass(mass,error);
	fPi0->Fill(mass);
	if( TMath::Abs( mass - kPi0Mass )<=0.03 ){
	  fNPi0s++;
	}
      }
    }

  
    if( fGamma ) PushBack( (TObject*) fGamma, kAliHLTDataTypeHistogram,fUID);
    
    if( fKShort ) PushBack( (TObject*) fKShort, kAliHLTDataTypeHistogram,fUID);
    
    if( fLambda ) PushBack( (TObject*) fLambda, kAliHLTDataTypeHistogram, fUID);
 
    if( fPi0 ) PushBack( (TObject*) fPi0, kAliHLTDataTypeHistogram, fUID);
    
    if( fAP ) PushBack( (TObject*) fAP, kAliHLTDataTypeHistogram, fUID);    

    if( fGammaXY ) PushBack( (TObject*) fGammaXY, kAliHLTDataTypeHistogram, fUID);
  }  
  if( fNPi0s>0 ){
    HLTInfo("Found %d Gammas, %d KShorts, %d Lambdas and %d Pi0's in %d events", fNGammas, fNKShorts, fNLambdas, fNPi0s, fNEvents );    
  }
  else HLTInfo("Found %d Gammas, %d KShorts and %d Lambdas in %d events", fNGammas, fNKShorts, fNLambdas, fNEvents );
  
  return 0;
}

int AliHLTV0HistoComponent::Configure(const char* arguments)
{
  
  int iResult=0;
  if (!arguments) return iResult;
  
  TString allArgs=arguments;
  TString argument;
  
  TObjArray* pTokens=allArgs.Tokenize(" ");
  int bMissingParam=0;

  if (pTokens) {
    for (int i=0; i<pTokens->GetEntries() && iResult>=0; i++) {
      argument=((TObjString*)pTokens->At(i))->GetString();
      if (argument.IsNull()) continue;

      if (argument.CompareTo("-cutsGamma")==0) {
	TString spar = "";	
	for( int j=0; j<8; j++ ){
	  if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
	  spar+=" ";
	  spar+=((TObjString*)pTokens->At(i))->GetString();
	  fGammaCuts[j]=((TObjString*)pTokens->At(i))->GetString().Atof();
	}
	if( !bMissingParam ){
	  HLTInfo("Gamma cuts are set to: %s", spar.Data());
	  continue;
	}
      } else if (argument.CompareTo("-cutsAP")==0) {
	TString spar = "";	
	for( int j=0; j<8; j++ ){
	  if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
	  spar+=" ";
	  spar+=((TObjString*)pTokens->At(i))->GetString();
	  fAPCuts[j]=((TObjString*)pTokens->At(i))->GetString().Atof();
	}
	if( !bMissingParam ){
	  HLTInfo("AP cuts are set to: %s", spar.Data());
	  continue;
	}
      }
      else if (argument.CompareTo("-cutsKs")==0) {
	TString spar = "";	
	for( int j=0; j<8; j++ ){
	  if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
	  spar+=" ";
	  spar+=((TObjString*)pTokens->At(i))->GetString();
	  fKsCuts[j]=((TObjString*)pTokens->At(i))->GetString().Atof();
	}
	if( !bMissingParam ){
	  HLTInfo("KShort cuts are set to: %s", spar.Data());
	  continue;
	}
      } else if (argument.CompareTo("-cutsLambda")==0) {
	TString spar = "";	
	for( int j=0; j<8; j++ ){
	  if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
	  spar+=" ";
	  spar+=((TObjString*)pTokens->At(i))->GetString();
	  fLambdaCuts[j]=((TObjString*)pTokens->At(i))->GetString().Atof();
	}
	if( !bMissingParam ){
	  HLTInfo("Lampda cuts are set to: %s", spar.Data());
	  continue;
	}
      }else {
	HLTError("unknown argument %s", argument.Data());
	iResult=-EINVAL;
	break;
      }
    }
    delete pTokens;
  }
  
  return iResult;
}

int AliHLTV0HistoComponent::Reconfigure(const char* cdbEntry, const char* chainId)
{
  // see header file for class documentation
  int iResult=0;
  const char* path="HLT/ConfigTPC/KryptonHistoComponent";
  const char* defaultNotify="";
  if (cdbEntry) {
    path=cdbEntry;
    defaultNotify=" (default)";
  }
  if (path) {
    HLTInfo("reconfigure from entry %s%s, chain id %s", path, defaultNotify,(chainId!=NULL && chainId[0]!=0)?chainId:"<none>");
    AliCDBEntry *pEntry = AliCDBManager::Instance()->Get(path/*,GetRunNo()*/);
    if (pEntry) {
      TObjString* pString=dynamic_cast<TObjString*>(pEntry->GetObject());
      if (pString) {
	HLTInfo("received configuration object string: \'%s\'", pString->GetString().Data());
	iResult=Configure(pString->GetString().Data());
      } else {
	HLTError("configuration object \"%s\" has wrong type, required TObjString", path);
      }
    } else {
      HLTError("can not fetch object \"%s\" from CDB", path);
    }
  }

  return iResult;
}
Commit	Line	Data
8125805f	1	//**************************************************************************
	2	//* This file is property of and copyright by the ALICE HLT Project *
	3	//* ALICE Experiment at CERN, All rights reserved. *
	4	//* *
08458675	5	//* Primary Authors: S.Gorbunov <sergey.gorbunov@kip.uni-heidelberg.de> *
8125805f	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	/** @file AliHLTV0HistoComponent.cxx
08458675	18	@author Sergey Gorbunov
8125805f	19	@brief Component for ploting charge in clusters
	20	*/
	21
	22	#if __GNUC__>= 3
	23	using namespace std;
	24	#endif
	25
	26	#include "AliHLTV0HistoComponent.h"
	27	#include "AliCDBEntry.h"
	28	#include "AliCDBManager.h"
	29	#include <TFile.h>
	30	#include <TString.h>
	31	#include "TObjString.h"
	32	#include "TObjArray.h"
	33	#include "AliKFParticle.h"
	34	#include "AliKFVertex.h"
	35	#include "TH1F.h"
	36	#include "TH2F.h"
c430b9f9	37	#include "TSystem.h"
8125805f	38	#include "AliESDEvent.h"
	39	#include "AliESDtrack.h"
	40	#include "AliESDv0.h"
	41	#include "AliHLTMessage.h"
e419c1ae	42	#include "TTimeStamp.h"
8125805f	43
	44	//#include "AliHLTTPC.h"
	45	//#include <stdlib.h>
	46	//#include <cerrno>
	47
	48	/** ROOT macro for the implementation of ROOT specific class methods */
	49	ClassImp(AliHLTV0HistoComponent)
	50
766aafea	51	AliHLTV0HistoComponent::AliHLTV0HistoComponent() :
e419c1ae	52	fUID(0),
01ce7f55	53	fGamma(0),
8125805f	54	fKShort(0),
8125805f	55	fLambda(0),
5c9f2410	56	fPi0(0),
8125805f	57	fAP(0),
01ce7f55	58	fGammaXY(0),
4d5ee3db	59	fNEvents(0),
01ce7f55	60	fNGammas(0),
8125805f	61	fNKShorts(0),
5c9f2410	62	fNLambdas(0),
5c9f2410	63	fNPi0s(0)
8125805f	64	{
	65	// see header file for class documentation
	66	// or
	67	// refer to README to build package
	68	// or
	69	// visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
	70
	71	}
	72
	73	AliHLTV0HistoComponent::~AliHLTV0HistoComponent()
	74	{
	75	// see header file for class documentation
	76	}
	77
	78	// Public functions to implement AliHLTComponent's interface.
	79	// These functions are required for the registration process
	80
	81	const char* AliHLTV0HistoComponent::GetComponentID()
	82	{
	83	// see header file for class documentation
	84
	85	return "V0Histo";
	86	}
	87
	88	void AliHLTV0HistoComponent::GetInputDataTypes(AliHLTComponentDataTypeList& list)
	89	{
	90	// see header file for class documentation
	91	list.clear();
	92	list.push_back( kAliHLTDataTypeESDObject\|kAliHLTDataOriginOut );
	93	}
	94
	95	AliHLTComponentDataType AliHLTV0HistoComponent::GetOutputDataType()
	96	{
	97	// see header file for class documentation
08458675	98	return kAliHLTDataTypeHistogram \| kAliHLTDataOriginOut;
8125805f	99	}
	100
	101	void AliHLTV0HistoComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier )
	102	{
	103	// see header file for class documentation
	104	// XXX TODO: Find more realistic values.
	105	constBase = 80000;
2bf7a9f5	106	inputMultiplier = 0;
8125805f	107	}
	108
	109	AliHLTComponent* AliHLTV0HistoComponent::Spawn()
	110	{
	111	// see header file for class documentation
	112	return new AliHLTV0HistoComponent;
	113	}
	114
766aafea	115	int AliHLTV0HistoComponent::DoInit( int argc, const char** argv ) {
8125805f	116	// init
e419c1ae	117
	118	fUID = 0;
	119
d0e01ea8	120	fGamma = new TH1F("hGamma","HLT: #gamma inv mass",50,-.06,.2);
01ce7f55	121	fGamma->SetFillColor(kGreen);
	122	fGamma->SetStats(0);
	123	fKShort = new TH1F("hKShort","HLT: K_{s}^{0} inv mass",80,0.4,.6);
	124	fKShort->SetFillColor(kGreen);
8125805f	125	fKShort->SetStats(0);
d0e01ea8	126	fLambda = new TH1F("hLambda","HLT: #Lambda^{0} inv mass",50,1.0,1.36);
01ce7f55	127	fLambda->SetFillColor(kGreen);
8125805f	128	fLambda->SetStats(0);
5c9f2410	129
d0e01ea8	130	fPi0 = new TH1F("hPi0","HLT: #Pi^{0} inv mass",50,0.0,0.38);
5c9f2410	131	fPi0->SetFillColor(kGreen);
	132	fPi0->SetStats(0);
	133
2bf7a9f5	134	fAP = new TH2F("hAP","HLT: Armenteros-Podolanski plot",60,-1.,1.,60,-0.02,0.3);
8125805f	135	fAP->SetMarkerStyle(8);
	136	fAP->SetMarkerSize(0.4);
	137	fAP->SetYTitle("p_{t}[GeV]");
	138	fAP->SetXTitle("(p^{+}_{L}-p^{-}_{L})/(p^{+}_{L}+p^{-}_{L})");
	139	fAP->SetStats(0);
d3792c06	140	fAP->SetOption("CONT4 Z");
01ce7f55	141
	142	fGammaXY = new TH2F("hGammaXY","HLT: #gamma conversions",100,-100.,100.,100,-100.,100.);
	143	fGammaXY->SetMarkerStyle(8);
	144	fGammaXY->SetMarkerSize(0.4);
	145	fGammaXY->SetYTitle("X[cm]");
	146	fGammaXY->SetXTitle("Y[cm]");
	147	fGammaXY->SetStats(0);
	148	fGammaXY->SetOption("COLZ");
	149
4d5ee3db	150	fNEvents =0;
01ce7f55	151	fNGammas = 0;
8125805f	152	fNKShorts = 0;
8125805f	153	fNLambdas = 0;
5c9f2410	154	fNPi0s = 0;
8125805f	155
608cfbda	156	// cuts:
	157	// [0] == 0 --- N clusters on each daughter track
	158	// [1] == 2.5 --- (daughter-primVtx)/sigma >= cut
	159	// [2] == 3.5 --- (v0 - primVtx)/sigma <= cut
	160	// [3] == 3.0 --- (decay_length)/sigma >= cut
	161	// [4] == 0.0 --- (decay_length)[cm] >= cut
	162	// [5] == 300.0 --- (v0 radius)[cm] <= cut
	163	// [6] == 3.5 --- (v0 mass - true value)/sigma <= cut (for identification)
	164	// [7] == 0.05 --- (v0 mass - true value) <= cut (for identification)
	165
	166	fGammaCuts[0] = 0;
	167	fGammaCuts[1] = 2.5;
	168	fGammaCuts[2] = 3.5;
	169	fGammaCuts[3] = 3.0;
	170	fGammaCuts[4] = 0.0;
	171	fGammaCuts[5] = 300.0;
d3792c06	172	fGammaCuts[6] = 3.5;
608cfbda	173	fGammaCuts[7] = 0.05;
608cfbda	174
d3792c06	175	fAPCuts[0] = 60;
608cfbda	176	fAPCuts[1] = 2.5;
	177	fAPCuts[2] = 3.5;
	178	fAPCuts[3] = 3.0;
	179	fAPCuts[4] = 0.0;
	180	fAPCuts[5] = 50.0;
	181	fAPCuts[6] = 4.0;
	182	fAPCuts[7] = 0.05;
	183
	184	fKsCuts[0] = 60;
	185	fKsCuts[1] = 2.5;
	186	fKsCuts[2] = 3.5;
	187	fKsCuts[3] = 3.0;
	188	fKsCuts[4] = 1.5;
	189	fKsCuts[5] = 50.0;
73a33d2e	190	fKsCuts[6] = 4.0;
d3792c06	191	fKsCuts[7] = 0.015;
608cfbda	192
d0e01ea8	193	fLambdaCuts[0] = 60; // [0] == 60 --- N clusters on each daughter track
	194	fLambdaCuts[1] = 3.0; // [1] == 3.0 --- (daughter-primVtx)/sigma >= cut
	195	fLambdaCuts[2] = 3.0; // [2] == 3.0 --- (v0 - primVtx)/sigma <= cut
	196	fLambdaCuts[3] = 3.5; // [3] == 3.5 --- (decay_length)/sigma >= cut
	197	fLambdaCuts[4] = 4.0; // [4] == 0.0 --- (decay_length)[cm] >= cut
	198	fLambdaCuts[5] = 50.0; // [5] == 300.0 --- (v0 radius)[cm] <= cut
	199	fLambdaCuts[6] = 4.0; // [6] == 3.5 --- (v0 mass - true value)/sigma <= cut (for identification)
	200	fLambdaCuts[7] = 0.03; // [7] == 0.05 --- (v0 mass - true value) <= cut (for identification)
608cfbda	201
8125805f	202	int iResult=0;
	203	TString configuration="";
	204	TString argument="";
	205	for (int i=0; i<argc && iResult>=0; i++) {
	206	argument=argv[i];
	207	if (!configuration.IsNull()) configuration+=" ";
	208	configuration+=argument;
	209	}
	210
	211	if (!configuration.IsNull()) {
	212	iResult=Configure(configuration.Data());
	213	}
	214
	215	return iResult;
	216	}
	217
	218	int AliHLTV0HistoComponent::DoDeinit()
	219	{
	220	// see header file for class documentation
01ce7f55	221	delete fGamma;
8125805f	222	delete fKShort;
	223	delete fLambda;
	224	delete fAP;
01ce7f55	225	delete fGammaXY;
e419c1ae	226	fUID = 0;
8125805f	227	return 0;
	228	}
	229
e419c1ae	230	int AliHLTV0HistoComponent::DoEvent(const AliHLTComponentEventData& evtData, AliHLTComponentTriggerData& /trigData/)
8125805f	231	{
	232
	233	if ( GetFirstInputBlock( kAliHLTDataTypeSOR ) \|\| GetFirstInputBlock( kAliHLTDataTypeEOR ) )
	234	return 0;
	235
e419c1ae	236	if( fUID == 0 ){
	237	TTimeStamp t;
	238	fUID = ( gSystem->GetPid() + t.GetNanoSec())*10 + evtData.fEventID;
	239	}
	240
4d5ee3db	241	fNEvents++;
4d5ee3db	242
8125805f	243	for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {
	244
	245	AliESDEvent event = dynamic_cast<AliESDEvent>(const_cast<TObject*>( iter ) );
	246	event->GetStdContent();
	247	Int_t nV0 = event->GetNumberOfV0s();
	248	AliKFParticle::SetField( event->GetMagneticField() );
01ce7f55	249
	250	const double kKsMass = 0.49767;
	251	const double kLambdaMass = 1.11568;
5c9f2410	252	const double kPi0Mass = 0.13498;
	253
	254	std::vector<AliKFParticle> vGammas;
01ce7f55	255
8125805f	256	for (Int_t iv=0; iv<nV0; iv++) {
8125805f	257
01ce7f55	258	AliESDtrack *t1=event->GetTrack( event->GetV0(iv)->GetNindex());
01ce7f55	259	AliESDtrack *t2=event->GetTrack( event->GetV0(iv)->GetPindex());
8125805f	260
c3c8af48	261	AliKFParticle kf1( *t1, 11 );
c3c8af48	262	AliKFParticle kf2( *t2, 11 );
8125805f	263
	264	AliKFVertex primVtx( *event->GetPrimaryVertexTracks() );
	265	double dev1 = kf1.GetDeviationFromVertex( primVtx );
	266	double dev2 = kf2.GetDeviationFromVertex( primVtx );
	267
01ce7f55	268	AliKFParticle v0( kf1, kf2 );
608cfbda	269	double devPrim = v0.GetDeviationFromVertex( primVtx );
01ce7f55	270	primVtx+=v0;
01ce7f55	271	v0.SetProductionVertex( primVtx );
8125805f	272
608cfbda	273	Double_t length, sigmaLength;
	274	if( v0.GetDecayLength( length, sigmaLength ) ) continue;
	275
	276	double dx = v0.GetX()-primVtx.GetX();
	277	double dy = v0.GetY()-primVtx.GetY();
	278	double r = sqrt(dxdx + dydy);
	279
2bf7a9f5	280	// AP plot
	281
	282	double pt=0, ap=0;
	283	{
	284	AliKFParticle kf01 = kf1, kf02 = kf2;
	285	kf01.SetProductionVertex(v0);
	286	kf02.SetProductionVertex(v0);
	287	kf01.TransportToProductionVertex();
	288	kf02.TransportToProductionVertex();
	289	double px1=kf01.GetPx(), py1=kf01.GetPy(), pz1=kf01.GetPz();
	290	double px2=kf02.GetPx(), py2=kf02.GetPy(), pz2=kf02.GetPz();
	291	double px = px1+px2, py = py1+py2, pz = pz1+pz2;
	292	double p = sqrt(pxpx+pypy+pz*pz);
	293	double l1 = (pxpx1 + pypy1 + pz*pz1)/p;
	294	double l2 = (pxpx2 + pypy2 + pz*pz2)/p;
	295	pt = sqrt(px1px1+py1py1+pz1pz1 - l1l1);
3d292756	296	ap = (l2-l1)/(l1+l2);
2bf7a9f5	297	}
	298
	299	if(
	300	t1->GetTPCNcls()>=fAPCuts[0]
	301	&& t2->GetTPCNcls()>=fAPCuts[0]
	302	&& dev1>=fAPCuts[1]
	303	&& dev2>=fAPCuts[1]
	304	&& devPrim <= fAPCuts[2]
	305	&& length >= fAPCuts[3]*sigmaLength
	306	&& length >= fAPCuts[4]
	307	&& r <= fAPCuts[5]
	308	){
	309	if( fAP ) fAP->Fill( ap, pt );
	310	}
	311
01ce7f55	312	// Gamma finder
608cfbda	313
01ce7f55	314	bool isGamma = 0;
608cfbda	315
	316	if(
	317	t1->GetTPCNcls()>=fGammaCuts[0]
	318	&& t2->GetTPCNcls()>=fGammaCuts[0]
	319	&& dev1>=fGammaCuts[1]
	320	&& dev2>=fGammaCuts[1]
	321	&& devPrim <= fGammaCuts[2]
	322	&& length >= fGammaCuts[3]*sigmaLength
	323	&& length >= fGammaCuts[4]
	324	&& r <= fGammaCuts[5]
	325	){
01ce7f55	326	double mass, error;
5c9f2410	327	v0.GetMass(mass,error);
01ce7f55	328	if( fGamma ) fGamma->Fill( mass );
01ce7f55	329
d3792c06	330	if( TMath::Abs(mass)<=fGammaCuts[6]*error \|\| TMath::Abs(mass)<=fGammaCuts[7] ){
01ce7f55	331	AliKFParticle gamma = v0;
01ce7f55	332	gamma.SetMassConstraint(0);
d3792c06	333	if( fGammaXY
	334	&& t1->GetTPCNcls()>=60
	335	&& t2->GetTPCNcls()>=60
	336	) fGammaXY->Fill(gamma.GetX(), gamma.GetY());
01ce7f55	337	isGamma = 1;
01ce7f55	338	fNGammas++;
5c9f2410	339	vGammas.push_back( gamma );
01ce7f55	340	}
01ce7f55	341	}
8125805f	342
01ce7f55	343	if( isGamma ) continue;
01ce7f55	344
01ce7f55	345
	346	// KShort finder
	347
	348	bool isKs = 0;
	349
608cfbda	350	if(
	351	t1->GetTPCNcls()>=fKsCuts[0]
	352	&& t2->GetTPCNcls()>=fKsCuts[0]
	353	&& dev1>=fKsCuts[1]
	354	&& dev2>=fKsCuts[1]
	355	&& devPrim <= fKsCuts[2]
	356	&& length >= fKsCuts[3]*sigmaLength
	357	&& length >= fKsCuts[4]
	358	&& r <= fKsCuts[5]
	359	){
	360
c3c8af48	361	AliKFParticle piN( *t1, 211 );
c3c8af48	362	AliKFParticle piP( *t2, 211 );
01ce7f55	363
	364	AliKFParticle Ks( piN, piP );
	365	Ks.SetProductionVertex( primVtx );
	366
	367	double mass, error;
	368	Ks.GetMass( mass, error);
73a33d2e	369	if( fKShort ) fKShort->Fill( mass );
d3792c06	370	if( TMath::Abs( mass - kKsMass )<=fKsCuts[6]*error \|\| TMath::Abs( mass - kKsMass )<=fKsCuts[7] ){
01ce7f55	371	isKs = 1;
	372	fNKShorts++;
	373	}
	374	}
	375
	376	if( isKs ) continue;
608cfbda	377
8125805f	378	// Lambda finder
01ce7f55	379
608cfbda	380	if(
	381	t1->GetTPCNcls()>=fLambdaCuts[0]
	382	&& t2->GetTPCNcls()>=fLambdaCuts[0]
	383	&& dev1>=fLambdaCuts[1]
	384	&& dev2>=fLambdaCuts[1]
	385	&& devPrim <= fLambdaCuts[2]
	386	&& length >= fLambdaCuts[3]*sigmaLength
	387	&& length >= fLambdaCuts[4]
	388	&& r <= fLambdaCuts[5]
d3792c06	389	&& TMath::Abs( ap )>.4
608cfbda	390	){
01ce7f55	391
8125805f	392	AliKFParticle kP, kpi;
01ce7f55	393	if( ap<0 ){
c3c8af48	394	kP = AliKFParticle( *t2, 2212 );
c3c8af48	395	kpi = AliKFParticle( *t1, 211 );
01ce7f55	396	} else {
c3c8af48	397	kP = AliKFParticle( *t1, 2212 );
c3c8af48	398	kpi = AliKFParticle( *t2, 211 );
8125805f	399	}
01ce7f55	400
	401	AliKFParticle lambda = AliKFParticle( kP, kpi );
	402	lambda.SetProductionVertex( primVtx );
	403	double mass, error;
	404	lambda.GetMass( mass, error);
	405	if( fLambda ) fLambda->Fill( mass );
d3792c06	406	if( TMath::Abs( mass - kLambdaMass )<=fLambdaCuts[6]*error \|\| TMath::Abs( mass - kLambdaMass )<=fLambdaCuts[7] ){
01ce7f55	407	fNLambdas++;
73a33d2e	408	}
8125805f	409	}
01ce7f55	410
8125805f	411	}// V0's
5c9f2410	412
	413
	414	// Pi0 finder
	415
	416	for(UInt_t g1=0;g1<vGammas.size();g1++){
	417	for(UInt_t g2=g1+1;g2<vGammas.size();g2++){
	418	AliKFParticle pi0(vGammas.at(g1),vGammas.at(g2));
	419	double mass, error;
	420	pi0.GetMass(mass,error);
	421	fPi0->Fill(mass);
8e635044	422	if( TMath::Abs( mass - kPi0Mass )<=0.03 ){
5c9f2410	423	fNPi0s++;
	424	}
	425	}
	426	}
	427
8125805f	428
e419c1ae	429	if( fGamma ) PushBack( (TObject*) fGamma, kAliHLTDataTypeHistogram,fUID);
01ce7f55	430
e419c1ae	431	if( fKShort ) PushBack( (TObject*) fKShort, kAliHLTDataTypeHistogram,fUID);
01ce7f55	432
e419c1ae	433	if( fLambda ) PushBack( (TObject*) fLambda, kAliHLTDataTypeHistogram, fUID);
5c9f2410	434
e419c1ae	435	if( fPi0 ) PushBack( (TObject*) fPi0, kAliHLTDataTypeHistogram, fUID);
01ce7f55	436
e419c1ae	437	if( fAP ) PushBack( (TObject*) fAP, kAliHLTDataTypeHistogram, fUID);
01ce7f55	438
e419c1ae	439	if( fGammaXY ) PushBack( (TObject*) fGammaXY, kAliHLTDataTypeHistogram, fUID);
8125805f	440	}
5c9f2410	441	if( fNPi0s>0 ){
d0e01ea8	442	HLTInfo("Found %d Gammas, %d KShorts, %d Lambdas and %d Pi0's in %d events", fNGammas, fNKShorts, fNLambdas, fNPi0s, fNEvents );
5c9f2410	443	}
d0e01ea8	444	else HLTInfo("Found %d Gammas, %d KShorts and %d Lambdas in %d events", fNGammas, fNKShorts, fNLambdas, fNEvents );
8125805f	445
	446	return 0;
	447	}
	448
	449	int AliHLTV0HistoComponent::Configure(const char* arguments)
	450	{
	451
	452	int iResult=0;
	453	if (!arguments) return iResult;
	454
	455	TString allArgs=arguments;
	456	TString argument;
	457
	458	TObjArray* pTokens=allArgs.Tokenize(" ");
608cfbda	459	int bMissingParam=0;
608cfbda	460
8125805f	461	if (pTokens) {
	462	for (int i=0; i<pTokens->GetEntries() && iResult>=0; i++) {
	463	argument=((TObjString*)pTokens->At(i))->GetString();
	464	if (argument.IsNull()) continue;
608cfbda	465
	466	if (argument.CompareTo("-cutsGamma")==0) {
	467	TString spar = "";
	468	for( int j=0; j<8; j++ ){
	469	if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
	470	spar+=" ";
	471	spar+=((TObjString*)pTokens->At(i))->GetString();
	472	fGammaCuts[j]=((TObjString*)pTokens->At(i))->GetString().Atof();
	473	}
	474	if( !bMissingParam ){
	475	HLTInfo("Gamma cuts are set to: %s", spar.Data());
	476	continue;
	477	}
	478	} else if (argument.CompareTo("-cutsAP")==0) {
	479	TString spar = "";
	480	for( int j=0; j<8; j++ ){
	481	if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
	482	spar+=" ";
	483	spar+=((TObjString*)pTokens->At(i))->GetString();
	484	fAPCuts[j]=((TObjString*)pTokens->At(i))->GetString().Atof();
	485	}
	486	if( !bMissingParam ){
	487	HLTInfo("AP cuts are set to: %s", spar.Data());
	488	continue;
	489	}
	490	}
	491	else if (argument.CompareTo("-cutsKs")==0) {
	492	TString spar = "";
	493	for( int j=0; j<8; j++ ){
	494	if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
	495	spar+=" ";
	496	spar+=((TObjString*)pTokens->At(i))->GetString();
	497	fKsCuts[j]=((TObjString*)pTokens->At(i))->GetString().Atof();
	498	}
	499	if( !bMissingParam ){
	500	HLTInfo("KShort cuts are set to: %s", spar.Data());
	501	continue;
	502	}
	503	} else if (argument.CompareTo("-cutsLambda")==0) {
	504	TString spar = "";
	505	for( int j=0; j<8; j++ ){
	506	if ((bMissingParam=(++i>=pTokens->GetEntries()))) break;
	507	spar+=" ";
	508	spar+=((TObjString*)pTokens->At(i))->GetString();
	509	fLambdaCuts[j]=((TObjString*)pTokens->At(i))->GetString().Atof();
	510	}
	511	if( !bMissingParam ){
	512	HLTInfo("Lampda cuts are set to: %s", spar.Data());
	513	continue;
	514	}
	515	}else {
8125805f	516	HLTError("unknown argument %s", argument.Data());
	517	iResult=-EINVAL;
	518	break;
	519	}
	520	}
	521	delete pTokens;
	522	}
	523
	524	return iResult;
	525	}
	526
	527	int AliHLTV0HistoComponent::Reconfigure(const char* cdbEntry, const char* chainId)
	528	{
	529	// see header file for class documentation
	530	int iResult=0;
	531	const char* path="HLT/ConfigTPC/KryptonHistoComponent";
	532	const char* defaultNotify="";
	533	if (cdbEntry) {
	534	path=cdbEntry;
	535	defaultNotify=" (default)";
	536	}
	537	if (path) {
	538	HLTInfo("reconfigure from entry %s%s, chain id %s", path, defaultNotify,(chainId!=NULL && chainId[0]!=0)?chainId:"<none>");
	539	AliCDBEntry pEntry = AliCDBManager::Instance()->Get(path/,GetRunNo()*/);
	540	if (pEntry) {
	541	TObjString* pString=dynamic_cast<TObjString*>(pEntry->GetObject());
	542	if (pString) {
	543	HLTInfo("received configuration object string: \'%s\'", pString->GetString().Data());
	544	iResult=Configure(pString->GetString().Data());
	545	} else {
	546	HLTError("configuration object \"%s\" has wrong type, required TObjString", path);
	547	}
	548	} else {
	549	HLTError("can not fetch object \"%s\" from CDB", path);
	550	}
	551	}
	552
	553	return iResult;
	554	}