* about the suitability of this software for any purpose. It is *
* provided "as is" without express or implied warranty. *
**************************************************************************/
-//#include <TBranch.h>
-//#include <TCanvas.h>
-//#include <TChain.h>
-//#include <TFile.h>
-//#include <TH1F.h>
-//#include <TH1I.h>
-//#include <TH2F.h>
+
#include <TList.h>
-//#include <TLorentzVector.h>
#include <TMath.h>
#include <TObjArray.h>
#include <TObject.h>
-//#include <TProfile.h>
-//#include <TRandom.h>
-//#include <TSystem.h>
-//#include <TTree.h>
#include <TVector3.h>
#include "AliAnalyseLeadingTrackUE.h"
-//#include "AliAnalysisTask.h"
-//#include "AliAnalysisHelperJetTasks.h"
-//#include "AliAnalysisManager.h"
#include "AliAODEvent.h"
-//#include "AliAODHandler.h"
-//#include "AliAODJet.h"
#include "AliAODMCParticle.h"
#include "AliAODTrack.h"
#include "AliESDEvent.h"
#include "AliESDtrack.h"
#include "AliESDtrackCuts.h"
-//#include "AliGenPythiaEventHeader.h"
#include "AliInputEventHandler.h"
-//#include "AliKFVertex.h"
-//#include "AliLog.h"
#include "AliMCEvent.h"
#include "AliVParticle.h"
#include "AliAODMCHeader.h"
+#include "TFormula.h"
#include "AliAnalysisManager.h"
#include "AliMCEventHandler.h"
#include "AliStack.h"
#include "AliPIDResponse.h"
+#include "AliHelperPID.h"
////////////////////////////////////////////////
TObject(),
fDebug(0),
fFilterBit(16),
+ fTrackStatus(0),
fOnlyHadrons(kFALSE),
+ fCheckMotherPDG(kTRUE),
fTrackEtaCut(0.8),
+ fTrackEtaCutMin(-1.),
fTrackPtMin(0),
fEventSelection(AliVEvent::kMB|AliVEvent::kUserDefined),
+ fDCAXYCut(0),
+ fSharedClusterCut(-1),
+ fCrossedRowsCut(-1),
+ fFoundFractionCut(-1),
fEsdTrackCuts(0x0),
fEsdTrackCutsExtra1(0x0),
- fEsdTrackCutsExtra2(0x0)
+ fEsdTrackCutsExtra2(0x0),
+ fHelperPID(0x0),
+ fEventCounter(0)
{
// constructor
}
fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
fEsdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kOff);
+ fEsdTrackCuts->SetMaxFractionSharedTPCClusters(0.4);
// Add SPD requirement
fEsdTrackCutsExtra1 = new AliESDtrackCuts("SPD", "Require 1 cluster in SPD");
AliVParticle* part = ParticleWithCuts( obj, ipart );
if (!part) continue;
// Accept leading-tracks in a limited pseudo-rapidity range
- if( TMath::Abs(part->Eta()) > fTrackEtaCut ) continue;
+ if( TMath::Abs(part->Eta()) > fTrackEtaCut || TMath::Abs(part->Eta()) < fTrackEtaCutMin ) continue;
tracks->AddLast( part );
}
// Order tracks by pT
Printf("WARNING: No mother found for particle %d:", part->GetLabel());
continue;
}
-
+
+// Printf("%d %d %d", i, part->GetLabel(), mother->GetLabel());
if (mother->GetLabel() >= maxLabel)
{
-// Printf("Removing %d with label %d", i, part->GetLabel()); part->Dump();
+// Printf("Removing %d with label %d", i, part->GetLabel()); ((AliMCParticle*)part)->Particle()->Print(); ((AliMCParticle*)mother)->Particle()->Print();
TObject* object = tracks->RemoveAt(i);
if (tracks->IsOwner())
delete object;
}
//-------------------------------------------------------------------
-TObjArray* AliAnalyseLeadingTrackUE::GetAcceptedParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts)
+TObjArray* AliAnalyseLeadingTrackUE::GetAcceptedParticles(TObject* obj, TObject* arrayMC, Bool_t onlyprimaries, Int_t particleSpecies, Bool_t useEtaPtCuts, Bool_t speciesOnTracks)
{
// Returns an array of particles that pass the cuts, if arrayMC is given each reconstructed particle is replaced by its corresponding MC particles, depending on the parameter onlyprimaries only for primaries
// particleSpecies: -1 all particles are returned
// 0 (pions) 1 (kaons) 2 (protons) 3 (others) particles
-
+ // speciesOnTracks if kFALSE, particleSpecies is only applied on the matched MC particle (not on the track itself)
+
Int_t nTracks = NParticles(obj);
TObjArray* tracks = new TObjArray;
// Loop over tracks or jets
for (Int_t ipart=0; ipart<nTracks; ++ipart) {
- AliVParticle* part = ParticleWithCuts( obj, ipart, onlyprimaries, particleSpecies );
+ AliVParticle* part = ParticleWithCuts( obj, ipart, onlyprimaries, (speciesOnTracks) ? particleSpecies : -1);
if (!part) continue;
if (useEtaPtCuts)
- if (TMath::Abs(part->Eta()) > fTrackEtaCut || part->Pt() < fTrackPtMin)
+ if (TMath::Abs(part->Eta()) > fTrackEtaCut || TMath::Abs(part->Eta()) < fTrackEtaCutMin || part->Pt() < fTrackPtMin)
{
if (hasOwnership)
delete part;
continue;
}
+
+// Printf("%p %p %d Accepted %d %f %f %f", obj, arrayMC, particleSpecies, ipart, part->Eta(), part->Phi(), part->Pt());
if (arrayMC) {
Int_t label = part->GetLabel();
if (!partReconstructed) continue;
if (useEtaPtCuts)
- if (TMath::Abs(partReconstructed->Eta()) > fTrackEtaCut || partReconstructed->Pt() < fTrackPtMin)
+ if (TMath::Abs(partReconstructed->Eta()) > fTrackEtaCut || TMath::Abs(partReconstructed->Eta()) < fTrackEtaCutMin || partReconstructed->Pt() < fTrackPtMin)
{
if (hasOwnership)
delete partReconstructed;
if (particleSpecies != -1) {
// find the primary mother
AliVParticle* mother = part;
- while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
- {
+ if(fCheckMotherPDG) {
+ while (!((AliAODMCParticle*)mother)->IsPhysicalPrimary())
+ {
if (((AliAODMCParticle*)mother)->GetMother() < 0)
{
mother = 0;
mother = (AliVParticle*) arrayMC->At(((AliAODMCParticle*)mother)->GetMother());
if (!mother)
break;
+ }
}
-
if (mother)
{
Int_t pdgCode = ((AliAODMCParticle*)mother)->GetPdgCode();
// find the primary mother
AliMCParticle* mother = (AliMCParticle*) part;
// Printf("");
- while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
- {
+ if(fCheckMotherPDG) {
+ while (!mcEvent->IsPhysicalPrimary(mother->GetLabel()))
+ {
// Printf("pdg = %d; mother = %d", mother->PdgCode(), mother->GetMother());
if (mother->GetMother() < 0)
{
mother = (AliMCParticle*) mcEvent->GetTrack(mother->GetMother());
if (!mother)
break;
+ }
}
-
if (mother)
{
Int_t pdgCode = mother->PdgCode();
}else if (obj->InheritsFrom("AliAODEvent")){ // RECO AOD TRACKS
AliAODEvent *aodEvent = static_cast<AliAODEvent*>(obj);
part = aodEvent->GetTrack(ipart);
+
// track selection cuts
if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) ) return 0;
- //if ( !(((AliAODTrack*)part)->TestFilterBit(fFilterBit)) && !(((AliAODTrack*)part)->TestFilterBit(32)) ) return 0;
- // only primary candidates
- //if ( ((AliAODTrack*)part)->IsPrimaryCandidate() )return 0;
+ if (fTrackStatus != 0 && !CheckTrack(part)) return 0;
+
+ // DCA XY
+ if (fDCAXYCut)
+ {
+ const AliVVertex* vertex = aodEvent->GetPrimaryVertex();
+ if (!vertex)
+ return 0;
+
+ Double_t pos[2];
+ Double_t covar[3];
+ AliAODTrack* clone = (AliAODTrack*) part->Clone();
+ Bool_t success = clone->PropagateToDCA(vertex, aodEvent->GetHeader()->GetMagneticField(), 3, pos, covar);
+ delete clone;
+ if (!success)
+ return 0;
+
+// Printf("%f", ((AliAODTrack*)part)->DCA());
+// Printf("%f", pos[0]);
+ if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(part->Pt()))
+ return 0;
+ }
+
+ if (fSharedClusterCut >= 0)
+ {
+ Double_t frac = Double_t(((AliAODTrack*)part)->GetTPCnclsS()) / Double_t(((AliAODTrack*)part)->GetTPCncls());
+ if (frac > fSharedClusterCut)
+ return 0;
+ }
+
+ if (fCrossedRowsCut >= 0)
+ {
+ if (((AliAODTrack*) part)->GetTPCNCrossedRows() < fCrossedRowsCut)
+ return 0;
+ }
+
+ if (fFoundFractionCut >= 0)
+ {
+ UInt_t findableClusters = ((AliAODTrack*) part)->GetTPCNclsF();
+ if (findableClusters == 0)
+ return 0;
+ if (((Double_t) ((AliAODTrack*) part)->GetTPCNCrossedRows() / findableClusters) < fFoundFractionCut)
+ return 0;
+ }
+
// eventually only hadrons
if (fOnlyHadrons){
Bool_t isHadron = ((AliAODTrack*)part)->GetMostProbablePID()==AliAODTrack::kPion ||
if (!isHadron) return 0;
}
- if (particleSpecies != -1 && GetParticleSpecies((AliVTrack*) part) != particleSpecies) return 0;
+ if (particleSpecies != -1 && fHelperPID->GetParticleSpecies((AliVTrack*) part,kTRUE) != particleSpecies) return 0;
}else if (obj->InheritsFrom("AliESDEvent")){ // RECO ESD TRACKS
AliESDEvent *esdEvent = static_cast<AliESDEvent*>(obj);
part = esdEvent->GetTrack(ipart);
if (!part)return 0;
+
// track selection cuts
-
if (!( ApplyCuts(part)) )
- return 0;
+ return 0;
+ if (fTrackStatus != 0 && !CheckTrack(part)) return 0;
+
if (fFilterBit == 128 || fFilterBit == 256 || fFilterBit == 512 || fFilterBit == 1024)
{
// create TPC only tracks constrained to the SPD vertex
if (!isHadron) return 0;
}
*/
- if (particleSpecies != -1 && GetParticleSpecies((AliVTrack*) part) != particleSpecies) return 0; // If it is -1 you take all the particles
+ if (particleSpecies != -1 && fHelperPID->GetParticleSpecies((AliVTrack*) part,kTRUE) != particleSpecies) return 0; // If it is -1 you take all the particles
}else {
if (fDebug > 1) AliFatal(" Analysis type not defined !!! ");
// only charged
if (!part->Charge())return 0;
+ part->SetUniqueID(fEventCounter * 100000 + ipart);
return part;
}
TVector3 partVect(part->Px(), part->Py(), part->Pz());
Int_t region = 0;
- if( TMath::Abs(partVect.Eta()) > fTrackEtaCut ) continue;
+ if( TMath::Abs(partVect.Eta()) > fTrackEtaCut || TMath::Abs(partVect.Eta()) < fTrackEtaCutMin) continue;
// transverse regions
if (leadVect.DeltaPhi(partVect) < -k60rad && leadVect.DeltaPhi(partVect) > -k120rad )region = -1; //left
if (leadVect.DeltaPhi(partVect) > k60rad && leadVect.DeltaPhi(partVect) < k120rad ) region = 1; //right
if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
// Select a quality vertex by number of tracks?
- if( nTracksPrim < ntracks || TMath::Abs(zVertex) > zed ) {
+ if( nTracksPrim < ntracks || TMath::Abs(zVertex) >= zed ) {
if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
return kFALSE;
}
if (obj->InheritsFrom("AliMCEvent"))
{
- if (TMath::Abs(((AliMCEvent*) obj)->GetPrimaryVertex()->GetZ()) > zed)
+ if (TMath::Abs(((AliMCEvent*) obj)->GetPrimaryVertex()->GetZ()) >= zed)
{
if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
return kFALSE;
if (obj->InheritsFrom("AliAODMCHeader"))
{
- if (TMath::Abs(((AliAODMCHeader*) obj)->GetVtxZ()) > zed)
+ if (TMath::Abs(((AliAODMCHeader*) obj)->GetVtxZ()) >= zed)
{
if (fDebug > 1) AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
return kFALSE;
if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return kFALSE;
// Select a quality vertex by number of tracks?
- if( nTracksPrim < ntracks || TMath::Abs(zVertex) > zed ) {
+ if( nTracksPrim < ntracks || TMath::Abs(zVertex) >= zed ) {
if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
return kFALSE;
}
return kTRUE;
}
+//____________________________________________________________________
-Int_t AliAnalyseLeadingTrackUE::GetParticleSpecies(AliVTrack * trk)
+Bool_t AliAnalyseLeadingTrackUE::CheckTrack(AliVParticle * part)
{
- // return PID according to detectors
- // Get PID response object, if needed
- // FIXME: THIS IS UGLY!! THE CODE IS COPIED AND ADAPTED FROM PWGLF/SPECTRA/PiKaPr/TestAOD/AliSpectraBothPID.cxx
- // TO BE REPLACED BY THE PROPER CLASS ASAP
- // MF 10/12/12
-
-
- enum BothPIDType_t {kNSigmaTPC,kNSigmaTOF,kNSigmaTPCTOF};
-
- enum BothParticleSpecies_t
- {
- kSpPion = 0, kSpKaon, kSpProton,
- kNSpecies,
- kSpUndefined,
- }; // Particle species used in plotting
-
- UInt_t fPIDType = kNSigmaTPCTOF;
-
- // Check TOF matching & PID
- // FIXME: this should go in the track selection!!
- UInt_t status;
- status=trk->GetStatus();
- if((status&AliAODTrack::kTOFout)==0 || (status&AliAODTrack::kTIME)==0){//kTOFout and kTIME
- return kSpUndefined;
- }
-
-
-
- // guess the particle based on the smaller nsigma
-
- // rec[kSpPion]=false;
- // rec[kSpKaon]=false;
- // rec[kSpProton]=false;
-
- static AliPIDResponse * fPIDResponse = 0;
-
- if(!fPIDResponse)
- {
- AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
- AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
- fPIDResponse = inputHandler->GetPIDResponse();
- }
-
- if(!fPIDResponse)
- {
- AliFatal("Cannot get pid response");
- return 0;
- }
-
-
- // Compute nsigma for each hypthesis
- AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk);
-
- const Double_t fshiftTPC = 0, fshiftTOF=0;
- const Double_t fNSigmaPID = 1; // FIXME
- // --- TPC
- Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton)+fshiftTPC);
- Double_t nsigmaTPCkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)+fshiftTPC);
- Double_t nsigmaTPCkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)+fshiftTPC);
-
-
- Double_t nsigmaTOFkProton=100.0,nsigmaTOFkKaon=100.0,nsigmaTOFkPion=100.0;
-
- Double_t nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
- Double_t nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
- Double_t nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
- if(fPIDType==kNSigmaTOF)
- {
- nsigmaTPCTOFkProton = 100.0;
- nsigmaTPCTOFkKaon = 100.0;
- nsigmaTPCTOFkPion = 100.0;
-
- }
-
-
- nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton)+fshiftTOF);
- nsigmaTOFkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)+fshiftTOF);
- nsigmaTOFkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)+fshiftTOF);
- // the TOF info is used in combined
- nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2.);
- nsigmaTPCTOFkKaon = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2.);
- nsigmaTPCTOFkPion = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2.);
-
- // select the nsigma to be used for the actual PID
- Double_t nsigmaPion=100;
- Double_t nsigmaKaon=100;
- Double_t nsigmaProton=100;
-
- switch (fPIDType)
- {
- case kNSigmaTPC:
- nsigmaProton = nsigmaTPCkProton;
- nsigmaKaon = nsigmaTPCkKaon ;
- nsigmaPion = nsigmaTPCkPion ;
- break;
- case kNSigmaTOF:
- nsigmaProton = nsigmaTOFkProton;
- nsigmaKaon = nsigmaTOFkKaon ;
- nsigmaPion = nsigmaTOFkPion ;
- break;
- case kNSigmaTPCTOF:
- nsigmaProton = nsigmaTPCTOFkProton;
- nsigmaKaon = nsigmaTPCTOFkKaon ;
- nsigmaPion = nsigmaTPCTOFkPion ;
- break;
- }
-
- // if(nsigmaPion < fNSigmaPID)
- // rec[kSpPion]=true;
- // if(nsigmaKaon < fNSigmaPID)
- // rec[kSpKaon]=true;
- // if(nsigmaProton < fNSigmaPID)
- // rec[kSpProton]=true;
-
- if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton ))
- return kSpUndefined;//if is the default value for the three
- if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID))
- {
- return kSpKaon;
- }
- if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID))
- {
- return kSpPion;
- }
- if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID))
- {
- return kSpProton;
- }
- // else, return undefined
- return kSpUndefined;
+ // check if the track status flags are set
+ UInt_t status=((AliVTrack*)part)->GetStatus();
+ if ((status & fTrackStatus) == fTrackStatus)
+ return kTRUE;
+ return kFALSE;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff