fHFEpid->AddDetector("TOF", 0);
fHFEpid->AddDetector("TPC", 1);
cout<<endl<<" ---> TPC and TOF added to the PID"<<endl;
- fHFEpid->ConfigureTPCrejection();
+ fHFEpid->ConfigureTOF();
+ fHFEpid->ConfigureTPCdefaultCut();
fHFEpid->InitializePID();
}
fDePID->SetHasMCData(HasMCData());
fDePID->AddDetector("TPC", 0);
fDePID->AddDetector("TOF", 1);
- fDePID->ConfigureTPCrejection();
+ fDePID->ConfigureTOF();
+ fDePID->ConfigureTPCdefaultCut();
fDePID->InitializePID(); // Only restrictions to TPC allowed
#include "AliMCEvent.h"
#include "AliMCEventHandler.h"
#include "AliMCParticle.h"
+#include "AliMultiplicity.h"
#include "AliPID.h"
#include "AliStack.h"
#include "AliTriggerAnalysis.h"
, fQAlevel(0)
, fPlugins(0)
, fFillSignalOnly(kTRUE)
+ , fFillNoCuts(kFALSE)
, fBackGroundFactorApply(kFALSE)
, fRemovePileUp(kFALSE)
, fIdentifiedAsPileUp(kFALSE)
, fHasSpecialTriggerSelection(kFALSE)
, fRejectKinkMother(kTRUE)
, fSpecialTrigger("NONE")
- , fCentralityF(99.0)
+ , fCentralityF(-1)
, fContributors(0.5)
, fWeightBackGround(0.)
, fVz(0.0)
, fQAlevel(0)
, fPlugins(0)
, fFillSignalOnly(kTRUE)
+ , fFillNoCuts(kFALSE)
, fBackGroundFactorApply(kFALSE)
, fRemovePileUp(kFALSE)
, fIdentifiedAsPileUp(kFALSE)
, fHasSpecialTriggerSelection(kFALSE)
, fRejectKinkMother(kTRUE)
, fSpecialTrigger("NONE")
- , fCentralityF(99.0)
+ , fCentralityF(-1)
, fContributors(0.5)
, fWeightBackGround(0.)
, fVz(0.0)
, fQAlevel(0)
, fPlugins(0)
, fFillSignalOnly(ref.fFillSignalOnly)
+ , fFillNoCuts(ref.fFillNoCuts)
, fBackGroundFactorApply(ref.fBackGroundFactorApply)
, fRemovePileUp(ref.fRemovePileUp)
, fIdentifiedAsPileUp(ref.fIdentifiedAsPileUp)
target.fQAlevel = fQAlevel;
target.fPlugins = fPlugins;
target.fFillSignalOnly = fFillSignalOnly;
+ target.fFillNoCuts = fFillNoCuts;
target.fBackGroundFactorApply = fBackGroundFactorApply;
target.fRemovePileUp = fRemovePileUp;
target.fIdentifiedAsPileUp = fIdentifiedAsPileUp;
// First Part: Make QA histograms
fQACollection = new AliHFEcollection("TaskQA", "QA histos from the Electron Task");
fQACollection->CreateTH1F("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100);
- fQACollection->CreateProfile("conr", "Electron PID contamination", 20, 0, 20);
- fQACollection->CreateTH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi());
- fQACollection->CreateTH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi());
fQACollection->CreateTH1F("nElectron", "Number of electrons", 100, 0, 100);
- fQACollection->CreateProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20);
- fQACollection->CreateProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20);
- fQACollection->CreateTH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20);
- fQACollection->CreateTH1F("mccharge", "MC Charge", 200, -100, 100);
fQACollection->CreateTH2F("radius", "Production Vertex", 100, 0.0, 5.0, 100, 0.0, 5.0);
- // Temporary histograms for TPC number of clusters for all signal tracks (MC true electrons) and for selected tracks (Markus Fasel)
- fQACollection->CreateTH2F("TPCclusters2_1_Signal", "TPCclusterInfo for findable clusters for 2 neighbors for signal tracks", 30, 0.1, 10., 162, 0., 161.);
- fQACollection->CreateTH2F("TPCclusters2_0_Signal", "TPCclusterInfo for the ratio for 2 neighbors for signal tracks", 30, 0.1, 10., 100, 0., 1.);
- fQACollection->CreateTH2F("TPCclusters2_1_Selected", "TPCclusterInfo for findable clusters for 2 neighbors for selected tracks", 30, 0.1, 10., 162, 0., 161.);
- fQACollection->CreateTH2F("TPCclusters2_0_Selected", "TPCclusterInfo for the ratio for 2 neighbors for selected tracks", 30, 0.1, 10., 110, 0., 1.1);
- fQACollection->CreateTH2F("TPCncls_Signal", "TPC Number of clusters for signal tracks", 30, 0.1, 10., 162, 0., 161.);
- fQACollection->CreateTH2F("TPCclr_Signal", "TPC cluster ratio for signal tracks", 30, 0.1, 10., 110, 0., 1.1);
- fQACollection->BinLogAxis("TPCclusters2_1_Signal", 0);
- fQACollection->BinLogAxis("TPCclusters2_0_Signal", 0);
- fQACollection->BinLogAxis("TPCclusters2_1_Selected", 0);
- fQACollection->BinLogAxis("TPCclusters2_0_Selected", 0);
- fQACollection->BinLogAxis("TPCncls_Signal", 0);
- fQACollection->BinLogAxis("TPCclr_Signal", 0);
// Temporary histograms for TRD efficiency maps (Markus Fasel)
fQACollection->CreateTH2F("TRDmapPosBefore", "Pos. charged tracks before TRD; #eta; #phi", 100, -1., 1., 180, 0., 2*TMath::Pi());
fQACollection->CreateTH2F("TRDmapNegBefore", "Neg. charged tracks before TRD; #eta; #phi", 100, -1., 1., 180, 0., 2*TMath::Pi());
}
// need the centrality for everything (MC also)
- fCentralityF = -100.0;
- if(!ReadCentrality()) fCentralityF = -100.0;
+ fCentralityF = -1;
+ if(!ReadCentrality()) fCentralityF = -1;
//printf("pass centrality\n");
//printf("Reading fCentralityF %f\n",fCentralityF);
fCFM->GetEventContainer()->Fill(eventContainer,AliHFEcuts::kEventStepGenerated);
Int_t nElectrons = 0;
if(IsESDanalysis()){
- if(!((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (TMath::Abs(fCentralityF+100.0) < 0.01))){ //kStepMCGeneratedZOutNoPileUpCentralityFine
+ if(!((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (fCentralityF < 0))){ //kStepMCGeneratedZOutNoPileUpCentralityFine
if (HasMCData() && IsQAOn(kMCqa)) {
AliDebug(2, "Running MC QA");
if(fTaggedTrackAnalysis) {
fTaggedTrackAnalysis->SetMagneticField(fESD->GetMagneticField());
fTaggedTrackAnalysis->SetCentrality(fCentralityF);
+ if(IsPbPb()) fTaggedTrackAnalysis->SetPbPb();
+ else fTaggedTrackAnalysis->SetPP();
}
// Do event Normalization
fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoPileUp);
//
- if(TMath::Abs(fCentralityF+100.0) < 0.01) return;
+ if(TMath::Abs(fCentralityF) < 0) return;
fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecCentralityOk);
//printf("In ProcessESD %f\n",fCentralityF);
Int_t nElectronCandidates = 0;
AliESDtrack *track = NULL, *htrack = NULL;
AliMCParticle *mctrack = NULL;
+ TParticle* mctrack4QA = NULL;
Int_t pid = 0;
Bool_t signal = kTRUE;
if(HasMCData()){
// Check if it is electrons near the vertex
if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
+ mctrack4QA = mctrack->Particle();
if(fFillSignalOnly && !fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
else AliDebug(3, "Signal Electron");
// Cache new Track information inside the var manager
fVarManager->NewTrack(track, mctrack, fCentralityF, -1, signal);
- if(signal) {
- fVarManager->FillContainer(fContainer, "recTrackContReco", AliHFEcuts::kStepRecNoCut, kFALSE);
- fVarManager->FillContainer(fContainer, "recTrackContMC", AliHFEcuts::kStepRecNoCut, kTRUE);
- if((track->GetStatus() & AliESDtrack::kTPCout)
- && (TMath::Abs(track->Eta()) < 0.8)
- && (track->GetKinkIndex(0) == 0)){
- fQACollection->Fill("TPCclusters2_1_Signal", track->Pt(), track->GetTPCClusterInfo(2,1));
- fQACollection->Fill("TPCclusters2_0_Signal", track->Pt(), track->GetTPCNclsF() > 0 ? track->GetTPCClusterInfo(2,1)/track->GetTPCNclsF() : 0.);
- fQACollection->Fill("TPCncls_Signal", track->Pt(), track->GetTPCNcls());
- fQACollection->Fill("TPCclr_Signal", track->Pt(), track->GetTPCNclsF() > 0 ? static_cast<Double_t>(track->GetTPCNcls())/static_cast<Double_t>(track->GetTPCNclsF()) : 0.);
+ if(fFillNoCuts) {
+ if(signal || !fFillSignalOnly){
+ fVarManager->FillContainer(fContainer, "recTrackContReco", AliHFEcuts::kStepRecNoCut, kFALSE);
+ fVarManager->FillContainer(fContainer, "recTrackContMC", AliHFEcuts::kStepRecNoCut, kTRUE);
}
}
// RecKine: ITSTPC cuts
if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
- // Check TRD criterions (outside the correction framework)
- if(track->GetTRDncls()){
- fQACollection->Fill("chi2TRD", track->GetTRDchi2()/track->GetTRDncls());
- fQACollection->Fill("alpha_rec", track->GetAlpha());
- fQACollection->Fill("pidquality", container[0], track->GetTRDpidQuality());
- fQACollection->Fill("ntrdclusters", container[0], track->GetTRDncls());
- }
-
// RecPrim
if(fRejectKinkMother) {
hfetrack.SetRecTrack(track);
if(HasMCData()) hfetrack.SetMCTrack(mctrack);
hfetrack.SetCentrality(fCentralityF);
+ if(IsPbPb()) hfetrack.SetPbPb();
+ else hfetrack.SetPP();
fPID->SetVarManager(fVarManager);
if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;
nElectronCandidates++;
- fQACollection->Fill("TPCclusters2_1_Selected", track->Pt(), track->GetTPCClusterInfo(2,1));
- fQACollection->Fill("TPCclusters2_0_Selected", track->Pt(), track->GetTPCClusterInfo(2,0));
// Fill Histogram for Hadronic Background
if(HasMCData()){
if(signal) {
// Apply weight for background contamination
if(fBackGroundFactorApply==kTRUE) {
- if(IsPbPb()) fWeightBackGround = fBackGroundFactorArray[(Int_t)fCentralityF]->Eval(TMath::Abs(track->P()));
+ if(IsPbPb() && fCentralityF >= 0) fWeightBackGround = fBackGroundFactorArray[fCentralityF]->Eval(TMath::Abs(track->P()));
else fWeightBackGround = fBackGroundFactorArray[0]->Eval(TMath::Abs(track->P())); // pp case
if(fWeightBackGround < 0.0) fWeightBackGround = 0.0;
if(fSecVtx->Process(track) && signal) {
fVarManager->FillContainer(fContainer, "recTrackContSecvtxReco", AliHFEcuts::kStepHFEcutsSecvtx, kFALSE);
fVarManager->FillContainer(fContainer, "recTrackContSecvtxMC", AliHFEcuts::kStepHFEcutsSecvtx, kTRUE);
- bTagged=kTRUE;
+ bTagged=kTRUE;
}
}
AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
if(!fAOD){
AliError("AOD Event required for AOD Analysis");
- return;
+ return;
}
//
}
if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, track)) return kFALSE;
- fQACollection->Fill("mccharge", signalContainer[3]);
fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGenerated, kFALSE);
signalContainer[4] = 0;
if(fSignalCuts->IsSelected(track)){
signalContainer[5] = 2;
}
fQACollection->Fill("SignalToBackgroundMC", signalContainer);
- fQACollection->Fill("alpha_sim", track->Phi() - TMath::Pi());
// Step GeneratedZOutNoPileUp
- if((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (TMath::Abs(fCentralityF+100.0) < 0.01)) return kFALSE;
+ if((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (fCentralityF < 0)) return kFALSE;
fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGeneratedZOutNoPileUpCentralityFine, kFALSE);
//printf("In ProcessMCtrack %f\n",fCentralityF);
//
// Recover the centrality of the event from ESD or AOD
//
- if(IsAODanalysis()){
-
- AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
- if(!fAOD){
- AliError("AOD Event required for AOD Analysis");
- return kFALSE;
- }
-
- if(IsPbPb()) {
- // Centrality
- AliCentrality *aodCentrality = fAOD->GetCentrality();
- Float_t fCentralityFtemp = aodCentrality->GetCentralityPercentile("V0M");
-
- if( fCentralityFtemp >= 0. && fCentralityFtemp < 10.) fCentralityF = 0;
- else if ( fCentralityFtemp >= 10. && fCentralityFtemp < 20.) fCentralityF = 1;
- else if ( fCentralityFtemp >= 20. && fCentralityFtemp < 30.) fCentralityF = 2;
- else if ( fCentralityFtemp >= 30. && fCentralityFtemp < 40.) fCentralityF = 3;
- else if ( fCentralityFtemp >= 40. && fCentralityFtemp < 50.) fCentralityF = 4;
- else if ( fCentralityFtemp >= 50. && fCentralityFtemp < 60.) fCentralityF = 5;
- else if ( fCentralityFtemp >= 60. && fCentralityFtemp < 90.) fCentralityF = 6;
- else if ( fCentralityFtemp >= 90. && fCentralityFtemp <= 100.) fCentralityF = 7;
- //else if ( fCentralityF_temp >= 90. && fCentralityF_temp < 95.) fCentralityF = 8;
- //else if ( fCentralityF_temp >= 95. && fCentralityF_temp < 90.) fCentralityF = 9;
- //else if ( fCentralityF_temp >= 90. && fCentralityF_temp <=100.) fCentralityF = 10;
- else return kFALSE;
+ Int_t bin = -1;
+ if(IsPbPb()) {
+ // Centrality
+ AliCentrality *centrality = fInputEvent->GetCentrality();
+ Float_t fCentralityFtemp = centrality->GetCentralityPercentile("V0M");
+ Float_t centralityLimits[12] = {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
+ for(Int_t ibin = 0; ibin < 11; ibin++){
+ if(fCentralityFtemp >= centralityLimits[ibin] && fCentralityFtemp < centralityLimits[ibin+1]){
+ bin = ibin;
+ break;
+ }
+ }
+ if(bin == -1) bin = 11; // Overflow
+ } else {
+ // PP: Tracklet multiplicity, use common definition
+ Int_t itsMultiplicity = GetITSMultiplicity(fInputEvent);
+ Int_t multiplicityLimits[8] = {0, 1, 9, 17, 25, 36, 60, 500};
+ for(Int_t ibin = 0; ibin < 7; ibin++){
+ if(itsMultiplicity >= multiplicityLimits[ibin] && itsMultiplicity < multiplicityLimits[ibin + 1]){
+ bin = ibin;
+ break;
+ }
+ }
+ if(bin == -1) bin = 7; // Overflow
+ }
+ fCentralityF = bin;
+ AliDebug(2, Form("Centrality class %d\n", fCentralityF));
- // contributors
- fContributors = 0.5;
- Int_t contributorstemp = 0;
- const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
- if(vtxAOD) contributorstemp = vtxAOD->GetNContributors();
-
- //printf("PbPb contributors_temp %d\n",contributors_temp);
-
- if( contributorstemp <= 0) fContributors = 0.5;
- else fContributors = 1.5;
-
-
-
- }
- else {
- fCentralityF = 0;
- Int_t centralityFtemp = 0;
- const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
- if(vtxAOD) centralityFtemp = vtxAOD->GetNContributors();
-
- //printf("pp centralityF_temp %d\n",centralityF_temp);
-
- if( centralityFtemp <= 0) fCentralityF = 0;
- else if ( centralityFtemp > 0 && centralityFtemp < 2) fCentralityF = 1;
- else if ( centralityFtemp >= 2 && centralityFtemp < 3) fCentralityF = 2;
- else if ( centralityFtemp >= 3 && centralityFtemp < 4) fCentralityF = 3;
- else if ( centralityFtemp >= 4 && centralityFtemp < 5) fCentralityF = 4;
- else if ( centralityFtemp >= 5 && centralityFtemp < 10) fCentralityF = 5;
- else if ( centralityFtemp >= 10 && centralityFtemp < 20) fCentralityF = 6;
- else if ( centralityFtemp >= 20 && centralityFtemp < 30) fCentralityF = 7;
- else if ( centralityFtemp >= 30 && centralityFtemp < 40) fCentralityF = 8;
- else if ( centralityFtemp >= 40 && centralityFtemp < 50) fCentralityF = 9;
- else if ( centralityFtemp >= 50) fCentralityF = 10;
-
- }
-
- return kTRUE;
-
-
- } else {
-
- AliDebug(3, "Processing ESD Centrality");
- AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
- if(!fESD){
- AliError("ESD Event required for ESD Analysis");
- return kFALSE;
- }
- const char* type = fESD->GetBeamType();
-
- if (strstr(type,"Pb-Pb")) {
-
- // Centrality
- AliCentrality *esdCentrality = fESD->GetCentrality();
- Float_t fCentralityFtemp = esdCentrality->GetCentralityPercentile("V0M");
- //printf("PbPb fCentralityF_temp %f\n",fCentralityF_temp);
-
- if( fCentralityFtemp >= 0. && fCentralityFtemp < 10.) fCentralityF = 0;
- else if ( fCentralityFtemp >= 10. && fCentralityFtemp < 20.) fCentralityF = 1;
- else if ( fCentralityFtemp >= 20. && fCentralityFtemp < 30.) fCentralityF = 2;
- else if ( fCentralityFtemp >= 30. && fCentralityFtemp < 40.) fCentralityF = 3;
- else if ( fCentralityFtemp >= 40. && fCentralityFtemp < 50.) fCentralityF = 4;
- else if ( fCentralityFtemp >= 50. && fCentralityFtemp < 60.) fCentralityF = 5;
- else if ( fCentralityFtemp >= 60. && fCentralityFtemp < 80.) fCentralityF = 6;
- else if ( fCentralityFtemp >= 80. && fCentralityFtemp < 90.) fCentralityF = 7;
- else if ( fCentralityFtemp >= 90. && fCentralityFtemp <= 100.) fCentralityF = 8;
- //else if ( fCentralityF_temp >= 80. && fCentralityF_temp < 90.) fCentralityF = 9;
- //else if ( fCentralityF_temp >= 90. && fCentralityF_temp <=100.) fCentralityF = 10;
- else return kFALSE;
-
- // Float_t fCentralityF_temp10 = esdCentrality->GetCentralityClass10("V0M");
- // printf("PbPb fCentralityF_temp %f %f %f \n",fCentralityF_temp, fCentralityF_temp10, fCentralityF);
-
- // contributors
- fContributors = 0.5;
- Int_t contributorstemp = 0;
- const AliESDVertex *vtxESD = fESD->GetPrimaryVertexSPD();
- if(vtxESD) contributorstemp = vtxESD->GetNContributors();
-
- //printf("PbPb contributors_temp %d\n",contributors_temp);
-
- if( contributorstemp <= 0) fContributors = 0.5;
- else fContributors = 1.5;
-
- return kTRUE;
-
- }
-
-
- if (strstr(type,"p-p")) {
- fCentralityF = 0;
- Int_t centralityFtemp = 0;
- const AliESDVertex *vtxESD = fESD->GetPrimaryVertexTracks();
- if(vtxESD && vtxESD->GetStatus()) centralityFtemp = vtxESD->GetNContributors();
-
- //printf("pp centralityF_temp %d\n",centralityF_temp);
-
- if( centralityFtemp <= 0) fCentralityF = 0;
- else if ( centralityFtemp > 0 && centralityFtemp < 2) fCentralityF = 1;
- else if ( centralityFtemp >= 2 && centralityFtemp < 3) fCentralityF = 2;
- else if ( centralityFtemp >= 3 && centralityFtemp < 4) fCentralityF = 3;
- else if ( centralityFtemp >= 4 && centralityFtemp < 5) fCentralityF = 4;
- else if ( centralityFtemp >= 5 && centralityFtemp < 10) fCentralityF = 5;
- else if ( centralityFtemp >= 10 && centralityFtemp < 20) fCentralityF = 6;
- else if ( centralityFtemp >= 20 && centralityFtemp < 30) fCentralityF = 7;
- else if ( centralityFtemp >= 30 && centralityFtemp < 40) fCentralityF = 8;
- else if ( centralityFtemp >= 40 && centralityFtemp < 50) fCentralityF = 9;
- else if ( centralityFtemp >= 50) fCentralityF = 10;
-
- return kTRUE;
-
- }
-
- return kFALSE;
-
- //printf("centrality %f\n",fCentralityF);
-
- }
+ // contributors, to be outsourced
+ const AliVVertex *vtx;
+ if(IsAODanalysis()){
+ AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
+ if(!fAOD){
+ AliError("AOD Event required for AOD Analysis");
+ return kFALSE;
+ }
+ vtx = fAOD->GetPrimaryVertex();
+ } else {
+ AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
+ if(!fESD){
+ AliError("ESD Event required for ESD Analysis");
+ return kFALSE;
+ }
+ vtx = fESD->GetPrimaryVertexSPD();
+ }
+ if(!vtx){
+ fContributors = 0.5;
+ return kFALSE;
+ }
+ else {
+ Int_t contributorstemp = vtx->GetNContributors();
+ if( contributorstemp <= 0) fContributors = 0.5;
+ else fContributors = 1.5;
+ }
+ return kTRUE;
+}
- //printf("centrality %f\n",fCentralityF);
+//___________________________________________________
+Int_t AliAnalysisTaskHFE::GetITSMultiplicity(AliVEvent *ev){
+ //
+ // Definition of the Multiplicity according to the JPSI group (F. Kramer)
+ //
+ Int_t nTracklets = 0;
+ Int_t nAcc = 0;
+ Double_t etaRange = 1.6;
+
+ if (ev->IsA() == AliAODEvent::Class()) {
+ AliAODTracklets *tracklets = ((AliAODEvent*)ev)->GetTracklets();
+ nTracklets = tracklets->GetNumberOfTracklets();
+ for (Int_t nn = 0; nn < nTracklets; nn++) {
+ Double_t theta = tracklets->GetTheta(nn);
+ Double_t eta = -TMath::Log(TMath::Tan(theta/2.0));
+ if (TMath::Abs(eta) < etaRange) nAcc++;
+ }
+ } else if (ev->IsA() == AliESDEvent::Class()) {
+ nTracklets = ((AliESDEvent*)ev)->GetMultiplicity()->GetNumberOfTracklets();
+ for (Int_t nn = 0; nn < nTracklets; nn++) {
+ Double_t eta = ((AliESDEvent*)ev)->GetMultiplicity()->GetEta(nn);
+ if (TMath::Abs(eta) < etaRange) nAcc++;
+ }
+ } else return -1;
+ return nAcc;
}
+
//___________________________________________________
void AliAnalysisTaskHFE::RejectionPileUpVertexRangeEventCut() {
//
AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
if(!fAOD){
AliError("AOD Event required for AOD Analysis");
- return;
+ return;
}
// PileUp
if(fRemovePileUp && fAOD->IsPileupFromSPD()) fIdentifiedAsPileUp = kTRUE;
AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
if(!fESD){
AliError("ESD Event required for ESD Analysis");
- return;
+ return;
}
// PileUp
fIdentifiedAsPileUp = kFALSE;
void SwitchOnPlugin(Int_t plug);
void SetHasMCData(Bool_t hasMC = kTRUE) { SetBit(kHasMCdata, hasMC); };
void SetFillSignalOnly(Bool_t signalOnly) { fFillSignalOnly = signalOnly; }
+ void SetFillNoCuts(Bool_t fillNoCuts) { fFillNoCuts = fillNoCuts; }
void SetRemovePileUp(Bool_t removePileUp) { fRemovePileUp = removePileUp; }
void SetPIDPreselect(AliHFEpid * const cuts) { fPIDpreselect = cuts; };
void SetAODAnalysis() { SetBit(kAODanalysis, kTRUE); };
void ProcessMC();
void ProcessESD();
void ProcessAOD();
+ Int_t GetITSMultiplicity(AliVEvent *ev);
Bool_t PreSelectTrack(AliESDtrack *track) const;
Bool_t ProcessMCtrack(AliVParticle *track);
Bool_t ProcessCutStep(Int_t cutStep, AliVParticle *track);
ULong_t fQAlevel; // QA level
UShort_t fPlugins; // Enabled Plugins
Bool_t fFillSignalOnly; // Fill container only with MC Signal Tracks
+ Bool_t fFillNoCuts; // Fill container before any cut
Bool_t fBackGroundFactorApply; // Apply Background Function Subtraction
Bool_t fRemovePileUp; // Remove Pile Up
Bool_t fIdentifiedAsPileUp; // Identified as pile-up
Bool_t fHasSpecialTriggerSelection; // Select special triggered events
Bool_t fRejectKinkMother; // Reject Kink Mother
TString fSpecialTrigger; // Special trigger selection
- Float_t fCentralityF; // Centrality
+ Int_t fCentralityF; // Centrality
Float_t fContributors; // Contributors
Double_t fWeightBackGround; // weight background function
Double_t fVz; // z position of the primary vertex
, fOutput(NULL)
, fEvents(NULL)
, fNNref(NULL)
+ , fTRDTotalChargeInSlice0(kFALSE)
{
//
// Default Constructor
, fOutput(NULL)
, fEvents(NULL)
, fNNref(NULL)
+ , fTRDTotalChargeInSlice0(kFALSE)
{
//
// Default Constructor
// Initialize PID QA
//
fOutput = new TList;
+ fOutput->SetOwner();
// Counter for number of events
fOutput->Add((fEvents = new TH1I("nEvents", "NumberOfEvents", 1, 1, 2)));
fPIDqa = new AliHFEpidQA;
+ if(fTRDTotalChargeInSlice0) fPIDqa->SetTRDTotalChargeInSlice0();
if(HasV0pidQA()) fPIDqa->SetV0pidQA();
if(HasRecalculateTRDpid()) fPIDqa->SetRecalculateTRDpid();
if(fNNref) fPIDqa->SetNNref(fNNref);
Bool_t HasRecalculateTRDpid() const { return TestBit(kRecalculateTRDpid); };
void SetV0pidQA(Bool_t v0pidQA = kTRUE) { SetBit(kV0pidQA, v0pidQA); };
void SetRecalculateTRDpid(Bool_t recal = kTRUE) { SetBit(kRecalculateTRDpid, recal); };
+ void SetTRDTotalChargeInSlice0() { fTRDTotalChargeInSlice0 = kTRUE; }
void SetNNref(TFile *f) { fNNref = f; };
TList *fOutput; //! Container for output histos
TH1 *fEvents; //! Number of Events
TFile *fNNref; // reference file for NN
+ Bool_t fTRDTotalChargeInSlice0; // Fix for Foreware/Backward compatibility
ClassDef(AliAnalysisTaskHFEpidQA, 1)
};
#include "AliLog.h"
#include "AliExternalTrackParam.h"
-#include "AliHFEcollection.h"
-
#include "AliHFEV0cuts.h"
ClassImp(AliHFEV0cuts)
, fCurrentV0id(0)
, fPdaughterPDG(0)
, fNdaughterPDG(0)
+ , fDestBits(0)
{
//
//
// destructor
//
- if (fQA) delete fQA;
- if (fQAmc) delete fQAmc;
+
+ if (fQA && TESTBIT(fDestBits, kBitQA)) delete fQA;
+ if (fQAmc && TESTBIT(fDestBits, kBitQAmc)) delete fQAmc;
}
//________________________________________________________________
, fCurrentV0id(0)
, fPdaughterPDG(0)
, fNdaughterPDG(0)
+ , fDestBits(0)
{
//
// Copy constructor
// [1] jus before the cut on given variable was applied, but after all the previous cuts
//
+ memset(&fDestBits, 0, sizeof(UInt_t));
+ // now set the first two bits to 1
+ SETBIT(fDestBits, kBitQA);
+ SETBIT(fDestBits, kBitQAmc);
+
+
fQA = new AliHFEcollection("fQA", name);
fQAmc = new AliHFEcollection("fQAmc", name);
}
}
- // loose cuts first
- //if(LooseRejectK0(v0) || LooseRejectLambda(v0)) return kFALSE;
-
AliVTrack* daughter[2];
Int_t pIndex = 0, nIndex = 0;
if(CheckSigns(v0)){
// possible new cuts
//
// separation cut at the entrance to the TPC
- const Double_t cutSeparation = 0.; // ORG 3.0 cm
+ const Double_t cutSeparation = 1.; // ORG 3.0 cm
fQAmc->Fill("h_Mass_L_as_K0", v0->GetEffMass(2, 0));
}
}
- // loose cuts first
- //if(LooseRejectK0(v0) || LooseRejectGamma(v0)) return kFALSE;
-
+
const Double_t cL0mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass(); // PDG lambda mass
AliVTrack* daughter[2];
fQAmc->Fill("h_alambda_p_B", iP);
}
//
- if(4 == TMath::Abs(fCurrentV0id)){
- fQAmc->Fill("h_ProtonL_P_S", 5, p[ixMC[0]]);
- fQAmc->Fill("h_PionL_P_S", 5, p[ixMC[1]]);
- }
- else if(-2 != fCurrentV0id){
- fQAmc->Fill("h_ProtonL_P_B", 5, p[ixMC[0]]);
- fQAmc->Fill("h_PionL_P_B", 5, p[ixMC[1]]);
+ if(fMCEvent){
+ if(4 == TMath::Abs(fCurrentV0id)){
+ fQAmc->Fill("h_ProtonL_P_S", 5, p[ixMC[0]]);
+ fQAmc->Fill("h_PionL_P_S", 5, p[ixMC[1]]);
+ }
+ else if(-2 != fCurrentV0id){
+ fQAmc->Fill("h_ProtonL_P_B", 5, p[ixMC[0]]);
+ fQAmc->Fill("h_PionL_P_B", 5, p[ixMC[1]]);
+ }
}
-
return kTRUE;
}
//________________________________________________________________
-Double_t AliHFEV0cuts::OpenAngle(AliESDv0 *v0) const {
+Double_t AliHFEV0cuts::OpenAngle(AliESDv0 const *v0){
//
// Opening angle between two daughter tracks
//
return TMath::Abs(openAngle);
}
//________________________________________________________________
-Double_t AliHFEV0cuts::PsiPair(AliESDv0 *v0) {
+Double_t AliHFEV0cuts::PsiPair(const AliESDv0 *v0) {
//
// Angle between daughter momentum plane and plane
//
return psiPair;
}
//________________________________________________________________
-AliKFParticle *AliHFEV0cuts::CreateMotherParticle(AliVTrack* const pdaughter, AliVTrack* const ndaughter, Int_t pspec, Int_t nspec){
+AliKFParticle *AliHFEV0cuts::CreateMotherParticle(AliVTrack const *pdaughter, AliVTrack const *ndaughter, Int_t pspec, Int_t nspec){
//
// Creates a mother particle
//
return m;
}
-//_________________________________________________
-Bool_t AliHFEV0cuts::LooseRejectK0(AliESDv0 * const v0) const {
- //
- // Reject K0 based on loose cuts
- //
- Double_t mass = v0->GetEffMass(AliPID::kPion, AliPID::kPion);
- if(mass > 0.494 && mass < 0.501) return kTRUE;
- return kFALSE;
-}
-
-//_________________________________________________
-Bool_t AliHFEV0cuts::LooseRejectLambda(AliESDv0 * const v0) const {
- //
- // Reject Lambda based on loose cuts
- //
- Double_t mass1 = v0->GetEffMass(AliPID::kPion, AliPID::kProton);
- Double_t mass2 = v0->GetEffMass(AliPID::kProton, AliPID::kPion);
-
- if(mass1 > 1.1 && mass1 < 1.12) return kTRUE;
- if(mass2 > 1.1 && mass2 < 1.12) return kTRUE;
- return kFALSE;
-}
-
-//_________________________________________________
-Bool_t AliHFEV0cuts::LooseRejectGamma(AliESDv0 * const v0) const {
- //
- // Reject Lambda based on loose cuts
- //
-
- Double_t mass = v0->GetEffMass(AliPID::kElectron, AliPID::kElectron);
-
- if(mass < 0.02) return kTRUE;
- return kFALSE;
-}
//___________________________________________________________________
-void AliHFEV0cuts::Armenteros(AliESDv0 *v0, Float_t val[2]){
+void AliHFEV0cuts::Armenteros(const AliESDv0 *v0, Float_t val[2]){
//
// computes the Armenteros variables for given V0
// fills the histogram
}
//___________________________________________________________________
-Bool_t AliHFEV0cuts::CheckSigns(AliESDv0* const v0){
+Bool_t AliHFEV0cuts::CheckSigns(AliESDv0 const *v0){
//
// check wheter the sign was correctly applied to
// V0 daughter tracks
class AliHFEV0cuts : public TObject {
public:
+ enum{
+ kBitQA = 1,
+ kBitQAmc = 2
+ };
enum{ // Reconstructed V0
kUndef = 0,
kRecoGamma = 1,
void SetInputEvent(AliVEvent* const e) { fInputEvent = e; };
void SetPrimaryVertex(AliKFVertex* const v) { fPrimaryVertex = v; };
- TList* GetList() { return fQA->GetList(); };
- TList* GetListMC() { return fQAmc->GetList(); };
+ TList* GetList() {
+ CLRBIT(fDestBits, kBitQA);
+ return fQA->GetList();
+ };
+ TList* GetListMC() {
+ CLRBIT(fDestBits, kBitQAmc);
+ return fQAmc->GetList();
+ };
Bool_t TrackCutsCommon(AliESDtrack* track);
Bool_t V0CutsCommon(AliESDv0 *v0);
Bool_t K0Cuts(AliESDv0 *v0);
Bool_t LambdaCuts(AliESDv0 *v0, Bool_t &isLambda);
- Bool_t LooseRejectK0(AliESDv0 * const v0) const;
- Bool_t LooseRejectLambda(AliESDv0 * const v0) const;
- Bool_t LooseRejectGamma(AliESDv0 * const v0) const;
-
- void Armenteros(AliESDv0 *v0, Float_t val[2]);
+ void Armenteros(const AliESDv0 *v0, Float_t val[2]);
- Double_t OpenAngle(AliESDv0 *v0) const;//opening angle between V0 daughters; close to zero for conversions
- Double_t PsiPair(AliESDv0 *v0);
+ Double_t OpenAngle(AliESDv0 const *v0);//opening angle between V0 daughters; close to zero for conversions
+ Double_t PsiPair(const AliESDv0 *v0);
- Bool_t CheckSigns(AliESDv0* const v0);
+ Bool_t CheckSigns(AliESDv0 const *v0);
Bool_t GetConvPosXY(AliESDtrack * const ptrack, AliESDtrack * const ntrack, Double_t convpos[2]);
Bool_t GetHelixCenter(AliESDtrack * const track, Double_t b,Int_t charge, Double_t center[2]);
void SetCurrentV0id(Int_t id) { fCurrentV0id = id; };
void SetDaughtersID(Int_t d[2]) {fPdaughterPDG = d[0]; fNdaughterPDG = d[1]; };
- AliKFParticle *CreateMotherParticle(AliVTrack* const pdaughter, AliVTrack* const ndaughter, Int_t pspec, Int_t nspec);
+ AliKFParticle *CreateMotherParticle(AliVTrack const *pdaughter, AliVTrack const *ndaughter, Int_t pspec, Int_t nspec);
private:
void Copy(TObject &ref) const;
Int_t fPdaughterPDG; // MC id of the positive daugeter
Int_t fNdaughterPDG; // MC id of the negative daugeter
-
+ UInt_t fDestBits; // status bits for destructor
+
ClassDef(AliHFEV0cuts, 1)
};
#include "AliHFEV0pid.h"
ClassImp(AliHFEV0pid)
-//____________________________________________________________
AliHFEV0pid::AliHFEV0pid():
- TObject()
+ TNamed()
, fInputEvent(NULL)
, fNtracks(0)
, fMCEvent(NULL)
, fQA(NULL)
, fV0cuts(NULL)
, fOutput(NULL)
+ , fDestBits(0)
+
{
//
// Default constructor
//
+}
+//____________________________________________________________
+AliHFEV0pid::AliHFEV0pid(const char *name):
+ TNamed(name, "")
+ , fInputEvent(NULL)
+ , fNtracks(0)
+ , fMCEvent(NULL)
+ , fMCon(kFALSE)
+ , fPrimaryVertex(NULL)
+ , fElectrons(NULL)
+ , fPionsK0(NULL)
+ , fPionsL(NULL)
+ , fKaons(NULL)
+ , fProtons(NULL)
+ , fGammas(NULL)
+ , fK0s(NULL)
+ , fLambdas(NULL)
+ , fAntiLambdas(NULL)
+ , fIndices(NULL)
+ , fQA(NULL)
+ , fV0cuts(NULL)
+ , fOutput(NULL)
+ , fDestBits(0)
+{
+ //
+ // Standard constructor
+ //
fElectrons = new TObjArray();
fPionsK0 = new TObjArray();
fPionsL = new TObjArray();
fIndices = new AliHFEV0pidTrackIndex();
}
-
//____________________________________________________________
+
AliHFEV0pid::~AliHFEV0pid(){
//
// Destructor
if(fAntiLambdas) delete fAntiLambdas;
if(fIndices) delete fIndices;
- if(fQA) delete fQA;
if(fV0cuts) delete fV0cuts;
- if(fOutput) delete fOutput;
+
+ if(TESTBIT(fDestBits, 1)){
+ if(fQA) delete fQA;
+ if(fOutput) delete fOutput;
+ }
}
//____________________________________________________________
// Initialize QA histograms
//
+ memset(&fDestBits, 0, sizeof(fDestBits));
+ SETBIT(fDestBits, 1);
+
fOutput = new TList();
+ fOutput->SetOwner();
fV0cuts = new AliHFEV0cuts();
fV0cuts->Init("V0cuts");
}
//____________________________________________________________
-void AliHFEV0pid::AliHFEV0pidTrackIndex::Flush(){
+const void AliHFEV0pid::AliHFEV0pidTrackIndex::Flush(){
//
// Reset containers
//
// Getter for V0 PID QA histograms
//
+ CLRBIT(fDestBits, 1);
+
TList *tmp = fV0cuts->GetList();
tmp->SetName("V0cuts");
fOutput->Add(tmp);
#ifndef ALIHFEV0PID_H
#define ALIHFEV0PID_H
-#ifndef ROOT_TObject
-#include <TObject.h>
+#ifndef ROOT_TNamed
+#include <TNamed.h>
#endif
class TObjArray;
class AliHFEV0cuts;
class AliHFEcollection;
-class AliHFEV0pid : public TObject{
+class AliHFEV0pid : public TNamed{
public:
AliHFEV0pid();
+ AliHFEV0pid(const char *name);
~AliHFEV0pid();
void Process(AliVEvent * const inputEvent);
Int_t GetNumberOfPionsL() const { return fNPionsL; };
Int_t GetNumberOfKaons() const { return fNKaons; };
Int_t GetNumberOfProtons() const { return fNProtons; };
- void Flush();
+ const void Flush();
private:
AliHFEV0pidTrackIndex(const AliHFEV0pidTrackIndex &ref);
AliHFEV0cuts *fV0cuts; // separate class for studying and applying the V0 cuts
TList *fOutput; // collection list
+ UInt_t fDestBits; // logical bits for destructor
+
ClassDef(AliHFEV0pid, 1) // V0 PID Class
};
AliHFEV0pidMC::AliHFEV0pidMC():
fMC(0x0)
, fColl(NULL)
- , fV0cuts(NULL)
+ , fDestBits(0)
{
//
// default constructor
// destructor
//
if (fColl) delete fColl;
- //if (fV0cuts) delete fV0cuts;
}
//____________________________________________________________
void AliHFEV0pidMC::Init(){
//
// initialize objects
//
+
+ memset(&fDestBits, 0, sizeof(UInt_t));
+ SETBIT(fDestBits, 1);
fColl = new AliHFEcollection("V0pidMC", "MC based V0 benchmarking");
// QA
//
// Get QA histograms
//
+
+ CLRBIT(fDestBits, 1);
if(fColl)
return fColl->GetList();
return NULL;
class AliMCEvent;
-class AliHFEV0cuts;
class AliHFEcollection;
class AliHFEV0pidMC : public TObject {
AliMCEvent* fMC; // MC event
AliHFEcollection* fColl; // Histogram collection
- AliHFEV0cuts* fV0cuts; // V0 cut class
+
+ UInt_t fDestBits; // status bits for the destructor
ClassDef(AliHFEV0pidMC, 1) // QA class for V0 PID
};
// Matus Kalisky <matus.kalisky@cern.ch>
//
-//#include <iostream>
-
#include <TH1F.h>
#include <TH2F.h>
#include <THnSparse.h>
#include <TMath.h>
#include "AliLog.h"
+
#include "AliHFEcollection.h"
using namespace std;
// Clone List Content
target.fList = new THashList();
+ target.fList->SetOwner();
for(Int_t ien = 0; ien < fList->GetEntries(); ien++)
target.fList->Add(fList->At(ien)->Clone());
}
}
}
}
+//____________________________________________________________________
+void AliHFEcollection::Print(Option_t *) const{
+ //
+ // Print content of the collection
+ //
+ TIter histIter(fList);
+ TObject *o = NULL;
+ Int_t nHistos = 0;
+ printf("Collection %s\n", GetName());
+ printf("Content of the collection:\n=========================================\n");
+ while((o = histIter())){
+ printf("Histo %s, Type %s\n", o->GetName(), o->IsA()->GetName());
+ nHistos++;
+ }
+ printf("Number of histos in the collection: %d\n", nHistos);
+ printf("\n");
+}
+
Long64_t Merge(TCollection *list);
+ virtual void Print(Option_t *) const;
// Get functions
TList* GetList() const { return fList; }
SetHFElectronTRDCuts();
SetHFElectronDcaCuts();
+ // Connect the event cuts
+ SetEventCutList(kEventStepGenerated);
+ SetEventCutList(kEventStepReconstructed);
+
+
}
//__________________________________________________________________
void AliHFEcuts::SetParticleGenCutList(){
//
// Initialize Particle Cuts for Monte Carlo Tracks
- // Production Vertex: < 1cm (Beampipe)
+ // Production Vertex Radius: < 3cm
// Particle Species: Electrons
- // Eta: < 0.9 (TRD-TOF acceptance)
+ // Eta: < 0.8
//
TObjArray *mcCuts = new TObjArray;
// Checks the cuts without using the correction framework manager
//
AliDebug(2, "Called\n");
- TString stepnames[kNcutStepsMCTrack + kNcutStepsRecTrack + kNcutStepsDETrack + kNcutStepsSecvtxTrack + 1] = {"fPartGenCuts","fPartEvCutPileupZ","fPartEvCut","fPartAccCuts","fPartRecNoCuts","fPartRecKineITSTPCCuts", "fPartPrimCuts", "fPartHFECutsITS","fPartHFECutsTRD","fPartHFECutsDca", "fPartHFECutsSecvtx"};
+ TString stepnames[kNcutStepsMCTrack + kNcutStepsRecTrack + kNcutStepsDETrack + kNcutStepsSecvtxTrack + 1] = {"fPartGenCuts","fPartEvCutPileupZ","fPartEvCut","fPartAccCuts","fPartRecNoCuts","fPartRecKineITSTPCCuts", "fPartPrimCuts", "fPartHFECutsITS","fPartHFECutsTOF","fPartHFECutsTRD","fPartHFECutsDca", "fPartHFECutsSecvtx"};
AliDebug(2, Form("Doing cut %s", stepnames[step].Data()));
TObjArray *cuts = dynamic_cast<TObjArray *>(fCutList->FindObject(stepnames[step].Data()));
if(!cuts) return kTRUE;
}
return status;
}
+
+
+//__________________________________________________________________
+Bool_t AliHFEcuts::CheckEventCuts(const char*namestep, TObject *o){
+ //
+ // Checks the cuts without using the correction framework manager
+ //
+ AliDebug(2, "Called\n");
+ TObjArray *cuts = dynamic_cast<TObjArray *>(fCutList->FindObject(namestep));
+ if(!cuts) return kTRUE;
+ TIter it(cuts);
+ AliCFCutBase *mycut;
+ Bool_t status = kTRUE;
+ while((mycut = dynamic_cast<AliCFCutBase *>(it()))){
+ status &= mycut->IsSelected(o);
+ }
+ return status;
+}
void Initialize();
Bool_t CheckParticleCuts(UInt_t step, TObject *o);
+ Bool_t CheckEventCuts(const char*namestep, TObject *o);
TList *GetQAhistograms() const { return fHistQA; }
#include <THnSparse.h>
#include <TString.h>
+#include <TMath.h>
+#include "AliESDInputHandler.h"
+//#include "AliVCluster.h"
+//#include "AliVCaloCells.h"
+//#include "AliVEvent.h"
+#include "AliMagF.h"
+
#include "AliLog.h"
#include "AliPID.h"
#include "AliVParticle.h"
-#include "AliVTrack.h"
-#include "AliESDtrack.h"
+//#include "AliVTrack.h"
+//#include "AliESDtrack.h"
#include "AliHFEcollection.h"
#include "AliHFEpid.h"
#include "AliHFEpidBase.h"
#include "AliHFEtools.h"
#include "AliHFEemcalPIDqa.h"
-ClassImp(AliHFEpidEMCAL)
+ClassImp(AliHFEemcalPIDqa)
//_________________________________________________________
AliHFEemcalPIDqa::AliHFEemcalPIDqa():
const Double_t kTPCSigMax = 140.;
// 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, Centrality)
- Int_t nBins[4] = {20, 20, 400, kCentralityBins};
- Double_t min[4] = {kMinP, kMinP, kTPCSigMim, 0};
- Double_t max[4] = {kMaxP, kMaxP, kTPCSigMax, 11.};
- fHistos->CreateTHnSparse("EMCAL_TPCdedx", "EMCAL signal; species; p [GeV/c]; pT [GeV/c] ; TPC signal [a.u.]; Centrality", 4, nBins, min, max);
-
+ Int_t nBins[6] = {AliPID::kSPECIES + 1, 20, 20, 400, kCentralityBins, 2};
+ Double_t min[6] = {-1, kMinP, kMinP, kTPCSigMim, 0, 0.};
+ Double_t max[6] = {AliPID::kSPECIES, kMaxP, kMaxP, kTPCSigMax, 11., 2.};
+ fHistos->CreateTHnSparse("EMCAL_TPCdedx", "EMCAL signal; species; p [GeV/c]; pT [GeV/c] ; TPC signal [a.u.]; Centrality; PID Step", 6, nBins, min, max);
+
+ //2nd histogram: EMCAL signal - E/p and Rmatch
+ Int_t nBins2[6] = {AliPID::kSPECIES + 1, 40, 40, 1000, 250, 2};
+ Double_t min2[6] = {-1, kMinP, kMinP, 0, 0, 0.};
+ Double_t max2[6] = {AliPID::kSPECIES, kMaxP, kMaxP, 10, 1., 2.};
+ fHistos->CreateTHnSparse("EMCAL_Signal", "EMCAL true signal; species; p [GeV/c]; pT [GeV/c] ; E/p; Rmatch; PID Step", 6, nBins2, min2, max2);
+
}
+
+
+
//_________________________________________________________
-void AliHFEemcalPIDqa::ProcessTrack(const AliHFEpidObject *track,AliHFEdetPIDqa::EStep_t /*step*/){
+void AliHFEemcalPIDqa::ProcessTrack(const AliHFEpidObject *track,AliHFEdetPIDqa::EStep_t step){
//
// Fill TPC histograms
//
//AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
Float_t centrality = track->GetCentrality();
+ //const AliVTrack *vtrack = dynamic_cast<const AliVTrack *>(track->GetRecTrack());
+ //const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(vtrack);
const AliESDtrack *esdtrack = static_cast<const AliESDtrack *>(track->GetRecTrack());
- Double_t contentSignal[4];
- contentSignal[0] = track->GetRecTrack()->P();
- contentSignal[1] = track->GetRecTrack()->Pt();
- contentSignal[2] = esdtrack->GetTPCsignal(); //?
- contentSignal[3] = centrality;
+ Int_t species = track->GetAbInitioPID();
+ if(species >= AliPID::kSPECIES) species = -1;
+
+ Double_t contentSignal[6];
+ contentSignal[0] = species;
+ contentSignal[1] = track->GetRecTrack()->P();
+ contentSignal[2] = track->GetRecTrack()->Pt();
+ contentSignal[3] = esdtrack->GetTPCsignal(); //?
+ contentSignal[4] = centrality;
+ contentSignal[5] = step == AliHFEdetPIDqa::kBeforePID ? 0. : 1.;
+
+ TVector3 emcsignal = MomentumEnergyMatchV2(esdtrack);
+
+
+ Double_t contentSignal2[6];
+ contentSignal2[0] = species;
+ contentSignal2[1] = track->GetRecTrack()->P();
+ contentSignal2[2] = track->GetRecTrack()->Pt();
+ contentSignal2[3] = emcsignal(0);
+ contentSignal2[4] = emcsignal(1);
+ contentSignal2[5] = contentSignal[5];
+
//printf("ProcessTrack ; Print Content %g; %g; %g; %g \n",contentSignal[0],contentSignal[1],contentSignal[2],contentSignal[3]);
fHistos->Fill("EMCAL_TPCdedx", contentSignal);
+ fHistos->Fill("EMCAL_Signal", contentSignal2);
}
//_________________________________________________________
return dynamic_cast<TH1 *>(histo);
}
+
+//___________________________________________________________________________
+TVector3 AliHFEemcalPIDqa::MomentumEnergyMatchV2(const AliESDtrack *esdtrack) const
+{ // Returns e/p & Rmatch
+
+ Float_t clsPos[3];
+ Double_t trkPos[3];
+ Double_t matchclsE = 9999.9;
+ TVector3 refVec(-9999,-9999,-9999);
+
+ AliESDEvent *evt = esdtrack->GetESDEvent();
+
+ //Int_t icl = esdtrack->GetEMCALcluster();
+ Int_t icl = (const_cast<AliESDtrack *>(esdtrack))->GetEMCALcluster();
+
+ AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
+ if(!cluster->IsEMCAL()) return refVec;
+
+ cluster->GetPosition(clsPos);
+ esdtrack->GetXYZ(trkPos);
+
+ TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
+ TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
+
+
+ Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi(); // track cluster matching
+ Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
+
+ double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
+
+ matchclsE = cluster->E();
+
+ //double feop = -9999.9;
+ //if(matchclsE<9999)
+ double feop = matchclsE/esdtrack->P();
+
+ // if(feop!=-9999.9)printf("%f\n",feop) ;
+ TVector3 emcsignal(feop,rmatch,0);
+
+ return emcsignal;
+
+}
+
+
+
+
+/*
+//___________________________________________________________________________
+TVector3 AliHFEemcalPIDqa::MomentumEnergyMatchV1(const AliESDtrack *esdtrack) const
+{ // Returns e/p & Rmatch
+
+ Float_t clsPos[3];
+ Double_t trkPos[3];
+ Double_t Rmatch;
+ Double_t matchclsE = 9999.9;
+ TVector3 refVec(-9999,-9999,-9999);
+
+ AliESDEvent *evt = esdtrack->GetESDEvent();
+ Double_t magF = evt->GetMagneticField();
+ Double_t magSign = 1.0;
+ if(magF<0)magSign = -1.0;
+ //printf("magF ; %g ; %g \n", magF,magSign);
+
+ if (!TGeoGlobalMagField::Instance()->GetField()) {
+ printf("Loading field map...\n");
+ //AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
+ AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG); // for 10d
+ TGeoGlobalMagField::Instance()->SetField(field);
+ }
+
+ AliEMCALTrack *emctrack = new AliEMCALTrack(*esdtrack);
+ Double_t fieldB[3];
+ emctrack->GetBxByBz(fieldB);
+ //printf("%g %g %g \n", fieldB[0], fieldB[1], fieldB[2]);
+
+ for(Int_t icl=0; icl<evt->GetNumberOfCaloClusters(); icl++)
+ {
+
+ AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
+ if(!cluster->IsEMCAL()) return refVec;
+
+ cluster->GetPosition(clsPos);
+ if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) return refVec;
+ emctrack->GetXYZ(trkPos);
+
+ TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
+ TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
+
+ Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi(); // track cluster matching
+ Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta(); // track cluster matching
+
+ double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
+
+ if(rmatch<0.02)
+ {
+ matchclsE = cluster->E();
+ Rmatch = rmatch;
+ }
+ }
+ delete emctrack;
+
+ //double feop = -9999.9;
+ //if(matchclsE<9999)
+ double feop = matchclsE/esdtrack->P();
+
+ // if(feop!=-9999.9)printf("%f\n",feop) ;
+ TVector3 emcsignal(feop, Rmatch, 0);
+
+ return emcsignal;
+
+}
+*/
+
class AliESDtrack;
class AliAODTrack;
+class TVector3;
+
+
class AliHFEemcalPIDqa : public AliHFEdetPIDqa{
public:
AliHFEemcalPIDqa();
TH1 *GetHistogram(const char *name);
AliHFEcollection *GetHistoCollection() const { return fHistos; }
+ TVector3 MomentumEnergyMatchV2(const AliESDtrack *esdtrack) const;
+ //TVector3 MomentumEnergyMatchV1(const AliESDtrack *esdtrack) const;
+
+
protected:
void ProcessESDtrack(const AliESDtrack *track, AliHFEdetPIDqa::EStep_t step, Int_t species);
void ProcessAODtrack(const AliAODTrack *track, AliHFEdetPIDqa::EStep_t step, Int_t species);
if(fCheck && !(statusL0 || statusL1))
SETBIT(survivedCut, kPixelITS);
break;
+ case kExclusiveSecond:
+ AliDebug(2, "Exlusive second");
+ if(fCheck){ // Cut out tracks which pass a dead ITS Layer 0
+ if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0) && statusL0)
+ SETBIT(survivedCut, kPixelITS);
+ } else {
+ if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0))
+ SETBIT(survivedCut, kPixelITS);
+ }
+ break;
default:
AliDebug(2, "None");
break;
kSecond = 1,
kBoth = 2,
kNone = 3,
- kAny = 4
+ kAny = 4,
+ kExclusiveSecond = 5
} ITSPixel_t;
typedef enum{
kFound = 0,
//_______________________________________________________________________________________________
void AliHFEmcQA::SetBackgroundWeightFactor(Double_t *elecBackgroundFactor, Double_t *binLimit)
{
+ //
+ // copy background weighting factors into data member
+ //
memcpy(fElecBackgroundFactor,elecBackgroundFactor,sizeof(Double_t) * kElecBgSpecies * kBgPtBins);
memcpy(fBinLimit,binLimit,sizeof(Double_t) * (kBgPtBins+1));
/*
}
//____________________________________________________________
-Bool_t AliHFEpid::IsSelected(AliHFEpidObject *track, AliHFEcontainer *cont, const Char_t *contname, AliHFEpidQAmanager *pidqa){
+Bool_t AliHFEpid::IsSelected(const AliHFEpidObject * const track, AliHFEcontainer *cont, const Char_t *contname, AliHFEpidQAmanager *pidqa){
//
// Select Tracks
//
}
//____________________________________________________________
-void AliHFEpid::ConfigureTPCrejection(const char *lowerCutParam, Double_t * const params, Float_t upperTPCCut, Float_t TOFCut){
+void AliHFEpid::ConfigureTOF(Float_t TOFCut){
//
- // Combined TPC-TOF PID
- // if no function parameterizaion is given, then the default one (exponential) is chosen
+ // Set Number of sigmas for TOF PID
//
- if(HasMCData()) AliInfo("Configuring TPC for MC\n");
- AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fDetectorPID[kTOFpid]);
if(tofpid) tofpid->SetTOFnSigma(TOFCut);
+}
+
+//____________________________________________________________
+void AliHFEpid::ConfigureTPCcentralityCut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
+ //
+ // Cofigure centrality dependent cut function for TPC PID
+ //
+ ConfigureTPCcut(centralityBin, lowerCutParam, params, upperTPCCut);
+}
+//____________________________________________________________
+void AliHFEpid::ConfigureTPCdefaultCut(const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
+ //
+ // Cofigure default cut function for TPC PID
+ //
+ ConfigureTPCcut(-1, lowerCutParam, params, upperTPCCut);
+}
+
+//____________________________________________________________
+void AliHFEpid::ConfigureTPCcut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
+ //
+ // Cofigure cut function for TPC PID
+ // if no function parameterizaion is given, then the default one (exponential) is chosen
+ //
+
+ if(HasMCData()) AliInfo("Configuring TPC for MC\n");
+ AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
//TF1 *upperCut = new TF1("upperCut", "[0] * TMath::Exp([1]*x)", 0, 20);
- TF1 *upperCut = new TF1("upperCut", "[0]", 0, 20); // Use constant upper cut
- TF1 *lowerCut = new TF1("lowerCut", lowerCutParam == NULL ? "[0] * TMath::Exp([1]*x) + [2]": lowerCutParam, 0, 20);
+ TF1 *upperCut = new TF1(Form("upperCut%s", centralityBin < 0 ? "Default" : Form("Bin%d", centralityBin)), "[0]", 0, 20); // Use constant upper cut
+ TF1 *lowerCut = new TF1(Form("lowerCut%s", centralityBin < 0 ? "Default" : Form("Bin%d", centralityBin)), lowerCutParam == NULL ? "[0] * TMath::Exp([1]*x) + [2]": lowerCutParam, 0, 20);
upperCut->SetParameter(0, upperTPCCut); // pp
-// upperCut->SetParameter(0, 3.5); //PbPb
-// printf("upper %f \n",upperCut->GetParameter(0));
- //upperCut->SetParameter(0, 2.7);
- //upperCut->SetParameter(1, -0.4357);
if(params){
- for(Int_t ipar = 0; ipar < lowerCut->GetNpar(); ipar++) lowerCut->SetParameter(ipar, params[ipar]);
+ for(Int_t ipar = 0; ipar < lowerCut->GetNpar(); ipar++)
+ {
+ lowerCut->SetParameter(ipar, params[ipar]);
+ // printf("printout %i %s %f \n", centralityBin, lowerCutParam, params[ipar]);
+ }
} else {
// Set default parameterization
- if(HasMCData()) lowerCut->SetParameter(0, -2.5);
- else lowerCut->SetParameter(0, -4.03); //pp
-// else lowerCut->SetParameter(0, -3.83); //PbPb
-
-
- //else lowerCut->SetParameter(0, -3.71769);
- //else lowerCut->SetParameter(0, -3.7);
-
+ if(HasMCData()) lowerCut->SetParameter(0, -2.5);
+ else lowerCut->SetParameter(0, -4.03); //pp
lowerCut->SetParameter(1, -0.22); // pp
-// lowerCut->SetParameter(1,-0.36 ); // PbPb
- //lowerCut->SetParameter(1, -0.40263);
- //lowerCut->SetParameter(1, -0.8);
if(HasMCData()) lowerCut->SetParameter(2, -2.2);
else lowerCut->SetParameter(2, 0.92); //pp
-// else lowerCut->SetParameter(2, 0.27); //PbPb
- // else lowerCut->SetParameter(2, 0.92); //pp
-// printf("lower0 %f \n",lowerCut->GetParameter(0));
-// printf("lower1 %f \n",lowerCut->GetParameter(1));
-// printf("lower2 %f \n",lowerCut->GetParameter(2));
- //else lowerCut->SetParameter(2, 0.267857);
- //else lowerCut->SetParameter(2, -0.35);
}
if(tpcpid){
tpcpid->SetTPCnSigma(2);
- tpcpid->SetUpperSigmaCut(upperCut);
- tpcpid->SetLowerSigmaCut(lowerCut);
+ if(centralityBin < 0){
+ tpcpid->SetUpperSigmaCutDefault(upperCut);
+ tpcpid->SetLowerSigmaCutDefault(lowerCut);
+ } else {
+ tpcpid->SetUpperSigmaCutCentrality(upperCut, centralityBin);
+ tpcpid->SetLowerSigmaCutCentrality(lowerCut, centralityBin);
+ }
}
AddCommonObject(upperCut);
AddCommonObject(lowerCut);
}
-//____________________________________________________________
-void AliHFEpid::ConfigureTPCstrategyParis(){
- //
- // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
- //
- AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
- if(pid){
- pid->SetTPCnSigma(2);
- pid->SetRejectParticle(AliPID::kProton, 0., -3., 10., 3.);
- pid->SetRejectParticle(AliPID::kKaon, 0., -3., 10., 3.);
- }
-}
-
//____________________________________________________________
void AliHFEpid::PrintStatus() const {
//
~AliHFEpid();
Bool_t InitializePID();
- Bool_t IsSelected(AliHFEpidObject *track, AliHFEcontainer *cont = NULL, const Char_t *contname = "trackContainer", AliHFEpidQAmanager *qa = NULL);
+ Bool_t IsSelected(const AliHFEpidObject * const track, AliHFEcontainer *cont = NULL, const Char_t *contname = "trackContainer", AliHFEpidQAmanager *qa = NULL);
Bool_t HasMCData() const { return TestBit(kHasMCData); };
else return fgkDetectorName[kNdetectorPID];
}
//-----Configure PID detectors with predefined stettings------
+ void ConfigureTOF(Float_t TOFcut = 3.);
void ConfigureTPCasymmetric(Double_t pmin = 0.1, Double_t pmax = 20., Double_t sigmamin = -0.2, Double_t sigmamax = 5.);
void ConfigureTPCrejectionSimple();
- void ConfigureTPCrejection(const char *lowerCutParam = NULL, Double_t * const params = NULL, Float_t upperTPCCut=3.0, Float_t TOFCut=3.0);
- void ConfigureTPCstrategyParis();
+ void ConfigureTPCcentralityCut(Int_t centralityBin, const char *lowerCutParam = NULL, const Double_t * const params = NULL, Float_t upperTPCCut=3.0);
+ void ConfigureTPCdefaultCut(const char *lowerCutParam = NULL, const Double_t * const params = NULL, Float_t upperTPCCut=3.0);
//------------------------------------------------------------
protected:
Bool_t MakePidTpcTof(AliHFEpidObject *track);
-
private:
enum{
kHasMCData = BIT(14)
return det < kNdetectorPID ? TESTBIT(fEnabledDetectors, det): kFALSE;
}
//--------------------------------------------------
+ void ConfigureTPCcut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut);
static const Char_t *fgkDetectorName[kNdetectorPID + 1]; // PID Detector Names
AliHFEpidBase *fDetectorPID[kNdetectorPID]; // Detector PID classes
fkRecTrack(NULL),
fAnalysisType(kESDanalysis),
fAbInitioPID(-1),
- fCentrality(99.)
+ fCentrality(99),
+ fIsPbPb(kFALSE) // Default: pp
{
}
AliHFEpidObject(const AliHFEpidObject &ref):
fkRecTrack(ref.fkRecTrack),
fAnalysisType(ref.fAnalysisType),
fAbInitioPID(ref.fAbInitioPID),
- fCentrality(ref.fCentrality)
+ fCentrality(ref.fCentrality),
+ fIsPbPb(ref.fIsPbPb)
{
}
AliHFEpidObject &operator=(const AliHFEpidObject &ref);
void SetMCTrack(const AliVParticle * mcTrack);
void SetAnalysisType(AnalysisType_t type) { fAnalysisType = type; }
void SetAbInitioPID(Int_t abInitioPID) { fAbInitioPID = abInitioPID; }
- void SetCentrality(Float_t centrality) { fCentrality = centrality; }
+ void SetCentrality(Int_t centrality) { fCentrality = centrality; }
+ void SetPbPb() { fIsPbPb = kTRUE; }
+ void SetPP() { fIsPbPb = kFALSE; }
const AliVParticle *GetRecTrack() const { return fkRecTrack; }
Int_t GetAbInitioPID() const { return fAbInitioPID; }
- Float_t GetCentrality() const { return fCentrality; }
+ Int_t GetCentrality() const { return fCentrality; }
Bool_t IsAODanalysis() const { return fAnalysisType == static_cast<UChar_t>(kAODanalysis); }
Bool_t IsESDanalysis() const { return fAnalysisType == static_cast<UChar_t>(kESDanalysis); }
+ Bool_t IsPbPb() const { return fIsPbPb; }
private:
- const AliVParticle *fkRecTrack; // Reconstructed track
- UChar_t fAnalysisType; // Analysis Mode (ESD or AOD)
- Int_t fAbInitioPID; // AbInitio PID
- Float_t fCentrality; // Centrality Information
+ const AliVParticle *fkRecTrack; // Reconstructed track
+ UChar_t fAnalysisType; // Analysis Mode (ESD or AOD)
+ Int_t fAbInitioPID; // AbInitio PID
+ Int_t fCentrality; // Centrality Information
+ Bool_t fIsPbPb; // Collision type
};
class AliHFEpidBase : public TNamed{
#include "AliESDpid.h"
#include "AliHFEdetPIDqa.h"
-//#include "AliHFEpidEMCAL.h"
+#include "AliHFEpidEMCAL.h"
#include "AliHFEpidQAmanager.h"
-//#include "AliEMCALRecoUtils.h"
-//#include "AliEMCALGeometry.h"
-#include "AliVCluster.h"
-#include "AliVCaloCells.h"
-#include "AliVEvent.h"
+#include "AliHFEemcalPIDqa.h"
+//#include "AliVCluster.h"
+//#include "AliVCaloCells.h"
+//#include "AliVEvent.h"
#include "AliLog.h"
-#include "AliEMCALPIDUtils.h"
#include "AliPID.h"
-#include "AliESDEvent.h"
-#include "AliESDtrack.h"
-//#include "AliEMCALTrack.h"
-//#include "AliEMCALTracker.h"
+//#include "AliESDEvent.h"
+//#include "AliESDtrack.h"
#include "AliHFEpidEMCAL.h"
ClassImp(AliHFEpidEMCAL)
}
//___________________________________________________________________
-Int_t AliHFEpidEMCAL::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager */*pidqa*/) const
+Int_t AliHFEpidEMCAL::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const
{ //Function to a return a code indicating whether or not an electron candidate is selected
//
//
AliDebug(2, "PID object available");
// EMCal not fESDpid (s.s Feb. 11)
- /*
- AliHFEpidObject::AnalysisType_t anaType = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
- if(!((dynamic_cast<const AliVTrack *>(track->GetRecTrack()))->GetStatus() & AliESDtrack::kEMCALpid)) return 0;
+
+ //AliHFEpidObject::AnalysisType_t anaType = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
+ if(!((dynamic_cast<const AliVTrack *>(track->GetRecTrack()))->GetStatus() & AliESDtrack::kEMCALmatch)) return 0;
AliDebug(2, "Track Has EMCAL PID");
- */
- //if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kBeforePID);
+
+ //if(pidqa)
+ pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kBeforePID);
// not QA for EMCal, will be added (s.s Feb. 11)
- Double_t eop = 0.;//MomentumEnergyMatch(track->GetRecTrack()); // get eop (What is GetRecTrack ?)
+ Double_t eop = MomentumEnergyMatchV2(track->GetRecTrack()); // get eop (What is GetRecTrack ?)
AliDebug(2, Form("Energy - Momentum Matching e/p : %f", eop));
Int_t pdg = 0;
if(eop>feopMim && eop<feopMax){
pdg = 11;
- //if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kAfterPID);
+ //if(pidqa)
+ pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kAfterPID);
}
else
{
pdg = 211; // EMCal doesn't separate pi,k.p by e/p. return pion code as same as TRD
}
- printf("eID %g ; %d \n",eop,pdg);
+ AliDebug(1, Form("eID %g ; %d \n",eop,pdg));
return pdg;
}
+
+
+//___________________________________________________________________________
+Double_t AliHFEpidEMCAL::MomentumEnergyMatchV2(const AliVParticle *const track) const
+{ // Returns e/p & Rmatch
+
+ Double_t matchclsE = 9999.9;
+ double feop = -9999.9;
+
+ const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
+ AliESDEvent *evt = esdtrack->GetESDEvent();
+
+ Int_t icl = (const_cast<AliESDtrack *>(esdtrack))->GetEMCALcluster();
+
+ AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
+ if(!cluster->IsEMCAL()) {return feop;}
+
+ matchclsE = cluster->E();
+
+
+ if(matchclsE<9999.0) feop = matchclsE/esdtrack->P();
+
+ return feop;
+
+}
+
+
/*
-Double_t AliHFEpidEMCAL::MomentumEnergyMatch(const AliVParticle *track) const
+//____________________________________________________________________________________
+Double_t AliHFEpidEMCAL::MomentumEnergyMatchV1(const AliVParticle *track) const
{ // Returns e/p if an electron is matched
Float_t clsPos[3];
class AliPID;
class AliHFEpidQAmanager;
-//class AliEMCALGeoUtils;
-//class AliEMCALGeometry;
-//class AliEMCALRecoUtils;
-/*
-class AliEMCALGeoUtils;
-class AliEMCALGeometry;
-#include "AliEMCALGeoParams.h"
-class AliEMCALRecoUtils;
-class AliEMCALPIDUtils;
-*/
-
class AliHFEpidEMCAL : public AliHFEpidBase{
public:
virtual Bool_t InitializePID();
virtual Int_t IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *piqa) const;
- //Double_t MomentumEnergyMatch(const AliVParticle *track) const;
+ //Double_t MomentumEnergyMatchV1(const AliVParticle *track) const;
+ Double_t MomentumEnergyMatchV2(const AliVParticle *track) const;
protected:
void Copy(TObject &ref) const;
AliHFEpidESD();
AliHFEpidESD(const AliHFEpidESD &ref);
AliHFEpidESD &operator=(const AliHFEpidESD &ref);
+ virtual ~AliHFEpidESD();
virtual void InitializePID();
virtual Int_t IsSelected(const AliVParticle *track);
fTRDpidQA(NULL),
fOutput(NULL),
fESDpid(NULL),
- fNNref(NULL)
+ fNNref(NULL),
+ fTRDTotalChargeInSlice0(kFALSE)
{
//
// Default constructor
fTRDpidQA(NULL),
fOutput(NULL),
fESDpid(NULL),
- fNNref(NULL)
+ fNNref(NULL),
+ fTRDTotalChargeInSlice0(ref.fTRDTotalChargeInSlice0)
{
//
// Copy constructor
AliHFEpidQA &target = dynamic_cast<AliHFEpidQA &>(o);
target.fMC = fMC;
+ target.fTRDTotalChargeInSlice0 = fTRDTotalChargeInSlice0;
if(target.fESDpid) delete target.fESDpid;
target.fESDpid = new AliESDpid;
fV0pidMC->Init();
fTRDpidQA = new AliHFEtrdPIDqa;
+ if(fTRDTotalChargeInSlice0) fTRDpidQA->SetTotalChargeInSlice0();
fTRDpidQA->Init();
fOutput = new AliHFEcollection("pidQA", "PID QA output");
// fTRDpidResponse->MakePID(track, pidProbs);
}
//___________________________________________________________________
-Float_t AliHFEpidQA::TOFbeta(AliESDtrack * const track) const {
+Float_t AliHFEpidQA::TOFbeta(const AliESDtrack * const track) const {
// computes the TOF beta
Double_t l = track->GetIntegratedLength(); // cm
Double_t t = track->GetTOFsignal();
void SetMCEvent(AliMCEvent * const mc) { fMC = mc; };
void SetV0pidQA(Bool_t v0pidQA = kTRUE) { SetBit(kV0pidQA, v0pidQA); };
void SetRecalculateTRDpid(Bool_t recal = kTRUE) { SetBit(kRecalculateTRDpid, recal); };
+ void SetTRDTotalChargeInSlice0() { fTRDTotalChargeInSlice0 = kTRUE; }
void SetESDpid(AliESDpid* const pid) { fESDpid = pid; }
- Float_t TOFbeta(AliESDtrack* const track) const;
+ Float_t TOFbeta(const AliESDtrack* const track) const;
void CheckEvent();
void SetNNref(TFile *f) { fNNref = f; };
private:
TFile *fNNref; // reference file for NN pid
TMultiLayerPerceptron *fNet[11]; // reference networks
+ Bool_t fTRDTotalChargeInSlice0; // Fix for Foreward/Backward compatibility
ClassDef(AliHFEpidQA, 1) // PID QA tool
};
// add a list here
AliHFEpidBase()
, fLineCrossingsEnabled(0)
- , fUpperSigmaCut(NULL)
- , fLowerSigmaCut(NULL)
- , fElectronMeanCorrection(NULL)
+ , fkElectronMeanCorrection(NULL)
+ , fHasCutModel(kFALSE)
, fNsigmaTPC(3)
, fRejectionEnabled(0)
- , fPID(NULL)
{
//
// default constructor
//
+
+ memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
+ memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
+
memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
memset(fPAsigCut, 0, sizeof(Float_t) * 2);
// add a list here
AliHFEpidBase(name)
, fLineCrossingsEnabled(0)
- , fUpperSigmaCut(NULL)
- , fLowerSigmaCut(NULL)
- , fElectronMeanCorrection(NULL)
+ , fkElectronMeanCorrection(NULL)
+ , fHasCutModel(kFALSE)
, fNsigmaTPC(3)
, fRejectionEnabled(0)
- , fPID(NULL)
{
//
// default constructor
//
//
+ memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
+ memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
+
memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
memset(fPAsigCut, 0, sizeof(Float_t) * 2);
memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
- fPID = new AliPID;
}
//___________________________________________________________________
AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
AliHFEpidBase("")
, fLineCrossingsEnabled(0)
- , fUpperSigmaCut(NULL)
- , fLowerSigmaCut(NULL)
- , fElectronMeanCorrection(NULL)
+ , fkElectronMeanCorrection(NULL)
+ , fHasCutModel(ref.fHasCutModel)
, fNsigmaTPC(2)
, fRejectionEnabled(0)
- , fPID(NULL)
{
//
// Copy constructor
AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
target.fLineCrossingsEnabled = fLineCrossingsEnabled;
- target.fUpperSigmaCut = fUpperSigmaCut;
- target.fLowerSigmaCut = fLowerSigmaCut;
- target.fElectronMeanCorrection = fElectronMeanCorrection;
+ target.fkElectronMeanCorrection = fkElectronMeanCorrection;
+ target.fHasCutModel = fHasCutModel;
target.fNsigmaTPC = fNsigmaTPC;
target.fRejectionEnabled = fRejectionEnabled;
- target.fPID = new AliPID(*fPID);
+
+ memcpy(target.fkUpperSigmaCut, fkUpperSigmaCut, sizeof(const TF1 *) * 12);
+ memcpy(target.fkLowerSigmaCut, fkLowerSigmaCut, sizeof(const TF1 *) * 12);
+
memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
//
// Destructor
//
- if(fPID) delete fPID;
}
//___________________________________________________________________
// Check if we have an asymmetric sigma model set
Int_t pdg = 0;
- if(fUpperSigmaCut || fLowerSigmaCut){
- pdg = CutSigmaModel(track->GetRecTrack(), anatype) ? 11 : 0;
+ if(fHasCutModel){
+ pdg = CutSigmaModel(track) ? 11 : 0;
} else {
// Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
Float_t p = 0.;
}
//___________________________________________________________________
-Bool_t AliHFEpidTPC::CutSigmaModel(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
+Bool_t AliHFEpidTPC::CutSigmaModel(const AliHFEpidObject * const track) const {
//
// N SigmaCut using parametrization of the cuts
//
Bool_t isSelected = kTRUE;
- Float_t nsigma = NumberOfSigmas(track, AliPID::kElectron, anaType);
- Double_t p = GetP(track, anaType);
- if(fUpperSigmaCut && nsigma > fUpperSigmaCut->Eval(p)) isSelected = kFALSE;
- if(fLowerSigmaCut && nsigma < fLowerSigmaCut->Eval(p)) isSelected = kFALSE;
+ AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
+ Float_t nsigma = NumberOfSigmas(track->GetRecTrack(), AliPID::kElectron, anatype);
+ Double_t p = GetP(track->GetRecTrack(), anatype);
+ Int_t centrality = track->IsPbPb() ? track->GetCentrality() + 1 : 0;
+ AliDebug(2, Form("Centrality: %d\n", centrality));
+ const TF1 *cutfunction;
+ if((cutfunction = fkUpperSigmaCut[centrality]) && nsigma > cutfunction->Eval(p)) isSelected = kFALSE;
+ if((cutfunction = fkLowerSigmaCut[centrality]) && nsigma < cutfunction->Eval(p)) isSelected = kFALSE;
return isSelected;
}
if(aodtrack && fAODpid) nSigmas = fAODpid->NumberOfSigmasTPC(aodtrack, species);
}
// Correct for the mean o
- if(fElectronMeanCorrection)
- nSigmas -= fElectronMeanCorrection->Eval(GetP(track, anaType));
+ if(fkElectronMeanCorrection)
+ nSigmas -= fkElectronMeanCorrection->Eval(GetP(track, anaType));
return nSigmas;
}
void AddTPCdEdxLineCrossing(Int_t species, Double_t sigma);
Bool_t HasAsymmetricSigmaCut() const { return TestBit(kAsymmetricSigmaCut);}
Bool_t HasParticleRejection() const { return TestBit(kRejection); }
- void SetElectronMeanCorrection(TF1 *electronLineCorrection) { fElectronMeanCorrection = electronLineCorrection; }
+ void SetElectronMeanCorrection(const TF1 * const electronLineCorrection) { fkElectronMeanCorrection = electronLineCorrection; }
void SetTPCnSigma(Short_t nSigma) { fNsigmaTPC = nSigma; };
inline void SetAsymmetricTPCsigmaCut(Float_t pmin, Float_t pmax, Float_t sigmaMin, Float_t sigmaMax);
inline void SetRejectParticle(Int_t species, Float_t pmin, Float_t sigmaMin, Float_t pmax, Float_t sigmaMax);
- void SetUpperSigmaCut(TF1 * const model) { fUpperSigmaCut = model; }
- void SetLowerSigmaCut(TF1 * const model) { fLowerSigmaCut = model; }
+ void SetUpperSigmaCutDefault(const TF1 * const model) { fkUpperSigmaCut[0] = model; fHasCutModel = kTRUE; }
+ void SetUpperSigmaCutCentrality(const TF1 * const model, Int_t centralityBin) { if(centralityBin < 11) fkUpperSigmaCut[centralityBin+1] = model; fHasCutModel = kTRUE; }
+ void SetLowerSigmaCutDefault(const TF1 * const model) { fkLowerSigmaCut[0] = model; fHasCutModel = kTRUE; }
+ void SetLowerSigmaCutCentrality(const TF1 * const model, Int_t centralityBin) { if(centralityBin < 11) fkLowerSigmaCut[centralityBin+1] = model; fHasCutModel = kTRUE; }
Double_t NumberOfSigmas(const AliVParticle *track, AliPID::EParticleType species, AliHFEpidObject::AnalysisType_t anatype) const;
Double_t GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const;
void Copy(TObject &o) const;
Int_t Reject(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const;
- Bool_t CutSigmaModel(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const;
+ Bool_t CutSigmaModel(const AliHFEpidObject *anaType) const;
private:
enum{
};
Double_t fLineCrossingSigma[AliPID::kSPECIES]; // with of the exclusion point
UChar_t fLineCrossingsEnabled; // Bitmap showing which line crossing is set
- TF1 *fUpperSigmaCut; // Upper Sigma Cut
- TF1 *fLowerSigmaCut; // Lower Sigma Cut
- TF1 *fElectronMeanCorrection; // Correct the mean of the electron line position as function of the momentum
+ const TF1 *fkUpperSigmaCut[12]; // Upper Sigma Cut
+ const TF1 *fkLowerSigmaCut[12]; // Lower Sigma Cut
+ const TF1 *fkElectronMeanCorrection; // Correct the mean of the electron line position as function of the momentum
+ Bool_t fHasCutModel; // Has cut model functions
Float_t fPAsigCut[2]; // Momentum region where to perform asymmetric sigma cut
Float_t fNAsigmaTPC[2]; // Asymmetric TPC Sigma band
Short_t fNsigmaTPC; // TPC sigma band
Float_t fRejection[4*AliPID::kSPECIES]; // All informations for Particle Rejection, order pmin, sigmin, pmax, sigmax
UChar_t fRejectionEnabled; // Bitmap for enabled particle rejection
- AliPID *fPID; //! PID Object
ClassDef(AliHFEpidTPC, 1) // TPC Electron ID class
};
, fMinP(1.)
, fElectronEfficiency(0.91)
, fPIDMethod(kNN)
+ , fTotalChargeInSlice0(kFALSE)
{
//
// default constructor
, fMinP(1.)
, fElectronEfficiency(0.91)
, fPIDMethod(kNN)
+ , fTotalChargeInSlice0(kFALSE)
{
//
// default constructor
, fMinP(ref.fMinP)
, fElectronEfficiency(ref.fElectronEfficiency)
, fPIDMethod(ref.fPIDMethod)
+ , fTotalChargeInSlice0(ref.fTotalChargeInSlice0)
{
//
// Copy constructor
target.SetUseDefaultParameters(defaultParameters);
target.fMinP = fMinP;
target.fPIDMethod = fPIDMethod;
+ target.fTotalChargeInSlice0 = fTotalChargeInSlice0;
target.fElectronEfficiency = fElectronEfficiency;
memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
AliHFEpidBase::Copy(ref);
Double_t charge = 0.;
if(anaType == AliHFEpidObject::kESDanalysis){
const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
- if(esdtrack)
- for(Int_t islice = 0; islice < esdtrack->GetNumberOfTRDslices(); islice++) charge += esdtrack->GetTRDslice(static_cast<UInt_t>(layer), islice);
+ if(esdtrack){
+ // Distinction between old and new reconstruction: in the new reconstruction, the total charge is stored in slice 0, slices 1 to 8 are used for the slices for
+ // the neural network.
+ if(fTotalChargeInSlice0)
+ charge = esdtrack->GetTRDslice(static_cast<UInt_t>(layer), 0);
+ else
+ for(Int_t islice = 0; islice < esdtrack->GetNumberOfTRDslices(); islice++) charge += esdtrack->GetTRDslice(static_cast<UInt_t>(layer), islice);
+ }
} else {
const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
AliAODPid *aoddetpid = aodtrack ? aodtrack->GetDetPid() : NULL;
- if(aoddetpid)
- for(Int_t islice = 0; islice < aoddetpid->GetTRDnSlices(); islice++) charge += aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices() + islice];
+ if(aoddetpid){
+ if(fTotalChargeInSlice0)
+ charge = aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices()];
+ else
+ for(Int_t islice = 0; islice < aoddetpid->GetTRDnSlices(); islice++) charge += aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices() + islice];
+ }
}
return charge;
}
Int_t indices[48];
Int_t icnt = 0;
for(Int_t idet = 0; idet < 6; idet++)
- for(Int_t islice = 0; islice < kSlicePerLayer; islice++){
+ for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0 ; islice < kSlicePerLayer; islice++){
AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
trdSlices[icnt++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
Int_t indices[kLayers*kSlicesLow];
Int_t cntLowTime=0, cntRemaining = 0;
for(Int_t idet = 0; idet < 6; idet++)
- for(Int_t islice = 0; islice < kSlicesLow+kSlicesHigh; islice++){
+ for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0; islice < kSlicesLow+kSlicesHigh; islice++){
if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
if(islice < kSlicesLow){
AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
Bool_t IsCalculateTRDSignals() const { return TestBit(kTRDsignals); }
Bool_t IsRenormalizeElPi() const { return TestBit(kTRDrenormalize); }
void SetPIDMethod(PIDMethodTRD_t method) { fPIDMethod = method; };
+ void SetTotalChargeInSlice0() { fTotalChargeInSlice0 = kTRUE; }
void SetRenormalizeElPi(Bool_t doRenorm = kTRUE) { if(doRenorm) SetBit(kTRDrenormalize, kTRUE); else SetBit(kTRDrenormalize, kFALSE);}
void SetElectronEfficiency(Double_t electronEfficiency) { fElectronEfficiency = electronEfficiency; }
void SetThresholdParameters(Double_t electronEff, Double_t *params);
Double_t fElectronEfficiency; // Cut on electron efficiency
PIDMethodTRD_t fPIDMethod; // PID Method: 2D Likelihood or Neural Network
Double_t fThreshParams[kThreshParams]; // Threshold parametrisation
+ Bool_t fTotalChargeInSlice0; // Flag for foreward/backward compatibility for the TRD total charge
ClassDef(AliHFEpidTRD, 1) // TRD electron ID class
};
#endif
}
//_______________________________________________________________________________________________
-Int_t AliHFEsecVtx::GetPDG(AliVTrack *track)
+Int_t AliHFEsecVtx::GetPDG(const AliVTrack *track)
{
//
// get KF particle input pdg for mass hypothesis
fSecVtxList->AddAt(fSecvtxQA,1);
for(Int_t ivar = 0; ivar < nDimPair; ivar++)
- delete binEdgesPair[ivar];
+ delete[] binEdgesPair[ivar];
for(Int_t ivar = 0; ivar < nDimSecvtx; ivar++)
- delete binEdgesSecvtx[ivar];
+ delete[] binEdgesSecvtx[ivar];
}
//____________________________________________________________
fSecVtxList->AddAt(new TH1F(hname+"ebgElec", "pT of e", iBin[0],binEdges[0]), step+4);
fSecVtxList->AddAt(new TH1F(hname+"hcontaminElec", "pT of e", iBin[0],binEdges[0]), step+5);
fSecVtxList->AddAt(new TH1F(hname+"elseElec", "pT of e", iBin[0],binEdges[0]), step+6);
+ delete[] binEdges[0];
}
//____________________________________________________________
Int_t GetPairOriginAOD(AliAODTrack* track1, AliAODTrack* track2); // return pair origin as a pdg code
Int_t GetPairCode(const AliVTrack* const track1, const AliVTrack* const track2); // return corresponding pair code to pdg code
Int_t GetElectronSource(Int_t mclabel); // return origin of the electron
- Int_t GetPDG(AliVTrack *track); // return pdg
+ Int_t GetPDG(const AliVTrack *track); // return pdg
void GetESDPID(const AliESDtrack *track, Int_t &recpid, Double_t &recprob); //return esd pid likelihood
void GetPrimaryCondition();
void RecalcPrimvtx(Int_t nkftrk, const Int_t * const, const AliKFParticle * const); //recalculate primary vertex
#include <TParticle.h>
#include <TString.h>
+#include "AliAODTrack.h"
#include "AliAODMCParticle.h"
+#include "AliESDtrack.h"
#include "AliLog.h"
#include "AliMCEvent.h"
#include "AliMCParticle.h"
return 0;
}
- TString sourcetype = track->IsA()->GetName();
+ TClass *tracktype;
const AliVParticle *mctrack = NULL;
TParticle *mcpart = NULL;
- if(!sourcetype.CompareTo("AliESDtrack") || !sourcetype.CompareTo("AliAODTrack")){
+ if((tracktype = track->IsA()) == AliESDtrack::Class() || tracktype == AliAODTrack::Class()){
mctrack = fMC->GetTrack(TMath::Abs(track->GetLabel()));
} else mctrack = track;
if(!mctrack) return 0;
- TString mctype = mctrack->IsA()->GetName();
Int_t eSource = 0;
- if(!mctype.CompareTo("AliMCParticle")){
+ if(mctrack->IsA() == AliMCParticle::Class()){
const AliMCParticle *esdmc = dynamic_cast<const AliMCParticle *>(mctrack);
if(esdmc){
mcpart = esdmc->Particle();
fInclusiveSpectrum(kTRUE),
fDumpToFile(kFALSE),
fUnSetCorrelatedErrors(kTRUE),
+ fSetSmoothing(kFALSE),
fIPanaHadronBgSubtract(kFALSE),
fIPanaCharmBgSubtract(kFALSE),
fIPanaConversionBgSubtract(kFALSE),
else if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
else if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD - 1))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
+
+ if(!mccorrelation) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
}
else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
//else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
if(fBeamType==1) {
- TCanvas * ccorrected_allspectra = new TCanvas("correctedallspectra","correctedallspectra",1000,700);
- ccorrected_allspectra->Divide(2,1);
+ TCanvas * ccorrectedallspectra = new TCanvas("correctedallspectra","correctedallspectra",1000,700);
+ ccorrectedallspectra->Divide(2,1);
TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
titleenamednormalized += "[";
Int_t nbEvents = 0;
for(Int_t k = fLowBoundaryCentralityBinAtTheEnd[binc]; k < fHighBoundaryCentralityBinAtTheEnd[binc]; k++) {
- //printf("Number of events %d in the bin %d added!!!\n",fNEvents[k],k);
+ printf("Number of events %d in the bin %d added!!!\n",fNEvents[k],k);
nbEvents += fNEvents[k];
}
- //Double_t lowedgega = cenaxisa->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
- //Double_t upedgega = cenaxisa->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
- //printf("Bin Low edge %f, up edge %f for a\n",lowedgega,upedgega);
- //Double_t lowedgegb = cenaxisb->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
- //Double_t upedgegb = cenaxisb->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
- //printf("Bin Low edge %f, up edge %f for b\n",lowedgegb,upedgegb);
+ Double_t lowedgega = cenaxisa->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
+ Double_t upedgega = cenaxisa->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
+ printf("Bin Low edge %f, up edge %f for a\n",lowedgega,upedgega);
+ Double_t lowedgegb = cenaxisb->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
+ Double_t upedgegb = cenaxisb->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
+ printf("Bin Low edge %f, up edge %f for b\n",lowedgegb,upedgegb);
cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
TCanvas * ccorrectedcheck = new TCanvas((const char*) titleec,(const char*) titleec,1000,700);
legcorrectedud->AddEntry(aftersun,"Corrected","p");
legcorrectedud->AddEntry(aftersdn,"Alltogether","p");
legcorrectedud->Draw("same");
- ccorrected_allspectra->cd(1);
+ ccorrectedallspectra->cd(1);
gPad->SetLogy();
TH1D *aftersunn = (TH1D *) aftersun->Clone();
aftersunn->SetMarkerStyle(stylee[binc]);
if(binc==0) aftersunn->Draw("AP");
else aftersunn->Draw("P");
legtotal->AddEntry(aftersunn,(const char*) titlee,"p");
- ccorrected_allspectra->cd(2);
+ ccorrectedallspectra->cd(2);
gPad->SetLogy();
TH1D *aftersdnn = (TH1D *) aftersdn->Clone();
aftersdnn->SetMarkerStyle(stylee[binc]);
ratiocorrectedbinc->Draw();
}
- ccorrected_allspectra->cd(1);
+ ccorrectedallspectra->cd(1);
legtotal->Draw("same");
- ccorrected_allspectra->cd(2);
+ ccorrectedallspectra->cd(2);
legtotalg->Draw("same");
cenaxisa->SetRange(0,nbbin);
//unfolding.SetUseCorrelatedErrors();
if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
- unfolding.UseSmoothing();
+ if(fSetSmoothing) unfolding.UseSmoothing();
unfolding.Unfold();
// Results
dN = input->GetBinError(ibin);
// New point
- nCorr = 0.5 * 1/1.6 * 1./(Double_t)(fNEvents[i]) * 1/(2. * TMath::Pi() * p) * n;
+ nCorr = 0.5 * 1./1.6 * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n;
errdN = 1./(2. * TMath::Pi() * p);
errdp = 1./(2. * TMath::Pi() * p*p) * n;
- dNcorr = 0.5 * 1/1.6 * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
+ dNcorr = 0.5 * 1./1.6 * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
spectrumNormalized->SetPoint(point, p, nCorr);
spectrumNormalized->SetPointError(point, dp, dNcorr);
dN = input->GetBinError(ibin);
// New point
- nCorr = 0.5 * 1/1.6 * 1./(Double_t)(normalization) * 1/(2. * TMath::Pi() * p) * n;
+ nCorr = 0.5 * 1./1.6 * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n;
errdN = 1./(2. * TMath::Pi() * p);
errdp = 1./(2. * TMath::Pi() * p*p) * n;
- dNcorr = 0.5 * 1/1.6 * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
+ dNcorr = 0.5 * 1./1.6 * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
spectrumNormalized->SetPoint(point, p, nCorr);
spectrumNormalized->SetPointError(point, dp, dNcorr);
void SetNbDimensions(Int_t nbDimensions) { fNbDimensions = nbDimensions; };
void SetChargeChoosen(Int_t chargechoosen) {fChargeChoosen = chargechoosen; };
void SetUnSetCorrelatedErrors(Bool_t unsetcorrelatederrors) {fUnSetCorrelatedErrors = unsetcorrelatederrors;};
+ void SetSmoothing(Bool_t setSmoothing) {fSetSmoothing = setSmoothing;};
void SetNCentralityBinAtTheEnd(Int_t nCentralityBinAtTheEnd) {fNCentralityBinAtTheEnd = nCentralityBinAtTheEnd; };
void SetLowHighBoundaryCentralityBinAtTheEnd(Int_t low, Int_t high, Int_t i) { fLowBoundaryCentralityBinAtTheEnd[i] = low; fHighBoundaryCentralityBinAtTheEnd[i] = high;};
Bool_t fDumpToFile; // Write Result in a file
Bool_t fUnSetCorrelatedErrors; // Unset correlated errors
+ Bool_t fSetSmoothing; // Set smoothing
Bool_t fIPanaHadronBgSubtract; // Hadron background subtraction
Bool_t fIPanaCharmBgSubtract; // Charm background subtraction
, fCuts(NULL)
, fCFM(NULL)
, fQAhistos(NULL)
- , fCentralityF(0.)
+ , fCentralityF(0)
, fClean(kFALSE)
, fMagneticField(0.0)
, fVariablesTRD(kFALSE)
+ , fIsPbPb(kFALSE)
{
//
// Dummy constructor
, fCuts(NULL)
, fCFM(NULL)
, fQAhistos(NULL)
- , fCentralityF(0.)
+ , fCentralityF(0)
, fClean(kFALSE)
, fMagneticField(0.0)
, fVariablesTRD(kFALSE)
+ , fIsPbPb(kFALSE)
{
//
// Default constructor
, fClean(ref.fClean)
, fMagneticField(ref.fMagneticField)
, fVariablesTRD(ref.fVariablesTRD)
+ , fIsPbPb(ref.fIsPbPb)
{
//
// Copy constructor
fClean = ref.fClean;
fMagneticField = ref.fMagneticField;
fVariablesTRD = ref.fVariablesTRD;
+ fIsPbPb = ref.fIsPbPb;
if(ref.fContainer) InitContainer();
hfetrack.SetRecTrack(track);
hfetrack.SetAbInitioPID(abinitioPID);
hfetrack.SetCentrality(fCentralityF);
+ if(fIsPbPb) hfetrack.SetPbPb();
+ else hfetrack.SetPP();
fPID->SetVarManager(fVarManager);
fPID->IsSelected(&hfetrack, fContainer, "taggedTrackContainer", fPIDqa);
}
TList * GetCutQA() const;
AliHFEcollection * GetQAcollection() const { return fQAhistos; }
Bool_t GetClean() const { return fClean; };
+ Bool_t IsPbPb() const { return fIsPbPb; }
+ Bool_t IsPP() const { return !fIsPbPb; }
Double_t GetMagneticField() const { return fMagneticField; };
AliHFEvarManager *GetVarManager() const { return fVarManager; }
void SetCuts(AliHFEcuts *cuts);
void SetPID(AliHFEpid *pid);
- void SetCentrality(Float_t centrality) { fCentralityF = centrality; };
+ void SetCentrality(Int_t centrality) { fCentralityF = centrality; };
+ void SetPbPb(){ fIsPbPb = kTRUE; }
+ void SetPP() { fIsPbPb = kFALSE; }
void SetClean(Bool_t clean) { fClean = clean; };
void SetMagneticField(Double_t magneticField) { fMagneticField = magneticField; };
void SetVariablesTRD(Bool_t variablesTRD) { fVariablesTRD = variablesTRD; };
AliHFEcuts *fCuts; // Single track cuts
AliCFManager *fCFM; // CF Manager used for the track filtering
AliHFEcollection *fQAhistos; // QA histos
- Float_t fCentralityF; // Centrality
+ Int_t fCentralityF; // Centrality
Bool_t fClean; // Clean
Double_t fMagneticField; // Magnetic field
- Bool_t fVariablesTRD; // Use phi angle at the first plane of the TRD
+ Bool_t fVariablesTRD; // Use phi angle at the first plane of the TRD
+ Bool_t fIsPbPb; // Analysis Type: pp or PbPb
ClassDef(AliHFEtaggedTrackAnalysis, 0)
};
AliHFEtpcPIDqa::AliHFEtpcPIDqa():
AliHFEdetPIDqa()
, fHistos(NULL)
+ , fBrowseCentrality(-1)
{
//
// Dummy constructor
AliHFEtpcPIDqa::AliHFEtpcPIDqa(const char* name):
AliHFEdetPIDqa(name, "QA for TPC")
, fHistos(NULL)
+ , fBrowseCentrality(-1)
{
//
// Default constructor
AliHFEtpcPIDqa::AliHFEtpcPIDqa(const AliHFEtpcPIDqa &o):
AliHFEdetPIDqa(o)
, fHistos(NULL)
+ , fBrowseCentrality(o.fBrowseCentrality)
{
//
// Copy constructor
target.fHistos = NULL;
}
if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
+ target.fBrowseCentrality = fBrowseCentrality;
}
//_________________________________________________________
TH2 *hptr = NULL;
for(Int_t ispec = 0; ispec < AliPID::kSPECIES+1; ispec++){
for(Int_t istep = 0; istep < 2; istep++){
- hptr = MakeSpectrumdEdx(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
- hptr->SetName(Form("hTPCdEdx%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
+ hptr = MakeSpectrumdEdx(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], fBrowseCentrality);
+ hptr->SetName(Form("hTPCdEdx%s%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After" , fBrowseCentrality == -1 ? "MinBias" : Form("Cent%d", fBrowseCentrality)));
listdEdx->Add(hptr);
- hptr = MakeSpectrumNSigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
- hptr->SetName(Form("hTPCnsigma%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
+ hptr = MakeSpectrumNSigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], fBrowseCentrality);
+ hptr->SetName(Form("hTPCnsigma%s%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After", fBrowseCentrality == -1 ? "MinBias" : Form("Cent%d", fBrowseCentrality)));
listNsigma->Add(hptr);
}
}
}
//_________________________________________________________
-TH2 *AliHFEtpcPIDqa::MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t istep, Int_t species){
+TH2 *AliHFEtpcPIDqa::MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
//
// Plot the Spectrum
//
hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
if(species >= 0 && species < AliPID::kSPECIES)
hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
- TH2 *hTmp = hSignal->Projection(2,1);
TString hname = Form("hTPCsignal%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after"),
htitle = Form("TPC dE/dx Spectrum %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
if(species > -1){
hname += AliPID::ParticleName(species);
htitle += Form(" for %ss", AliPID::ParticleName(species));
}
+ TString centname, centtitle;
+ Bool_t hasCentralityInfo = kTRUE;
+ if(centralityClass > -1){
+ if(hSignal->GetNdimensions() < 5){
+ AliError("Centrality Information not available");
+ centname = centtitle = "MinBias";
+ hasCentralityInfo= kFALSE;
+ } else {
+ // Project centrality classes
+ // -1 is Min. Bias
+ hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
+ centname = Form("Cent%d", centralityClass);
+ centtitle = Form(" Centrality %d", centralityClass);
+ }
+ } else {
+ centname = centtitle = "MinBias";
+ hasCentralityInfo= kFALSE;
+ }
+ hname += centtitle;
+ htitle += centtitle;
+
+ TH2 *hTmp = hSignal->Projection(2,1);
hTmp->SetName(hname.Data());
hTmp->SetTitle(htitle.Data());
hTmp->SetStats(kFALSE);
hTmp->GetYaxis()->SetTitle("TPC signal [a.u.]");
hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
+ if(hasCentralityInfo) hSignal->GetAxis(4)->SetRange(0, hSignal->GetAxis(4)->GetNbins());
return hTmp;
}
//_________________________________________________________
-TH2 *AliHFEtpcPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species){
+TH2 *AliHFEtpcPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
//
// Plot the Spectrum
//
hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
if(species >= 0 && species < AliPID::kSPECIES)
hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
- TH2 *hTmp = hSignal->Projection(2,1);
TString hname = Form("hTPCsigma%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after"),
htitle = Form("TPC dE/dx Spectrum[#sigma] %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
if(species > -1){
hname += AliPID::ParticleName(species);
htitle += Form(" for %ss", AliPID::ParticleName(species));
}
+ TString centname, centtitle;
+ Bool_t hasCentralityInfo = kTRUE;
+ if(centralityClass > -1){
+ if(hSignal->GetNdimensions() < 5){
+ AliError("Centrality Information not available");
+ centname = centtitle = "MinBias";
+ hasCentralityInfo= kFALSE;
+ } else {
+ // Project centrality classes
+ // -1 is Min. Bias
+ hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
+ centname = Form("Cent%d", centralityClass);
+ centtitle = Form(" Centrality %d", centralityClass);
+ }
+ } else {
+ centname = centtitle = "MinBias";
+ hasCentralityInfo= kFALSE;
+ }
+ hname += centtitle;
+ htitle += centtitle;
+
+ TH2 *hTmp = hSignal->Projection(2,1);
hTmp->SetName(hname.Data());
hTmp->SetTitle(htitle.Data());
hTmp->SetStats(kFALSE);
hTmp->GetYaxis()->SetTitle("TPC dE/dx - <dE/dx>|_{el} [#sigma]");
hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
+ if(hasCentralityInfo) hSignal->GetAxis(4)->SetRange(0, hSignal->GetAxis(4)->GetNbins());
return hTmp;
}
virtual void Initialize();
virtual void ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step);
+ void SetBrowseCentrality(Int_t browseCentrality) { browseCentrality < 11 && browseCentrality >= -1 ? fBrowseCentrality = browseCentrality : -1;} // *MENU*
+
AliHFEcollection *GetHistograms() const { return fHistos; }
- TH2 *MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t step, Int_t species = -1);
- TH2 *MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t step, Int_t species = -1);
+ TH2 *MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t step, Int_t species = -1, Int_t centralityClass = -1);
+ TH2 *MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t step, Int_t species = -1, Int_t centralityClass = -1);
protected:
Double_t GetTPCsignal(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype);
private:
AliHFEcollection *fHistos; // Container for Histograms
+ Int_t fBrowseCentrality; // Centrality Class for Browser
ClassDef(AliHFEtpcPIDqa, 1);
};
}
//__________________________________________________________________
-void AliHFEtrackFilter::SetRecEvent(AliVEvent *rec){
+void AliHFEtrackFilter::SetRecEvent(const AliVEvent *rec){
//
// Publish MC event to the single cut steps
//
// Base Functionality
void AddCutStep(AliHFEcutStep *cuts);
void GenerateCutSteps();
- void SetRecEvent(AliVEvent *rec);
+ void SetRecEvent(const AliVEvent *rec);
AliHFEcutStep *GetCutStep(Int_t istep);
AliHFEcutStep *GetCutStep(const Char_t *name);
void FilterTracks(const AliESDEvent *const esd);
fProtonEfficiencies(NULL),
fKaonEfficiencies(NULL),
fThresholds(NULL),
- fShowMessage(kFALSE)
+ fShowMessage(kFALSE),
+ fTotalChargeInSlice0(kFALSE)
{
//
// Default Constructor
fProtonEfficiencies(NULL),
fKaonEfficiencies(NULL),
fThresholds(NULL),
- fShowMessage(kFALSE)
+ fShowMessage(kFALSE),
+ fTotalChargeInSlice0(kFALSE)
{
//
// Main Constructor
fProtonEfficiencies(NULL),
fKaonEfficiencies(NULL),
fThresholds(NULL),
- fShowMessage(kFALSE)
+ fShowMessage(kFALSE),
+ fTotalChargeInSlice0(ref.fTotalChargeInSlice0)
{
//
// Copy constructor
AliHFEtrdPIDqa &target = dynamic_cast<AliHFEtrdPIDqa &>(ref);
target.fTRDpid = fTRDpid;
- target.fHistos = dynamic_cast<AliHFEcollection *>(fHistos->Clone());
+ target.fHistos = dynamic_cast<AliHFEcollection *>(fHistos->Clone());
+ target.fTotalChargeInSlice0 = fTotalChargeInSlice0;
}
//__________________________________________________________________
fHistos->CreateTHnSparse("fTRDtruncMean","TRD TruncatedMean studies", kQuantitiesTruncMean, nbins, binMin, binMax);
fHistos->BinLogAxis("fTRDtruncMean", kP);
- fHistos->CreateTH2F("fTRDslicesPions","TRD dEdx per slice for Pions", 8, 0, 8, 500, 0, 2000);
- fHistos->CreateTH2F("fTRDslicesElectrons","TRD dEdx per slice for Electrons", 8, 0, 8, 500, 0, 2000);
+ fHistos->CreateTH2F("fTRDslicesPions","TRD dEdx per slice for Pions", 8, 0, 8, 2000, 0, 8000);
+ fHistos->CreateTH2F("fTRDslicesElectrons","TRD dEdx per slice for Electrons", 8, 0, 8, 2000, 0, 8000);
}
quantitiesdEdx[kP] = track->GetTRDmomentum(iplane);
dEdxSum = 0.;
for(Int_t islice = 0; islice < nSlices; islice++){
- qSlice = track->GetTRDslice(iplane, islice);
+ if(fTotalChargeInSlice0 && islice >= 7) break;
+ qSlice = track->GetTRDslice(iplane, fTotalChargeInSlice0 ? islice + 1 : islice); // hack by mfasel: For data with the new reconstruction, slice 0 is used to store the total charge, the total number of slices is 7 instead of 8
if(qSlice > 1e-1){
// cut out 0 slices
nSlicesNonZero++;
}
}
quantitiesdEdx[kNonZeroSlices] = nSlicesNonZero;
- quantitiesdEdx[kdEdx] = dEdxSum;
+ quantitiesdEdx[kdEdx] = fTotalChargeInSlice0 ? track->GetTRDslice(iplane, 0) : dEdxSum; // hack by mfasel: In the new reconstruction, the total charge is stored in the first slice, in the old reconstruction it has to be calculated from the slice charges.
if(dEdxSum) ntrackletsNonZero++;
// Fill dEdx histogram
if(dEdxSum > 1e-1) fHistos->Fill("fQAdEdx", quantitiesdEdx); // Cut out 0 entries
if(!fPionEfficiencies){
fPionEfficiencies = new TList;
fPionEfficiencies->SetName("pionEfficiencies");
+ fPionEfficiencies->SetOwner();
}
if(!fProtonEfficiencies){
fProtonEfficiencies = new TList;
fProtonEfficiencies->SetName("protonEfficiencies");
+ fProtonEfficiencies->SetOwner();
}
if(!fThresholds){
fThresholds = new TList;
fThresholds->SetName("thresholds");
+ fThresholds->SetOwner();
}
for(Int_t itr = 4; itr <= 6; itr++){
}
//__________________________________________________________________
-void AliHFEtrdPIDqa::SaveThresholdParameters(const Char_t *name){
+void AliHFEtrdPIDqa::SaveThresholdParameters(const Char_t *name, Double_t lowerLimit, Double_t upperLimit){
//
// Fit the threshold histograms with the given parametrisation
// and store the TF1 in the file
threshhist = dynamic_cast<TGraph *>(lHistos->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
if(!threshhist) continue;
- threshparam = MakeThresholds(threshhist);
+ threshparam = MakeThresholds(threshhist, lowerLimit, upperLimit);
threshparam->SetName(Form("thresh_%d_%d", itracklet, static_cast<Int_t>(fgkElectronEff[ieff] * 100)));
lFormulas->Add(threshparam);
}
hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kNTracklets)->GetNbins());
// Prepare List for output
- TList *listPions = new TList; listPions->SetName(Form("%dTracklets", nTracklets));
- TList *listProtons = new TList; listProtons->SetName(Form("%dTracklets", nTracklets));
- TList *listThresholds = new TList; listThresholds->SetName(Form("%dTracklets", nTracklets));
+ TList *listPions = new TList; listPions->SetName(Form("%dTracklets", nTracklets)); listPions->SetOwner();
+ TList *listProtons = new TList; listProtons->SetName(Form("%dTracklets", nTracklets)); listProtons->SetOwner();
+ TList *listThresholds = new TList; listThresholds->SetName(Form("%dTracklets", nTracklets)); listThresholds->SetOwner();
fPionEfficiencies->Add(listPions);
fProtonEfficiencies->Add(listProtons);
fThresholds->Add(listThresholds);
// Optionally do Fit
if(doFit){
- threshfit = MakeThresholds(tr);
+ threshfit = MakeThresholds(tr, pmin, pmax);
threshfit->SetLineColor(kBlack);
threshfit->Draw("same");
}
}
//__________________________________________________________________
-TF1 *AliHFEtrdPIDqa::MakeThresholds(TGraph *threshist){
+TF1 *AliHFEtrdPIDqa::MakeThresholds(TGraph *threshist, Double_t lowerLimit, Double_t upperLimit){
//
// Create TF1 containing the threshold parametrisation
//
TF1 *threshparam = new TF1("thresh", "1-[0]-[1]*x-[2]*TMath::Exp(-[3]*x)", 0.1, 10);
- threshist->Fit(threshparam, "NE", "", 0.1, 3.5);
+ threshist->Fit(threshparam, "NE", "", lowerLimit, upperLimit);
return threshparam;
}
void Init();
void FinishAnalysis();
void ShowMessages() { fShowMessage = kTRUE; }
+ void SetTotalChargeInSlice0() { fTotalChargeInSlice0 = kTRUE; }
void StoreResults(const Char_t *filename = "HFEtrdPIDqa.root");
- void SaveThresholdParameters(const Char_t * filename = "TRD.Thresholds.root");
+ void SaveThresholdParameters(const Char_t * filename = "TRD.Thresholds.root", Double_t lowerLimit = 0.5, Double_t upperLimit = 3.5);
void DrawTracklet(Int_t tracklet, Double_t pmin = 0., Double_t pmax = 0., Bool_t doFit = kFALSE);
void ClearLists();
void AnalyseNTracklets(Int_t nTracklets);
Int_t GetThresholdBin(const TH1 * const input, Double_t efficiency);
Bool_t CalculateEfficiency(const TH1 * const input, Int_t threshbin, Double_t *params);
- TF1 *MakeThresholds(TGraph *input);
+ TF1 *MakeThresholds(TGraph *input, Double_t lowerLimit, Double_t upperLimit);
void CreateLikelihoodHistogram();
void CreateQAHistogram();
TList *fThresholds; //! List for Threshold Graphs
Bool_t fShowMessage; // Display debug messages
+ Bool_t fTotalChargeInSlice0; // Flag for Foreward/Backward compatibility in TRD total charge calculation
ClassDef(AliHFEtrdPIDqa, 3) // QA class for TRD PID
};
listTPCnsigma->SetOwner();
TH2 *hptr = NULL;
- for(Int_t ispec = 0; ispec < 4; ispec++){
- for(Int_t istep = 0; istep < 2; istep++){
- hptr = MakeTRDspectrumTM(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
- if(hptr){
- hptr->SetName(Form("hTRDtm%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
- listTM->Add(hptr);
+ TList *lctm, *lclike, *lccharge, *lcTPCsig;
+ for(Int_t icent = -1; icent < 11; icent++){
+ lctm = new TList;
+ lctm->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
+ lctm->SetOwner();
+ listTM->Add(lctm);
+ lclike = new TList;
+ lclike->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
+ lclike->SetOwner();
+ listLike->Add(lclike);
+ lccharge = new TList;
+ lccharge->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
+ lccharge->SetOwner();
+ listCharge->Add(lccharge);
+ lcTPCsig = new TList;
+ lcTPCsig->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
+ lcTPCsig->SetOwner();
+ listTPCnsigma->Add(lcTPCsig);
+ for(Int_t ispec = 0; ispec < 4; ispec++){
+ for(Int_t istep = 0; istep < 2; istep++){
+ hptr = MakeTRDspectrumTM(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
+ if(hptr){
+ hptr->SetName(Form("hTRDtm%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
+ lctm->Add(hptr);
+ }
+ hptr = MakeTRDlikelihoodDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
+ hptr->SetName(Form("hTRDlike%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
+ lclike->Add(hptr);
+ hptr = MakeTRDchargeDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
+ hptr->SetName(Form("hTRDcharge%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
+ lccharge->Add(hptr);
+ hptr = MakeTPCspectrumNsigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
+ hptr->SetName(Form("hTPCspectrum%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
+ lcTPCsig->Add(hptr);
}
- hptr = MakeTRDlikelihoodDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
- hptr->SetName(Form("hTRDlike%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
- listLike->Add(hptr);
- hptr = MakeTRDchargeDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
- hptr->SetName(Form("hTRDcharge%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
- listCharge->Add(hptr);
- hptr = MakeTPCspectrumNsigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
- hptr->SetName(Form("hTPCspectrum%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
- listTPCnsigma->Add(hptr);
}
}
-
+
b->Add(listTM, "Projections Truncated Mean");
b->Add(listLike, "Projections Likelihood distribution");
b->Add(listCharge, "Projections Tracklet Charge");
Int_t nBinsTPCSigma[5] = {kPIDbins, kPbins, tpcSigmaBins, kSteps, kCentralityBins};
Double_t minTPCSigma[5] = {kMinPID, kMinP, -12., 0., 0.};
Double_t maxTPCSigma[5] = {kMaxPID, kMaxP, 12., 2., 11.};
- fHistos->CreateTHnSparse("hTPCsigma", "TPC sigma; species p [GeV/c]; TPC dEdx - <dE/dx>|_{el} [#sigma]; selection step", 4, nBinsTPCSigma, minTPCSigma, maxTPCSigma);
+ fHistos->CreateTHnSparse("hTPCsigma", "TPC sigma; species p [GeV/c]; TPC dEdx - <dE/dx>|_{el} [#sigma]; selection step", 5, nBinsTPCSigma, minTPCSigma, maxTPCSigma);
fHistos->Sumw2("hTPCsigma");
// Create Monitoring histogram for the Likelihood distribution
Int_t nBinsTRDlike[5] = {kPIDbins, kPbins, trdLikelihoodBins, kSteps, kCentralityBins};
Double_t minTRDlike[5] = {kMinPID, kMinP, 0., 0., 0.};
Double_t maxTRDlike[5] = {kMaxPID, kMaxP, 1., 2., 11.};
- fHistos->CreateTHnSparse("hTRDlikelihood", "TRD Likelihood Distribution; species; p [GeV/c]; TRD electron Likelihood; selection step", 4, nBinsTRDlike, minTRDlike, maxTRDlike);
+ fHistos->CreateTHnSparse("hTRDlikelihood", "TRD Likelihood Distribution; species; p [GeV/c]; TRD electron Likelihood; selection step", 5, nBinsTRDlike, minTRDlike, maxTRDlike);
fHistos->Sumw2("hTRDlikelihood");
// Create Monitoring histogram for the TRD total charge
const Int_t kTRDchargeBins = 1000;
Int_t nBinsTRDcharge[5] = {kPIDbins, kPbins, kTRDchargeBins, kSteps, kCentralityBins};
Double_t minTRDcharge[5] = {kMinPID, kMinP, 0., 0., 0.};
Double_t maxTRDcharge[5] = {kMaxPID, kMaxP, 100000., 2., 11.};
- fHistos->CreateTHnSparse("hTRDcharge", "Total TRD charge; species; p [GeV/c]; TRD charge [a.u.]; selection step", 4, nBinsTRDcharge, minTRDcharge, maxTRDcharge);
+ fHistos->CreateTHnSparse("hTRDcharge", "Total TRD charge; species; p [GeV/c]; TRD charge [a.u.]; selection step", 5, nBinsTRDcharge, minTRDcharge, maxTRDcharge);
fHistos->Sumw2("hTRDcharge");
// Monitoring of the TRD truncated mean according to version 1
const Int_t kTRDtmBins = 1000;
Int_t nBinsTRDtm[5] = {kPIDbins, kPbins, kTRDtmBins, kSteps, kCentralityBins};
Double_t minTRDtm[5] = {kMinPID, kMinP, 0., 0., 0.};
Double_t maxTRDtm[5] = {kMaxPID, kMaxP, 20000., 2., 11.};
- fHistos->CreateTHnSparse("hTRDtruncatedMean", "TRD truncated Mean; species; p [GeV/c]; TRD signal [a.u.]; selection step", 4, nBinsTRDtm, minTRDtm, maxTRDtm);
+ fHistos->CreateTHnSparse("hTRDtruncatedMean", "TRD truncated Mean; species; p [GeV/c]; TRD signal [a.u.]; selection step", 5, nBinsTRDtm, minTRDtm, maxTRDtm);
fHistos->Sumw2("hTRDtruncatedMean");
}
}
for(UInt_t ily = 0; ily < 6; ily++){
container[2] = trdpid ? trdpid->GetChargeLayer(track->GetRecTrack(), ily, anatype) : 0;
+ if(container[2] < 1e-3) continue; // Filter out 0 entries
fHistos->Fill("hTRDcharge", container);
}
}
//_________________________________________________________
-TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species){
+TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
//
// Get the TPC control histogram for the TRD selection step (either before or after PID)
//
// cut on species (if available)
histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
}
+ TString centname, centtitle;
+ Bool_t hasCentralityInfo = kTRUE;
+ if(centralityClass > -1){
+ if(histo->GetNdimensions() < 5){
+ AliError("Centrality Information not available");
+ centname = centtitle = "MinBias";
+ hasCentralityInfo = kFALSE;
+ } else {
+ // Project centrality classes
+ // -1 is Min. Bias
+ histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
+ centname = Form("Cent%d", centralityClass);
+ centtitle = Form("Centrality %d", centralityClass);
+ }
+ } else {
+ histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
+ centname = centtitle = "MinBias";
+ hasCentralityInfo = kTRUE;
+ }
histo->GetAxis(3)->SetRange(step + 1, step + 1);
TH2 *hSpec = histo->Projection(2, 1);
TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
- TString histname = Form("hSigmaTPC%s%s", specID.Data(), stepname.Data());
- TString histtitle = Form("TPC Sigma for %s %s PID", speciesname.Data(), stepname.Data());
+ TString histname = Form("hSigmaTPC%s%s%s", specID.Data(), stepname.Data(), centname.Data());
+ TString histtitle = Form("TPC Sigma for %s %s PID %s", speciesname.Data(), stepname.Data(), centtitle.Data());
hSpec->SetName(histname.Data());
hSpec->SetTitle(histtitle.Data());
// Unset range on the original histogram
histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
- histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
+ histo->GetAxis(3)->SetRange(0, histo->GetAxis(3)->GetNbins());
+ if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
return hSpec;
}
//_________________________________________________________
-TH2 *AliHFEtrdPIDqaV1::MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species){
+TH2 *AliHFEtrdPIDqaV1::MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
//
// Get the TPC control histogram for the TRD selection step (either before or after PID)
//
// cut on species (if available)
histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
}
+ TString centname, centtitle;
+ Bool_t hasCentralityInfo = kTRUE;
+ if(centralityClass > -1){
+ if(histo->GetNdimensions() < 5){
+ AliError("Centrality Information not available");
+ centname = centtitle = "MinBias";
+ hasCentralityInfo= kFALSE;
+ } else {
+ // Project centrality classes
+ // -1 is Min. Bias
+ histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
+ centname = Form("Cent%d", centralityClass);
+ centtitle = Form("Centrality %d", centralityClass);
+ }
+ } else {
+ histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
+ centname = centtitle = "MinBias";
+ hasCentralityInfo= kFALSE;
+ }
histo->GetAxis(3)->SetRange(step + 1, step + 1);
TH2 *hSpec = histo->Projection(2, 1);
TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
- TString histname = Form("hTMTRD%s%s", specID.Data(), stepname.Data());
- TString histtitle = Form("TRD Truncated Mean for %s %s PID", speciesname.Data(), stepname.Data());
+ TString histname = Form("hTMTRD%s%s%s", specID.Data(), stepname.Data(), centname.Data());
+ TString histtitle = Form("TRD Truncated Mean for %s %s PID %s", speciesname.Data(), stepname.Data(), centtitle.Data());
hSpec->SetName(histname.Data());
hSpec->SetTitle(histtitle.Data());
// Unset range on the original histogram
histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
+ if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
return hSpec;
}
//_________________________________________________________
-TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
+TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
//
// Make Histogram for TRD Likelihood distribution
//
// cut on species (if available)
histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
}
+ TString centname, centtitle;
+ Bool_t hasCentralityInfo = kTRUE;
+ if(centralityClass > -1){
+ if(histo->GetNdimensions() < 5){
+ AliError("Centrality Information not available");
+ centname = centtitle = "MinBias";
+ hasCentralityInfo= kFALSE;
+ } else {
+ // Project centrality classes
+ // -1 is Min. Bias
+ histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
+ centname = Form("Cent%d", centralityClass);
+ centtitle = Form("Centrality %d", centralityClass);
+ }
+ } else {
+ histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
+ centname = centtitle = "MinBias";
+ hasCentralityInfo= kTRUE;
+ }
histo->GetAxis(3)->SetRange(step + 1, step + 1);
TH2 *hSpec = histo->Projection(2, 1);
TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
- TString histname = Form("hLikeElTRD%s%s", specID.Data(), stepname.Data());
- TString histtitle = Form("TRD electron Likelihood for %s %s PID", speciesname.Data(), stepname.Data());
+ TString histname = Form("hLikeElTRD%s%s%s", specID.Data(), stepname.Data(),centname.Data());
+ TString histtitle = Form("TRD electron Likelihood for %s %s PID %s", speciesname.Data(), stepname.Data(),centtitle.Data());
hSpec->SetName(histname.Data());
hSpec->SetTitle(histtitle.Data());
// Unset range on the original histogram
histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
+ if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
return hSpec;
}
//_________________________________________________________
-TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
+TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
//
// Make Histogram for TRD Likelihood distribution
//
// cut on species (if available)
histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
}
+ TString centname, centtitle;
+ Bool_t hasCentralityInfo = kTRUE;
+ if(centralityClass > -1){
+ if(histo->GetNdimensions() < 5){
+ AliError("Centrality Information not available");
+ centname = centtitle = "MinBias";
+ hasCentralityInfo= kFALSE;
+ } else {
+ // Project centrality classes
+ // -1 is Min. Bias
+ histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
+ centname = Form("Cent%d", centralityClass);
+ centtitle = Form("Centrality %d", centralityClass);
+ }
+ } else {
+ histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
+ centname = centtitle = "MinBias";
+ hasCentralityInfo= kTRUE;
+ }
histo->GetAxis(3)->SetRange(step + 1, step + 1);
TH2 *hSpec = histo->Projection(2, 1);
TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
- TString histname = Form("hChargeTRD%s%s", specID.Data(), stepname.Data());
- TString histtitle = Form("TRD total charge for %s %s PID", speciesname.Data(), stepname.Data());
+ TString histname = Form("hChargeTRD%s%s%s", specID.Data(), stepname.Data(), centname.Data());
+ TString histtitle = Form("TRD total charge for %s %s PID %s", speciesname.Data(), stepname.Data(), centtitle.Data());
hSpec->SetName(histname.Data());
hSpec->SetTitle(histtitle.Data());
// Unset range on the original histogram
histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
+ if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
return hSpec;
}
virtual void Initialize();
virtual void ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step);
- TH2 *MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species = -1);
- TH2 *MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species = -1);
- TH2 *MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species = -1);
- TH2 *MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species = -1);
+ AliHFEcollection *GetListOfHistograms() const { return fHistos; }
+
+ TH2 *MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species = -1, Int_t centralityClass = -1);
+ TH2 *MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species = -1, Int_t centralityClass = -1);
+ TH2 *MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species = -1, Int_t centralityClass = -1);
+ TH2 *MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species = -1, Int_t centralityClass = -1);
protected:
void ProcessESDtrack(const AliESDtrack *esdtrack, AliHFEdetPIDqa::EStep_t step, Int_t species);
void ProcessAODtrack(const AliAODTrack *aodtrack, AliHFEdetPIDqa::EStep_t step, Int_t species);
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff