#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()
}
Int_t AliAnalysisEtSelectorEmcal::Init(const AliESDEvent* event)
-{ // Init
-
+{ // Init
+
AliAnalysisEtSelector::Init(event);
Printf("Initializing selector for run: %d", event->GetRunNumber());
fInitialized = kTRUE;
}
TRefArray* AliAnalysisEtSelectorEmcal::GetClusters()
-{ // Get clusters
-
- if(!fClusterArray) fClusterArray = new TRefArray;
-
- if(fClusterArray)
- {
- fEvent->GetEMCALClusters(fClusterArray);
- }
- else
- {
- Printf("Could not initialize cluster array");
- }
- return fClusterArray;
+{ // 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
{
- return cl.E() > fCuts->GetReconstructedEmcalClusterEnergyCut();
+ 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 p.Energy() > fCuts->GetReconstructedEmcalClusterEnergyCut();
+ 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& ) const
+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) const
+Bool_t AliAnalysisEtSelectorEmcal::CutGeometricalAcceptance(const TParticle& part)
{
- return TMath::Abs(part.Eta()) < fCuts->GetGeometryEmcalEtaAccCut()
- && part.Phi() < fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180.
- && part.Phi() > fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.;
+ 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) const
+Bool_t AliAnalysisEtSelectorEmcal::CutGeometricalAcceptance(const AliVTrack& part)
{
- return TMath::Abs(part.Eta()) < fCuts->GetGeometryEmcalEtaAccCut()
- && part.Phi() < fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180.
- && part.Phi() > fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.;
+ 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 ?
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff