[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisEtSelectorEmcal.cxx

//_________________________________________________________________________
//  Utility Class for transverse energy studies
//  Selection class for EMCAL
//
//*-- Authors: Oystein Djuvsland (Bergen)
//_________________________________________________________________________


#include "AliAnalysisEtSelectorEmcal.h"
#include "AliAnalysisEtCuts.h"
#include "AliEMCALTrack.h"
#include "TParticle.h"
#include <iostream>
#include "AliStack.h"

AliAnalysisEtSelectorEmcal::AliAnalysisEtSelectorEmcal(AliAnalysisEtCuts* cuts):AliAnalysisEtSelector(cuts)
{

}
AliAnalysisEtSelectorEmcal::AliAnalysisEtSelectorEmcal():AliAnalysisEtSelector()
{
}

AliAnalysisEtSelectorEmcal::~AliAnalysisEtSelectorEmcal()
{

}

void AliAnalysisEtSelectorEmcal::Init()
{
    AliAnalysisEtSelector::Init();
}

Int_t AliAnalysisEtSelectorEmcal::Init(const AliESDEvent* event)
{   // Init

    AliAnalysisEtSelector::Init(event);
    Printf("Initializing selector for run: %d", event->GetRunNumber());
    fInitialized = kTRUE;
    return 0;
}

TRefArray* AliAnalysisEtSelectorEmcal::GetClusters()
{   // Get clusters

    if(!fClusterArray) fClusterArray = new TRefArray;

    if(fClusterArray)
    {
        fEvent->GetEMCALClusters(fClusterArray);
    }
    else
    {
        Printf("Could not initialize cluster array");
    }
    return fClusterArray;
}

Bool_t AliAnalysisEtSelectorEmcal::PassMinEnergyCut(const AliESDCaloCluster& cl) const
{
  Float_t pos[3];
  cl.GetPosition(pos);
  TVector3 cp(pos);
  return TMath::Sin(cp.Theta())*cl.E() > fCuts->GetReconstructedEmcalClusterEnergyCut();
}

Bool_t AliAnalysisEtSelectorEmcal::PassMinEnergyCut(const TParticle& p) const
{
  return TMath::Sin(p.Theta())*p.Energy() > fCuts->GetReconstructedEmcalClusterEnergyCut();
}
Bool_t AliAnalysisEtSelectorEmcal::PassMinEnergyCut(Double_t e) const
{
  return e > fCuts->GetReconstructedEmcalClusterEnergyCut();
}

Bool_t AliAnalysisEtSelectorEmcal::PassDistanceToBadChannelCut(const AliESDCaloCluster& ) const
{
    return kTRUE;
}

Bool_t AliAnalysisEtSelectorEmcal::PassTrackMatchingCut(const AliESDCaloCluster& cluster) const
{

  Int_t nTracksMatched = cluster.GetNTracksMatched();
  if(nTracksMatched == 0){
    return kTRUE;
  }
  
  Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
  if(trackMatchedIndex < 0){
    return kTRUE;
  }
  AliVParticle *track = fEvent->GetTrack(trackMatchedIndex);
  if(track->Pt()<0.5) return kTRUE;//Track matches below about 500 MeV are mostly random.  It takes ~460 MeV to reach the EMCal

  //Float_t recoE = cluster.E();
  //Float_t pos[3];

//     //cluster.GetPosition(pos);
//     Int_t trackMatchIdx = cluster.GetTrackMatchedIndex();
//     //Double_t distance = 9999.0;
//     if(trackMatchIdx>-1)
//     {
//       return kTRUE;
//       //distance = CalcTrackClusterDistance(pos, &trackMatchIdx);
//     }


    return kFALSE;
    //return distance > fCuts->GetEmcalTrackDistanceCut();
}

Bool_t AliAnalysisEtSelectorEmcal::CutGeometricalAcceptance(const TParticle& part)
{
  float myphi = part.Phi();
  myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
    return TMath::Abs(part.Eta()) < fCuts->GetGeometryEmcalEtaAccCut()
    && myphi < fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180.
    && myphi > fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.;
}

Bool_t AliAnalysisEtSelectorEmcal::CutGeometricalAcceptance(const AliVTrack& part)
{
  float myphi = part.Phi();
  myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
  return TMath::Abs(part.Eta()) < fCuts->GetGeometryEmcalEtaAccCut()
    && myphi < fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180.
    && myphi > fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.;
}

Bool_t AliAnalysisEtSelectorEmcal::CutGeometricalAcceptance(const AliESDCaloCluster& cluster)
{
  Float_t pos[3];
  cluster.GetPosition(pos);
  TVector3 cp(pos);
  float myphi = cp.Phi();
  myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
  return TMath::Abs(cp.Eta()) < fCuts->GetGeometryEmcalEtaAccCut()
    && myphi < fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180.
    && myphi > fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.;
}


Double_t
AliAnalysisEtSelectorEmcal::CalcTrackClusterDistance(const Float_t clsPos[3],Int_t *trkMatchId) const
{   // calculate distance between cluster and closest track

    Double_t trkPos[3] = {0,0,0};

    Int_t bestTrkMatchId = -1;
    Double_t distance = 9999; // init to a big number

    Double_t dist = 0;
    Double_t distX = 0, distY = 0, distZ = 0;

    for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
        AliESDtrack *track = static_cast<AliESDtrack*>( fEvent->GetTrack(iTrack) );
        if (!track) {
            AliError(Form("ERROR: Could not get track %d", iTrack));
            continue;
        }

        // check for approx. eta and phi range before we propagate..
        // TBD

        AliEMCALTrack emctrack(*track);
        if (!emctrack.PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) {
            continue;
        }
        emctrack.GetXYZ(trkPos);

        distX = clsPos[0]-trkPos[0];
        distY = clsPos[1]-trkPos[1];
        distZ = clsPos[2]-trkPos[2];
        dist = TMath::Sqrt(distX*distX + distY*distY + distZ*distZ);

        if (dist < distance) {
            distance = dist;
            bestTrkMatchId = iTrack;
        }
    } // iTrack

    // printf("CalcTrackClusterDistance: bestTrkMatch %d origTrkMatch %d distance %f\n", bestTrkMatchId, *trkMatchId, distance);
    *trkMatchId = bestTrkMatchId;
    return distance;
}

UInt_t AliAnalysisEtSelectorEmcal::GetLabel(const AliESDCaloCluster *cluster, AliStack *stack){
  UInt_t correctLabel = TMath::Abs(cluster->GetLabel());
  correctLabel = GetFirstMotherNotFromDetectorCover(correctLabel,*stack);

  //Now we want to see if this particle is really just something that converted in the cover of the detector and if so, override the label
//   if(!stack){return 0;}
//   if( correctLabel > 0 && stack->IsSecondaryFromMaterial(correctLabel)){//if this is flagged as a secondary then we look to see where it really came from
//     TParticle *hitParticle = stack->Particle(correctLabel);
//     if(hitParticle){
//       Bool_t partVtxSecondary = (TMath::Sqrt(hitParticle->Vx()*hitParticle->Vx() + hitParticle->Vy()*hitParticle->Vy()) >420);
//       if(partVtxSecondary){//at this point we have something which converted near the detector.  Let's find the mother particle
// 	UInt_t mothIdx = stack->Particle(correctLabel)->GetMother(0);
// 	if(mothIdx>0){
// 	  TParticle *mother = stack->Particle(mothIdx);
// 	  if(mother){
// 	    if(AliAnalysisEtSelector::CutGeometricalAcceptance(*mother)){//and the mother is in the acceptance
// 	      if( !(mother->GetPdgCode()==fgPi0Code)){//some of these are decays that just happen this far out
// 		cout<<"I am declaring that "<<hitParticle->GetName()<<" with a vertex of "<< TMath::Sqrt(hitParticle->Vx()*hitParticle->Vx() + hitParticle->Vy()*hitParticle->Vy()) <<" is actually "<<mother->GetName()<<endl;
// 		cout<<"ID check "<<mothIdx<<" vs "<<mother->GetUniqueID()<<endl;
// 		//so now we know that the particle originated near the cover and within the acceptance of the detector
// 		return mothIdx;
// 	      }
// 	    }
// 	  }
// 	}
//       }

//     }
//   }
  return  correctLabel ; // DS: should this line be inside n>0 check, and return another value if n<=0 ? 

}
Commit	Line	Data
69fe2599	1	//_________________________________________________________________________
	2	// Utility Class for transverse energy studies
	3	// Selection class for EMCAL
	4	//
	5	//*-- Authors: Oystein Djuvsland (Bergen)
	6	//_________________________________________________________________________
	7
	8
	9	#include "AliAnalysisEtSelectorEmcal.h"
	10	#include "AliAnalysisEtCuts.h"
f20da103	11	#include "AliEMCALTrack.h"
69fe2599	12	#include "TParticle.h"
6844c491	13	#include <iostream>
6844c491	14	#include "AliStack.h"
69fe2599	15
7f03b597	16	AliAnalysisEtSelectorEmcal::AliAnalysisEtSelectorEmcal(AliAnalysisEtCuts* cuts):AliAnalysisEtSelector(cuts)
69fe2599	17	{
69fe2599	18
c31562f7	19	}
	20	AliAnalysisEtSelectorEmcal::AliAnalysisEtSelectorEmcal():AliAnalysisEtSelector()
	21	{
69fe2599	22	}
69fe2599	23
7f03b597	24	AliAnalysisEtSelectorEmcal::~AliAnalysisEtSelectorEmcal()
69fe2599	25	{
	26
	27	}
	28
	29	void AliAnalysisEtSelectorEmcal::Init()
	30	{
	31	AliAnalysisEtSelector::Init();
	32	}
	33
	34	Int_t AliAnalysisEtSelectorEmcal::Init(const AliESDEvent* event)
f20da103	35	{ // Init
f20da103	36
69fe2599	37	AliAnalysisEtSelector::Init(event);
	38	Printf("Initializing selector for run: %d", event->GetRunNumber());
	39	fInitialized = kTRUE;
	40	return 0;
	41	}
	42
	43	TRefArray* AliAnalysisEtSelectorEmcal::GetClusters()
f20da103	44	{ // Get clusters
	45
	46	if(!fClusterArray) fClusterArray = new TRefArray;
	47
	48	if(fClusterArray)
	49	{
	50	fEvent->GetEMCALClusters(fClusterArray);
	51	}
	52	else
	53	{
	54	Printf("Could not initialize cluster array");
	55	}
	56	return fClusterArray;
69fe2599	57	}
69fe2599	58
86e7d5db	59	Bool_t AliAnalysisEtSelectorEmcal::PassMinEnergyCut(const AliESDCaloCluster& cl) const
69fe2599	60	{
cfe23ff0	61	Float_t pos[3];
	62	cl.GetPosition(pos);
	63	TVector3 cp(pos);
	64	return TMath::Sin(cp.Theta())*cl.E() > fCuts->GetReconstructedEmcalClusterEnergyCut();
69fe2599	65	}
69fe2599	66
86e7d5db	67	Bool_t AliAnalysisEtSelectorEmcal::PassMinEnergyCut(const TParticle& p) const
69fe2599	68	{
cfe23ff0	69	return TMath::Sin(p.Theta())*p.Energy() > fCuts->GetReconstructedEmcalClusterEnergyCut();
69fe2599	70	}
02d47689	71	Bool_t AliAnalysisEtSelectorEmcal::PassMinEnergyCut(Double_t e) const
	72	{
	73	return e > fCuts->GetReconstructedEmcalClusterEnergyCut();
	74	}
69fe2599	75
86e7d5db	76	Bool_t AliAnalysisEtSelectorEmcal::PassDistanceToBadChannelCut(const AliESDCaloCluster& ) const
69fe2599	77	{
7f03b597	78	return kTRUE;
69fe2599	79	}
69fe2599	80
f20da103	81	Bool_t AliAnalysisEtSelectorEmcal::PassTrackMatchingCut(const AliESDCaloCluster& cluster) const
69fe2599	82	{
f20da103	83
fe4603e7	84	Int_t nTracksMatched = cluster.GetNTracksMatched();
	85	if(nTracksMatched == 0){
	86	return kTRUE;
	87	}
	88
	89	Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
	90	if(trackMatchedIndex < 0){
	91	return kTRUE;
	92	}
ea3ce170	93	AliVParticle *track = fEvent->GetTrack(trackMatchedIndex);
	94	if(track->Pt()<0.5) return kTRUE;//Track matches below about 500 MeV are mostly random. It takes ~460 MeV to reach the EMCal
	95
	96	//Float_t recoE = cluster.E();
fe4603e7	97	//Float_t pos[3];
	98
	99	// //cluster.GetPosition(pos);
	100	// Int_t trackMatchIdx = cluster.GetTrackMatchedIndex();
	101	// //Double_t distance = 9999.0;
	102	// if(trackMatchIdx>-1)
	103	// {
	104	// return kTRUE;
	105	// //distance = CalcTrackClusterDistance(pos, &trackMatchIdx);
	106	// }
	107
fe4603e7	108
	109	return kFALSE;
	110	//return distance > fCuts->GetEmcalTrackDistanceCut();
69fe2599	111	}
69fe2599	112
43dd5a38	113	Bool_t AliAnalysisEtSelectorEmcal::CutGeometricalAcceptance(const TParticle& part)
69fe2599	114	{
43dd5a38	115	float myphi = part.Phi();
43dd5a38	116	myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
f20da103	117	return TMath::Abs(part.Eta()) < fCuts->GetGeometryEmcalEtaAccCut()
43dd5a38	118	&& myphi < fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180.
43dd5a38	119	&& myphi > fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.;
69fe2599	120	}
69fe2599	121
43dd5a38	122	Bool_t AliAnalysisEtSelectorEmcal::CutGeometricalAcceptance(const AliVTrack& part)
69fe2599	123	{
43dd5a38	124	float myphi = part.Phi();
	125	myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
	126	return TMath::Abs(part.Eta()) < fCuts->GetGeometryEmcalEtaAccCut()
	127	&& myphi < fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180.
	128	&& myphi > fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.;
69fe2599	129	}
69fe2599	130
43dd5a38	131	Bool_t AliAnalysisEtSelectorEmcal::CutGeometricalAcceptance(const AliESDCaloCluster& cluster)
31c813d5	132	{
	133	Float_t pos[3];
	134	cluster.GetPosition(pos);
	135	TVector3 cp(pos);
43dd5a38	136	float myphi = cp.Phi();
43dd5a38	137	myphi = AliAnalysisEtSelector::ShiftAngle(myphi);
31c813d5	138	return TMath::Abs(cp.Eta()) < fCuts->GetGeometryEmcalEtaAccCut()
43dd5a38	139	&& myphi < fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180.
43dd5a38	140	&& myphi > fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.;
31c813d5	141	}
31c813d5	142
69fe2599	143
f20da103	144	Double_t
	145	AliAnalysisEtSelectorEmcal::CalcTrackClusterDistance(const Float_t clsPos[3],Int_t *trkMatchId) const
	146	{ // calculate distance between cluster and closest track
	147
	148	Double_t trkPos[3] = {0,0,0};
	149
	150	Int_t bestTrkMatchId = -1;
	151	Double_t distance = 9999; // init to a big number
	152
	153	Double_t dist = 0;
	154	Double_t distX = 0, distY = 0, distZ = 0;
	155
	156	for (Int_t iTrack = 0; iTrack < fEvent->GetNumberOfTracks(); iTrack++) {
	157	AliESDtrack track = static_cast<AliESDtrack>( fEvent->GetTrack(iTrack) );
	158	if (!track) {
	159	AliError(Form("ERROR: Could not get track %d", iTrack));
	160	continue;
	161	}
	162
	163	// check for approx. eta and phi range before we propagate..
	164	// TBD
	165
	166	AliEMCALTrack emctrack(*track);
	167	if (!emctrack.PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) {
	168	continue;
	169	}
	170	emctrack.GetXYZ(trkPos);
	171
	172	distX = clsPos[0]-trkPos[0];
	173	distY = clsPos[1]-trkPos[1];
	174	distZ = clsPos[2]-trkPos[2];
	175	dist = TMath::Sqrt(distXdistX + distYdistY + distZ*distZ);
	176
	177	if (dist < distance) {
	178	distance = dist;
	179	bestTrkMatchId = iTrack;
	180	}
	181	} // iTrack
	182
	183	// printf("CalcTrackClusterDistance: bestTrkMatch %d origTrkMatch %d distance %f\n", bestTrkMatchId, *trkMatchId, distance);
	184	*trkMatchId = bestTrkMatchId;
	185	return distance;
	186	}
	187
