[u/mrichter/AliRoot.git] / HLT / trigger / AliHLTD0Trigger.cxx

// $Id: AliHLTD0Trigger.cxx 
//**************************************************************************
//* This file is property of and copyright by the ALICE HLT Project        * 
//* ALICE Experiment at CERN, All rights reserved.                         *
//*                                                                        *
//* Primary Authors: Gaute Ovrebekk                                        *
//*                  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   AliHLTD0Trigger.cxx
/// @author Gaute Ovrebekk
/// @date   2009-10-28
/// @brief  HLT trigger component for D0->Kpi

// see header file for class documentation
// or
// refer to README to build package
// or
// visit http://web.ift.uib.no/~kjeks/doc/alice-hlt

#include "AliHLTD0Trigger.h"
#include "AliESDEvent.h"
#include "AliESDv0.h"
#include "AliHLTTriggerDecision.h"
#include "AliHLTDomainEntry.h"
#include "AliHLTGlobalBarrelTrack.h"
#include "TObjArray.h"
#include "TObjString.h"
#include "TDatabasePDG.h"
#include "AliESDVertex.h"
#include "TH1F.h"
#include "AliHLTD0toKpi.h"
#include "AliAODVertex.h"
#include "AliESDVertex.h"
#include "AliAODRecoDecay.h"
#include "TMap.h"
#include "AliHLTD0Candidate.h"
#include "TFile.h"
#include "TClonesArray.h"

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

AliHLTD0Trigger::AliHLTD0Trigger()
  : AliHLTTrigger()
  , fPtMin(0.0)
  , fdca(0.0)
  , finvMass(0.0)     
  , fcosThetaStar(0.0)
  , fd0(0.0) 
  , fd0d0(0.0)
  , fcosPoint(0.0)
  , fplothisto(false)
  , fUseV0(false)
  , fD0PDG(TDatabasePDG::Instance()->GetParticle(421)->Mass())
  , fD0mass(NULL)
  , fD0pt(NULL)
  , fPos()
  , fNeg()
  , fd0calc(NULL)                  
  , ftwoTrackArray(NULL)
  , fTotalD0(0)
  , fuseKF(false)
  , fSendCandidate(false)
  , fCandidateTree(NULL)
  , fCandidateArray(NULL)
  , fWriteFile(false)
{
  
  // see header file for class documentation
  // or
  // refer to README to build package
  // or
  // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt

}

const char* AliHLTD0Trigger::fgkOCDBEntry="HLT/ConfigHLT/D0Trigger";

AliHLTD0Trigger::~AliHLTD0Trigger()
{
  // destructor 
}
  
const char* AliHLTD0Trigger::GetTriggerName() const
{
  // see header file for class documentation
  return "D0Trigger";
}

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

int AliHLTD0Trigger::DoTrigger()
{
  // -- Iterator over Data Blocks --
  //const AliHLTComponentBlockData* iter = NULL;

  if (!IsDataEvent()) return 0;

  HLTDebug("Used  Cuts: pT: %f , DCA: %f , InvMass: %f , CosThetaStar: %f , d0: %f , d0xd0: %f , CosPointingAngle: %f",fPtMin,fdca,finvMass,fcosThetaStar,fd0,fd0d0,fcosPoint);

  Int_t nD0=0;
  TString description;

  const AliESDVertex *Vertex=NULL;
  if(fCandidateArray){fCandidateArray->Clear();}
  
  for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {   
    if(fUseV0){
      nD0=RecV0(iter);
    }
    else{
       AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
       event->GetStdContent();
       Double_t field = event->GetMagneticField();
       Vertex = event->GetPrimaryVertexTracks();
       if(!Vertex || Vertex->GetNContributors()<2){
	 HLTDebug("Contributors in ESD vertex to low or not been set");
	 continue;
       }
       for(Int_t it=0;it<event->GetNumberOfTracks();it++){
	 SingleTrackSelect(event->GetTrack(it),Vertex,field);
       }
       RecD0(nD0,Vertex,field);       
    }
  }

  for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDVertex|kAliHLTDataOriginOut); 
	iter != NULL; iter = GetNextInputObject() ) { 
    Vertex = dynamic_cast<const AliESDVertex*>(iter);
    if(!Vertex){
      HLTError("Vertex object is corrupted");
      //iResult = -EINVAL;    
    }
  }  
   
  for ( const AliHLTComponentBlockData* iter = GetFirstInputBlock(kAliHLTDataTypeTrack|kAliHLTDataOriginITS); 
	iter != NULL; iter = GetNextInputBlock() ) { 
    vector<AliHLTGlobalBarrelTrack> tracksVector;
    AliHLTGlobalBarrelTrack::ConvertTrackDataArray(reinterpret_cast<const AliHLTTracksData*>(iter->fPtr), iter->fSize, tracksVector);
    
    Double_t field = GetBz();
    if (Vertex){
      for(UInt_t i=0;i<tracksVector.size();i++){
	SingleTrackSelect(&tracksVector[i],Vertex,field);
      }
      RecD0(nD0,Vertex,field);
    }    
  }
  
  fTotalD0+=nD0;

  ftwoTrackArray->Clear();
  fPos.clear();
  fNeg.clear();
     
  HLTDebug("Number of D0 found: %d",nD0);
  HLTDebug("Total Number of D0 found: %d",fTotalD0);
 
  if(fplothisto){
    PushBack( (TObject*) fD0mass, kAliHLTDataTypeHistogram,0);
    PushBack( (TObject*) fD0pt, kAliHLTDataTypeHistogram,0);
  }
  if(fSendCandidate){
    PushBack( (TObject*) fCandidateArray, kAliHLTDataTypeTObjArray,0);
  }
  if(fWriteFile && fCandidateTree!=NULL){
    fCandidateTree->Fill();
  }
  if(fCandidateArray){
    fCandidateArray->Clear(); 
  }

  //if (iResult>=0) {
  if (1) {
   
    if (nD0>=1) {
      description.Form("Event contains %d D0(s)", nD0);
      SetDescription(description.Data());
      // Enable the central detectors for readout.
      GetReadoutList().Enable(
			      AliHLTReadoutList::kITSSPD |
			      AliHLTReadoutList::kITSSDD |
			      AliHLTReadoutList::kITSSSD |
			      AliHLTReadoutList::kTPC |
			      AliHLTReadoutList::kTRD |
			      AliHLTReadoutList::kTOF |
			      AliHLTReadoutList::kHMPID |
			      AliHLTReadoutList::kPHOS
			      );
      // Add the available HLT information for readout too.
      GetTriggerDomain().Add("CLUSTERS", "TPC ");
      TriggerEvent(true);
      return 0;
    }
    description.Form("No D0");
  } else {
    description.Form("No input blocks found");
  }
  SetDescription(description.Data());
  TriggerEvent(false);
  //return iResult;
  return 0;
}

int AliHLTD0Trigger::DoInit(int argc, const char** argv)
{
  // component initialization

  fd0calc = new AliHLTD0toKpi();
  ftwoTrackArray = new TObjArray(2);
  // see header file for class documentation
  fD0mass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
  fD0pt = new TH1F("hPt","D^{0} Pt plot",20,0,20);
  // first configure the default
  int iResult=0;
  if (iResult>=0) iResult=ConfigureFromCDBTObjString(fgkOCDBEntry);

  // configure from the command line parameters if specified
  if (iResult>=0 && argc>0)
    iResult=ConfigureFromArgumentString(argc, argv);

  if (iResult<0) return iResult;

  if(fSendCandidate || fWriteFile){
    fCandidateArray = new TClonesArray("AliHLTD0Candidate",1000);
  }
  if(fWriteFile){
    fCandidateTree = new TTree("D0Candidates","Tree with D0 Candidates");
    fCandidateTree->Branch("Candidates","TClonesArray",&fCandidateArray,32000,99);
  }

  return iResult;
}

int AliHLTD0Trigger::DoDeinit()
{
  // see header file for class documentation
  if(fd0calc){delete fd0calc; fd0calc = NULL;}  
  if(fD0mass){delete fD0mass; fD0mass = NULL;}
  if(fD0pt){delete fD0pt; fD0pt=NULL;}
  if(ftwoTrackArray){delete ftwoTrackArray; ftwoTrackArray=NULL;}
  if(fWriteFile){
    TFile *mcFile = new TFile("D0_Candidates.root","RECREATE");
    fCandidateTree->Write();
    mcFile->Close();
  }
  if(fCandidateTree){delete fCandidateTree; fCandidateTree=NULL;}
  if(fCandidateArray){delete fCandidateArray; fCandidateArray=NULL;}
  return 0;
}

int AliHLTD0Trigger::Reconfigure(const char* cdbEntry, const char* /*chainId*/)
{
  // see header file for class documentation

  // configure from the specified antry or the default one
  const char* entry=cdbEntry;
  if (!entry || entry[0]==0) {
    entry=fgkOCDBEntry;
  }

  return ConfigureFromCDBTObjString(entry);
}

int AliHLTD0Trigger::ScanConfigurationArgument(int argc, const char** argv)
{
  // see header file for class documentation
  if (argc<=0) return 0;
  int i=0;
  TString argument=argv[i];
  // -minpt for decay
  if (argument.CompareTo("-pt")==0) {
    if (++i>=argc) return -EPROTO;
    argument=argv[i];
    fPtMin=argument.Atof();
    return 2;
  }    
  // minimum dca for decay tracks
  if (argument.CompareTo("-dca")==0) {
    if (++i>=argc) return -EPROTO;
    argument=argv[i];
    fdca=argument.Atof();
    return 2;
  }
  // inv. mass half width.
  if (argument.CompareTo("-invmass")==0) {
    if (++i>=argc) return -EPROTO;
    argument=argv[i];
    finvMass=argument.Atof();
    return 2;
  }

// cos theta for decay angle
  if (argument.CompareTo("-costhetastar")==0) {
    if (++i>=argc) return -EPROTO;
    argument=argv[i];
    fcosThetaStar=argument.Atof();
    return 2;
  }

  // impact parameter for decay
  if (argument.CompareTo("-d0")==0) {
    if (++i>=argc) return -EPROTO;
    argument=argv[i];
    fd0=argument.Atof();
    return 2;
  }
  // product of impact parameter
  if (argument.CompareTo("-d0d0")==0) {
    if (++i>=argc) return -EPROTO;
    argument=argv[i];
    fd0d0=argument.Atof();
    return 2;
  }
  // product of impact parameter
  if (argument.CompareTo("-cospoint")==0) {
    if (++i>=argc) return -EPROTO;
    argument=argv[i];
    fcosPoint=argument.Atof();
    return 2;
  }
  if (argument.CompareTo("-plothistogram")==0) {
    fplothisto=true;
    return 1;
  }
  if (argument.CompareTo("-usev0")==0) {
    fUseV0=true;
    return 1;
  }
  if (argument.CompareTo("-useKF")==0) {
    fuseKF=true;
    return 1;
  }
  if (argument.CompareTo("-send-candidates")==0) {
    fSendCandidate=true;
    return 1;
  }
  if (argument.CompareTo("-write-file")==0) {
    fWriteFile=true;
    return 1;
  }
 
  // unknown argument
  return -EINVAL;
}

void AliHLTD0Trigger::SingleTrackSelect(AliExternalTrackParam* t, const AliESDVertex* pv, Double_t field){
  // Offline uses || on the cuts of decay products
  if (!pv) return;

  Double_t pvpos[3];
  pv->GetXYZ(pvpos);

  if(t->Pt()<fPtMin){return;}
  if(TMath::Abs(t->GetD(pvpos[0],pvpos[1],field)) > fd0){return;}

  if(t->Charge()>0){
    fPos.push_back(t);
  }
  else{
    fNeg.push_back(t);
  }
}

void AliHLTD0Trigger::RecD0(Int_t& nD0,const AliESDVertex* pv, Double_t field){
  // Default way of reconstructing D0's. Both from ESD and Barreltracks
  Double_t D0,D0bar,xdummy,ydummy; 
  Double_t d0[2];
  Double_t svpos[3];
  Double_t pvpos[3];
  
  if(!pv){
    HLTError("No Vertex is set");
    return;
  }
  pv->GetXYZ(pvpos);
    
  for(UInt_t ip=0;ip<fPos.size();ip++){
    AliExternalTrackParam *tP=fPos[ip];
    for(UInt_t in=0;in<fNeg.size();in++){
      AliExternalTrackParam *tN=fNeg[in];
          
      tP->PropagateToDCA(pv,field,kVeryBig);  //do I need this??????
      tN->PropagateToDCA(pv,field,kVeryBig);
      
      Double_t dcatPtN = tP->GetDCA(tN,field,xdummy,ydummy);
      if(dcatPtN>fdca) { continue; }
      
      ftwoTrackArray->AddAt(tP,0);
      ftwoTrackArray->AddAt(tN,1);
      AliAODVertex *vertexp1n1 = fd0calc->ReconstructSecondaryVertex(ftwoTrackArray,field,pv,fuseKF);
      if(!vertexp1n1) { 
	ftwoTrackArray->Clear();
	continue; 
      }
      
      vertexp1n1->GetXYZ(svpos);
      
      tP->PropagateToDCA(vertexp1n1,field,kVeryBig); 
      tN->PropagateToDCA(vertexp1n1,field,kVeryBig);
      
      /*
      Double_t px[2],py[2],pz[2];
      Double_t momentum[3];
      tP->GetPxPyPz(momentum);
      px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
      tN->GetPxPyPz(momentum);
      px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];

      Short_t dummycharge=0;
      Double_t *dummyd0 = new Double_t[2];
      Int_t nprongs = 2;

      for(Int_t ipr=0;ipr<nprongs;ipr++) dummyd0[ipr]=0.;
      AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
      delete [] dummyd0; dummyd0=NULL;

      UInt_t pdg2[2],pdg2bar[2];
      Double_t mPDG,minv,minvbar;

      pdg2[0]=211; pdg2[1]=321; pdg2bar[0]=321; pdg2bar[1]=211;
      minv = rd->InvMass(nprongs,pdg2);
      minvbar = rd->InvMass(nprongs,pdg2bar);
      if(TMath::Abs(minv-fD0PDG)>finvMass && TMath::Abs(minv-fD0PDG)>finvMass) {continue; delete vertexp1n1; delete rd;}
      */
      if((TMath::Abs(fd0calc->InvMass(tN,tP)-fD0PDG)) > finvMass && TMath::Abs((fd0calc->InvMass(tP,tN))-fD0PDG) > finvMass){
	vertexp1n1->RemoveCovMatrix();
	delete vertexp1n1;
	ftwoTrackArray->Clear();
	continue;
      }
      
      fd0calc->cosThetaStar(tN,tP,D0,D0bar);
      if(TMath::Abs(D0) > fcosThetaStar && TMath::Abs(D0bar) > fcosThetaStar){
	vertexp1n1->RemoveCovMatrix();
	delete vertexp1n1;
	ftwoTrackArray->Clear();
	continue;
      }
      
      d0[0] = tP->GetD(pvpos[0],pvpos[1],field);
      d0[1] = tN->GetD(pvpos[0],pvpos[1],field);
      if((d0[0]*d0[1]) > fd0d0){
	vertexp1n1->RemoveCovMatrix();
	delete vertexp1n1;
	ftwoTrackArray->Clear();
	continue;
      }
      
      if(fd0calc->pointingAngle(tN,tP,pvpos,svpos) < fcosPoint){
	vertexp1n1->RemoveCovMatrix();
	delete vertexp1n1;
	ftwoTrackArray->Clear();	    
	continue;
      }
      
      if(fplothisto){
	//fD0mass->Fill(minv);
	//fD0mass->Fill(minvbar);
	fD0mass->Fill(fd0calc->InvMass(tN,tP));
	fD0mass->Fill(fd0calc->InvMass(tP,tN));
	fD0pt->Fill(fd0calc->Pt(tP,tN));
	/*
	if((fd0calc->InvMass(tN,tP) - fD0PDG) > finvMass){
	  fD0mass->Fill(fd0calc->InvMass(tN,tP));
	}
	else{
	  fD0mass->Fill(fd0calc->InvMass(tP,tN));
	  }
	*/
      }

      if(fSendCandidate && fCandidateArray){
	new(fCandidateArray->AddrAt(nD0)) AliHLTD0Candidate(fd0calc->InvMass(tN,tP),fd0calc->InvMass(tP,tN),fd0calc->Pt(tP,tN),tP->GetLabel(),tN->GetLabel(),tP->Pt(),tN->Pt());
      }

      nD0++;
      vertexp1n1->RemoveCovMatrix();
      delete vertexp1n1;
      ftwoTrackArray->Clear();
    }
  }
  ftwoTrackArray->Clear();
}

Int_t AliHLTD0Trigger::RecV0(const TObject* iter){
  // Reconstruction D0's from the V0 found by the V0 finder
  int nD0=0;
  Double_t d0[2];
  Double_t D0,D0bar; 
  Double_t svpos[3];
  Double_t pvpos[3];
  
  AliESDEvent *event = dynamic_cast<AliESDEvent*>(const_cast<TObject*>( iter ) );
  if (!event) {
    HLTError("input object is not of type AliESDEvent");
    return 0;
  }

  event->GetStdContent();
  Int_t nV0 = event->GetNumberOfV0s();
  Double_t field = event->GetMagneticField();
  const AliESDVertex* pv = event->GetPrimaryVertexTracks();
  pv->GetXYZ(pvpos);

  for (Int_t iv=0; iv<nV0; iv++) {
      
      AliESDtrack *tN=event->GetTrack( event->GetV0(iv)->GetNindex());
      AliESDtrack *tP=event->GetTrack( event->GetV0(iv)->GetPindex());      
      
      if(tN->Pt()<fPtMin && tP->Pt()<fPtMin){continue;} //||????
      
      d0[0] = tP->GetD(pvpos[0],pvpos[1],field);
      d0[1] = tN->GetD(pvpos[0],pvpos[1],field);
      
      if(d0[0]>fd0 && d0[0]>fd0){continue;} // make sure < or>
      
      event->GetV0(iv)->GetXYZ(svpos[0],svpos[1],svpos[2]);
      
      if(!tN->PropagateTo(svpos[0],field) && !tP->PropagateTo(svpos[0],field)){
      	HLTInfo("Tracks could not be propagated to secondary vertex");
      	continue;
      }
      
      Double_t tmp1, tmp2;
      if(tN->GetDCA(tP,field,tmp1,tmp2) > fdca){continue;}
      
      if((fd0calc->InvMass(tN,tP) - fD0PDG) > finvMass && (fd0calc->InvMass(tP,tN) - fD0PDG) > finvMass){continue;}
      fd0calc->cosThetaStar(tN,tP,D0,D0bar);
      if(D0 > fcosThetaStar && D0bar > fcosThetaStar){continue;}
      if((d0[0]*d0[1]) > fd0d0){continue;}
      if(fd0calc->pointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
      
      nD0++;
      if((fd0calc->InvMass(tN,tP) - fD0PDG) > finvMass){
	fD0mass->Fill(fd0calc->InvMass(tN,tP));
      }
      else{
	fD0mass->Fill(fd0calc->InvMass(tP,tN));
      }
  }
  return nD0;
  
}
void AliHLTD0Trigger::GetOCDBObjectDescription( TMap* const targetMap)
{
  // Get a list of OCDB object description.
  if (!targetMap) return;
  targetMap->Add(new TObjString("HLT/ConfigHLT/D0Trigger"),new TObjString("Object with default cuts for D0 reconstruction" ));
}
Commit	Line	Data
5f4502cc	1	// $Id: AliHLTD0Trigger.cxx
	2	//**************************************************************************
	3	//* This file is property of and copyright by the ALICE HLT Project *
	4	//* ALICE Experiment at CERN, All rights reserved. *
	5	//* *
	6	//* Primary Authors: Gaute Ovrebekk *
	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 AliHLTD0Trigger.cxx
	19	/// @author Gaute Ovrebekk
	20	/// @date 2009-10-28
	21	/// @brief HLT trigger component for D0->Kpi
	22
	23	// see header file for class documentation
	24	// or
	25	// refer to README to build package
	26	// or
	27	// visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
	28
	29	#include "AliHLTD0Trigger.h"
	30	#include "AliESDEvent.h"
5f4502cc	31	#include "AliESDv0.h"
	32	#include "AliHLTTriggerDecision.h"
	33	#include "AliHLTDomainEntry.h"
	34	#include "AliHLTGlobalBarrelTrack.h"
	35	#include "TObjArray.h"
	36	#include "TObjString.h"
	37	#include "TDatabasePDG.h"
	38	#include "AliESDVertex.h"
	39	#include "TH1F.h"
	40	#include "AliHLTD0toKpi.h"
9929f8f5	41	#include "AliAODVertex.h"
629b904b	42	#include "AliESDVertex.h"
3a626a12	43	#include "AliAODRecoDecay.h"
10d746a6	44	#include "TMap.h"
0f4240c3	45	#include "AliHLTD0Candidate.h"
	46	#include "TFile.h"
	47	#include "TClonesArray.h"
5f4502cc	48
	49	/** ROOT macro for the implementation of ROOT specific class methods */
	50	ClassImp(AliHLTD0Trigger)
	51
	52	AliHLTD0Trigger::AliHLTD0Trigger()
	53	: AliHLTTrigger()
	54	, fPtMin(0.0)
	55	, fdca(0.0)
	56	, finvMass(0.0)
	57	, fcosThetaStar(0.0)
	58	, fd0(0.0)
	59	, fd0d0(0.0)
	60	, fcosPoint(0.0)
9929f8f5	61	, fplothisto(false)
9929f8f5	62	, fUseV0(false)
052f676e	63	, fD0PDG(TDatabasePDG::Instance()->GetParticle(421)->Mass())
5f4502cc	64	, fD0mass(NULL)
252bba71	65	, fD0pt(NULL)
9929f8f5	66	, fPos()
	67	, fNeg()
	68	, fd0calc(NULL)
	69	, ftwoTrackArray(NULL)
	70	, fTotalD0(0)
15f6cee9	71	, fuseKF(false)
0f4240c3	72	, fSendCandidate(false)
	73	, fCandidateTree(NULL)
	74	, fCandidateArray(NULL)
	75	, fWriteFile(false)
5f4502cc	76	{
9929f8f5	77
5f4502cc	78	// see header file for class documentation
	79	// or
	80	// refer to README to build package
	81	// or
	82	// visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
	83
	84	}
	85
	86	const char* AliHLTD0Trigger::fgkOCDBEntry="HLT/ConfigHLT/D0Trigger";
	87
	88	AliHLTD0Trigger::~AliHLTD0Trigger()
	89	{
0f4240c3	90	// destructor
5f4502cc	91	}
9929f8f5	92
5f4502cc	93	const char* AliHLTD0Trigger::GetTriggerName() const
	94	{
	95	// see header file for class documentation
	96	return "D0Trigger";
	97	}
	98
	99	AliHLTComponent* AliHLTD0Trigger::Spawn()
	100	{
	101	// see header file for class documentation
	102	return new AliHLTD0Trigger;
	103	}
	104
	105	int AliHLTD0Trigger::DoTrigger()
	106	{
9929f8f5	107	// -- Iterator over Data Blocks --
9929f8f5	108	//const AliHLTComponentBlockData* iter = NULL;
5f4502cc	109
9929f8f5	110	if (!IsDataEvent()) return 0;
5f4502cc	111
9f96f527	112	HLTDebug("Used Cuts: pT: %f , DCA: %f , InvMass: %f , CosThetaStar: %f , d0: %f , d0xd0: %f , CosPointingAngle: %f",fPtMin,fdca,finvMass,fcosThetaStar,fd0,fd0d0,fcosPoint);
9f96f527	113
9929f8f5	114	Int_t nD0=0;
9929f8f5	115	TString description;
5f4502cc	116
a7cf6e02	117	const AliESDVertex *Vertex=NULL;
0f4240c3	118	if(fCandidateArray){fCandidateArray->Clear();}
0f4240c3	119
9929f8f5	120	for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDObject); iter != NULL; iter = GetNextInputObject() ) {
	121	if(fUseV0){
	122	nD0=RecV0(iter);
	123	}
	124	else{
629b904b	125	AliESDEvent event = dynamic_cast<AliESDEvent>(const_cast<TObject*>( iter ) );
629b904b	126	event->GetStdContent();
14978adc	127	Double_t field = event->GetMagneticField();
a7cf6e02	128	Vertex = event->GetPrimaryVertexTracks();
a7cf6e02	129	if(!Vertex \|\| Vertex->GetNContributors()<2){
9f96f527	130	HLTDebug("Contributors in ESD vertex to low or not been set");
3a626a12	131	continue;
3a626a12	132	}
629b904b	133	for(Int_t it=0;it<event->GetNumberOfTracks();it++){
14978adc	134	SingleTrackSelect(event->GetTrack(it),Vertex,field);
629b904b	135	}
0f4240c3	136	RecD0(nD0,Vertex,field);
5f4502cc	137	}
5f4502cc	138	}
9929f8f5	139
3a626a12	140	for ( const TObject *iter = GetFirstInputObject(kAliHLTDataTypeESDVertex\|kAliHLTDataOriginOut);
629b904b	141	iter != NULL; iter = GetNextInputObject() ) {
a7cf6e02	142	Vertex = dynamic_cast<const AliESDVertex*>(iter);
a7cf6e02	143	if(!Vertex){
3a626a12	144	HLTError("Vertex object is corrupted");
629b904b	145	//iResult = -EINVAL;
	146	}
	147	}
	148
	149	for ( const AliHLTComponentBlockData* iter = GetFirstInputBlock(kAliHLTDataTypeTrack\|kAliHLTDataOriginITS);
	150	iter != NULL; iter = GetNextInputBlock() ) {
	151	vector<AliHLTGlobalBarrelTrack> tracksVector;
	152	AliHLTGlobalBarrelTrack::ConvertTrackDataArray(reinterpret_cast<const AliHLTTracksData*>(iter->fPtr), iter->fSize, tracksVector);
	153
14978adc	154	Double_t field = GetBz();
a7cf6e02	155	if (Vertex){
a7cf6e02	156	for(UInt_t i=0;i<tracksVector.size();i++){
14978adc	157	SingleTrackSelect(&tracksVector[i],Vertex,field);
a7cf6e02	158	}
0f4240c3	159	RecD0(nD0,Vertex,field);
a7cf6e02	160	}
	161	}
	162
9929f8f5	163	fTotalD0+=nD0;
2a743e38	164
9929f8f5	165	ftwoTrackArray->Clear();
	166	fPos.clear();
	167	fNeg.clear();
629b904b	168
2a743e38	169	HLTDebug("Number of D0 found: %d",nD0);
2a743e38	170	HLTDebug("Total Number of D0 found: %d",fTotalD0);
0f4240c3	171
252bba71	172	if(fplothisto){
	173	PushBack( (TObject*) fD0mass, kAliHLTDataTypeHistogram,0);
	174	PushBack( (TObject*) fD0pt, kAliHLTDataTypeHistogram,0);
	175	}
0f4240c3	176	if(fSendCandidate){
	177	PushBack( (TObject*) fCandidateArray, kAliHLTDataTypeTObjArray,0);
	178	}
	179	if(fWriteFile && fCandidateTree!=NULL){
	180	fCandidateTree->Fill();
	181	}
	182	if(fCandidateArray){
	183	fCandidateArray->Clear();
	184	}
	185
9929f8f5	186	//if (iResult>=0) {
9929f8f5	187	if (1) {
5f4502cc	188
	189	if (nD0>=1) {
	190	description.Form("Event contains %d D0(s)", nD0);
	191	SetDescription(description.Data());
	192	// Enable the central detectors for readout.
	193	GetReadoutList().Enable(
	194	AliHLTReadoutList::kITSSPD \|
	195	AliHLTReadoutList::kITSSDD \|
	196	AliHLTReadoutList::kITSSSD \|
	197	AliHLTReadoutList::kTPC \|
	198	AliHLTReadoutList::kTRD \|
	199	AliHLTReadoutList::kTOF \|
	200	AliHLTReadoutList::kHMPID \|
	201	AliHLTReadoutList::kPHOS
	202	);
	203	// Add the available HLT information for readout too.
	204	GetTriggerDomain().Add("CLUSTERS", "TPC ");
	205	TriggerEvent(true);
	206	return 0;
	207	}
	208	description.Form("No D0");
	209	} else {
	210	description.Form("No input blocks found");
	211	}
	212	SetDescription(description.Data());
	213	TriggerEvent(false);
9929f8f5	214	//return iResult;
9929f8f5	215	return 0;
5f4502cc	216	}
	217
	218	int AliHLTD0Trigger::DoInit(int argc, const char** argv)
	219	{
0f4240c3	220	// component initialization
0f4240c3	221
9929f8f5	222	fd0calc = new AliHLTD0toKpi();
9929f8f5	223	ftwoTrackArray = new TObjArray(2);
5f4502cc	224	// see header file for class documentation
9929f8f5	225	fD0mass = new TH1F("hMass","D^{0} mass plot",100,1.7,2);
252bba71	226	fD0pt = new TH1F("hPt","D^{0} Pt plot",20,0,20);
5f4502cc	227	// first configure the default
	228	int iResult=0;
	229	if (iResult>=0) iResult=ConfigureFromCDBTObjString(fgkOCDBEntry);
	230
	231	// configure from the command line parameters if specified
	232	if (iResult>=0 && argc>0)
	233	iResult=ConfigureFromArgumentString(argc, argv);
0f4240c3	234
	235	if (iResult<0) return iResult;
	236
	237	if(fSendCandidate \|\| fWriteFile){
	238	fCandidateArray = new TClonesArray("AliHLTD0Candidate",1000);
	239	}
	240	if(fWriteFile){
	241	fCandidateTree = new TTree("D0Candidates","Tree with D0 Candidates");
	242	fCandidateTree->Branch("Candidates","TClonesArray",&fCandidateArray,32000,99);
	243	}
	244
5f4502cc	245	return iResult;
	246	}
	247
	248	int AliHLTD0Trigger::DoDeinit()
	249	{
	250	// see header file for class documentation
a7cf6e02	251	if(fd0calc){delete fd0calc; fd0calc = NULL;}
a7cf6e02	252	if(fD0mass){delete fD0mass; fD0mass = NULL;}
925bbd08	253	if(fD0pt){delete fD0pt; fD0pt=NULL;}
925bbd08	254	if(ftwoTrackArray){delete ftwoTrackArray; ftwoTrackArray=NULL;}
0f4240c3	255	if(fWriteFile){
	256	TFile *mcFile = new TFile("D0_Candidates.root","RECREATE");
	257	fCandidateTree->Write();
	258	mcFile->Close();
	259	}
	260	if(fCandidateTree){delete fCandidateTree; fCandidateTree=NULL;}
	261	if(fCandidateArray){delete fCandidateArray; fCandidateArray=NULL;}
5f4502cc	262	return 0;
	263	}
	264
	265	int AliHLTD0Trigger::Reconfigure(const char* cdbEntry, const char* /chainId/)
	266	{
	267	// see header file for class documentation
	268
	269	// configure from the specified antry or the default one
	270	const char* entry=cdbEntry;
	271	if (!entry \|\| entry[0]==0) {
	272	entry=fgkOCDBEntry;
	273	}
	274
	275	return ConfigureFromCDBTObjString(entry);
	276	}
	277
	278	int AliHLTD0Trigger::ScanConfigurationArgument(int argc, const char** argv)
	279	{
	280	// see header file for class documentation
	281	if (argc<=0) return 0;
	282	int i=0;
	283	TString argument=argv[i];
	284	// -minpt for decay
	285	if (argument.CompareTo("-pt")==0) {
	286	if (++i>=argc) return -EPROTO;
	287	argument=argv[i];
	288	fPtMin=argument.Atof();
	289	return 2;
	290	}
	291	// minimum dca for decay tracks
	292	if (argument.CompareTo("-dca")==0) {
	293	if (++i>=argc) return -EPROTO;
	294	argument=argv[i];
	295	fdca=argument.Atof();
	296	return 2;
	297	}
	298	// inv. mass half width.
	299	if (argument.CompareTo("-invmass")==0) {
	300	if (++i>=argc) return -EPROTO;
	301	argument=argv[i];
	302	finvMass=argument.Atof();
	303	return 2;
	304	}
	305
	306	// cos theta for decay angle
	307	if (argument.CompareTo("-costhetastar")==0) {
	308	if (++i>=argc) return -EPROTO;
	309	argument=argv[i];
	310	fcosThetaStar=argument.Atof();
	311	return 2;
	312	}
	313
	314	// impact parameter for decay
	315	if (argument.CompareTo("-d0")==0) {
	316	if (++i>=argc) return -EPROTO;
	317	argument=argv[i];
	318	fd0=argument.Atof();
	319	return 2;
	320	}
	321	// product of impact parameter
	322	if (argument.CompareTo("-d0d0")==0) {
	323	if (++i>=argc) return -EPROTO;
	324	argument=argv[i];
	325	fd0d0=argument.Atof();
326	return 2;
327	}
328	// product of impact parameter
329	if (argument.CompareTo("-cospoint")==0) {
330	if (++i>=argc) return -EPROTO;
331	argument=argv[i];
332	fcosPoint=argument.Atof();
333	return 2;
334	}
335	if (argument.CompareTo("-plothistogram")==0) {
5f4502cc	336	fplothisto=true;
9929f8f5	337	return 1;
	338	}
	339	if (argument.CompareTo("-usev0")==0) {
	340	fUseV0=true;
	341	return 1;
5f4502cc	342	}
15f6cee9	343	if (argument.CompareTo("-useKF")==0) {
	344	fuseKF=true;
	345	return 1;
	346	}
0f4240c3	347	if (argument.CompareTo("-send-candidates")==0) {
	348	fSendCandidate=true;
	349	return 1;
	350	}
	351	if (argument.CompareTo("-write-file")==0) {
	352	fWriteFile=true;
	353	return 1;
	354	}
	355
5f4502cc	356	// unknown argument
	357	return -EINVAL;
	358	}
