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

#include "AliHLTTriggerFastJet.h"
#include "AliESDEvent.h"
#include "AliESDtrackCuts.h"
#include "AliESDCaloCluster.h"
#include "AliHLTCaloClusterReader.h"
#include "AliHLTCaloClusterDataStruct.h"
#include "AliESDVZERO.h"
#include "AliHLTTriggerDecision.h"
#include "AliHLTDomainEntry.h"
#include "TLorentzVector.h"
#include "TRefArray.h"
#include "TString.h"
#include "TMap.h"
#include "TObjArray.h"
#include "TArrayI.h"

#include "AliFJWrapper.h"


ClassImp(AliHLTTriggerFastJet)

AliHLTTriggerFastJet::AliHLTTriggerFastJet() :
AliHLTTrigger(),
  fEThreshold(0.0),
  fDetector("EMCAL"),
  fFastJetWrapper(NULL),
  EsdTrackCuts(NULL),
  fOCDBEntry("HLT/ConfigHLT/EmcalJetTrigger")
{
  fFastJetWrapper = new AliFJWrapper("FastJet","FastJet");

}
//_____________________________________________________________

AliHLTTriggerFastJet::~AliHLTTriggerFastJet() {

  if (fFastJetWrapper)
    delete fFastJetWrapper;
  fFastJetWrapper = NULL;

  if (EsdTrackCuts)
    delete EsdTrackCuts;
  EsdTrackCuts = NULL;
  
}
//_____________________________________________________________

int AliHLTTriggerFastJet::DoInit(int argc, const char** argv) {
  
  int iResult = ConfigureFromCDBTObjString(fOCDBEntry);
  
  if ( iResult>=0 && argc>0 ) {
    iResult = ConfigureFromArgumentString(argc, argv);
    HLTImportant("Trigger threshold set from argument string :  %.02f GeV:", fEThreshold );
  } else if ( iResult>=0 ) {
    HLTImportant("Trigger threshold set from OCDB database entry : %.02f Gev:", fEThreshold );
  }
  return iResult;
}
//______________________________________________________________

int AliHLTTriggerFastJet::DoDeInit() {
  
  return 0;

}
//______________________________________________________________


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


const char* AliHLTTriggerFastJet::GetTriggerName() const {
  // see header file for class documentation
  return "EmcalJetTrigger";
}

Int_t AliHLTTriggerFastJet::DoTrigger() {

  Int_t iResult = 0;

  if ( GetFirstInputBlock( kAliHLTDataTypeSOR ) || GetFirstInputBlock( kAliHLTDataTypeEOR ) )
    return 0;

  const TObject *obj = GetFirstInputObject( kAliHLTAllDataTypes, "AliESDEvent" );
  AliESDEvent   *esd = dynamic_cast<AliESDEvent*>(const_cast<TObject*>(obj));
  
  // check for MatchedTrack to avoid double counting
  

  if ( esd != NULL ) {
    esd->GetStdContent();
    
    //Double_t vertex[3] = {0,0,0};
    //esd->GetVertex()->GetXYZ(vertex);
    //cout << "Vertex " << vertex[0] << vertex[1] << vertex [2] << endl;
    //TLorentzVector gamma;
    
    // -- add VZERO
    // uncomment for usage (to avoid warning)    
    
    /*
    AliESDVZERO *v0  = esd->GetVZEROData();
    
    Float_t MtotVOA = v0->GetMTotV0A();
    Float_t MtotVOC = v0->GetMTotV0A();
    */
    Int_t nTracks = esd->GetNumberOfTracks();

    TObjArray *tracks = new TObjArray(nTracks);
    tracks->SetOwner(1);
    

    EsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
    EsdTrackCuts->SetMinNClustersTPC(70);

    // add the tracks to FastJet Wrapper                                                          

    for ( Int_t i = 0; i < nTracks; i++ ) {

      const AliESDVertex *vtxTPC = esd->GetPrimaryVertexTPC();
      if (!vtxTPC) continue;
      cout << "Got primary vertex " << endl;

      AliESDtrack *esdtrack = esd->GetTrack(i);
      if (!esdtrack) continue;
      tracks->Add(esdtrack);
      cout << "Got esd track " << endl;

      AliESDtrack *tpctrack = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdtrack->GetID());
      if (!tpctrack) continue;
      cout << "Got tpc track " << endl;

      AliExternalTrackParam exParam;
      Bool_t relate = tpctrack->RelateToVertexTPC(vtxTPC,esd->GetMagneticField(),kVeryBig,&exParam);
      if (!relate) {
	delete tpctrack;
	continue;
      }

      tpctrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
      cout << "Set param to tpctrack " << endl;

      fFastJetWrapper->AddInputVector(tpctrack->Px(),tpctrack->Py(),tpctrack->Pz(),tpctrack->P());
      cout << "Added track with P: " << tpctrack->P() << endl;
    }


//    // add the Calo Clusters to FastJet Wrapper
//    TRefArray *caloClustersArr = new TRefArray();
//    esd->GetEMCALClusters(caloClustersArr);
//    const Int_t nEMCALClusters = caloClustersArr->GetEntries();
//
//    for ( UInt_t ic = 0; ic < nEMCALClusters; ic++ )
//      {
//	AliESDCaloCluster *cluster = dynamic_cast<AliESDCaloCluster*>( caloClustersArr->At(ic) );
//        Double_t clusterEn = cluster->E();
//	Double_t subE = 0;
//
//        TArrayI *matched = cluster->GetTracksMatched();
//	for ( UInt_t jk = 0; jk < matched->GetSize(); jk++ )
//          {
//            if ( matched->At(jk) > -1 && matched->At(jk) < esd ->GetNumberOfTracks() )
//              {
//                AliESDtrack *track = dynamic_cast<AliESDtrack*>( tracks->At( (matched->At(jk)) ) );
//                if ( !track ) continue;
//                subE += track->P();
//              }
//          }
//        clusterEn -= subE;
//        if ( clusterEn < 0)
//          clusterEn = 0;
//        cluster->SetE( clusterEn );
//        if ( cluster->E() == 0 ) continue;
//
//        cluster->GetMomentum( gamma, vertex );
//        fFastJetWrapper->AddInputVector( gamma.Px(), gamma.Py(), gamma.Pz(), clusterEn );
//        gamma.Clear();
//      }
//  }
    for ( Int_t j = 0; j< esd->GetNumberOfCaloClusters(); j++) {
      //AliESDCaloCluster *cluster = static_cast<AliESDCaloCluster*>(esd->GetCaloCluster(j));
      AliVCluster *cluster = dynamic_cast<AliVCluster*>(esd->GetCaloCluster(j));
      if ( cluster->IsEMCAL() ) {
        Int_t trackindex = cluster->GetTrackMatchedIndex();
   	if (trackindex<0) continue;
   	AliESDtrack *track = esd->GetTrack(trackindex);
   	if (!track) continue;
   	fFastJetWrapper->AddInputVector(track->Px(),track->Py(),track->Pz(),cluster->E());
   	cout << "Added cluster with E: " << cluster->E() << endl;
      }
    }
  }
  
  fFastJetWrapper->SetupStrategyfromOpt("Best");
  fFastJetWrapper->SetupAlgorithmfromOpt("antikt");
  fFastJetWrapper->SetupSchemefromOpt("BIpt");
  fFastJetWrapper->SetupAreaTypefromOpt("active");
  fFastJetWrapper->SetNRepeats(1);
  fFastJetWrapper->SetGhostArea(0.01);
  fFastJetWrapper->SetMaxRap(0.9);
  fFastJetWrapper->SetR(0.4);
  
  double median = 0.0;
  double sigma  = 0.0;

  fFastJetWrapper->Run();
  fFastJetWrapper->GetMedianAndSigma(median,sigma);

  if ( TriggerOnJet(fFastJetWrapper->GetSubtractedJetsPts(median)) ) {
    return iResult;
  }
  
  TString description;
  description.Form(" No jets with energy >  %.02f GeV found! ", fDetector.Data(), fEThreshold);
  SetDescription(description.Data());
  TriggerEvent(kFALSE);
  
  return iResult;
}

//______________________________________________________________

template <class T>
Bool_t AliHLTTriggerFastJet::TriggerOnJet(T Jet) {

  for ( unsigned int ij = 0; ij<Jet.size(); ij++ ) {
    if ( Jet[ij] > fEThreshold ) {
      TString description;
      description.Form(" Event contains at least one %s jet with energy > %.02f GeV! ", fDetector.Data(), fEThreshold);
      SetDescription(description.Data());
      
      // Enable the detectors for readout
      GetReadoutList().Enable(AliHLTReadoutList::kEMCAL |
			      AliHLTReadoutList::kTPC);
      
      // Add the available HLT info for readout
      GetTriggerDomain().Add(kAliHLTAnyDataTypeID, fDetector.Data());
      
      // Set trigger desicion
      TriggerEvent(kTRUE);
      
      return(kTRUE);
    }
  }
  
  return kFALSE;

}
//______________________________________________________________

int AliHLTTriggerFastJet::Reconfigure(const char* cdbEntry, const char* /*chainId*/) {
  
  const char* entry = cdbEntry;
  if ( !entry || entry[0]==0 ) entry=fOCDBEntry;

  return ConfigureFromCDBTObjString(entry);

}

//______________________________________________________________

int AliHLTTriggerFastJet::ScanConfigurationArgument(int argc, const char** argv) {

  if ( argc<=0 ) return 0;
  int i = 0;
  TString argument = argv[i];

  if ( argument.CompareTo("-energy") == 0 ) {
    if (++i>=argc) return -EPROTO;
    argument = argv[i];
    fEThreshold=argument.Atof();
    return 2;
  }
  
  return -EINVAL;

}
//______________________________________________________________

void AliHLTTriggerFastJet::GetOutputDataSize(unsigned long &constBase, double &inputMultiplier) {
  
  constBase = sizeof(AliHLTTriggerDecision) + sizeof(AliHLTDomainEntry)*14;
  inputMultiplier = 1;

}
//______________________________________________________________

void AliHLTTriggerFastJet::GetOCDBObjectDescription(TMap* const targetMap) {
  
  if ( !targetMap ) return;
  targetMap->Add(new TObjString(fOCDBEntry),
		 new TObjString(Form("%s threshold trigger OCDB object", fDetector.Data()) ) 
		 );

}
//_______________________________________________________________
Commit	Line	Data
7634d323	1	#include "AliHLTTriggerFastJet.h"
	2	#include "AliESDEvent.h"
	3	#include "AliESDtrackCuts.h"
	4	#include "AliESDCaloCluster.h"
917b0861	5	#include "AliHLTCaloClusterReader.h"
917b0861	6	#include "AliHLTCaloClusterDataStruct.h"
7634d323	7	#include "AliESDVZERO.h"
	8	#include "AliHLTTriggerDecision.h"
	9	#include "AliHLTDomainEntry.h"
	10	#include "TLorentzVector.h"
	11	#include "TRefArray.h"
	12	#include "TString.h"
	13	#include "TMap.h"
917b0861	14	#include "TObjArray.h"
917b0861	15	#include "TArrayI.h"
7634d323	16
7634d323	17	#include "AliFJWrapper.h"
	18
	19
	20	ClassImp(AliHLTTriggerFastJet)
	21