6844c491	188	UInt_t AliAnalysisEtSelectorEmcal::GetLabel(const AliESDCaloCluster cluster, AliStack stack){
	189	UInt_t correctLabel = TMath::Abs(cluster->GetLabel());
	190	correctLabel = GetFirstMotherNotFromDetectorCover(correctLabel,*stack);
	191
	192	//Now we want to see if this particle is really just something that converted in the cover of the detector and if so, override the label
	193	// if(!stack){return 0;}
	194	// if( correctLabel > 0 && stack->IsSecondaryFromMaterial(correctLabel)){//if this is flagged as a secondary then we look to see where it really came from
	195	// TParticle *hitParticle = stack->Particle(correctLabel);
	196	// if(hitParticle){
	197	// Bool_t partVtxSecondary = (TMath::Sqrt(hitParticle->Vx()hitParticle->Vx() + hitParticle->Vy()hitParticle->Vy()) >420);
	198	// if(partVtxSecondary){//at this point we have something which converted near the detector. Let's find the mother particle
	199	// UInt_t mothIdx = stack->Particle(correctLabel)->GetMother(0);
	200	// if(mothIdx>0){
	201	// TParticle *mother = stack->Particle(mothIdx);
	202	// if(mother){
	203	// if(AliAnalysisEtSelector::CutGeometricalAcceptance(*mother)){//and the mother is in the acceptance
	204	// if( !(mother->GetPdgCode()==fgPi0Code)){//some of these are decays that just happen this far out
	205	// cout<<"I am declaring that "<<hitParticle->GetName()<<" with a vertex of "<< TMath::Sqrt(hitParticle->Vx()hitParticle->Vx() + hitParticle->Vy()hitParticle->Vy()) <<" is actually "<<mother->GetName()<<endl;
	206	// cout<<"ID check "<<mothIdx<<" vs "<<mother->GetUniqueID()<<endl;
	207	// //so now we know that the particle originated near the cover and within the acceptance of the detector
	208	// return mothIdx;
	209	// }
	210	// }
	211	// }
	212	// }
	213	// }
69fe2599	214
6844c491	215	// }
	216	// }
	217	return correctLabel ; // DS: should this line be inside n>0 check, and return another value if n<=0 ?
69fe2599	218
6844c491	219	}