9929f8f5	359
0f4240c3	360	void AliHLTD0Trigger::SingleTrackSelect(AliExternalTrackParam* t, const AliESDVertex* pv, Double_t field){
a7cf6e02	361	// Offline uses \|\| on the cuts of decay products
a7cf6e02	362	if (!pv) return;
8d9d674e	363
a7cf6e02	364	Double_t pvpos[3];
a7cf6e02	365	pv->GetXYZ(pvpos);
629b904b	366
9929f8f5	367	if(t->Pt()<fPtMin){return;}
14978adc	368	if(TMath::Abs(t->GetD(pvpos[0],pvpos[1],field)) > fd0){return;}
9929f8f5	369
	370	if(t->Charge()>0){
	371	fPos.push_back(t);
	372	}
	373	else{
	374	fNeg.push_back(t);
	375	}
	376	}
	377
0f4240c3	378	void AliHLTD0Trigger::RecD0(Int_t& nD0,const AliESDVertex* pv, Double_t field){
a7cf6e02	379	// Default way of reconstructing D0's. Both from ESD and Barreltracks
9929f8f5	380	Double_t D0,D0bar,xdummy,ydummy;
	381	Double_t d0[2];
	382	Double_t svpos[3];
	383	Double_t pvpos[3];
	384
a7cf6e02	385	if(!pv){
629b904b	386	HLTError("No Vertex is set");
2a743e38	387	return;
9929f8f5	388	}
a7cf6e02	389	pv->GetXYZ(pvpos);
629b904b	390
	391	for(UInt_t ip=0;ip<fPos.size();ip++){
	392	AliExternalTrackParam *tP=fPos[ip];
	393	for(UInt_t in=0;in<fNeg.size();in++){
	394	AliExternalTrackParam *tN=fNeg[in];
	395
14978adc	396	tP->PropagateToDCA(pv,field,kVeryBig); //do I need this??????
14978adc	397	tN->PropagateToDCA(pv,field,kVeryBig);
9929f8f5	398
14978adc	399	Double_t dcatPtN = tP->GetDCA(tN,field,xdummy,ydummy);
9929f8f5	400	if(dcatPtN>fdca) { continue; }
9929f8f5	401
629b904b	402	ftwoTrackArray->AddAt(tP,0);
629b904b	403	ftwoTrackArray->AddAt(tN,1);
14978adc	404	AliAODVertex *vertexp1n1 = fd0calc->ReconstructSecondaryVertex(ftwoTrackArray,field,pv,fuseKF);
9929f8f5	405	if(!vertexp1n1) {
	406	ftwoTrackArray->Clear();
	407	continue;
	408	}
	409
	410	vertexp1n1->GetXYZ(svpos);
	411
14978adc	412	tP->PropagateToDCA(vertexp1n1,field,kVeryBig);
14978adc	413	tN->PropagateToDCA(vertexp1n1,field,kVeryBig);
9929f8f5	414
3a626a12	415	/*
	416	Double_t px[2],py[2],pz[2];
	417	Double_t momentum[3];
	418	tP->GetPxPyPz(momentum);
	419	px[0] = momentum[0]; py[0] = momentum[1]; pz[0] = momentum[2];
	420	tN->GetPxPyPz(momentum);
	421	px[1] = momentum[0]; py[1] = momentum[1]; pz[1] = momentum[2];
	422
	423	Short_t dummycharge=0;
	424	Double_t *dummyd0 = new Double_t[2];
	425	Int_t nprongs = 2;
	426
	427	for(Int_t ipr=0;ipr<nprongs;ipr++) dummyd0[ipr]=0.;
	428	AliAODRecoDecay *rd = new AliAODRecoDecay(0x0,nprongs,dummycharge,px,py,pz,dummyd0);
	429	delete [] dummyd0; dummyd0=NULL;
	430
	431	UInt_t pdg2[2],pdg2bar[2];
	432	Double_t mPDG,minv,minvbar;
	433
	434	pdg2[0]=211; pdg2[1]=321; pdg2bar[0]=321; pdg2bar[1]=211;
	435	minv = rd->InvMass(nprongs,pdg2);
	436	minvbar = rd->InvMass(nprongs,pdg2bar);
052f676e	437	if(TMath::Abs(minv-fD0PDG)>finvMass && TMath::Abs(minv-fD0PDG)>finvMass) {continue; delete vertexp1n1; delete rd;}
3a626a12	438	*/
0f4240c3	439	if((TMath::Abs(fd0calc->InvMass(tN,tP)-fD0PDG)) > finvMass && TMath::Abs((fd0calc->InvMass(tP,tN))-fD0PDG) > finvMass){
	440	vertexp1n1->RemoveCovMatrix();
	441	delete vertexp1n1;
	442	ftwoTrackArray->Clear();
	443	continue;
	444	}
	445
9929f8f5	446	fd0calc->cosThetaStar(tN,tP,D0,D0bar);
0f4240c3	447	if(TMath::Abs(D0) > fcosThetaStar && TMath::Abs(D0bar) > fcosThetaStar){
	448	vertexp1n1->RemoveCovMatrix();
	449	delete vertexp1n1;
	450	ftwoTrackArray->Clear();
	451	continue;
	452	}
	453
14978adc	454	d0[0] = tP->GetD(pvpos[0],pvpos[1],field);
14978adc	455	d0[1] = tN->GetD(pvpos[0],pvpos[1],field);
0f4240c3	456	if((d0[0]*d0[1]) > fd0d0){
	457	vertexp1n1->RemoveCovMatrix();
	458	delete vertexp1n1;
	459	ftwoTrackArray->Clear();
	460	continue;
	461	}
	462
	463	if(fd0calc->pointingAngle(tN,tP,pvpos,svpos) < fcosPoint){
	464	vertexp1n1->RemoveCovMatrix();
	465	delete vertexp1n1;
	466	ftwoTrackArray->Clear();
	467	continue;
	468	}
9929f8f5	469
9929f8f5	470	if(fplothisto){
3a626a12	471	//fD0mass->Fill(minv);
	472	//fD0mass->Fill(minvbar);
	473	fD0mass->Fill(fd0calc->InvMass(tN,tP));
	474	fD0mass->Fill(fd0calc->InvMass(tP,tN));
252bba71	475	fD0pt->Fill(fd0calc->Pt(tP,tN));
3a626a12	476	/*
052f676e	477	if((fd0calc->InvMass(tN,tP) - fD0PDG) > finvMass){
9929f8f5	478	fD0mass->Fill(fd0calc->InvMass(tN,tP));
	479	}
	480	else{
	481	fD0mass->Fill(fd0calc->InvMass(tP,tN));
3a626a12	482	}
3a626a12	483	*/
9929f8f5	484	}
0f4240c3	485
	486	if(fSendCandidate && fCandidateArray){
	487	new(fCandidateArray->AddrAt(nD0)) AliHLTD0Candidate(fd0calc->InvMass(tN,tP),fd0calc->InvMass(tP,tN),fd0calc->Pt(tP,tN),tP->GetLabel(),tN->GetLabel(),tP->Pt(),tN->Pt());
9f96f527	488	}
2a743e38	489
9929f8f5	490	nD0++;
052f676e	491	vertexp1n1->RemoveCovMatrix();
9929f8f5	492	delete vertexp1n1;
0f4240c3	493	ftwoTrackArray->Clear();
9929f8f5	494	}
9929f8f5	495	}
0f4240c3	496	ftwoTrackArray->Clear();
9929f8f5	497	}
	498
	499	Int_t AliHLTD0Trigger::RecV0(const TObject* iter){
a7cf6e02	500	// Reconstruction D0's from the V0 found by the V0 finder
9929f8f5	501	int nD0=0;
	502	Double_t d0[2];
	503	Double_t D0,D0bar;
	504	Double_t svpos[3];
	505	Double_t pvpos[3];
	506
	507	AliESDEvent event = dynamic_cast<AliESDEvent>(const_cast<TObject*>( iter ) );
55f4e3b5	508	if (!event) {
	509	HLTError("input object is not of type AliESDEvent");
	510	return 0;
	511	}
	512
9929f8f5	513	event->GetStdContent();
	514	Int_t nV0 = event->GetNumberOfV0s();
	515	Double_t field = event->GetMagneticField();
	516	const AliESDVertex* pv = event->GetPrimaryVertexTracks();
	517	pv->GetXYZ(pvpos);
	518
	519	for (Int_t iv=0; iv<nV0; iv++) {
	520
	521	AliESDtrack *tN=event->GetTrack( event->GetV0(iv)->GetNindex());
	522	AliESDtrack *tP=event->GetTrack( event->GetV0(iv)->GetPindex());
	523
	524	if(tN->Pt()<fPtMin && tP->Pt()<fPtMin){continue;} //\|\|????
	525
	526	d0[0] = tP->GetD(pvpos[0],pvpos[1],field);
	527	d0[1] = tN->GetD(pvpos[0],pvpos[1],field);
	528
	529	if(d0[0]>fd0 && d0[0]>fd0){continue;} // make sure < or>
	530
	531	event->GetV0(iv)->GetXYZ(svpos[0],svpos[1],svpos[2]);
	532
	533	if(!tN->PropagateTo(svpos[0],field) && !tP->PropagateTo(svpos[0],field)){
	534	HLTInfo("Tracks could not be propagated to secondary vertex");
	535	continue;
	536	}
	537
	538	Double_t tmp1, tmp2;
	539	if(tN->GetDCA(tP,field,tmp1,tmp2) > fdca){continue;}
	540
052f676e	541	if((fd0calc->InvMass(tN,tP) - fD0PDG) > finvMass && (fd0calc->InvMass(tP,tN) - fD0PDG) > finvMass){continue;}
9929f8f5	542	fd0calc->cosThetaStar(tN,tP,D0,D0bar);
	543	if(D0 > fcosThetaStar && D0bar > fcosThetaStar){continue;}
	544	if((d0[0]*d0[1]) > fd0d0){continue;}
	545	if(fd0calc->pointingAngle(tN,tP,pvpos,svpos) < fcosPoint){continue;}
	546
	547	nD0++;
052f676e	548	if((fd0calc->InvMass(tN,tP) - fD0PDG) > finvMass){
9929f8f5	549	fD0mass->Fill(fd0calc->InvMass(tN,tP));
	550	}
	551	else{
	552	fD0mass->Fill(fd0calc->InvMass(tP,tN));
	553	}
	554	}
	555	return nD0;
	556
	557	}
10d746a6	558	void AliHLTD0Trigger::GetOCDBObjectDescription( TMap* const targetMap)
	559	{
	560	// Get a list of OCDB object description.
	561	if (!targetMap) return;
	562	targetMap->Add(new TObjString("HLT/ConfigHLT/D0Trigger"),new TObjString("Object with default cuts for D0 reconstruction" ));
	563	}