917b0861	22	AliHLTTriggerFastJet::AliHLTTriggerFastJet() :
7634d323	23	AliHLTTrigger(),
7634d323	24	fEThreshold(0.0),
917b0861	25	fDetector("EMCAL"),
7634d323	26	fFastJetWrapper(NULL),
917b0861	27	EsdTrackCuts(NULL),
917b0861	28	fOCDBEntry("HLT/ConfigHLT/EmcalJetTrigger")
7634d323	29	{
	30	fFastJetWrapper = new AliFJWrapper("FastJet","FastJet");
	31
7634d323	32	}
	33	//_____________________________________________________________
	34
	35	AliHLTTriggerFastJet::~AliHLTTriggerFastJet() {
	36
	37	if (fFastJetWrapper)
	38	delete fFastJetWrapper;
	39	fFastJetWrapper = NULL;
	40
917b0861	41	if (EsdTrackCuts)
	42	delete EsdTrackCuts;
	43	EsdTrackCuts = NULL;
	44
7634d323	45	}
	46	//_____________________________________________________________
	47
	48	int AliHLTTriggerFastJet::DoInit(int argc, const char** argv) {
	49
	50	int iResult = ConfigureFromCDBTObjString(fOCDBEntry);
	51
	52	if ( iResult>=0 && argc>0 ) {
	53	iResult = ConfigureFromArgumentString(argc, argv);
	54	HLTImportant("Trigger threshold set from argument string : %.02f GeV:", fEThreshold );
	55	} else if ( iResult>=0 ) {
	56	HLTImportant("Trigger threshold set from OCDB database entry : %.02f Gev:", fEThreshold );
	57	}
	58	return iResult;
	59	}
	60	//______________________________________________________________
	61
	62	int AliHLTTriggerFastJet::DoDeInit() {
	63
	64	return 0;
	65
	66	}
	67	//______________________________________________________________
	68
7eefc52d	69
	70
	71	AliHLTComponent* AliHLTTriggerFastJet::Spawn() {
	72	// see header file for class documentation
	73	return new AliHLTTriggerFastJet;
	74	}
	75
	76
	77	const char* AliHLTTriggerFastJet::GetTriggerName() const {
	78	// see header file for class documentation
	79	return "EmcalJetTrigger";
	80	}
	81
