[u/mrichter/AliRoot.git] / PWG / EMCAL / AliEmcalEsdTrackFilterTask.cxx

// $Id$
//
// Task to filter Esd tracks and propagate to Emcal surface.
//
// Author: C.Loizides

#include "AliEmcalEsdTrackFilterTask.h"
#include <TClonesArray.h>
#include <TRandom3.h>
#include <TGeoGlobalMagField.h>
#include <AliAnalysisManager.h>
#include <AliEMCALRecoUtils.h>
#include <AliESDEvent.h>
#include <AliESDtrackCuts.h>
#include <AliMagF.h>
#include <AliTrackerBase.h>


ClassImp(AliEmcalEsdTrackFilterTask)

//________________________________________________________________________
AliEmcalEsdTrackFilterTask::AliEmcalEsdTrackFilterTask() : 
  AliAnalysisTaskSE("AliEmcalEsdTrackFilterTask"),
  fEsdTrackCuts(0),
  fDoSpdVtxCon(0),
  fHybridTrackCuts(0),
  fTracksName(),
  fIncludeNoITS(kTRUE),
  fDoPropagation(kFALSE),
  fDist(440),
  fTrackEfficiency(1),
  fEsdEv(0),
  fTracks(0)
{
  // Constructor.
}

//________________________________________________________________________
AliEmcalEsdTrackFilterTask::AliEmcalEsdTrackFilterTask(const char *name) : 
  AliAnalysisTaskSE(name),
  fEsdTrackCuts(0),
  fDoSpdVtxCon(0),
  fHybridTrackCuts(0),
  fTracksName("EsdTracksOut"),
  fIncludeNoITS(kTRUE),
  fDoPropagation(kFALSE),
  fDist(440),
  fTrackEfficiency(1),
  fEsdEv(0),
  fTracks(0)
{
  // Constructor.

  if (!name)
    return;

  SetName(name);

  fBranchNames = "ESD:AliESDHeader.,AliESDRun.,SPDVertex.,Tracks";
}

//________________________________________________________________________
AliEmcalEsdTrackFilterTask::~AliEmcalEsdTrackFilterTask()
{
  //Destructor

  delete fEsdTrackCuts;
}

//________________________________________________________________________
void AliEmcalEsdTrackFilterTask::UserCreateOutputObjects()
{
  // Create histograms.

  fTracks = new TClonesArray("AliESDtrack");
  fTracks->SetName(fTracksName);

  if (fDoSpdVtxCon) {
    if (!fEsdTrackCuts) {
      AliInfo("No track cuts given, creating default (standard only TPC) cuts");
      fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
      fEsdTrackCuts->SetPtRange(0.15,1e3);
    } 
  } else {
    AliWarning("No track cuts given, but maybe this is indeed intended?");
  }
}

//________________________________________________________________________
void AliEmcalEsdTrackFilterTask::UserExec(Option_t *) 
{
  // Main loop, called for each event.

  fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
  if (!fEsdEv) {
    AliError("Task works only on ESD events, returning");
    return;
  }

  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
  if (!am) {
    AliError("Manager zero, returning");
    return;
  }

  // add tracks to event if not yet there
  fTracks->Delete();
  if (!(InputEvent()->FindListObject(fTracksName)))
    InputEvent()->AddObject(fTracks);

  if (!fHybridTrackCuts) { // constrain TPC tracks to SPD vertex if fDoSpdVtxCon==kTRUE
    am->LoadBranch("AliESDRun.");
    am->LoadBranch("AliESDHeader.");
    am->LoadBranch("Tracks");

    if (fDoSpdVtxCon) {
      if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
        fEsdEv->InitMagneticField();
      }
      am->LoadBranch("SPDVertex.");
      const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
      if (!vtxSPD) {
        AliError("No SPD vertex, returning");
        return;
      }
      Int_t ntr = fEsdEv->GetNumberOfTracks();
      for (Int_t i=0, ntrnew=0; i<ntr; ++i) {
        AliESDtrack *etrack = fEsdEv->GetTrack(i);
        if (!etrack)
          continue;

	if (fTrackEfficiency < 1) {
	  Double_t r = gRandom->Rndm();
	  if (fTrackEfficiency < r)
	    continue;
	}

        if (!fEsdTrackCuts->AcceptTrack(etrack))
          continue;

        AliESDtrack *ntrack = AliESDtrackCuts::GetTPCOnlyTrack(fEsdEv,etrack->GetID());
        if (!ntrack)
          continue;
        if (ntrack->Pt()<=0) {
          delete ntrack;
          continue;
        }
        Double_t bfield[3] = {0,0,0};
        ntrack->GetBxByBz(bfield);
        AliExternalTrackParam exParam;
        Bool_t relate = ntrack->RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&exParam);
        if (!relate) {
          delete ntrack;
          continue;
        }
        // set the constraint parameters to the track
        ntrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
        if (ntrack->Pt()<=0) {
          delete ntrack;
          continue;
        }
	if (fDoPropagation) 	
	  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(ntrack,fDist);
        new ((*fTracks)[ntrnew++]) AliESDtrack(*ntrack);
        delete ntrack;
      }
    } else { /* no spd vtx constraint */
      Int_t ntr = fEsdEv->GetNumberOfTracks();
      for (Int_t i=0, ntrnew=0; i<ntr; ++i) {
        AliESDtrack *etrack = fEsdEv->GetTrack(i);
        if (!etrack)
          continue;
	if (fTrackEfficiency < 1) {
	  Double_t r = gRandom->Rndm();
	  if (fTrackEfficiency < r)
	    continue;
	}
	if ((fEsdTrackCuts!=0) && !fEsdTrackCuts->AcceptTrack(etrack))
          continue;
        AliESDtrack *ntrack = new ((*fTracks)[ntrnew++]) AliESDtrack(*etrack);
	if (fDoPropagation) 	
	  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(ntrack,fDist);
      }
    }

  } else { // use hybrid track cuts

    am->LoadBranch("Tracks");
    Int_t ntr = fEsdEv->GetNumberOfTracks();
    for (Int_t i=0, ntrnew=0; i<ntr; ++i) {
      AliESDtrack *etrack = fEsdEv->GetTrack(i);
      if (!etrack) 
	continue;

      if (fTrackEfficiency < 1) {
	Double_t r = gRandom->Rndm();
	if (fTrackEfficiency < r)
	  continue;
      }

      if (fEsdTrackCuts->AcceptTrack(etrack)) {
        AliESDtrack *newTrack = new ((*fTracks)[ntrnew]) AliESDtrack(*etrack);
	if (fDoPropagation) 
	  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(newTrack,fDist);
        newTrack->SetBit(BIT(22),0); 
        newTrack->SetBit(BIT(23),0);
        ++ntrnew;
      } else if (fHybridTrackCuts->AcceptTrack(etrack)) {
	if (!etrack->GetConstrainedParam())
	  continue;
	UInt_t status = etrack->GetStatus();
	if (!fIncludeNoITS && ((status&AliESDtrack::kITSrefit)==0))
	  continue;
	AliESDtrack *newTrack = new ((*fTracks)[ntrnew]) AliESDtrack(*etrack);
	if (fDoPropagation) 	
	  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(newTrack,fDist);
	const AliExternalTrackParam* constrainParam = etrack->GetConstrainedParam();
	newTrack->Set(constrainParam->GetX(),
		      constrainParam->GetAlpha(),
		      constrainParam->GetParameter(),
		      constrainParam->GetCovariance());
	if ((status&AliESDtrack::kITSrefit)==0) {
	  newTrack->SetBit(BIT(22),0); //type 2
	  newTrack->SetBit(BIT(23),1);
	} else {
	  newTrack->SetBit(BIT(22),1); //type 1
	  newTrack->SetBit(BIT(23),0);
	}
	++ntrnew;
      }
    }
  }
}
Commit	Line	Data
cc2b5fb5	1	// $Id$
b0a53615	2	//
cdd4e0e2	3	// Task to filter Esd tracks and propagate to Emcal surface.
b0a53615	4	//
cd231d42	5	// Author: C.Loizides
cc2b5fb5	6
33a8056d	7	#include "AliEmcalEsdTrackFilterTask.h"
cc2b5fb5	8	#include <TClonesArray.h>
d47bb504	9	#include <TRandom3.h>
cc2b5fb5	10	#include <TGeoGlobalMagField.h>
fe75ca27	11	#include <AliAnalysisManager.h>
33a8056d	12	#include <AliEMCALRecoUtils.h>
fe75ca27 CL	13	#include <AliESDEvent.h>
	14	#include <AliESDtrackCuts.h>
	15	#include <AliMagF.h>
	16	#include <AliTrackerBase.h>
cc2b5fb5	17
d47bb504	18
cdd4e0e2	19	ClassImp(AliEmcalEsdTrackFilterTask)
cc2b5fb5	20
cc2b5fb5	21	//________________________________________________________________________
cdd4e0e2 CL	22	AliEmcalEsdTrackFilterTask::AliEmcalEsdTrackFilterTask() :
cdd4e0e2 CL	23	AliAnalysisTaskSE("AliEmcalEsdTrackFilterTask"),
cc2b5fb5	24	fEsdTrackCuts(0),
4d1e3842	25	fDoSpdVtxCon(0),
57184cc6	26	fHybridTrackCuts(0),
cc2b5fb5	27	fTracksName(),
6231f841	28	fIncludeNoITS(kTRUE),
35fb503f	29	fDoPropagation(kFALSE),
af2f024d	30	fDist(440),
d47bb504	31	fTrackEfficiency(1),
cc2b5fb5	32	fEsdEv(0),
	33	fTracks(0)
	34	{
a93ef37c	35	// Constructor.
cc2b5fb5	36	}
	37
	38	//________________________________________________________________________
cdd4e0e2	39	AliEmcalEsdTrackFilterTask::AliEmcalEsdTrackFilterTask(const char *name) :
cc2b5fb5	40	AliAnalysisTaskSE(name),
cc2b5fb5	41	fEsdTrackCuts(0),
4d1e3842	42	fDoSpdVtxCon(0),
57184cc6	43	fHybridTrackCuts(0),
cdd4e0e2	44	fTracksName("EsdTracksOut"),
6231f841	45	fIncludeNoITS(kTRUE),
35fb503f	46	fDoPropagation(kFALSE),
af2f024d	47	fDist(440),
d47bb504	48	fTrackEfficiency(1),
cc2b5fb5	49	fEsdEv(0),
	50	fTracks(0)
	51	{
a93ef37c	52	// Constructor.
27fde6c7	53
	54	if (!name)
	55	return;
	56
	57	SetName(name);
	58
a93ef37c	59	fBranchNames = "ESD:AliESDHeader.,AliESDRun.,SPDVertex.,Tracks";
cc2b5fb5	60	}
	61
	62	//________________________________________________________________________
cdd4e0e2	63	AliEmcalEsdTrackFilterTask::~AliEmcalEsdTrackFilterTask()
cc2b5fb5	64	{
	65	//Destructor
	66
	67	delete fEsdTrackCuts;
	68	}
	69
	70	//________________________________________________________________________
cdd4e0e2	71	void AliEmcalEsdTrackFilterTask::UserCreateOutputObjects()
cc2b5fb5	72	{
a93ef37c	73	// Create histograms.
cc2b5fb5	74
	75	fTracks = new TClonesArray("AliESDtrack");
	76	fTracks->SetName(fTracksName);
	77
04410c6f	78	if (fDoSpdVtxCon) {
	79	if (!fEsdTrackCuts) {
	80	AliInfo("No track cuts given, creating default (standard only TPC) cuts");
	81	fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
	82	fEsdTrackCuts->SetPtRange(0.15,1e3);
	83	}
	84	} else {
	85	AliWarning("No track cuts given, but maybe this is indeed intended?");
cc2b5fb5	86	}
	87	}
	88
	89	//________________________________________________________________________
cdd4e0e2	90	void AliEmcalEsdTrackFilterTask::UserExec(Option_t *)
cc2b5fb5	91	{
	92	// Main loop, called for each event.
	93
	94	fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
	95	if (!fEsdEv) {
	96	AliError("Task works only on ESD events, returning");
	97	return;
	98	}
	99
	100	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
	101	if (!am) {
	102	AliError("Manager zero, returning");
	103	return;
	104	}
	105
	106	// add tracks to event if not yet there
71c1dd76	107	fTracks->Delete();
cc2b5fb5	108	if (!(InputEvent()->FindListObject(fTracksName)))
	109	InputEvent()->AddObject(fTracks);
	110
ce2a8490	111	if (!fHybridTrackCuts) { // constrain TPC tracks to SPD vertex if fDoSpdVtxCon==kTRUE
57184cc6	112	am->LoadBranch("AliESDRun.");
57184cc6	113	am->LoadBranch("AliESDHeader.");
57184cc6	114	am->LoadBranch("Tracks");
4d1e3842	115
	116	if (fDoSpdVtxCon) {
	117	if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
	118	fEsdEv->InitMagneticField();
57184cc6	119	}
4d1e3842	120	am->LoadBranch("SPDVertex.");
	121	const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
	122	if (!vtxSPD) {
	123	AliError("No SPD vertex, returning");
	124	return;
57184cc6	125	}
4d1e3842	126	Int_t ntr = fEsdEv->GetNumberOfTracks();
	127	for (Int_t i=0, ntrnew=0; i<ntr; ++i) {
	128	AliESDtrack *etrack = fEsdEv->GetTrack(i);
	129	if (!etrack)
	130	continue;
d47bb504	131
	132	if (fTrackEfficiency < 1) {
	133	Double_t r = gRandom->Rndm();
	134	if (fTrackEfficiency < r)
	135	continue;
	136	}
	137
4d1e3842	138	if (!fEsdTrackCuts->AcceptTrack(etrack))
4d1e3842	139	continue;
d47bb504	140
4d1e3842	141	AliESDtrack *ntrack = AliESDtrackCuts::GetTPCOnlyTrack(fEsdEv,etrack->GetID());
	142	if (!ntrack)
	143	continue;
	144	if (ntrack->Pt()<=0) {
	145	delete ntrack;
	146	continue;
	147	}
	148	Double_t bfield[3] = {0,0,0};
	149	ntrack->GetBxByBz(bfield);
	150	AliExternalTrackParam exParam;
	151	Bool_t relate = ntrack->RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&exParam);
	152	if (!relate) {
	153	delete ntrack;
	154	continue;
	155	}
	156	// set the constraint parameters to the track
	157	ntrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
	158	if (ntrack->Pt()<=0) {
	159	delete ntrack;
	160	continue;
	161	}
04410c6f	162	if (fDoPropagation)
04410c6f	163	AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(ntrack,fDist);
4d1e3842	164	new ((fTracks)[ntrnew++]) AliESDtrack(ntrack);
57184cc6	165	delete ntrack;
57184cc6	166	}
4d1e3842	167	} else { /* no spd vtx constraint */
	168	Int_t ntr = fEsdEv->GetNumberOfTracks();
	169	for (Int_t i=0, ntrnew=0; i<ntr; ++i) {
	170	AliESDtrack *etrack = fEsdEv->GetTrack(i);
	171	if (!etrack)
	172	continue;
d47bb504	173	if (fTrackEfficiency < 1) {
	174	Double_t r = gRandom->Rndm();
	175	if (fTrackEfficiency < r)
	176	continue;
	177	}
04410c6f	178	if ((fEsdTrackCuts!=0) && !fEsdTrackCuts->AcceptTrack(etrack))
4d1e3842	179	continue;
04410c6f	180	AliESDtrack ntrack = new ((fTracks)[ntrnew++]) AliESDtrack(*etrack);
	181	if (fDoPropagation)
	182	AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(ntrack,fDist);
4d1e3842	183	}
cc2b5fb5	184	}
cc2b5fb5	185
10d1f756	186	} else { // use hybrid track cuts
57184cc6	187
	188	am->LoadBranch("Tracks");
	189	Int_t ntr = fEsdEv->GetNumberOfTracks();
	190	for (Int_t i=0, ntrnew=0; i<ntr; ++i) {
	191	AliESDtrack *etrack = fEsdEv->GetTrack(i);
35fb503f CL	192	if (!etrack)
35fb503f CL	193	continue;
d47bb504	194
	195	if (fTrackEfficiency < 1) {
	196	Double_t r = gRandom->Rndm();
	197	if (fTrackEfficiency < r)
	198	continue;
	199	}
	200
57184cc6	201	if (fEsdTrackCuts->AcceptTrack(etrack)) {
b3bef4aa	202	AliESDtrack newTrack = new ((fTracks)[ntrnew]) AliESDtrack(*etrack);
cdd4e0e2 CL	203	if (fDoPropagation)
cdd4e0e2 CL	204	AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(newTrack,fDist);
32c97abe	205	newTrack->SetBit(BIT(22),0);
32c97abe	206	newTrack->SetBit(BIT(23),0);
57184cc6	207	++ntrnew;
57184cc6	208	} else if (fHybridTrackCuts->AcceptTrack(etrack)) {
3f4c638b	209	if (!etrack->GetConstrainedParam())
3f4c638b	210	continue;
6231f841	211	UInt_t status = etrack->GetStatus();
3f4c638b	212	if (!fIncludeNoITS && ((status&AliESDtrack::kITSrefit)==0))
	213	continue;
	214	AliESDtrack newTrack = new ((fTracks)[ntrnew]) AliESDtrack(*etrack);
	215	if (fDoPropagation)
	216	AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(newTrack,fDist);
	217	const AliExternalTrackParam* constrainParam = etrack->GetConstrainedParam();
	218	newTrack->Set(constrainParam->GetX(),
	219	constrainParam->GetAlpha(),
	220	constrainParam->GetParameter(),
	221	constrainParam->GetCovariance());
	222	if ((status&AliESDtrack::kITSrefit)==0) {
	223	newTrack->SetBit(BIT(22),0); //type 2
	224	newTrack->SetBit(BIT(23),1);
	225	} else {
	226	newTrack->SetBit(BIT(22),1); //type 1
	227	newTrack->SetBit(BIT(23),0);
af2f024d	228	}
3f4c638b	229	++ntrnew;
57184cc6	230	}
57184cc6	231	}
cc2b5fb5	232	}
cc2b5fb5	233	}