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

// $Id$
//**************************************************************************
//* 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
*/

#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>

using namespace std;

/** 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
  for( int i=0; i<8; i++){
    fGammaCuts[i] = 0;
    fKsCuts[i] = 0;
    fLambdaCuts[i] = 0;
    fAPCuts[i] = 0;
  }
}

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