7634d323	82	Int_t AliHLTTriggerFastJet::DoTrigger() {
	83
	84	Int_t iResult = 0;
	85
	86	if ( GetFirstInputBlock( kAliHLTDataTypeSOR ) \|\| GetFirstInputBlock( kAliHLTDataTypeEOR ) )
	87	return 0;
	88
	89	const TObject *obj = GetFirstInputObject( kAliHLTAllDataTypes, "AliESDEvent" );
	90	AliESDEvent esd = dynamic_cast<AliESDEvent>(const_cast<TObject*>(obj));
	91
	92	// check for MatchedTrack to avoid double counting
	93
	94
	95	if ( esd != NULL ) {
	96	esd->GetStdContent();
	97
917b0861	98	//Double_t vertex[3] = {0,0,0};
	99	//esd->GetVertex()->GetXYZ(vertex);
	100	//cout << "Vertex " << vertex[0] << vertex[1] << vertex [2] << endl;
	101	//TLorentzVector gamma;
7634d323	102
7634d323	103	// -- add VZERO
7eefc52d	104	// uncomment for usage (to avoid warning)
	105
	106	/*
7634d323	107	AliESDVZERO *v0 = esd->GetVZEROData();
7eefc52d	108
7634d323	109	Float_t MtotVOA = v0->GetMTotV0A();
7634d323	110	Float_t MtotVOC = v0->GetMTotV0A();
7eefc52d	111	*/
917b0861	112	Int_t nTracks = esd->GetNumberOfTracks();
	113
	114	TObjArray *tracks = new TObjArray(nTracks);
	115	tracks->SetOwner(1);
	116
	117
	118	EsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
	119	EsdTrackCuts->SetMinNClustersTPC(70);
	120
	121	// add the tracks to FastJet Wrapper
	122
	123	for ( Int_t i = 0; i < nTracks; i++ ) {
	124
	125	const AliESDVertex *vtxTPC = esd->GetPrimaryVertexTPC();
	126	if (!vtxTPC) continue;
	127	cout << "Got primary vertex " << endl;
7634d323	128
7634d323	129	AliESDtrack *esdtrack = esd->GetTrack(i);
917b0861	130	if (!esdtrack) continue;
	131	tracks->Add(esdtrack);
	132	cout << "Got esd track " << endl;
	133
	134	AliESDtrack tpctrack = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent>(esd),esdtrack->GetID());
	135	if (!tpctrack) continue;
	136	cout << "Got tpc track " << endl;
	137
	138	AliExternalTrackParam exParam;
	139	Bool_t relate = tpctrack->RelateToVertexTPC(vtxTPC,esd->GetMagneticField(),kVeryBig,&exParam);
	140	if (!relate) {
	141	delete tpctrack;
	142	continue;
	143	}
	144
	145	tpctrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
	146	cout << "Set param to tpctrack " << endl;
	147
	148	fFastJetWrapper->AddInputVector(tpctrack->Px(),tpctrack->Py(),tpctrack->Pz(),tpctrack->P());
	149	cout << "Added track with P: " << tpctrack->P() << endl;
7634d323	150	}
917b0861	151
	152
	153	// // add the Calo Clusters to FastJet Wrapper
	154	// TRefArray *caloClustersArr = new TRefArray();
	155	// esd->GetEMCALClusters(caloClustersArr);
	156	// const Int_t nEMCALClusters = caloClustersArr->GetEntries();
	157	//
	158	// for ( UInt_t ic = 0; ic < nEMCALClusters; ic++ )
	159	// {
	160	// AliESDCaloCluster cluster = dynamic_cast<AliESDCaloCluster>( caloClustersArr->At(ic) );
	161	// Double_t clusterEn = cluster->E();
	162	// Double_t subE = 0;
	163	//
	164	// TArrayI *matched = cluster->GetTracksMatched();
	165	// for ( UInt_t jk = 0; jk < matched->GetSize(); jk++ )
	166	// {
	167	// if ( matched->At(jk) > -1 && matched->At(jk) < esd ->GetNumberOfTracks() )
	168	// {
	169	// AliESDtrack track = dynamic_cast<AliESDtrack>( tracks->At( (matched->At(jk)) ) );
	170	// if ( !track ) continue;
	171	// subE += track->P();
	172	// }
	173	// }
	174	// clusterEn -= subE;
	175	// if ( clusterEn < 0)
	176	// clusterEn = 0;
	177	// cluster->SetE( clusterEn );
	178	// if ( cluster->E() == 0 ) continue;
	179	//
	180	// cluster->GetMomentum( gamma, vertex );
	181	// fFastJetWrapper->AddInputVector( gamma.Px(), gamma.Py(), gamma.Pz(), clusterEn );
	182	// gamma.Clear();
	183	// }
	184	// }
	185	for ( Int_t j = 0; j< esd->GetNumberOfCaloClusters(); j++) {
	186	//AliESDCaloCluster cluster = static_cast<AliESDCaloCluster>(esd->GetCaloCluster(j));
	187	AliVCluster cluster = dynamic_cast<AliVCluster>(esd->GetCaloCluster(j));
7634d323	188	if ( cluster->IsEMCAL() ) {
917b0861	189	Int_t trackindex = cluster->GetTrackMatchedIndex();
	190	if (trackindex<0) continue;
	191	AliESDtrack *track = esd->GetTrack(trackindex);
	192	if (!track) continue;
	193	fFastJetWrapper->AddInputVector(track->Px(),track->Py(),track->Pz(),cluster->E());
	194	cout << "Added cluster with E: " << cluster->E() << endl;
7634d323	195	}
7634d323	196	}
7634d323	197	}
7634d323	198
ab9926fa	199	fFastJetWrapper->SetupStrategyfromOpt("Best");
	200	fFastJetWrapper->SetupAlgorithmfromOpt("antikt");
	201	fFastJetWrapper->SetupSchemefromOpt("BIpt");
	202	fFastJetWrapper->SetupAreaTypefromOpt("active");
7634d323	203	fFastJetWrapper->SetNRepeats(1);
	204	fFastJetWrapper->SetGhostArea(0.01);
	205	fFastJetWrapper->SetMaxRap(0.9);
	206	fFastJetWrapper->SetR(0.4);
	207
	208	double median = 0.0;
	209	double sigma = 0.0;
	210
	211	fFastJetWrapper->Run();
	212	fFastJetWrapper->GetMedianAndSigma(median,sigma);
	213
917b0861	214	if ( TriggerOnJet(fFastJetWrapper->GetSubtractedJetsPts(median)) ) {
	215	return iResult;
	216	}
7634d323	217
7634d323	218	TString description;
917b0861	219	description.Form(" No jets with energy > %.02f GeV found! ", fDetector.Data(), fEThreshold);
7634d323	220	SetDescription(description.Data());
	221	TriggerEvent(kFALSE);
	222
	223	return iResult;
7634d323	224	}
917b0861	225
7634d323	226	//______________________________________________________________
	227
	228	template <class T>
	229	Bool_t AliHLTTriggerFastJet::TriggerOnJet(T Jet) {
	230
	231	for ( unsigned int ij = 0; ij<Jet.size(); ij++ ) {
	232	if ( Jet[ij] > fEThreshold ) {
	233	TString description;
	234	description.Form(" Event contains at least one %s jet with energy > %.02f GeV! ", fDetector.Data(), fEThreshold);
	235	SetDescription(description.Data());
	236
	237	// Enable the detectors for readout
	238	GetReadoutList().Enable(AliHLTReadoutList::kEMCAL \|
	239	AliHLTReadoutList::kTPC);
	240
	241	// Add the available HLT info for readout
	242	GetTriggerDomain().Add(kAliHLTAnyDataTypeID, fDetector.Data());
	243
	244	// Set trigger desicion
	245	TriggerEvent(kTRUE);
	246
	247	return(kTRUE);
	248	}
	249	}
	250
	251	return kFALSE;
	252
	253	}
	254	//______________________________________________________________
	255
	256	int AliHLTTriggerFastJet::Reconfigure(const char* cdbEntry, const char* /chainId/) {
	257
	258	const char* entry = cdbEntry;
	259	if ( !entry \|\| entry[0]==0 ) entry=fOCDBEntry;
	260
	261	return ConfigureFromCDBTObjString(entry);
	262
	263	}
	264
	265	//______________________________________________________________
	266
	267	int AliHLTTriggerFastJet::ScanConfigurationArgument(int argc, const char** argv) {
	268
	269	if ( argc<=0 ) return 0;
	270	int i = 0;
	271	TString argument = argv[i];
	272
	273	if ( argument.CompareTo("-energy") == 0 ) {
	274	if (++i>=argc) return -EPROTO;
	275	argument = argv[i];
	276	fEThreshold=argument.Atof();
	277	return 2;
	278	}
	279
	280	return -EINVAL;
	281
	282	}
	283	//______________________________________________________________
	284
	285	void AliHLTTriggerFastJet::GetOutputDataSize(unsigned long &constBase, double &inputMultiplier) {
	286
	287	constBase = sizeof(AliHLTTriggerDecision) + sizeof(AliHLTDomainEntry)*14;
	288	inputMultiplier = 1;
	289
290	}
291	//______________________________________________________________
292
293	void AliHLTTriggerFastJet::GetOCDBObjectDescription(TMap* const targetMap) {
294
295	if ( !targetMap ) return;
296	targetMap->Add(new TObjString(fOCDBEntry),
297	new TObjString(Form("%s threshold trigger OCDB object", fDetector.Data()) )
298	);
299
300	}
301	//_______________________________________________________________