// Matus Kalisky <matus.kalisky@cern.ch>
//
+#include <sstream>
+
#include <TH1F.h>
#include <TList.h>
//
if (fOutput) delete fOutput;
- if (fColl) delete fColl;
- if (fCollMC) delete fCollMC;
+ //if (fColl) delete fColl;
+ //if (fCollMC) delete fCollMC;
if (fEvents) delete fEvents;
}
//__________________________________________________________
// the V0s tagged by the V0 tender
//
- const char * type[4] = {"Gamma", "K0", "Lambda", "ALambda"};
- char name[256] = "";
+ const TString type[4] = {"Gamma", "K0", "Lambda", "ALambda"};
+ TString name;
Int_t nV0s = fInputEvent->GetNumberOfV0s();
for(Int_t i=0; i<nV0s; ++i){
if(!esdV0) continue;
if(!esdV0->GetOnFlyStatus()) continue; // Take only V0s from the On-the-fly v0 finder
Int_t pid = GetTenderPidV0(esdV0);
- if(pid < 0) continue;
+ //if(pid < 0) continue;
+ pid = 2;
fColl->Fill("h_NumberOf_V0s", pid);
- sprintf(name, "h_%s_pt", type[pid]);
Float_t pT = esdV0->Pt();
+ name = "h_" + type[pid] + "_pt";
fColl->Fill(name, pT);
- Float_t mass = MassV0(esdV0, pid);
- sprintf(name, "h_%s_mass", type[pid]);
+ Float_t mass = MassV0(esdV0, pid);
+ name = "h_" + type[pid] + "_mass";
+ //printf(" -D: name: %s \n", name.Data());
fColl->Fill(name, mass);
}
// produce some check plots for V0 tender tagged single tracks
//
- const char * type[3] = {"Electron", "Pion", "Proton"};
- char name[256] = "";
+ const TString type[3] = {"Electron", "Pion", "Proton"};
+ TString name;
Int_t nTracks = fInputEvent->GetNumberOfTracks();
for(Int_t i=0; i<nTracks; ++i){
if(pid < 0) continue;
fColl->Fill("h_NumberOfDaughters", pid*1.0);
Float_t pT = track->Pt();
- sprintf(name, "h_%s_pt", type[pid]);
+ name = "h_" + type[pid] + "_pt";
fColl->Fill(name, pT);
}
//
const Int_t pid2pdg[4] = {22, 310, 3122, -3122};
- const char * type[4] = {"Gamma", "K0", "Lambda", "ALambda"};
- char name[256] = "";
+ const TString type[4] = {"Gamma", "K0", "Lambda", "ALambda"};
+ TString name;
Int_t nV0s = fInputEvent->GetNumberOfV0s();
Int_t nTracks = fInputEvent->GetNumberOfTracks();
Int_t pdg = m->PdgCode();
if((lPm == lNm) && (pdg == pid2pdg[pid])) id = kTRUE;
- if(id) sprintf(name, "h_%s_pt_S", type[pid]);
- else sprintf(name, "h_%s_pt_B", type[pid]);
+ if(id) name = "h_" + type[pid] + "_pt_S";
+ else name = "h_" + type[pid] + "_pt_B";
Float_t pT = esdV0->Pt();
fCollMC->Fill(name, pT);
- if(id) sprintf(name, "h_%s_mass_S", type[pid]);
- else sprintf(name, "h_%s_mass_B", type[pid]);
+ if(id) name = "h_" + type[pid] + "_mass_S";
+ else name = "h_" + type[pid] + "_mass_B";
Float_t mass = MassV0(esdV0, pid);
fCollMC->Fill(name, mass);
//
const Int_t pid2pdg [3] = {11, 211, 2212};
- const char * type[3] = {"Electron", "Pion", "Proton"};
- char name[256] = "";
+ const TString type[3] = {"Electron", "Pion", "Proton"};
+ TString name;
Int_t nTracks = fInputEvent->GetNumberOfTracks();
if(!mcp) continue;
Int_t pdg = TMath::Abs(mcp->PdgCode());
if(pdg == pid2pdg[pid]) id = kTRUE;
- if(id) sprintf(name, "h_%s_pt_S", type[pid]);
- else sprintf(name, "h_%s_pt_B", type[pid]);
+ if(id) name = "h_" + type[pid] + "_pt_S";
+ else name = "h_" + type[pid] + "_pt_B";
fCollMC->Fill(name, pT);
}
}
* provided "as is" without express or implied warranty. *
**************************************************************************/
//
-// Task fir checking the performance of the V0 tender
+// Task for checking the performance of the V0 tender
//
//
// Authors:
, fCollMC(0x0)
, fV0cuts(0x0)
, fPrimaryVertex(0x0)
- , fpdgV0(0)
- , fpdgP(0)
- , fpdgN(0)
+ , fpdgV0(-1)
+ , fpdgP(-1)
+ , fpdgN(-1)
, fEvents(0x0)
{
//
//
if (fOutput) delete fOutput;
- if (fColl) delete fColl;
- if (fCollMC) delete fCollMC;
+ //if (fColl) delete fColl;
+ //if (fCollMC) delete fCollMC;
if (fV0cuts) delete fV0cuts;
if (fPrimaryVertex) delete fPrimaryVertex;
if (fEvents) delete fEvents;
fCollMC->CreateTH1F("h_Electron_pt_B", "MC-B: p_{T} spectrum of tagged electrons; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
fCollMC->CreateTH1F("h_Pion_pt_B", "MC-B: p_{T} spectrum of tagged pions; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
fCollMC->CreateTH1F("h_Proton_pt_B", "MC-B: p_{T} spectrum of tagged protons; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
-
+ // background study - MC only
+ fCollMC->CreateTH1Fvector1(6, "h_Gamma_Bg", "MC gamma bg. - different sources; p_{T} (GeV/c), counts", 20, 0.1, 20, 0);
+ fCollMC->CreateTH1Fvector1(6, "h_K0_Bg", "MC K0 bg. - different sources; p_{T} (GeV/c), counts", 20, 0.1, 20, 0);
+ fCollMC->CreateTH1Fvector1(6, "h_Lambda_Bg", "MC gamma bg. - different sources; p_{T} (GeV/c), counts", 20, 0.1, 20, 0);
+ fCollMC->CreateTH1Fvector1(6, "h_ALambda_Bg", "MC gamma bg. - different sources; p_{T} (GeV/c), counts", 20, 0.1, 20, 0);
+
TList *tmp = fColl->GetList();
tmp->SetName(fColl->GetName());
if(fMCEvent && !mcHandler->TreeK() ) return;
if(fMCEvent && !mcHandler->TreeTR() ) return;
- fPrimaryVertex = new AliKFVertex(*(fInputEvent->GetPrimaryVertex()));
- if(!fPrimaryVertex) return;
-
+ AliKFVertex *primaryVertex = new AliKFVertex(*(fInputEvent->GetPrimaryVertex()));
+ if(!primaryVertex) return;
+
+ fV0cuts->SetPrimaryVertex(primaryVertex);
fV0cuts->SetEvent(fInputEvent);
- fV0cuts->SetPrimaryVertex(fPrimaryVertex);
-
+
ProcessV0s();
+
+ if (primaryVertex) delete primaryVertex;
fEvents->Fill(1.1);
PostData(1, fEvents);
PostData(2, fOutput);
// the V0s tagged by the V0 tender
//
- const char * type[4] = {"Gamma", "K0", "Lambda", "ALambda"};
- char name[256] = "";
+ const TString type[4] = {"Gamma", "K0", "Lambda", "ALambda"};
+ TString name;
Int_t nV0s = fInputEvent->GetNumberOfV0s();
for(Int_t i=0; i<nV0s; ++i){
if(!esdV0) continue;
if(!esdV0->GetOnFlyStatus()) continue; // Take only V0s from the On-the-fly v0 finder
+ ResetPDGcodes();
+
// New standalone V0 selection software
if( ! (fV0cuts->ProcessV0(esdV0, fpdgV0, fpdgP, fpdgN)) ) return;
if(pid <= -1) continue;
fColl->Fill("h_NumberOf_V0s", pid);
- sprintf(name, "h_%s_pt", type[pid]);
+ name = "h_" + type[pid] + "_pt";
Float_t pT = esdV0->Pt();
fColl->Fill(name, pT);
Float_t mass = MassV0(esdV0, pid);
- sprintf(name, "h_%s_mass", type[pid]);
+ name = "h_" + type[pid] + "_mass";
fColl->Fill(name, mass);
ProcessDaughters(esdV0);
if(fMCEvent){
ProcessV0sMC(esdV0);
ProcessDaughtersMC(esdV0);
+ ProcessBackground(esdV0);
}
}
// produce some check plots for V0 tender tagged single tracks
//
- const char * type[3] = {"Electron", "Pion", "Proton"};
- char name[256] = "";
+ const TString type[3] = {"Electron", "Pion", "Proton"};
+ TString name;
if(!v0) return;
if(pid < 0) continue;
fColl->Fill("h_NumberOfDaughters", pid*1.0);
Float_t pT = d[i]->Pt();
- sprintf(name, "h_%s_pt", type[pid]);
+ name = "h_" + type[pid] + "_pt";
fColl->Fill(name, pT);
-
}
}
//__________________________________________________________
// check all V0tender selected V0 on their true identity
//
- const char * type[4] = {"Gamma", "K0", "Lambda", "ALambda"};
- char name[256] = "";
-
+ const TString type[4] = {"Gamma", "K0", "Lambda", "ALambda"};
+ TString name;
Int_t pid = PDGtoPIDv0(fpdgV0);
if(pid < 0) return;
if((lPm == lNm) && (pdg == fpdgV0)) id = kTRUE;
}
- if(id) sprintf(name, "h_%s_pt_S", type[pid]);
- else sprintf(name, "h_%s_pt_B", type[pid]);
+ if(id) name = "h_" + type[pid] + "_pt_S";
+ else name = "h_" + type[pid] + "_pt_B";
Float_t pT = v0->Pt();
fCollMC->Fill(name, pT);
- if(id) sprintf(name, "h_%s_mass_S", type[pid]);
- else sprintf(name, "h_%s_mass_B", type[pid]);
+ if(id) name = "h_" + type[pid] + "_mass_S";
+ else name = "h_" + type[pid] + "_mass_B";
Float_t mass = MassV0(v0, pid);
fCollMC->Fill(name, mass);
// from primary vertex or pion dalitz deca or true gamma conversion) !!!!
//
- const char * type[3] = {"Electron", "Pion", "Proton"};
- char name[256] = "";
+ const TString type[3] = {"Electron", "Pion", "Proton"};
+ TString name;
if(!v0) return;
Int_t pdgMC = mcp->PdgCode();
//printf(" * pid: %i, pdg: %i, pdgMC: %i \n", pid, pdg, pdgMC);
if(pdgMC == pdg) id = kTRUE;
- if(id) sprintf(name, "h_%s_pt_S", type[pid]);
- else sprintf(name, "h_%s_pt_B", type[pid]);
+ if(id) name = "h_" + type[pid] + "_pt_S";
+ else name = "h_" + type[pid] + "_pt_B";
fCollMC->Fill(name, pT);
}
+}
+//__________________________________________________________
+void AliAnalysisTaskCheckV0tenderII::ProcessBackground(AliESDv0 * const v0){
+ //
+ // look at the compostition and the properties of the
+ // miss-identified V0
+ //
+
+ if(!v0) return;
+ Float_t pt = v0->Pt();
+ const TString type[4] = {"Gamma", "K0", "Lambda", "ALambda"};
+ TString name;
+ // daughter tracks
+ const Int_t nTracks = fInputEvent->GetNumberOfTracks();
+ AliESDtrack *d[2];
+ Int_t index[2] = {-1, -1};
+ index[0] = v0->GetPindex();
+ index[1] = v0->GetNindex();
+ if(index[0] < 0 || index[1] < 0) return;
+ if(index[0] >= nTracks || index[1] >= nTracks) return;
+ d[0] = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(index[0]));
+ d[1] = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(index[1]));
+ if(!d[0] || !d[1]) return;
+
+ // daughter MC particles
+ Int_t label[2] = {d[0]->GetLabel(), d[1]->GetLabel()};
+ if(label[0] < 0 || label[1] < 0) return;
+ AliMCParticle *dmc[2] = {0x0, 0x0};
+ dmc[0] = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(label[0]));
+ dmc[1] = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(label[1]));
+ if(!dmc[0] || !dmc[1]) return;
+
+ Bool_t primary[2] = {dmc[0]->Particle()->IsPrimary(), dmc[1]->Particle()->IsPrimary()};
+
+ // mother MC particles
+ Int_t labelM[2] = {dmc[0]->GetMother(), dmc[1]->GetMother()};
+ AliMCParticle *mmc[2] = {0x0, 0x0};
+ if(!primary[0]) mmc[0] = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(labelM[0]));
+ if(!primary[1]) mmc[1] = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(labelM[1]));
+
+ // if particle is not primary it MUST have a mother particle
+ if(!primary[0] && !mmc[0]) return;
+ if(!primary[1] && !mmc[1]) return;
+
+ Int_t pdgM[2] = {0, 0};
+ if(mmc[0]) pdgM[0] = mmc[0]->PdgCode();
+ if(mmc[1]) pdgM[1] = mmc[1]->PdgCode();
+
+ // select only miss-identified V0s
+ if(labelM[0] == labelM[1] && pdgM[0] == fpdgV0) return;
+
+ // indices for V0 background histograms
+ // 0 - random match
+ // 1 - gamma
+ // 2 - K0
+ // 3 - lambda
+ // 4 - a-lambda
+ // 5 - other dacay or interaction
+
+ Int_t ixM = -1;
+
+ Int_t typeV0 = (labelM[0] != labelM[1]) ? 0 : 1;
+ if(0 == typeV0) ixM = 0;
+ else{
+ ixM = PDGtoPIDv0(pdgM[0]) + 1;
+ if(ixM < 0) ixM = 5;
+ }
+ name = "h_" + type[PDGtoPIDv0(fpdgV0)] + "_Bg";
+ fCollMC->Fill(name, ixM, pt);
+
+ // now look at the daughter tracks
+
+
}
//__________________________________________________________
Float_t AliAnalysisTaskCheckV0tenderII::MassV0(AliESDv0 * const v0, Int_t id){
return -1;
}
+//__________________________________________________________
+void AliAnalysisTaskCheckV0tenderII::ResetPDGcodes(){
+ //
+ // reset the PDG codes of the V0 and daughter
+ // particles to default values
+ //
+
+ fpdgV0 = -1;
+ fpdgP = -1;
+ fpdgN = -1;
+
+}
* provided "as is" without express or implied warranty. *
**************************************************************************/
//
-// Task for PID QA
-// Using AliHFEpidQA and AliHFEMCpidQA
+// Task for perfomance studies of V0 selection code
// More information can be found in the source file
//
#ifndef ALIANALYSISTASKCHECKV0TENDERII_H
void ProcessDaughters(AliESDv0 * const v0);
void ProcessV0sMC(AliESDv0 * const v0);
void ProcessDaughtersMC(AliESDv0 * const v0);
+ void ProcessBackground(AliESDv0 * const v0);
private:
AliAnalysisTaskCheckV0tenderII(const AliAnalysisTaskCheckV0tenderII &ref);
Int_t PDGtoPIDv0(Int_t pdgV0) const;
Int_t PDGtoPID(Int_t pdg) const;
-
+
+ void ResetPDGcodes();
+
TList *fOutput; //! Container for output histos
AliHFEcollection *fColl; //! collection of Data output
AliHFEcollection *fCollMC; //! collection of MC output
AliESDv0KineCuts *fV0cuts; //! standalone V0 selection class
AliKFVertex *fPrimaryVertex; //! primary vertex of the current event
- Int_t fpdgV0; // PDG code of teh reconstructed V0
+ Int_t fpdgV0; // PDG code of the reconstructed V0
Int_t fpdgP; // PDG code of the positive daughter (sign corrected)
Int_t fpdgN; // PDG code of the negative daughter (sign coreccted)
, fCuts(NULL)
, fTaggedTrackCuts(NULL)
, fCleanTaggedTrack(kFALSE)
+ , fVariablesTRDTaggedTrack(kFALSE)
, fCutspreselect(NULL)
, fSecVtx(NULL)
, fElecBackGround(NULL)
, fCuts(NULL)
, fTaggedTrackCuts(NULL)
, fCleanTaggedTrack(kFALSE)
+ , fVariablesTRDTaggedTrack(kFALSE)
, fCutspreselect(NULL)
, fSecVtx(NULL)
, fElecBackGround(NULL)
, fCuts(NULL)
, fTaggedTrackCuts(NULL)
, fCleanTaggedTrack(ref.fCleanTaggedTrack)
+ , fVariablesTRDTaggedTrack(ref.fVariablesTRDTaggedTrack)
, fCutspreselect(NULL)
, fSecVtx(NULL)
, fElecBackGround(NULL)
target.fCuts = fCuts;
target.fTaggedTrackCuts = fTaggedTrackCuts;
target.fCleanTaggedTrack = fCleanTaggedTrack;
+ target.fVariablesTRDTaggedTrack = fVariablesTRDTaggedTrack;
target.fCutspreselect = fCutspreselect;
target.fSecVtx = fSecVtx;
target.fElecBackGround = fElecBackGround;
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->BinLogAxis("TPCclusters2_1_Signal", 0);
+ fQACollection->BinLogAxis("TPCclusters2_0_Signal", 0);
+ fQACollection->BinLogAxis("TPCclusters2_1_Selected", 0);
+ fQACollection->BinLogAxis("TPCclusters2_0_Selected", 0);
+
InitPIDperformanceQA();
InitContaminationQA();
fQA->Add(fQACollection->GetList());
fTaggedTrackAnalysis = new AliHFEtaggedTrackAnalysis;
fTaggedTrackAnalysis->SetCuts(fTaggedTrackCuts);
fTaggedTrackAnalysis->SetClean(fCleanTaggedTrack);
+ if(fPIDqa->HasHighResolutionHistos())
+ fTaggedTrackAnalysis->GetPIDqa()->SetHighResolutionHistos();
fTaggedTrackAnalysis->SetPID(fPID);
+ fTaggedTrackAnalysis->SetVariablesTRD(fVariablesTRDTaggedTrack);
fTaggedTrackAnalysis->InitContainer();
fOutput->Add(fTaggedTrackAnalysis->GetContainer());
fQA->Add(fTaggedTrackAnalysis->GetPIDQA());
}
fPID->SetESDpid(workingPID);
if(fPIDpreselect) fPIDpreselect->SetESDpid(workingPID);
-
+
ProcessESD();
}
// Done!!!
fMCQA->SetGenEventHeader(fMCEvent->GenEventHeader());
fMCQA->Init();
- Int_t nMCTracks = fMCEvent->Stack()->GetNtrack();
-
// loop over all tracks for decayed electrons
- for (Int_t igen = 0; igen < nMCTracks; igen++){
+ for (Int_t igen = 0; igen < fMCEvent->GetNumberOfTracks(); igen++){
TParticle* mcpart = fMCEvent->Stack()->Particle(igen);
+ if(!mcpart) continue;
fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kCharm);
fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kBeauty);
fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kCharm);
return;
}
+ // Set magnetic field if V0 task on
+ if(fTaggedTrackAnalysis) fTaggedTrackAnalysis->SetMagneticField(fESD->GetMagneticField());
+
// Do event Normalization
Double_t eventContainer[3];
eventContainer[0] = 0.0;
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.);
+ }
}
// RecKine: ITSTPC cuts
// RecPrim
+ if(track->GetKinkIndex(0) != 0) continue; // Quick and dirty fix to reject both kink mothers and daughters
if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue;
// HFEcuts: ITS layers cuts
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(fBackGroundFactorsFunction) {
- Double_t weightBackGround = fBackGroundFactorsFunction->Eval(TMath::Abs(track->P()));
- if(weightBackGround < 0.0) weightBackGround = 0.0;
+ Double_t weightBackGround = fBackGroundFactorsFunction->Eval(TMath::Abs(track->P()));
+ if(weightBackGround < 0.0) weightBackGround = 0.0;
else if(weightBackGround > 1.0) weightBackGround = 1.0;
// weightBackGround as special weight
fVarManager->FillContainer(fContainer, "hadronicBackground", 1, kFALSE, weightBackGround);
void SetHFECuts(AliHFEcuts * const cuts) { fCuts = cuts; };
void SetTaggedTrackCuts(AliHFEcuts * const cuts) { fTaggedTrackCuts = cuts; }
void SetCleanTaggedTrack(Bool_t clean) { fCleanTaggedTrack = clean; };
+ void SetVariablesTRDTaggedTrack(Bool_t variablesTRD) { fVariablesTRDTaggedTrack = variablesTRD; };
void SetHFECutsPreselect(AliHFEcuts * const cuts) { fCutspreselect = cuts; };
void SetHFEElecBackGround(AliHFEelecbackground * const elecBackGround) { fElecBackGround = elecBackGround; };
void SetQAOn(Int_t qaLevel) { SETBIT(fQAlevel, qaLevel); };
AliHFEcuts *fCuts; // Cut Collection
AliHFEcuts *fTaggedTrackCuts; // Cut Collection for V0 tagged tracks
Bool_t fCleanTaggedTrack; // Loose cleaning of the V0 tagged tracks electron
+ Bool_t fVariablesTRDTaggedTrack; // Take the variables at the TRD for the V0 tagged tracks electron
AliHFEcuts *fCutspreselect; // Cut Collection for pre-selected tracks
AliHFEsecVtx *fSecVtx; //! Secondary Vertex Analysis
AliHFEelecbackground *fElecBackGround;//! Background analysis
// process the selected V0 daughter tracks
//
- Char_t hname[256] = "";
- const Char_t *typeName[5] = {"Electron", "PionK0", "PionL", "Kaon", "Proton"};
+ TString hname;
+ const TString typeName[5] = {"Electron", "PionK0", "PionL", "Kaon", "Proton"};
const Int_t typePID[5] = {0, 2, 2, 3, 4};
if(!fMC) return kFALSE;
//Int_t pdgM = mcpM->PdgCode();
// all candidates
- sprintf(hname, "h_%s", typeName[type]);
+ hname = "h_" + typeName[type];
fColl->Fill(hname, p);
Int_t pidD = PDGtoPIDdaughter(pdgD);
// all misidentified candidates
- sprintf(hname, "h_mis_%s", typeName[type]);
+ hname = "h_mis_" + typeName[type];
if(typePID[type] != pidD){
fColl->Fill(hname, p);
}
// for every particle fill detailed information about the missidentified particles
- sprintf(hname, "h_tag_%s", typeName[type]);
+ hname = "h_tag_" + typeName[type];
Int_t aliPID = PDGtoAliPID(pdgD);
fColl->Fill(hname, aliPID, p);
this->~AliHFEcontainer(); // cleanup old object before creating the new onwe
TNamed::operator=(ref);
fContainers = new THashList();
+ fCorrelationMatrices = NULL;
fNVars = ref.fNVars;
- AliCFContainer *ctmp = NULL;
- for(Int_t ien = 0; ien < ref.fContainers->GetEntries(); ien++){
- ctmp = dynamic_cast<AliCFContainer *>(ref.fContainers->At(ien));
- fContainers->Add(new AliCFContainer(*ctmp));
- }
if(fNVars){
fVariables = new TObjArray(fNVars);
AliHFEvarInfo *vtmp = NULL;
if(vtmp) fVariables->AddAt(new AliHFEvarInfo(*vtmp), ivar);
}
} else {
+ // No varible defined, do not try to copy anything
fVariables = NULL;
+ return *this;
+ }
+
+ // Reference contains content, try copying also the containers and the correlation matrices
+ fContainers = new THashList();
+ AliCFContainer *ctmp = NULL;
+ for(Int_t ien = 0; ien < ref.fContainers->GetEntries(); ien++){
+ ctmp = dynamic_cast<AliCFContainer *>(ref.fContainers->At(ien));
+ if(ctmp) fContainers->Add(new AliCFContainer(*ctmp));
}
// Copy also correlation matrices
if(ref.fCorrelationMatrices){
AliDebug(1, Form("Cut Step %s: Number of cut objects: %d", GetName(), fCuts->GetEntriesFast()));
if(!fCuts->GetEntriesFast()) return kTRUE;
Bool_t isSelected = kTRUE;
+ AliAnalysisCuts *cuts = NULL;
for(Int_t iCut = 0; iCut < fCuts->GetEntriesFast(); iCut++){
- if(!(dynamic_cast<AliAnalysisCuts *>(fCuts->UncheckedAt(iCut)))->IsSelected(o)) isSelected = kFALSE;
+ cuts = dynamic_cast<AliAnalysisCuts *>(fCuts->UncheckedAt(iCut));
+ if(cuts && (!cuts->IsSelected(o))) isSelected = kFALSE;
}
AliDebug(1, Form("Accepted: %s", isSelected ? "yes" : "no"));
return isSelected;
AliHFEcuts::AliHFEcuts():
TNamed(),
fRequirements(0),
- fTPCiter1(kFALSE),
+ fTPCclusterDef(AliHFEextraCuts::kFound),
+ fTPCratioDef(AliHFEextraCuts::kFoundOverFindable),
fMinClustersTPC(0),
fMinClustersITS(0),
fMinTrackletsTRD(0),
fMinClusterRatioTPC(0.),
fSigmaToVtx(0.),
fVertexRangeZ(20.),
+ fTOFPIDStep(kFALSE),
fHistQA(0x0),
fCutList(0x0),
fDebugLevel(0)
AliHFEcuts::AliHFEcuts(const Char_t *name, const Char_t *title):
TNamed(name, title),
fRequirements(0),
- fTPCiter1(kFALSE),
+ fTPCclusterDef(AliHFEextraCuts::kFound),
+ fTPCratioDef(AliHFEextraCuts::kFoundOverFindable),
fMinClustersTPC(0),
fMinClustersITS(0),
fMinTrackletsTRD(0),
fMinClusterRatioTPC(0.),
fSigmaToVtx(0.),
fVertexRangeZ(20.),
+ fTOFPIDStep(kFALSE),
fHistQA(0x0),
fCutList(0x0),
fDebugLevel(0)
AliHFEcuts::AliHFEcuts(const AliHFEcuts &c):
TNamed(c),
fRequirements(c.fRequirements),
- fTPCiter1(c.fTPCiter1),
+ fTPCclusterDef(c.fTPCclusterDef),
+ fTPCratioDef(c.fTPCratioDef),
fMinClustersTPC(0),
fMinClustersITS(0),
fMinTrackletsTRD(0),
fMinClusterRatioTPC(0),
fSigmaToVtx(0),
fVertexRangeZ(20.),
+ fTOFPIDStep(kFALSE),
fHistQA(0x0),
fCutList(0x0),
fDebugLevel(0)
AliHFEcuts &target = dynamic_cast<AliHFEcuts &>(c);
target.fRequirements = fRequirements;
- target.fTPCiter1 = fTPCiter1;
+ target.fTPCclusterDef = fTPCclusterDef;
+ target.fTPCratioDef = fTPCratioDef;
target.fMinClustersTPC = fMinClustersTPC;
target.fMinClustersITS = fMinClustersITS;
target.fMinTrackletsTRD = fMinTrackletsTRD;
target.fMinClusterRatioTPC = fMinClusterRatioTPC;
target.fSigmaToVtx = fSigmaToVtx;
target.fVertexRangeZ = fVertexRangeZ;
+ target.fTOFPIDStep = fTOFPIDStep;
target.fDebugLevel = 0;
memcpy(target.fProdVtx, fProdVtx, sizeof(Double_t) * 4);
//trackQuality->SetMaxCovDiagonalElements(2., 2., 0.5, 0.5, 2);
AliHFEextraCuts *hfecuts = new AliHFEextraCuts("fCutsHFElectronGroupTPC","Extra cuts from the HFE group");
- if(fMinClusterRatioTPC > 0.) hfecuts->SetClusterRatioTPC(fMinClusterRatioTPC);
hfecuts->SetDebugLevel(fDebugLevel);
// Set the cut in the TPC number of clusters
- if(fTPCiter1){
- hfecuts->SetMinNClustersTPC(fMinClustersTPC);
- hfecuts->SetTPCIter1(kTRUE);
- }
- else
- trackQuality->SetMinNClusterTPC(fMinClustersTPC);
+ hfecuts->SetMinNClustersTPC(fMinClustersTPC, fTPCclusterDef);
+ hfecuts->SetClusterRatioTPC(fMinClusterRatioTPC, fTPCratioDef);
AliCFTrackKineCuts *kineCuts = new AliCFTrackKineCuts((Char_t *)"fCutsKineRec", (Char_t *)"REC Kine Cuts");
kineCuts->SetPtRange(fPtRange[0], fPtRange[1]);
AliDebug(2, "Called\n");
AliHFEextraCuts *hfecuts = new AliHFEextraCuts("fCutsHFElectronGroupTRD","Extra cuts from the HFE group");
if(fMinTrackletsTRD > 0.) hfecuts->SetMinTrackletsTRD(fMinTrackletsTRD);
+ if(fTOFPIDStep) hfecuts->SetTOFPID(kTRUE);
if(IsQAOn()) hfecuts->SetQAOn(fHistQA);
hfecuts->SetDebugLevel(fDebugLevel);
void SetPtRange(Double_t ptmin, Double_t ptmax){fPtRange[0] = ptmin; fPtRange[1] = ptmax;};
inline void SetProductionVertex(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax);
inline void SetSigmaToVertex(Double_t sig);
- void SetTPCiter1(Bool_t iter1) { fTPCiter1 = iter1; }
- void SetVertexRange(Double_t zrange){fVertexRangeZ = zrange;};
+ void SetTPCmodes(AliHFEextraCuts::ETPCclusterDef_t clusterDef, AliHFEextraCuts::ETPCclrDef_t ratioDef) {
+ fTPCclusterDef= clusterDef;
+ fTPCratioDef = ratioDef;
+ }
+ void SetVertexRange(Double_t zrange){fVertexRangeZ = zrange;};
+ void SetTOFPIDStep(Bool_t tofPidStep) {fTOFPIDStep = tofPidStep;};
inline void CreateStandardCuts();
static const Char_t* fgkUndefined; // Name for undefined (overflow)
ULong64_t fRequirements; // Bitmap for requirements
- Bool_t fTPCiter1; // TPC iter1
+ AliHFEextraCuts::ETPCclusterDef_t fTPCclusterDef; // TPC cluster definition
+ AliHFEextraCuts::ETPCclrDef_t fTPCratioDef; // TPC cluster ratio Definition
Double_t fDCAtoVtx[2]; // DCA to Vertex
Double_t fProdVtx[4]; // Production Vertex
Double_t fPtRange[2]; // pt range
Double_t fMinClusterRatioTPC; // Min. Ratio findable / found TPC clusters
Double_t fSigmaToVtx; // Sigma To Vertex
Double_t fVertexRangeZ; // Vertex Range reconstructed
+ Bool_t fTOFPIDStep; // TOF matching step efficiency
TList *fHistQA; //! QA Histograms
// Event Loop
// Filter track, fill Efficiency container
//
+ AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fInputEvent);
+ if(!esdevent){
+ AliError("ESD Event required");
+ return;
+ }
fEfficiency->NewEvent();
fFilter->SetRecEvent(fInputEvent);
if(fMCEvent){
AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
- if(!mcH->InitOk()) return;
- if(!mcH->TreeK()) return;
- if(!mcH->TreeTR()) return;
+ if(mcH &&(!mcH->InitOk())) return;
+ if(mcH &&(!mcH->TreeK())) return;
+ if(mcH &&(!mcH->TreeTR())) return;
fFilter->SetMC(fMCEvent);
FilterMC();
}
- fFilter->FilterTracks(dynamic_cast<AliESDEvent *>(fInputEvent));
+ fFilter->FilterTracks(esdevent);
TObjArray *tracks = fFilter->GetFilteredTracks();
TIterator *iter = tracks->MakeIterator();
fOutput->Fill("hNtracks", tracks->GetEntriesFast());
AliMCParticle *track = NULL;
for(Int_t itrack = 0; itrack < fMCEvent->GetNumberOfTracks(); itrack++){
track = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(itrack));
+ if(!track) continue;
if(!fMCcut->IsSelected(track)) continue;
cont[0] = track->Pt();
cont[1] = track->Eta();
cinvmass2D->cd(12);
if(invmassSplittedTrackosproj2D) invmassSplittedTrackosproj2D->Draw("lego");
-
- ////////////////////////
- // Cut efficiencies
- ////////////////////////
-
- THnSparseF *hsSparseMCe = dynamic_cast<THnSparseF *>(fList->FindObject("CutPassedMC"));
-
- if(hsSparseMCe) {
-
- // init histos
- TAxis *axissources = hsSparseMCe->GetAxis(2);
- Int_t nbsources = axissources->GetNbins();
- TAxis *axiscuts = hsSparseMCe->GetAxis(1);
- Int_t nbcuts = axiscuts->GetNbins();
- TH1D **histopassedcuts = new TH1D*[nbsources*nbcuts];
- Double_t *nbEntriesCuts = new Double_t[nbsources*nbcuts];
-
- //printf("Number of cuts %d\n",nbcuts);
-
- // canvas
- TCanvas * chsSparseMCeeff =new TCanvas("hsSparseMCeeffDebug","hsSparseMCeeffDebug",800,800);
- chsSparseMCeeff->Divide(3,1);
-
- // histos
- for(Int_t sourceid = 0; sourceid < nbsources; sourceid++) {
- hsSparseMCe->GetAxis(2)->SetRange(sourceid+1,sourceid+1);
- for(Int_t cut = 0; cut < nbcuts; cut++){
- hsSparseMCe->GetAxis(1)->SetRange(cut+1,cut+1);
- histopassedcuts[sourceid*nbcuts+cut] = hsSparseMCe->Projection(0);
- hsSparseMCe->GetAxis(1)->SetRange(1,hsSparseMCe->GetAxis(1)->GetNbins());
- }
- hsSparseMCe->GetAxis(2)->SetRange(1,hsSparseMCe->GetAxis(2)->GetNbins());
- }
-
- // calcul efficiencies
- ///////////////////////
- // histos
- for(Int_t sourceid = 0; sourceid < nbsources; sourceid++) {
- // Next is compared to the partner tracked
- for(Int_t cut = 2; cut < nbcuts; cut++){
- nbEntriesCuts[sourceid*nbcuts+cut] = histopassedcuts[sourceid*nbcuts+cut]->GetEntries();
- if(histopassedcuts[sourceid*nbcuts+1]->GetEntries() > 0.0) histopassedcuts[sourceid*nbcuts+cut]->Divide(histopassedcuts[sourceid*nbcuts+1]);
- }
- // First one is if the partner is tracked.
- nbEntriesCuts[sourceid*nbcuts+1] = histopassedcuts[sourceid*nbcuts+1]->GetEntries();
- if(histopassedcuts[sourceid*nbcuts]->GetEntries() > 0.0) histopassedcuts[sourceid*nbcuts+1]->Divide(histopassedcuts[sourceid*nbcuts]);
- // First one is input
- nbEntriesCuts[sourceid*nbcuts] = histopassedcuts[sourceid*nbcuts]->GetEntries();
- }
-
- /////////////
- // ratios
- ////////////
- for(Int_t sourceid = 0; sourceid < nbsources; sourceid++) {
- for(Int_t cut = 1; cut < nbcuts; cut++){
- if(nbEntriesCuts[sourceid*nbcuts] > 0.0) nbEntriesCuts[sourceid*nbcuts+cut] = nbEntriesCuts[sourceid*nbcuts+cut]/nbEntriesCuts[sourceid*nbcuts];
- }
- }
- TH1F *ratioHistoEntriesGamma = new TH1F("ratioHistoEntriesGamma","", nbcuts-1, 0.0, nbcuts-1.0);
- TH1F *ratioHistoEntriesPi0 = new TH1F("ratioHistoEntriesPi0","", nbcuts-1, 0.0, nbcuts-1.0);
- TH1F *ratioHistoEntriesC = new TH1F("ratioHistoEntriesC","", nbcuts-1, 0.0, nbcuts-1.0);
- for(Int_t k = 1; k < nbcuts; k++){
- ratioHistoEntriesGamma->SetBinContent(k,nbEntriesCuts[nbcuts+k]);
- ratioHistoEntriesPi0->SetBinContent(k,nbEntriesCuts[2*nbcuts+k]);
- ratioHistoEntriesC->SetBinContent(k,nbEntriesCuts[4*nbcuts+k]);
- }
- //
- TAxis *xAxisGamma = ratioHistoEntriesGamma->GetXaxis();
- xAxisGamma->SetBinLabel(1,"Partner tracked");
- xAxisGamma->SetBinLabel(2,"Opposite sign");
- xAxisGamma->SetBinLabel(3,"Single Track Cut");
- xAxisGamma->SetBinLabel(4,"Shared Clusters");
- xAxisGamma->SetBinLabel(5,"PID");
- xAxisGamma->SetBinLabel(6,"DCA");
- xAxisGamma->SetBinLabel(7,"Chi^{2}/Ndf");
- xAxisGamma->SetBinLabel(8,"Opening angle");
- xAxisGamma->SetBinLabel(9,"Invariant mass");
- //
- TAxis *xAxisPi0 = ratioHistoEntriesPi0->GetXaxis();
- xAxisPi0->SetBinLabel(1,"Partner tracked");
- xAxisPi0->SetBinLabel(2,"Opposite sign");
- xAxisPi0->SetBinLabel(3,"Single Track Cut");
- xAxisPi0->SetBinLabel(4,"Shared Clusters");
- xAxisPi0->SetBinLabel(5,"PID");
- xAxisPi0->SetBinLabel(6,"DCA");
- xAxisPi0->SetBinLabel(7,"Chi^{2}/Ndf");
- xAxisPi0->SetBinLabel(8,"Opening angle");
- xAxisPi0->SetBinLabel(9,"Invariant mass");
- //
- TAxis *xAxisC = ratioHistoEntriesC->GetXaxis();
- xAxisC->SetBinLabel(1,"Partner tracked");
- xAxisC->SetBinLabel(2,"Opposite sign");
- xAxisC->SetBinLabel(3,"Single Track Cut");
- xAxisC->SetBinLabel(4,"Shared Clusters");
- xAxisC->SetBinLabel(5,"PID");
- xAxisC->SetBinLabel(6,"DCA");
- xAxisC->SetBinLabel(7,"Chi^{2}/Ndf");
- xAxisC->SetBinLabel(8,"Opening angle");
- xAxisC->SetBinLabel(9,"Invariant mass");
- //
- TCanvas * cRatioHistoEntries =new TCanvas("cRatioHistoEntries","cRatioHistoEntries",800,800);
- cRatioHistoEntries->cd(1);
- ratioHistoEntriesGamma->SetStats(0);
- ratioHistoEntriesGamma->Draw();
- ratioHistoEntriesPi0->SetStats(0);
- ratioHistoEntriesPi0->Draw("same");
- ratioHistoEntriesC->SetStats(0);
- //ratioHistoEntriesC->Draw("same");
- TLegend *legEntries = new TLegend(0.4,0.6,0.89,0.89);
- legEntries->AddEntry(ratioHistoEntriesGamma,"#gamma","l");
- legEntries->AddEntry(ratioHistoEntriesPi0,"#pi^{0}","l");
- //legEntries->AddEntry(ratioHistoEntriesC,"c","p");
- legEntries->Draw("same");
-
- ////////////////////
- // plot Debug
- ///////////////////
- Int_t source = 1;
- chsSparseMCeeff->cd(1);
- histopassedcuts[source*nbcuts+0]->SetTitle("#gamma");
- histopassedcuts[source*nbcuts+1]->SetTitle("#gamma");
- histopassedcuts[source*nbcuts+2]->SetTitle("#gamma");
- histopassedcuts[source*nbcuts+3]->SetTitle("#gamma");
- histopassedcuts[source*nbcuts+4]->SetTitle("#gamma");
- histopassedcuts[source*nbcuts+5]->SetTitle("#gamma");
- histopassedcuts[source*nbcuts+6]->SetTitle("#gamma");
- histopassedcuts[source*nbcuts+7]->SetTitle("#gamma");
- histopassedcuts[source*nbcuts+8]->SetTitle("#gamma");
- histopassedcuts[source*nbcuts+9]->SetTitle("#gamma");
- //histopassedcuts[source*nbcuts+0]->SetStats(0);
- histopassedcuts[source*nbcuts+1]->SetStats(0);
- histopassedcuts[source*nbcuts+2]->SetStats(0);
- histopassedcuts[source*nbcuts+3]->SetStats(0);
- histopassedcuts[source*nbcuts+4]->SetStats(0);
- histopassedcuts[source*nbcuts+5]->SetStats(0);
- histopassedcuts[source*nbcuts+6]->SetStats(0);
- histopassedcuts[source*nbcuts+7]->SetStats(0);
- histopassedcuts[source*nbcuts+8]->SetStats(0);
- histopassedcuts[source*nbcuts+9]->SetStats(0);
- //histopassedcuts[source*nbcuts+0]->Draw();
- //histopassedcuts[source*nbcuts+1]->Draw("");
- histopassedcuts[source*nbcuts+2]->Draw();
- histopassedcuts[source*nbcuts+3]->Draw("same");
- //histopassedcuts[source*nbcuts+4]->Draw("same");
- histopassedcuts[source*nbcuts+5]->Draw("same");
- histopassedcuts[source*nbcuts+6]->Draw("same");
- //histopassedcuts[source*nbcuts+7]->Draw("same");
- histopassedcuts[source*nbcuts+8]->Draw("same");
- histopassedcuts[source*nbcuts+9]->Draw("same");
- TLegend *legb = new TLegend(0.4,0.6,0.89,0.89);
- //legb->AddEntry(histopassedcuts[source*nbcuts+0],"all","p");
- //legb->AddEntry(histopassedcuts[source*nbcuts+1],"Partner tracked","p");
- legb->AddEntry(histopassedcuts[source*nbcuts+2],"Opposite sign","p");
- legb->AddEntry(histopassedcuts[source*nbcuts+3],"SingleTrackPart","p");
- //legb->AddEntry(histopassedcuts[source*nbcuts+4],"SharedCluster","p");
- legb->AddEntry(histopassedcuts[source*nbcuts+5],"PID","p");
- legb->AddEntry(histopassedcuts[source*nbcuts+6],"DCA","p");
- //legb->AddEntry(histopassedcuts[source*nbcuts+7],"Chi2Ndf","p");
- legb->AddEntry(histopassedcuts[source*nbcuts+8],"OpeningAngle","p");
- legb->AddEntry(histopassedcuts[source*nbcuts+9],"InvMass","p");
- legb->Draw("same");
-
- source = 2;
- chsSparseMCeeff->cd(2);
- histopassedcuts[source*nbcuts+0]->SetTitle("#pi^{0}");
- histopassedcuts[source*nbcuts+1]->SetTitle("#pi^{0}");
- histopassedcuts[source*nbcuts+2]->SetTitle("#pi^{0}");
- histopassedcuts[source*nbcuts+3]->SetTitle("#pi^{0}");
- histopassedcuts[source*nbcuts+4]->SetTitle("#pi^{0}");
- histopassedcuts[source*nbcuts+5]->SetTitle("#pi^{0}");
- histopassedcuts[source*nbcuts+6]->SetTitle("#pi^{0}");
- histopassedcuts[source*nbcuts+7]->SetTitle("#pi^{0}");
- histopassedcuts[source*nbcuts+8]->SetTitle("#pi^{0}");
- histopassedcuts[source*nbcuts+9]->SetTitle("#pi^{0}");
- //histopassedcuts[source*nbcuts+0]->SetStats(0);
- histopassedcuts[source*nbcuts+1]->SetStats(0);
- histopassedcuts[source*nbcuts+2]->SetStats(0);
- histopassedcuts[source*nbcuts+3]->SetStats(0);
- histopassedcuts[source*nbcuts+4]->SetStats(0);
- histopassedcuts[source*nbcuts+5]->SetStats(0);
- histopassedcuts[source*nbcuts+6]->SetStats(0);
- histopassedcuts[source*nbcuts+7]->SetStats(0);
- histopassedcuts[source*nbcuts+8]->SetStats(0);
- histopassedcuts[source*nbcuts+9]->SetStats(0);
- //histopassedcuts[source*nbcuts+0]->Draw();
- //histopassedcuts[source*nbcuts+1]->Draw();
- histopassedcuts[source*nbcuts+2]->Draw();
- histopassedcuts[source*nbcuts+3]->Draw("same");
- //histopassedcuts[source*nbcuts+4]->Draw("same");
- histopassedcuts[source*nbcuts+5]->Draw("same");
- histopassedcuts[source*nbcuts+6]->Draw("same");
- //histopassedcuts[source*nbcuts+7]->Draw("same");
- histopassedcuts[source*nbcuts+8]->Draw("same");
- histopassedcuts[source*nbcuts+9]->Draw("same");
- TLegend *legc = new TLegend(0.4,0.6,0.89,0.89);
- //legc->AddEntry(histopassedcuts[source*nbcuts+0],"all","p");
- //legc->AddEntry(histopassedcuts[source*nbcuts+1],"Partner tracked","p");
- legc->AddEntry(histopassedcuts[source*nbcuts+2],"Opposite sign","p");
- legc->AddEntry(histopassedcuts[source*nbcuts+3],"SingleTrackPart","p");
- //legc->AddEntry(histopassedcuts[source*nbcuts+4],"SharedCluster","p");
- legc->AddEntry(histopassedcuts[source*nbcuts+5],"PID","p");
- legc->AddEntry(histopassedcuts[source*nbcuts+6],"DCA","p");
- //legc->AddEntry(histopassedcuts[source*nbcuts+7],"Chi2Ndf","p");
- legc->AddEntry(histopassedcuts[source*nbcuts+8],"OpeningAngle","p");
- legc->AddEntry(histopassedcuts[source*nbcuts+9],"InvMass","p");
- legc->Draw("same");
-
- source = 4;
- chsSparseMCeeff->cd(3);
- histopassedcuts[source*nbcuts+0]->SetTitle("C");
- histopassedcuts[source*nbcuts+1]->SetTitle("C");
- histopassedcuts[source*nbcuts+2]->SetTitle("C");
- histopassedcuts[source*nbcuts+3]->SetTitle("C");
- histopassedcuts[source*nbcuts+4]->SetTitle("C");
- histopassedcuts[source*nbcuts+5]->SetTitle("C");
- histopassedcuts[source*nbcuts+6]->SetTitle("C");
- histopassedcuts[source*nbcuts+7]->SetTitle("C");
- histopassedcuts[source*nbcuts+8]->SetTitle("C");
- histopassedcuts[source*nbcuts+9]->SetTitle("C");
- //histopassedcuts[source*nbcuts+0]->SetStats(0);
- histopassedcuts[source*nbcuts+1]->SetStats(0);
- histopassedcuts[source*nbcuts+2]->SetStats(0);
- histopassedcuts[source*nbcuts+3]->SetStats(0);
- histopassedcuts[source*nbcuts+4]->SetStats(0);
- histopassedcuts[source*nbcuts+5]->SetStats(0);
- histopassedcuts[source*nbcuts+6]->SetStats(0);
- histopassedcuts[source*nbcuts+7]->SetStats(0);
- histopassedcuts[source*nbcuts+8]->SetStats(0);
- histopassedcuts[source*nbcuts+9]->SetStats(0);
- //histopassedcuts[source*nbcuts+0]->Draw();
- //histopassedcuts[source*nbcuts+1]->Draw();
- histopassedcuts[source*nbcuts+2]->Draw();
- histopassedcuts[source*nbcuts+3]->Draw("same");
- //histopassedcuts[source*nbcuts+4]->Draw("same");
- histopassedcuts[source*nbcuts+5]->Draw("same");
- histopassedcuts[source*nbcuts+6]->Draw("same");
- //histopassedcuts[source*nbcuts+7]->Draw("same");
- histopassedcuts[source*nbcuts+8]->Draw("same");
- histopassedcuts[source*nbcuts+9]->Draw("same");
- TLegend *lege = new TLegend(0.4,0.6,0.89,0.89);
- //lege->AddEntry(histopassedcuts[source*nbcuts+0],"all","p");
- //lege->AddEntry(histopassedcuts[source*nbcuts+1],"Partner tracked","p");
- lege->AddEntry(histopassedcuts[source*nbcuts+2],"Opposite sign","p");
- lege->AddEntry(histopassedcuts[source*nbcuts+3],"SingleTrackPart","p");
- //lege->AddEntry(histopassedcuts[source*nbcuts+4],"SharedCluster","p");
- lege->AddEntry(histopassedcuts[source*nbcuts+5],"PID","p");
- lege->AddEntry(histopassedcuts[source*nbcuts+6],"DCA","p");
- //lege->AddEntry(histopassedcuts[source*nbcuts+7],"Chi2Ndf","p");
- lege->AddEntry(histopassedcuts[source*nbcuts+8],"OpeningAngle","p");
- lege->AddEntry(histopassedcuts[source*nbcuts+9],"InvMass","p");
- lege->Draw("same");
-
- //////////////////////
- // Input
- //////////////////////
-
- TCanvas * chsSparseMCein =new TCanvas("hsSparseMCeinput","hsSparseMCeinput",800,800);
- chsSparseMCein->cd(1);
- Double_t nbGamma = 0.0;
- source = 1;
- nbGamma = histopassedcuts[source*nbcuts+0]->GetEntries();
- histopassedcuts[source*nbcuts+0]->SetStats(0);
- histopassedcuts[source*nbcuts+0]->Draw();
- TLegend *leginput = new TLegend(0.4,0.6,0.89,0.89);
- leginput->AddEntry(histopassedcuts[source*nbcuts+0],"#gamma","p");
- Double_t nbPi0 = 0.0;
- source = 2;
- nbPi0 = histopassedcuts[source*nbcuts+0]->GetEntries();
- histopassedcuts[source*nbcuts+0]->SetStats(0);
- histopassedcuts[source*nbcuts+0]->Draw("same");
- leginput->AddEntry(histopassedcuts[source*nbcuts+0],"#pi^{0}","p");
- Double_t nbEta = 0.0;
- source = 3;
- nbEta = histopassedcuts[source*nbcuts+0]->GetEntries();
- histopassedcuts[source*nbcuts+0]->SetStats(0);
- histopassedcuts[source*nbcuts+0]->Draw("same");
- leginput->AddEntry(histopassedcuts[source*nbcuts+0],"#eta","p");
- Double_t nbC = 0.0;
- source = 4;
- nbC = histopassedcuts[source*nbcuts+0]->GetEntries();
- histopassedcuts[source*nbcuts+0]->SetStats(0);
- histopassedcuts[source*nbcuts+0]->Draw("same");
- leginput->AddEntry(histopassedcuts[source*nbcuts+0],"c","p");
- leginput->Draw("same");
-
- //printf("Gamma %f, pi^{0} %f and #eta %f, c %f\n",nbGamma,nbPi0,nbEta,nbC);
-
- //////////////////////
- // Tracked
- //////////////////////
-
- TCanvas * cTracked = new TCanvas("cTracked","cTracked",800,800);
- cTracked->cd(1);
- source = 1;
- histopassedcuts[source*nbcuts+1]->Draw();
- TLegend *legTracked = new TLegend(0.4,0.6,0.89,0.89);
- legTracked->AddEntry(histopassedcuts[source*nbcuts+1],"#gamma","p");
- source = 2;
- histopassedcuts[source*nbcuts+1]->Draw("same");
- legTracked->AddEntry(histopassedcuts[source*nbcuts+1],"#pi^{0}","p");
- legTracked->Draw("same");
-
- }
-
/////////////////////////////////////
// Data Radius and chi2Ndf if AliKF
////////////////////////////////////
legITS5->Draw("same");
- }
+ }
+ ////////////////////////
+ // Cut efficiencies
+ ////////////////////////
+
+ THnSparseF *hsSparseMCe = dynamic_cast<THnSparseF *>(fList->FindObject("CutPassedMC"));
+ if(!hsSparseMCe) return;
+
+ // init histos
+ TAxis *axissources = hsSparseMCe->GetAxis(2);
+ Int_t nbsources = axissources->GetNbins();
+ TAxis *axiscuts = hsSparseMCe->GetAxis(1);
+ Int_t nbcuts = axiscuts->GetNbins();
+ TH1D **histopassedcuts = new TH1D*[nbsources*nbcuts];
+ Double_t *nbEntriesCuts = new Double_t[nbsources*nbcuts];
+ for(Int_t k =0; k < nbsources*nbcuts; k++){
+ nbEntriesCuts[k] = 0.0;
+ histopassedcuts[k] = 0x0;
+ }
+
+ //printf("Number of cuts %d\n",nbcuts);
+
+ // canvas
+ TCanvas * chsSparseMCeeff =new TCanvas("hsSparseMCeeffDebug","hsSparseMCeeffDebug",800,800);
+ chsSparseMCeeff->Divide(3,1);
+
+ // histos
+ for(Int_t sourceid = 0; sourceid < nbsources; sourceid++) {
+ hsSparseMCe->GetAxis(2)->SetRange(sourceid+1,sourceid+1);
+ for(Int_t cut = 0; cut < nbcuts; cut++){
+ hsSparseMCe->GetAxis(1)->SetRange(cut+1,cut+1);
+ histopassedcuts[sourceid*nbcuts+cut] = hsSparseMCe->Projection(0);
+ hsSparseMCe->GetAxis(1)->SetRange(1,hsSparseMCe->GetAxis(1)->GetNbins());
+ }
+ hsSparseMCe->GetAxis(2)->SetRange(1,hsSparseMCe->GetAxis(2)->GetNbins());
+ }
+
+ // calcul efficiencies
+
+ // histos
+ for(Int_t sourceid = 0; sourceid < nbsources; sourceid++) {
+ // Next is compared to the partner tracked
+ for(Int_t cut = 2; cut < nbcuts; cut++){
+ nbEntriesCuts[sourceid*nbcuts+cut] = histopassedcuts[sourceid*nbcuts+cut]->GetEntries();
+ if(histopassedcuts[sourceid*nbcuts+1]->GetEntries() > 0.0) histopassedcuts[sourceid*nbcuts+cut]->Divide(histopassedcuts[sourceid*nbcuts+1]);
+ }
+ // First one is if the partner is tracked.
+ nbEntriesCuts[sourceid*nbcuts+1] = histopassedcuts[sourceid*nbcuts+1]->GetEntries();
+ if(histopassedcuts[sourceid*nbcuts]->GetEntries() > 0.0) histopassedcuts[sourceid*nbcuts+1]->Divide(histopassedcuts[sourceid*nbcuts]);
+ // First one is input
+ nbEntriesCuts[sourceid*nbcuts] = histopassedcuts[sourceid*nbcuts]->GetEntries();
+ }
+
+ // ratios
+ for(Int_t sourceid = 0; sourceid < nbsources; sourceid++) {
+ for(Int_t cut = 1; cut < nbcuts; cut++){
+ if(nbEntriesCuts[sourceid*nbcuts] > 0.0) nbEntriesCuts[sourceid*nbcuts+cut] = nbEntriesCuts[sourceid*nbcuts+cut]/nbEntriesCuts[sourceid*nbcuts];
+ }
+ }
+ TH1F *ratioHistoEntriesGamma = new TH1F("ratioHistoEntriesGamma","", nbcuts-1, 0.0, nbcuts-1.0);
+ TH1F *ratioHistoEntriesPi0 = new TH1F("ratioHistoEntriesPi0","", nbcuts-1, 0.0, nbcuts-1.0);
+ TH1F *ratioHistoEntriesC = new TH1F("ratioHistoEntriesC","", nbcuts-1, 0.0, nbcuts-1.0);
+ for(Int_t k = 1; k < nbcuts; k++){
+ if((nbcuts+k) < (nbsources*nbcuts)) ratioHistoEntriesGamma->SetBinContent(k,nbEntriesCuts[nbcuts+k]);
+ if((2*nbcuts+k) < (nbsources*nbcuts)) ratioHistoEntriesPi0->SetBinContent(k,nbEntriesCuts[2*nbcuts+k]);
+ if((4*nbcuts+k) < (nbsources*nbcuts)) ratioHistoEntriesC->SetBinContent(k,nbEntriesCuts[4*nbcuts+k]);
+ }
+ //
+ TAxis *xAxisGamma = ratioHistoEntriesGamma->GetXaxis();
+ xAxisGamma->SetBinLabel(1,"Partner tracked");
+ xAxisGamma->SetBinLabel(2,"Opposite sign");
+ xAxisGamma->SetBinLabel(3,"Single Track Cut");
+ xAxisGamma->SetBinLabel(4,"Shared Clusters");
+ xAxisGamma->SetBinLabel(5,"PID");
+ xAxisGamma->SetBinLabel(6,"DCA");
+ xAxisGamma->SetBinLabel(7,"Chi^{2}/Ndf");
+ xAxisGamma->SetBinLabel(8,"Opening angle");
+ xAxisGamma->SetBinLabel(9,"Invariant mass");
+ //
+ TAxis *xAxisPi0 = ratioHistoEntriesPi0->GetXaxis();
+ xAxisPi0->SetBinLabel(1,"Partner tracked");
+ xAxisPi0->SetBinLabel(2,"Opposite sign");
+ xAxisPi0->SetBinLabel(3,"Single Track Cut");
+ xAxisPi0->SetBinLabel(4,"Shared Clusters");
+ xAxisPi0->SetBinLabel(5,"PID");
+ xAxisPi0->SetBinLabel(6,"DCA");
+ xAxisPi0->SetBinLabel(7,"Chi^{2}/Ndf");
+ xAxisPi0->SetBinLabel(8,"Opening angle");
+ xAxisPi0->SetBinLabel(9,"Invariant mass");
+ //
+ TAxis *xAxisC = ratioHistoEntriesC->GetXaxis();
+ xAxisC->SetBinLabel(1,"Partner tracked");
+ xAxisC->SetBinLabel(2,"Opposite sign");
+ xAxisC->SetBinLabel(3,"Single Track Cut");
+ xAxisC->SetBinLabel(4,"Shared Clusters");
+ xAxisC->SetBinLabel(5,"PID");
+ xAxisC->SetBinLabel(6,"DCA");
+ xAxisC->SetBinLabel(7,"Chi^{2}/Ndf");
+ xAxisC->SetBinLabel(8,"Opening angle");
+ xAxisC->SetBinLabel(9,"Invariant mass");
+ //
+ TCanvas * cRatioHistoEntries =new TCanvas("cRatioHistoEntries","cRatioHistoEntries",800,800);
+ cRatioHistoEntries->cd(1);
+ ratioHistoEntriesGamma->SetStats(0);
+ ratioHistoEntriesGamma->Draw();
+ ratioHistoEntriesPi0->SetStats(0);
+ ratioHistoEntriesPi0->Draw("same");
+ ratioHistoEntriesC->SetStats(0);
+ //ratioHistoEntriesC->Draw("same");
+ TLegend *legEntries = new TLegend(0.4,0.6,0.89,0.89);
+ legEntries->AddEntry(ratioHistoEntriesGamma,"#gamma","l");
+ legEntries->AddEntry(ratioHistoEntriesPi0,"#pi^{0}","l");
+ //legEntries->AddEntry(ratioHistoEntriesC,"c","p");
+ legEntries->Draw("same");
+
+ // plot Debug
+ Int_t source = 1;
+ chsSparseMCeeff->cd(1);
+ if(((source*nbcuts+0)> (nbsources*nbcuts-1)) || ((source*nbcuts+1)> (nbsources*nbcuts-1)) || ((source*nbcuts+2)> (nbsources*nbcuts-1)) || ((source*nbcuts+3)> (nbsources*nbcuts-1)) || ((source*nbcuts+4)> (nbsources*nbcuts-1)) || ((source*nbcuts+5)> (nbsources*nbcuts-1)) || ((source*nbcuts+6)> (nbsources*nbcuts-1)) || ((source*nbcuts+7)> (nbsources*nbcuts-1)) || ((source*nbcuts+8)> (nbsources*nbcuts-1)) || ((source*nbcuts+9)> (nbsources*nbcuts-1))) return;
+ if((!histopassedcuts[source*nbcuts+0]) || (!histopassedcuts[source*nbcuts+1]) || (!histopassedcuts[source*nbcuts+2]) || (!histopassedcuts[source*nbcuts+3]) || (!histopassedcuts[source*nbcuts+4]) || (!histopassedcuts[source*nbcuts+5]) || (!histopassedcuts[source*nbcuts+6]) || (!histopassedcuts[source*nbcuts+7]) || (!histopassedcuts[source*nbcuts+8]) || (!histopassedcuts[source*nbcuts+9])) return;
+ histopassedcuts[source*nbcuts+0]->SetTitle("#gamma");
+ histopassedcuts[source*nbcuts+1]->SetTitle("#gamma");
+ histopassedcuts[source*nbcuts+2]->SetTitle("#gamma");
+ histopassedcuts[source*nbcuts+3]->SetTitle("#gamma");
+ histopassedcuts[source*nbcuts+4]->SetTitle("#gamma");
+ histopassedcuts[source*nbcuts+5]->SetTitle("#gamma");
+ histopassedcuts[source*nbcuts+6]->SetTitle("#gamma");
+ histopassedcuts[source*nbcuts+7]->SetTitle("#gamma");
+ histopassedcuts[source*nbcuts+8]->SetTitle("#gamma");
+ histopassedcuts[source*nbcuts+9]->SetTitle("#gamma");
+ //histopassedcuts[source*nbcuts+0]->SetStats(0);
+ histopassedcuts[source*nbcuts+1]->SetStats(0);
+ histopassedcuts[source*nbcuts+2]->SetStats(0);
+ histopassedcuts[source*nbcuts+3]->SetStats(0);
+ histopassedcuts[source*nbcuts+4]->SetStats(0);
+ histopassedcuts[source*nbcuts+5]->SetStats(0);
+ histopassedcuts[source*nbcuts+6]->SetStats(0);
+ histopassedcuts[source*nbcuts+7]->SetStats(0);
+ histopassedcuts[source*nbcuts+8]->SetStats(0);
+ histopassedcuts[source*nbcuts+9]->SetStats(0);
+ //histopassedcuts[source*nbcuts+0]->Draw();
+ //histopassedcuts[source*nbcuts+1]->Draw("");
+ histopassedcuts[source*nbcuts+2]->Draw();
+ histopassedcuts[source*nbcuts+3]->Draw("same");
+ //histopassedcuts[source*nbcuts+4]->Draw("same");
+ histopassedcuts[source*nbcuts+5]->Draw("same");
+ histopassedcuts[source*nbcuts+6]->Draw("same");
+ //histopassedcuts[source*nbcuts+7]->Draw("same");
+ histopassedcuts[source*nbcuts+8]->Draw("same");
+ histopassedcuts[source*nbcuts+9]->Draw("same");
+ TLegend *legb = new TLegend(0.4,0.6,0.89,0.89);
+ //legb->AddEntry(histopassedcuts[source*nbcuts+0],"all","p");
+ //legb->AddEntry(histopassedcuts[source*nbcuts+1],"Partner tracked","p");
+ legb->AddEntry(histopassedcuts[source*nbcuts+2],"Opposite sign","p");
+ legb->AddEntry(histopassedcuts[source*nbcuts+3],"SingleTrackPart","p");
+ //legb->AddEntry(histopassedcuts[source*nbcuts+4],"SharedCluster","p");
+ legb->AddEntry(histopassedcuts[source*nbcuts+5],"PID","p");
+ legb->AddEntry(histopassedcuts[source*nbcuts+6],"DCA","p");
+ //legb->AddEntry(histopassedcuts[source*nbcuts+7],"Chi2Ndf","p");
+ legb->AddEntry(histopassedcuts[source*nbcuts+8],"OpeningAngle","p");
+ legb->AddEntry(histopassedcuts[source*nbcuts+9],"InvMass","p");
+ legb->Draw("same");
+
+ source = 2;
+ if(((source*nbcuts+0)> (nbsources*nbcuts-1)) || ((source*nbcuts+1)> (nbsources*nbcuts-1)) || ((source*nbcuts+2)> (nbsources*nbcuts-1)) || ((source*nbcuts+3)> (nbsources*nbcuts-1)) || ((source*nbcuts+4)> (nbsources*nbcuts-1)) || ((source*nbcuts+5)> (nbsources*nbcuts-1)) || ((source*nbcuts+6)> (nbsources*nbcuts-1)) || ((source*nbcuts+7)> (nbsources*nbcuts-1)) || ((source*nbcuts+8)> (nbsources*nbcuts-1)) || ((source*nbcuts+9)> (nbsources*nbcuts-1))) return;
+ if((!histopassedcuts[source*nbcuts+0]) || (!histopassedcuts[source*nbcuts+1]) || (!histopassedcuts[source*nbcuts+2]) || (!histopassedcuts[source*nbcuts+3]) || (!histopassedcuts[source*nbcuts+4]) || (!histopassedcuts[source*nbcuts+5]) || (!histopassedcuts[source*nbcuts+6]) || (!histopassedcuts[source*nbcuts+7]) || (!histopassedcuts[source*nbcuts+8]) || (!histopassedcuts[source*nbcuts+9])) return;
+ chsSparseMCeeff->cd(2);
+ histopassedcuts[source*nbcuts+0]->SetTitle("#pi^{0}");
+ histopassedcuts[source*nbcuts+1]->SetTitle("#pi^{0}");
+ histopassedcuts[source*nbcuts+2]->SetTitle("#pi^{0}");
+ histopassedcuts[source*nbcuts+3]->SetTitle("#pi^{0}");
+ histopassedcuts[source*nbcuts+4]->SetTitle("#pi^{0}");
+ histopassedcuts[source*nbcuts+5]->SetTitle("#pi^{0}");
+ histopassedcuts[source*nbcuts+6]->SetTitle("#pi^{0}");
+ histopassedcuts[source*nbcuts+7]->SetTitle("#pi^{0}");
+ histopassedcuts[source*nbcuts+8]->SetTitle("#pi^{0}");
+ histopassedcuts[source*nbcuts+9]->SetTitle("#pi^{0}");
+ //histopassedcuts[source*nbcuts+0]->SetStats(0);
+ histopassedcuts[source*nbcuts+1]->SetStats(0);
+ histopassedcuts[source*nbcuts+2]->SetStats(0);
+ histopassedcuts[source*nbcuts+3]->SetStats(0);
+ histopassedcuts[source*nbcuts+4]->SetStats(0);
+ histopassedcuts[source*nbcuts+5]->SetStats(0);
+ histopassedcuts[source*nbcuts+6]->SetStats(0);
+ histopassedcuts[source*nbcuts+7]->SetStats(0);
+ histopassedcuts[source*nbcuts+8]->SetStats(0);
+ histopassedcuts[source*nbcuts+9]->SetStats(0);
+ //histopassedcuts[source*nbcuts+0]->Draw();
+ //histopassedcuts[source*nbcuts+1]->Draw();
+ histopassedcuts[source*nbcuts+2]->Draw();
+ histopassedcuts[source*nbcuts+3]->Draw("same");
+ //histopassedcuts[source*nbcuts+4]->Draw("same");
+ histopassedcuts[source*nbcuts+5]->Draw("same");
+ histopassedcuts[source*nbcuts+6]->Draw("same");
+ //histopassedcuts[source*nbcuts+7]->Draw("same");
+ histopassedcuts[source*nbcuts+8]->Draw("same");
+ histopassedcuts[source*nbcuts+9]->Draw("same");
+ TLegend *legc = new TLegend(0.4,0.6,0.89,0.89);
+ //legc->AddEntry(histopassedcuts[source*nbcuts+0],"all","p");
+ //legc->AddEntry(histopassedcuts[source*nbcuts+1],"Partner tracked","p");
+ legc->AddEntry(histopassedcuts[source*nbcuts+2],"Opposite sign","p");
+ legc->AddEntry(histopassedcuts[source*nbcuts+3],"SingleTrackPart","p");
+ //legc->AddEntry(histopassedcuts[source*nbcuts+4],"SharedCluster","p");
+ legc->AddEntry(histopassedcuts[source*nbcuts+5],"PID","p");
+ legc->AddEntry(histopassedcuts[source*nbcuts+6],"DCA","p");
+ //legc->AddEntry(histopassedcuts[source*nbcuts+7],"Chi2Ndf","p");
+ legc->AddEntry(histopassedcuts[source*nbcuts+8],"OpeningAngle","p");
+ legc->AddEntry(histopassedcuts[source*nbcuts+9],"InvMass","p");
+ legc->Draw("same");
+
+ source = 4;
+ if(((source*nbcuts+0)> (nbsources*nbcuts-1)) || ((source*nbcuts+1)> (nbsources*nbcuts-1)) || ((source*nbcuts+2)> (nbsources*nbcuts-1)) || ((source*nbcuts+3)> (nbsources*nbcuts-1)) || ((source*nbcuts+4)> (nbsources*nbcuts-1)) || ((source*nbcuts+5)> (nbsources*nbcuts-1)) || ((source*nbcuts+6)> (nbsources*nbcuts-1)) || ((source*nbcuts+7)> (nbsources*nbcuts-1)) || ((source*nbcuts+8)> (nbsources*nbcuts-1)) || ((source*nbcuts+9)> (nbsources*nbcuts-1))) return;
+ if((!histopassedcuts[source*nbcuts+0]) || (!histopassedcuts[source*nbcuts+1]) || (!histopassedcuts[source*nbcuts+2]) || (!histopassedcuts[source*nbcuts+3]) || (!histopassedcuts[source*nbcuts+4]) || (!histopassedcuts[source*nbcuts+5]) || (!histopassedcuts[source*nbcuts+6]) || (!histopassedcuts[source*nbcuts+7]) || (!histopassedcuts[source*nbcuts+8]) || (!histopassedcuts[source*nbcuts+9])) return;
+ chsSparseMCeeff->cd(3);
+ histopassedcuts[source*nbcuts+0]->SetTitle("C");
+ histopassedcuts[source*nbcuts+1]->SetTitle("C");
+ histopassedcuts[source*nbcuts+2]->SetTitle("C");
+ histopassedcuts[source*nbcuts+3]->SetTitle("C");
+ histopassedcuts[source*nbcuts+4]->SetTitle("C");
+ histopassedcuts[source*nbcuts+5]->SetTitle("C");
+ histopassedcuts[source*nbcuts+6]->SetTitle("C");
+ histopassedcuts[source*nbcuts+7]->SetTitle("C");
+ histopassedcuts[source*nbcuts+8]->SetTitle("C");
+ histopassedcuts[source*nbcuts+9]->SetTitle("C");
+ //histopassedcuts[source*nbcuts+0]->SetStats(0);
+ histopassedcuts[source*nbcuts+1]->SetStats(0);
+ histopassedcuts[source*nbcuts+2]->SetStats(0);
+ histopassedcuts[source*nbcuts+3]->SetStats(0);
+ histopassedcuts[source*nbcuts+4]->SetStats(0);
+ histopassedcuts[source*nbcuts+5]->SetStats(0);
+ histopassedcuts[source*nbcuts+6]->SetStats(0);
+ histopassedcuts[source*nbcuts+7]->SetStats(0);
+ histopassedcuts[source*nbcuts+8]->SetStats(0);
+ histopassedcuts[source*nbcuts+9]->SetStats(0);
+ //histopassedcuts[source*nbcuts+0]->Draw();
+ //histopassedcuts[source*nbcuts+1]->Draw();
+ histopassedcuts[source*nbcuts+2]->Draw();
+ histopassedcuts[source*nbcuts+3]->Draw("same");
+ //histopassedcuts[source*nbcuts+4]->Draw("same");
+ histopassedcuts[source*nbcuts+5]->Draw("same");
+ histopassedcuts[source*nbcuts+6]->Draw("same");
+ //histopassedcuts[source*nbcuts+7]->Draw("same");
+ histopassedcuts[source*nbcuts+8]->Draw("same");
+ histopassedcuts[source*nbcuts+9]->Draw("same");
+ TLegend *lege = new TLegend(0.4,0.6,0.89,0.89);
+ //lege->AddEntry(histopassedcuts[source*nbcuts+0],"all","p");
+ //lege->AddEntry(histopassedcuts[source*nbcuts+1],"Partner tracked","p");
+ lege->AddEntry(histopassedcuts[source*nbcuts+2],"Opposite sign","p");
+ lege->AddEntry(histopassedcuts[source*nbcuts+3],"SingleTrackPart","p");
+ //lege->AddEntry(histopassedcuts[source*nbcuts+4],"SharedCluster","p");
+ lege->AddEntry(histopassedcuts[source*nbcuts+5],"PID","p");
+ lege->AddEntry(histopassedcuts[source*nbcuts+6],"DCA","p");
+ //lege->AddEntry(histopassedcuts[source*nbcuts+7],"Chi2Ndf","p");
+ lege->AddEntry(histopassedcuts[source*nbcuts+8],"OpeningAngle","p");
+ lege->AddEntry(histopassedcuts[source*nbcuts+9],"InvMass","p");
+ lege->Draw("same");
+
+ //////////////////////
+ // Input
+ //////////////////////
+
+ TCanvas * chsSparseMCein =new TCanvas("hsSparseMCeinput","hsSparseMCeinput",800,800);
+ chsSparseMCein->cd(1);
+ Double_t nbGamma = 0.0;
+ source = 1;
+ if(((source*nbcuts+0)> (nbsources*nbcuts-1)) || (!histopassedcuts[source*nbcuts+0])) return;
+ nbGamma = histopassedcuts[source*nbcuts+0]->GetEntries();
+ histopassedcuts[source*nbcuts+0]->SetStats(0);
+ histopassedcuts[source*nbcuts+0]->Draw();
+ TLegend *leginput = new TLegend(0.4,0.6,0.89,0.89);
+ leginput->AddEntry(histopassedcuts[source*nbcuts+0],"#gamma","p");
+ Double_t nbPi0 = 0.0;
+ source = 2;
+ if(((source*nbcuts+0)> (nbsources*nbcuts-1)) || (!histopassedcuts[source*nbcuts+0])) return;
+ nbPi0 = histopassedcuts[source*nbcuts+0]->GetEntries();
+ histopassedcuts[source*nbcuts+0]->SetStats(0);
+ histopassedcuts[source*nbcuts+0]->Draw("same");
+ leginput->AddEntry(histopassedcuts[source*nbcuts+0],"#pi^{0}","p");
+ Double_t nbEta = 0.0;
+ source = 3;
+ if(((source*nbcuts+0)> (nbsources*nbcuts-1)) || (!histopassedcuts[source*nbcuts+0])) return;
+ nbEta = histopassedcuts[source*nbcuts+0]->GetEntries();
+ histopassedcuts[source*nbcuts+0]->SetStats(0);
+ histopassedcuts[source*nbcuts+0]->Draw("same");
+ leginput->AddEntry(histopassedcuts[source*nbcuts+0],"#eta","p");
+ Double_t nbC = 0.0;
+ source = 4;
+ if(((source*nbcuts+0)> (nbsources*nbcuts-1)) || (!histopassedcuts[source*nbcuts+0])) return;
+ nbC = histopassedcuts[source*nbcuts+0]->GetEntries();
+ histopassedcuts[source*nbcuts+0]->SetStats(0);
+ histopassedcuts[source*nbcuts+0]->Draw("same");
+ leginput->AddEntry(histopassedcuts[source*nbcuts+0],"c","p");
+ leginput->Draw("same");
+
+ //printf("Gamma %f, pi^{0} %f and #eta %f, c %f\n",nbGamma,nbPi0,nbEta,nbC);
+
+ //////////////////////
+ // Tracked
+ //////////////////////
+
+ TCanvas * cTracked = new TCanvas("cTracked","cTracked",800,800);
+ cTracked->cd(1);
+ source = 1;
+ if(((source*nbcuts+1)> (nbsources*nbcuts-1)) || (!histopassedcuts[source*nbcuts+1])) return;
+ histopassedcuts[source*nbcuts+1]->Draw();
+ TLegend *legTracked = new TLegend(0.4,0.6,0.89,0.89);
+ legTracked->AddEntry(histopassedcuts[source*nbcuts+1],"#gamma","p");
+ source = 2;
+ if(((source*nbcuts+1)> (nbsources*nbcuts-1)) || (!histopassedcuts[source*nbcuts+1])) return;
+ histopassedcuts[source*nbcuts+1]->Draw("same");
+ legTracked->AddEntry(histopassedcuts[source*nbcuts+1],"#pi^{0}","p");
+ legTracked->Draw("same");
+
+
}
//_____________________________________________________________________________
Double_t AliHFEelecbackground::BetheBlochElectronITS(const Double_t *x, const Double_t * /*par*/)
fEvent(NULL),
fCutCorrelation(0),
fRequirements(0),
- fTPCiter1(kFALSE),
fMinNClustersTPC(0),
fClusterRatioTPC(0.),
fMinTrackletsTRD(0),
fPixelITS(0),
+ fTOFpid(kFALSE),
+ fTPCclusterDef(0),
+ fTPCclusterRatioDef(0),
fCheck(kFALSE),
fQAlist(0x0) ,
fDebugLevel(0)
fEvent(c.fEvent),
fCutCorrelation(c.fCutCorrelation),
fRequirements(c.fRequirements),
- fTPCiter1(c.fTPCiter1),
fMinNClustersTPC(c.fMinNClustersTPC),
fClusterRatioTPC(c.fClusterRatioTPC),
fMinTrackletsTRD(c.fMinTrackletsTRD),
fPixelITS(c.fPixelITS),
+ fTOFpid(c.fTOFpid),
+ fTPCclusterDef(c.fTPCclusterDef),
+ fTPCclusterRatioDef(c.fTPCclusterRatioDef),
fCheck(c.fCheck),
fQAlist(0x0),
fDebugLevel(0)
fEvent = c.fEvent;
fCutCorrelation = c.fCutCorrelation;
fRequirements = c.fRequirements;
- fTPCiter1 = c.fTPCiter1;
fClusterRatioTPC = c.fClusterRatioTPC;
fMinNClustersTPC = c.fMinNClustersTPC;
fMinTrackletsTRD = c.fMinTrackletsTRD;
fPixelITS = c.fPixelITS;
+ fTPCclusterDef = c.fTPCclusterDef;
+ fTPCclusterRatioDef = c.fTPCclusterRatioDef;
+ fTOFpid = c.fTOFpid;
fCheck = c.fCheck;
fDebugLevel = c.fDebugLevel;
memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
ULong64_t survivedCut = 0; // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
if(IsQAOn()) FillQAhistosRec(track, kBeforeCuts);
// Apply cuts
- Float_t impactR, impactZ, ratioTPC;
+ Float_t impactR, impactZ;
Double_t hfeimpactR, hfeimpactnsigmaR;
Double_t hfeimpactRcut, hfeimpactnsigmaRcut;
+ Bool_t tofstep = kTRUE;
GetImpactParameters(track, impactR, impactZ);
if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
// Protection for PbPb
GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
}
- UInt_t nTPCf = GetTPCfindableClusters(track, fTPCiter1), nclsTPC = GetTPCncls(track, fTPCiter1);
+ UInt_t nclsTPC = GetTPCncls(track);
// printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
- ratioTPC = nTPCf > 0. ? static_cast<Float_t>(nclsTPC)/static_cast<Float_t>(nTPCf) : 1.;
+ Double_t ratioTPC = GetTPCclusterRatio(track);
UChar_t trdTracklets;
trdTracklets = GetTRDnTrackletsPID(track);
UChar_t itsPixel = track->GetITSClusterMap();
}
AliDebug(1, Form("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
}
- if(fRequirements == survivedCut){
+ if(fTOFpid){
+ // cut on TOF pid
+ tofstep = kFALSE;
+ if(track->GetStatus() & AliESDtrack::kTOFpid) tofstep = kTRUE;
+
+ }
+ if(fRequirements == survivedCut && tofstep){
//
// Track selected
//
const Int_t kNhistos = 6;
Float_t impactR, impactZ;
GetImpactParameters(track, impactR, impactZ);
- Int_t nTPCf = GetTPCfindableClusters(track, fTPCiter1), nclsTPC = GetTPCncls(track, fTPCiter1);
TH1 *htmp = NULL;
if((htmp = dynamic_cast<TH1F *>(fQAlist->At(0 + when * kNhistos)))) htmp->Fill(impactR);
if((htmp = dynamic_cast<TH1F *>(fQAlist->At(1 + when * kNhistos)))) htmp->Fill(impactZ);
// printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
- if((htmp = dynamic_cast<TH1F *>(fQAlist->At(2 + when * kNhistos)))) htmp->Fill(nTPCf > 0. ? static_cast<Float_t>(nclsTPC)/static_cast<Float_t>(nTPCf) : 1.);
+ if((htmp = dynamic_cast<TH1F *>(fQAlist->At(2 + when * kNhistos)))) htmp->Fill(GetTPCclusterRatio(track));
if((htmp = dynamic_cast<TH1F *>(fQAlist->At(3 + when * kNhistos)))) htmp->Fill(GetTRDnTrackletsPID(track));
- if((htmp = dynamic_cast<TH1F *>(fQAlist->At(4 + when * kNhistos)))) htmp->Fill(nclsTPC);
+ if((htmp = dynamic_cast<TH1F *>(fQAlist->At(4 + when * kNhistos)))) htmp->Fill(GetTPCncls(track));
UChar_t itsPixel = track->GetITSClusterMap();
TH1 *pixelHist = dynamic_cast<TH1F *>(fQAlist->At(5 + when * kNhistos));
//Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
fQAlist->AddAt(histo1D, 1 + icond * kNhistos);
histo1D->GetXaxis()->SetTitle("Impact Parameter");
histo1D->GetYaxis()->SetTitle("Number of Tracks");
- qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 10, 0, 1)), 2 + icond * kNhistos);
+ qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1)), 2 + icond * kNhistos);
fQAlist->AddAt(histo1D, 2 + icond * kNhistos);
histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
histo1D->GetYaxis()->SetTitle("Number of Tracks");
}
//______________________________________________________
-Int_t AliHFEextraCuts::GetTPCfindableClusters(AliVTrack *track, Bool_t iter1){
- //
- // Get Number of findable clusters in the TPC
- //
- AliDebug(1, Form("Using TPC clusters from iteration 1: %s", iter1 ? "Yes" : "No"));
- Int_t nClusters = 159; // in case no Information available consider all clusters findable
- if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
- AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
- if(esdtrack) nClusters = esdtrack->GetTPCNclsF();
- }
- return nClusters;
-}
-
-//______________________________________________________
-Int_t AliHFEextraCuts::GetTPCncls(AliVTrack *track, Bool_t iter1){
+UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
//
// Get Number of findable clusters in the TPC
//
- AliDebug(1, Form("Using TPC clusters from iteration 1: %s", iter1 ? "Yes" : "No"));
Int_t nClusters = 0; // in case no Information available consider all clusters findable
TString type = track->IsA()->GetName();
if(!type.CompareTo("AliESDtrack")){
AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
if(esdtrack){ // coverity
- if(iter1)
+ if(TESTBIT(fTPCclusterDef, kFoundIter1)){
nClusters = esdtrack->GetTPCNclsIter1();
- else
+ AliDebug(2, ("Using def kFoundIter1"));
+ } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
+ AliDebug(2, ("Using def kCrossedRows"));
+ nClusters = static_cast<UInt_t>(esdtrack->GetTPCClusterInfo(2,1));
+ } else{
+ AliDebug(2, ("Using def kFound"));
nClusters = esdtrack->GetTPCNcls();
+ }
}
}
else if(!type.CompareTo("AliAODTrack")){
AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
- const TBits &tpcmap = aodtrack->GetTPCClusterMap();
- for(UInt_t ibit = 0; ibit < tpcmap.GetNbits(); ibit++)
- if(tpcmap.TestBitNumber(ibit)) nClusters++;
+ if(aodtrack){
+ const TBits &tpcmap = aodtrack->GetTPCClusterMap();
+ for(UInt_t ibit = 0; ibit < tpcmap.GetNbits(); ibit++)
+ if(tpcmap.TestBitNumber(ibit)) nClusters++;
+ }
}
return nClusters;
}
+//______________________________________________________
+Double_t AliHFEextraCuts::GetTPCclusterRatio(AliVTrack *track){
+ //
+ // Get Ratio of found / findable clusters for different definitions
+ // Only implemented for ESD tracks
+ //
+ Double_t clusterRatio = 1.; // in case no Information available consider all clusters findable
+ TString type = track->IsA()->GetName();
+ if(!type.CompareTo("AliESDtrack")){
+ AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
+ if(esdtrack){ // coverity
+ if(TESTBIT(fTPCclusterRatioDef, kCROverFindable)){
+ AliDebug(2, "Using ratio def kCROverFindable");
+ clusterRatio = esdtrack->GetTPCNclsF() ? esdtrack->GetTPCClusterInfo(2,1)/static_cast<Double_t>(esdtrack->GetTPCNclsF()) : 1.; // crossed rows/findable
+ } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverCR)){
+ AliDebug(2, "Using ratio def kFoundOverCR");
+ clusterRatio = esdtrack->GetTPCClusterInfo(2,0); // found/crossed rows
+ } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindableIter1)){
+ AliDebug(2, "Using ratio def kFoundOverFindableIter1");
+ clusterRatio = esdtrack->GetTPCNclsFIter1() ? static_cast<Double_t>(esdtrack->GetTPCNclsIter1())/static_cast<Double_t>(esdtrack->GetTPCNclsFIter1()) : 1.; // crossed
+ } else {
+ AliDebug(2, "Using ratio def kFoundOverFindable");
+ clusterRatio = esdtrack->GetTPCNclsF() ? static_cast<Double_t>(esdtrack->GetTPCNcls())/static_cast<Double_t>(esdtrack->GetTPCNclsF()) : 1.; // found/findable
+ }
+ }
+ }
+ return clusterRatio;
+}
+
//______________________________________________________
void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
//
// Case AOD track: take different constructor
esdtrack = new AliESDtrack(track);
}
- if(esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dca, cov)){
+ if(esdtrack && esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dca, cov)){
// protection
dcaxy = dca[0];
if(cov[0]) dcansigmaxy = dcaxy/TMath::Sqrt(cov[0]);
kNone = 3,
kAny = 4
} ITSPixel_t;
+ typedef enum{
+ kFound = 0,
+ kFoundIter1 = 1,
+ kCrossedRows = 2
+ } ETPCclusterDef_t;
+ typedef enum{
+ kFoundOverFindable = 0,
+ kFoundOverFindableIter1 = 1,
+ kFoundOverCR = 2,
+ kCROverFindable = 3
+ } ETPCclrDef_t;
AliHFEextraCuts(const Char_t *name, const Char_t *title);
AliHFEextraCuts(const AliHFEextraCuts &c);
AliHFEextraCuts &operator=(const AliHFEextraCuts &c);
virtual Bool_t IsSelected(TList *) { return kTRUE; };
virtual void SetRecEventInfo(const TObject *event);
- inline void SetClusterRatioTPC(Double_t ratio);
+ inline void SetClusterRatioTPC(Double_t ratio, ETPCclrDef_t def);
inline void SetRequireITSpixel(ITSPixel_t pixel);
inline void SetMinImpactParamR(Double_t impactParam);
inline void SetMaxImpactParamR(Double_t impactParam);
inline void SetMinHFEImpactParamR();
inline void SetMinHFEImpactParamNsigmaR();
inline void SetMinTrackletsTRD(Int_t minTracklets);
- inline void SetMinNClustersTPC(Int_t minclusters);
- void SetTPCIter1(Bool_t tpcIter1) { fTPCiter1 = tpcIter1; }
+ inline void SetMinNClustersTPC(Int_t minclusters, ETPCclusterDef_t def);
+ void SetTOFPID(Bool_t tofPid) { fTOFpid = tofPid;}
void SetCheckITSstatus(Bool_t check) { fCheck = check; };
Bool_t GetCheckITSstatus() const { return fCheck; };
// Getter Functions for ESD/AOD compatible mode
Int_t GetTRDnTrackletsPID(AliVTrack *track);
Int_t GetITSstatus(AliVTrack *track, Int_t layer);
- Int_t GetTPCfindableClusters(AliVTrack *track, Bool_t iter1 = kFALSE);
- Int_t GetTPCncls(AliVTrack *track, Bool_t iter1 = kFALSE);
+ UInt_t GetTPCncls(AliVTrack *track);
+ Double_t GetTPCclusterRatio(AliVTrack *track);
void GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z);
void GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy);
void GetHFEImpactParameterCuts(AliVTrack *track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut);
kBeforeCuts =0,
kAfterCuts = 1
};
- AliVEvent *fEvent; //! working event
- ULong64_t fCutCorrelation; // Cut Correlation
- ULong64_t fRequirements; // Cut Requirements
- Bool_t fTPCiter1; // Tracking iteration from which the number of clusters is taken
- Float_t fImpactParamCut[4]; // Impact Parmameter Cut
- UInt_t fMinNClustersTPC; // Minimum TPC clusters cut
- Float_t fClusterRatioTPC; // Ratio of findable vs. found clusters in TPC
- UChar_t fMinTrackletsTRD; // Min. Number of Tracklets inside TRD
- UChar_t fPixelITS; // Cut on ITS Pixels
+ AliVEvent *fEvent; //! working event
+ ULong64_t fCutCorrelation; // Cut Correlation
+ ULong64_t fRequirements; // Cut Requirements
+ Float_t fImpactParamCut[4]; // Impact Parmameter Cut
+ UInt_t fMinNClustersTPC; // Minimum TPC clusters cut
+ Float_t fClusterRatioTPC; // Ratio of findable vs. found clusters in TPC
+ UChar_t fMinTrackletsTRD; // Min. Number of Tracklets inside TRD
+ UChar_t fPixelITS; // Cut on ITS Pixels
+ Bool_t fTOFpid; // TOF pid
+ UChar_t fTPCclusterDef; // TPC cluster definition Bitmap
+ UChar_t fTPCclusterRatioDef; // TPC cluster ratio definition Bitmap
Bool_t fCheck; // check
TList *fQAlist; //! Directory for QA histograms
};
//__________________________________________________________
-void AliHFEextraCuts::SetClusterRatioTPC(Double_t ratio) {
+void AliHFEextraCuts::SetClusterRatioTPC(Double_t ratio, ETPCclrDef_t def) {
SETBIT(fRequirements, kClusterRatioTPC);
+ SETBIT(fTPCclusterRatioDef, def);
fClusterRatioTPC = ratio;
}
}
//__________________________________________________________
-void AliHFEextraCuts::SetMinNClustersTPC(Int_t minClusters){
+void AliHFEextraCuts::SetMinNClustersTPC(Int_t minClusters, ETPCclusterDef_t tpcdef){
SETBIT(fRequirements, kMinNClustersTPC);
+ SETBIT(fTPCclusterDef, tpcdef);
fMinNClustersTPC = minClusters;
}
#endif
Double_t* binEdges[1];
binEdges[0] = AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
-
TString hname;
if(kquark == kOthers){
for (Int_t iqType = 0; iqType < 4; iqType++ ){
hname = hnopt+"Pt_"+kqTypeLabel[iqType];
- fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
- //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",500,0,50);
+ //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]); //mj to compare with FONLL
+ fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25);
hname = hnopt+"Y_"+kqTypeLabel[iqType];
fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
hname = hnopt+"Eta_"+kqTypeLabel[iqType];
hname = hnopt+"PdgCode_"+kqTypeLabel[iqType];
fHist[iq][iqType][icut].fPdgCode = new TH1F(hname,hname,20001,-10000.5,10000.5);
hname = hnopt+"Pt_"+kqTypeLabel[iqType];
- fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
- //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",500,0,50);
+ //fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",iBin[0],binEdges[0]);
+ fHist[iq][iqType][icut].fPt = new TH1F(hname,hname+";p_{T} (GeV/c)",60,0.25,30.25); // mj to compare with FONLL
hname = hnopt+"Y_"+kqTypeLabel[iqType];
fHist[iq][iqType][icut].fY = new TH1F(hname,hname,150,-7.5,7.5);
hname = hnopt+"Eta_"+kqTypeLabel[iqType];
if (icut == 0){
hname = hnopt+"Nq_"+kqTypeLabel[kQuark];
- fHistComm[iq][icut].fNq = new TH1F(hname,hname,10,-0.5,9.5);
+ fHistComm[iq][icut].fNq = new TH1F(hname,hname,50,-0.5,49.5);
hname = hnopt+"ProcessID_"+kqTypeLabel[kQuark];
fHistComm[iq][icut].fProcessID = new TH1F(hname,hname,21,-10.5,10.5);
}
#include "AliHFEpidTPC.h"
#include "AliHFEpidTRD.h"
#include "AliHFEpidTOF.h"
+#include "AliHFEpidEMCAL.h"
#include "AliHFEpidMC.h"
#include "AliHFEvarManager.h"
"TPCPID",
"TRDPID",
"TOFPID",
+ "EMCALPID",
"UndefinedPID"
};
fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPCPID");
fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRDPID");
fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOFPID");
+ fDetectorPID[kEMCALpid] = new AliHFEpidEMCAL("EMCALPID");
}
else if(!detector.CompareTo("TPC")) detectorID = kTPCpid;
else if(!detector.CompareTo("TRD")) detectorID = kTRDpid;
else if(!detector.CompareTo("TOF")) detectorID = kTOFpid;
+ else if(!detector.CompareTo("EMCAL")) detectorID = kEMCALpid;
else AliError("Detector not available");
if(detectorID == kUndefined) return;
AliDebug(2, Form("MC Information available, Filling container %s", mccontname.Data()));
if(fVarManager->IsSignalTrack()) {
fVarManager->FillContainerStepname(cont, mccontname.Data(), SortedDetectorName(idet), kTRUE);
- if(cont->GetCorrelationMatrix("correlationstepafterTOF")){
- TString tstept("TOFPID");
- if(!tstept.CompareTo(SortedDetectorName(idet))) {
- fVarManager->FillCorrelationMatrix(cont->GetCorrelationMatrix("correlationstepafterTOF"));
- //printf("Step %s\n",(const char*) SortedDetectorName(idet));
- }
- }
- }
+ if(cont->GetCorrelationMatrix("correlationstepafterTOF")){
+ TString tstept("TOFPID");
+ if(!tstept.CompareTo(SortedDetectorName(idet))) {
+ fVarManager->FillCorrelationMatrix(cont->GetCorrelationMatrix("correlationstepafterTOF"));
+ //printf("Step %s\n",(const char*) SortedDetectorName(idet));
+ }
+ }
+ }
}
// The PID will NOT fill the double counting information
}
}
//____________________________________________________________
-void AliHFEpid::ConfigureTPCrejection(){
+void AliHFEpid::ConfigureTPCrejection(const char *lowerCutParam, Double_t *params){
//
- // Combined TPC-TOF PID, combination is discribed in the funtion MakePidTpcTof
+ // Combined TPC-TOF PID
+ // if no function parameterizaion is given, then the default one (exponential) is chosen
//
- if(HasMCData()) printf("Configuring TPC for MC\n");
+ 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(3);
//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", "[0] * TMath::Exp([1]*x) + [2]", 0, 20);
+ TF1 *lowerCut = new TF1("lowerCut", lowerCutParam == NULL ? "[0] * TMath::Exp([1]*x) + [2]": lowerCutParam, 0, 20);
upperCut->SetParameter(0, 3.);
//upperCut->SetParameter(0, 2.7);
//upperCut->SetParameter(1, -0.4357);
-
- if(HasMCData()) lowerCut->SetParameter(0, -2.5);
- else lowerCut->SetParameter(0, -3.71769);
- //else lowerCut->SetParameter(0, -3.7);
-
- lowerCut->SetParameter(1, -0.40263);
- //lowerCut->SetParameter(1, -0.8);
-
- if(HasMCData()) lowerCut->SetParameter(2, -2.2);
- else lowerCut->SetParameter(2, 0.267857);
- //else lowerCut->SetParameter(2, -0.35);
-
+ if(params){
+ for(Int_t ipar = 0; ipar < lowerCut->GetNpar(); ipar++) lowerCut->SetParameter(ipar, params[ipar]);
+ } else {
+ // Set default parameterization
+ if(HasMCData()) lowerCut->SetParameter(0, -2.5);
+ else lowerCut->SetParameter(0, -3.71769);
+ //else lowerCut->SetParameter(0, -3.7);
+
+ lowerCut->SetParameter(1, -0.40263);
+ //lowerCut->SetParameter(1, -0.8);
+
+ if(HasMCData()) lowerCut->SetParameter(2, -2.2);
+ else lowerCut->SetParameter(2, 0.267857);
+ //else lowerCut->SetParameter(2, -0.35);
+ }
+
if(tpcpid){
tpcpid->SetTPCnSigma(2);
tpcpid->SetUpperSigmaCut(upperCut);
kTPCpid = 3,
kTRDpid = 4,
kTOFpid = 5,
- kNdetectorPID = 6
+ kEMCALpid = 6,
+ kNdetectorPID = 7
};
AliHFEpid();
AliHFEpid(const Char_t *name);
//-----Configure PID detectors with predefined stettings------
void ConfigureTPCasymmetric(Double_t pmin = 0.1, Double_t pmax = 20., Double_t sigmamin = -0.2, Double_t sigmamax = 5.);
void ConfigureTPCrejectionSimple();
- void ConfigureTPCrejection();
+ void ConfigureTPCrejection(const char *lowerCutParam = NULL, Double_t *params = NULL);
void ConfigureTPCstrategyParis();
//------------------------------------------------------------
while((o = trackIter())){
if(!TString(o->IsA()->GetName()).CompareTo("AliESDtrack")){
// work on local copy in order to not spoil others
- esdtrack = new AliESDtrack(*(dynamic_cast<AliESDtrack *>(o)));
+ esdtrack = new AliESDtrack(*(static_cast<AliESDtrack *>(o)));
if(!esdtrack) continue;
} else if(!TString(o->IsA()->GetName()).CompareTo("AliAODrack")){
// Bad hack: Fill ESD track with AOD information
if(!fMC) return;
AliDebug(3, Form("Doing Purity checks for species %d", species));
Int_t pdg = GetPDG(species);
- Char_t hname[256];
+ TString hname;
switch(species){
case AliHFEV0cuts::kRecoElectron:
- sprintf(hname, "purityElectron");
+ hname = "purityElectron";
break;
case AliHFEV0cuts::kRecoPionK0:
- sprintf(hname, "purityPionK0");
+ hname = "purityPionK0";
break;
case AliHFEV0cuts::kRecoPionL:
- sprintf(hname, "purityPionL");
+ hname = "purityPionL";
break;
case AliHFEV0cuts::kRecoProton:
- sprintf(hname, "purityProton");
+ hname = "purityProton";
break;
default: // non investigated species
AliDebug(3, "Species not investigated");
// pion and proton efficiency and the thresholds
//
Long_t status = 0;
- const Char_t *detname[4] = {"ITS", "TPC", "TRD", "TOF"};
- Char_t specname[256];
+ const TString detname[4] = {"ITS", "TPC", "TRD", "TOF"};
+ TString specname;
switch(species){
case AliHFEV0cuts::kRecoElectron:
- sprintf(specname, "Electron");
+ specname = "Electron";
break;
case AliHFEV0cuts::kRecoPionK0:
- sprintf(specname, "PionK0");
+ specname = "PionK0";
break;
case AliHFEV0cuts::kRecoPionL:
- sprintf(specname, "PionL");
+ specname = "PionL";
break;
case AliHFEV0cuts::kRecoProton:
- sprintf(specname, "Proton");
+ specname = "Proton";
break;
default:
AliDebug(2, Form("Species %d not investigated", species));
quantities[0] = pTPC;
Bool_t detFlagSet = kFALSE;
for(Int_t idet = 0; idet < 4; idet++){
- Char_t histname[256], histnameMC[256];
- sprintf(histname, "h%s_El_like_%s", detname[idet], specname);
- sprintf(histnameMC, "h%s_El_like_MC_%s", detname[idet], specname);
+ TString histname, histnameMC;
+ histname = "h" + detname[idet] = "_El_like_" + specname;
+ histnameMC = "h" + detname[idet] = "_El_like_MC_" + specname;
+
switch(idet){
case kITS: esdTrack->GetITSpid(pidProbs);
detFlagSet = status & AliESDtrack::kITSpid;
//
// Fill the PID response of different detectors to V0 daughter particles
//
- Char_t hname[256] = "";
- Char_t hname2[256] = "";
- Char_t hname3[256] = "";
+ TString hname, hname2, hname3;
- const Char_t *typeName[5] = {"Electron", "PionK0", "PionL", "Kaon", "Proton"};
+ const TString typeName[5] = {"Electron", "PionK0", "PionL", "Kaon", "Proton"};
const Int_t typePID[5] = {0, 2, 2, 3, 4};
// PID THnSparse
esdTrack->GetITSdEdxSamples(dEdxSamples);
Int_t nSamples = 0;
Double_t dEdxSum = 0.;
- sprintf(hname, "hITS_dEdx_%s", typeName[species]);
+ hname = "hITS_dEdx_" + typeName[species];
for(Int_t i=0; i<4; ++i){
if(dEdxSamples[i] > 0){
nSamples++;
fOutput->Fill("hITS_dEdx_nSamples", nSamples);
Double_t signal = esdTrack->GetITSsignal();
- sprintf(hname, "hITS_Signal_%s", typeName[species]);
+ hname = "hITS_Signal_" + typeName[species];
fOutput->Fill(hname, p, signal);
data[4] = signal;
// ITS number of signas
Double_t nsigma = fESDpid->NumberOfSigmasITS(esdTrack,(AliPID::EParticleType)typePID[species]);
- sprintf(hname, "hITS_nSigma_%s", typeName[species]);
+ hname = "hITS_nSigma_" + typeName[species];
fOutput->Fill(hname, p, nsigma);
-
// ITS PID response
Double_t itsPID[5] = {-1, -1, -1, -1, -1};
esdTrack->GetITSpid(itsPID);
Int_t ix = GetMaxPID(itsPID);
- sprintf(hname, "hITS_PID_p_%s", typeName[species]);
+ hname = "hITS_PID_" + typeName[species];
fOutput->Fill(hname, p, ix);
}//.. kITSpid
Double_t p = esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P();
// TPC dEdx
Double_t dEdx = esdTrack->GetTPCsignal();
- sprintf(hname, "hTPC_dEdx_%s", typeName[species]);
+ hname = "hTPC_dEdx_" + typeName[species];
fOutput->Fill(hname, p, dEdx);
data[6] = dEdx;
//TPC number of sigmas
Double_t nsigma = fESDpid->NumberOfSigmasTPC(esdTrack,(AliPID::EParticleType)typePID[species]);
- sprintf(hname, "hTPC_nSigma_%s", typeName[species]);
+ hname = "hTPC_nSigma_" + typeName[species];
fOutput->Fill(hname, p, nsigma);
data[7] = nsigma;
// TPC PID response
- sprintf(hname, "hTPC_PID_p_%s", typeName[species]);
+ hname = "hTPC_PID_p_" + typeName[species];
Double_t tpcPID[5] = {-1, -1, -1, -1, -1};
esdTrack->GetTPCpid(tpcPID);
Int_t ix = GetMaxPID(tpcPID);
// TRD number of tracklets
Int_t ntrk = esdTrack->GetTRDntrackletsPID();
- sprintf(hname, "hTRD_trk_%s", typeName[species]);
+ hname = "hTRD_trk_" + typeName[species];
fOutput->Fill(hname, p, ntrk);
data[8] = ntrk;
// TRD PID response
- sprintf(hname, "hTRD_PID_p_%s", typeName[species]);
+ hname = "hTRD_PID_p_" + typeName[species];
Double_t trdPID[5] = {-1., -1., -1., -1., -1};
esdTrack->GetTRDpid(trdPID);
Int_t ix = GetMaxPID(trdPID);
fOutput->Fill(hname, p, ix);
// TRD n clusters
Int_t ncls = esdTrack->GetTRDncls();
- sprintf(hname, "hTRD_cls_%s", typeName[species]);
+ hname = "hTRD_cls_" + typeName[species];
fOutput->Fill(hname, p, ncls);
data[9] = ncls;
// TRD - per tracklet - dEdx per, likelihood
- sprintf(hname, "hTRD_Nslc_%s", typeName[species]);
- sprintf(hname2, "hTRD_slc_%s", typeName[species]);
- sprintf(hname3, "hTRD_dEdx_%s", typeName[species]);
+ hname = "hTRD_Nslc_" + typeName[species];
+ hname2 = "hTRD_slc_" + typeName[species];
+ hname3 = "hTRD_dEdx_" + typeName[species];
Int_t nSlices = esdTrack->GetNumberOfTRDslices();
Double_t sumTot = 0.;
Int_t not0Tot = 0;
Double_t t0 = fESDpid->GetTOFResponse().GetStartTime(esdTrack->P());
//TOF beta
- sprintf(hname, "hTOF_beta_%s", typeName[species]);
+ hname = "hTOF_beta_" + typeName[species];
Float_t beta = TOFbeta(esdTrack);
fOutput->Fill(hname, p, beta);
fOutput->Fill("hTOF_beta_all", p, beta);
// TOF nSigma
Double_t nsigma = fESDpid->NumberOfSigmasTOF(esdTrack,(AliPID::EParticleType)typePID[species], t0);
- sprintf(hname, "hTOF_nSigma_%s", typeName[species]);
+ hname = "hTOF_nSigma_" + typeName[species];
fOutput->Fill(hname, p, nsigma);
if(beta > 0.97 && beta < 1.03){
data[11] = 1;
}
// TOF PID response
- sprintf(hname, "hTOF_PID_p_%s", typeName[species]);
+ hname = "hTOF_PID_p_" + typeName[species];
Double_t tofPID[5] = {-1., -1., -1., -1., -1};
esdTrack->GetTOFpid(tofPID);
Int_t ix = GetMaxPID(tofPID);
if((esd = dynamic_cast<AliESDEvent *>(fEvent))){
AliESDtrack *track = NULL, *partnerTrack = NULL;
while((hfetrack = dynamic_cast<AliHFEV0info *>(candidates()))){
- if(!hfetrack) continue;
track = dynamic_cast<AliESDtrack *>(hfetrack->GetTrack());
if(!track) continue;
partnerTrack = esd->GetTrack(hfetrack->GetPartnerID());
AliHFEpidObject::AnalysisType_t anaType = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
const AliVTrack *vtrack = dynamic_cast<const AliVTrack *>(track->GetRecTrack());
- if(!(vtrack && vtrack->GetStatus() & AliESDtrack::kTOFpid)) return 0;
+ if(!(vtrack && (vtrack->GetStatus() & AliESDtrack::kTOFpid) && !(vtrack->GetStatus() & AliESDtrack::kTOFmismatch))) return 0;
AliDebug(2, "Track Has TOF PID");
if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTOFpid, AliHFEdetPIDqa::kBeforePID);
(dynamic_cast<TH1F *>(fSecVtxList->At(step)))->Fill(track->Pt()); // electrons tagged
- if(HasMCData()){
+ if(HasMCData() && fMCQA){
if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) return;
mcpart = mctrack->Particle();
// Return PDG code of a particle itself
//
- if(!fMC){
- AliDebug(1, "No MC Event Available\n");
- return 0;
- }
- if(!fMCQA){
- AliDebug(1, "No MCQA Available\n");
- return 0;
- }
- if(!track){
- AliDebug(1, "Track not Available\n");
- return 0;
- }
+ if(!fMC){
+ AliDebug(1, "No MC Event Available\n");
+ return 0;
+ }
+ if(!fMCQA){
+ AliDebug(1, "No MCQA Available\n");
+ return 0;
+ }
+ if(!track){
+ AliDebug(1, "Track not Available\n");
+ return 0;
+ }
- TString sourcetype = track->IsA()->GetName();
- const AliVParticle *mctrack = NULL;
- TParticle *mcpart = NULL;
- if(!sourcetype.CompareTo("AliESDtrack") || !sourcetype.CompareTo("AliAODTrack")){
- mctrack = fMC->GetTrack(TMath::Abs(track->GetLabel()));
- } else mctrack = track;
- if(!mctrack) return 0;
+ TString sourcetype = track->IsA()->GetName();
+ const AliVParticle *mctrack = NULL;
+ TParticle *mcpart = NULL;
+ if(!sourcetype.CompareTo("AliESDtrack") || !sourcetype.CompareTo("AliAODTrack")){
+ 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")){
- const AliMCParticle *esdmc = dynamic_cast<const AliMCParticle *>(mctrack);
- mcpart = esdmc->Particle();
- eSource=fMCQA->GetElecSource(mcpart);
- } else {
- return -1;
- }
-
- return eSource;
+ TString mctype = mctrack->IsA()->GetName();
+ Int_t eSource = 0;
+ if(!mctype.CompareTo("AliMCParticle")){
+ const AliMCParticle *esdmc = dynamic_cast<const AliMCParticle *>(mctrack);
+ if(esdmc){
+ mcpart = esdmc->Particle();
+ eSource=fMCQA->GetElecSource(mcpart);
+ }
+ } else {
+ return -1;
+ }
+ return eSource;
}
#include "AliESDtrack.h"
#include "AliLog.h"
#include "AliPID.h"
+#include "AliExternalTrackParam.h"
#include "AliHFEcollection.h"
#include "AliHFEcontainer.h"
, fCFM(NULL)
, fQAhistos(NULL)
, fClean(kFALSE)
+ , fMagneticField(0.0)
+ , fVariablesTRD(kFALSE)
{
//
// Default constructor
, fCFM(ref.fCFM)
, fQAhistos(ref.fQAhistos)
, fClean(ref.fClean)
+ , fMagneticField(ref.fMagneticField)
+ , fVariablesTRD(ref.fVariablesTRD)
{
//
// Copy constructor
fCFM = ref.fCFM;
fQAhistos = ref.fQAhistos;
fClean = ref.fClean;
+ fMagneticField = ref.fMagneticField;
+ fVariablesTRD = ref.fVariablesTRD;
if(ref.fContainer) InitContainer();
// temporarily special QA
fQAhistos = new AliHFEcollection("taggedTrackQA", "Special QA for the TaggedTrackAnalysis");
fQAhistos->CreateTH2F("TPCclusters2_1", "TPCclusterInfo for findable clusters for 2 neighbors", 30, 0.1, 10., 162, 0., 161.);
- fQAhistos->CreateTH2F("TPCclusters2_0", "TPCclusterInfo for the ratio for 2 neighbors", 30, 0.1, 10., 100, 0., 1.);
+ fQAhistos->CreateTH2F("TPCclusters2_0", "TPCclusterInfo for the ratio for 2 neighbors", 30, 0.1, 10., 110, 0., 1.1);
fQAhistos->BinLogAxis("TPCclusters2_1", 0); // pt axis in logarithmic binning
fQAhistos->BinLogAxis("TPCclusters2_0", 0); // pt axis in logarithmic binning
}
//
// Filter tracks tagged by V0 PID class
//
+ //
fVarManager->NewTrack(track, NULL, 0., abinitioPID, kTRUE);
+
+
+
+ // Phi Angle
+ if(fVariablesTRD) {
+ AliESDtrack *esdtrackc = dynamic_cast<AliESDtrack *>(track);
+ if(esdtrackc) {
+
+ const AliExternalTrackParam *trueparam = NULL;
+ if(esdtrackc->GetOuterParam()) {
+ trueparam = esdtrackc->GetOuterParam();
+ fVarManager->NewTrack((AliVParticle *)trueparam, NULL, 0., abinitioPID, kTRUE);
+ }
+ else return;
+ }
+ }
+
+
+ // Try a loose cut to reject pion contamination
+ if(fClean) {
+ if(abinitioPID == AliPID::kElectron){
+ AliHFEpidTPC *pidTPC = (AliHFEpidTPC *) fPID->GetDetPID(AliHFEpid::kTPCpid);
+ if(pidTPC) {
+ Double_t numberOfSigmaTPC = pidTPC->NumberOfSigmas(track,AliPID::kElectron,AliHFEpidObject::kESDanalysis);
+ if(numberOfSigmaTPC < -5) return;
+ }
+ }
+ }
+ // temporarily monitoring of the number of TPC clusters
+ AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
+ if(esdtrack && abinitioPID == AliPID::kElectron){
+ if((esdtrack->GetITSClusterMap() & (BIT(0) | BIT(1))) && (TMath::Abs(esdtrack->Eta()) < 0.8)){ // Only select quasi-primary tracks
+ fQAhistos->Fill("TPCclusters2_1", track->Pt(), esdtrack->GetTPCClusterInfo(2,1));
+ fQAhistos->Fill("TPCclusters2_0", track->Pt(), esdtrack->GetTPCNclsF() > 0 ? esdtrack->GetTPCClusterInfo(2,1)/esdtrack->GetTPCNclsF() : 0.);
+ }
+ }
+
Int_t offset = AliHFEcuts::kStepRecKineITSTPC;
fVarManager->FillContainer(fCFM->GetParticleContainer(), 0); // Fill Container without filtering
Bool_t survived = kTRUE;
for(Int_t icut = AliHFEcuts::kStepRecKineITSTPC; icut <= AliHFEcuts::kStepHFEcutsTRD; icut++){
AliDebug(2, Form("Checking cut %d for species %s", icut + AliHFEcuts::kNcutStepsMCTrack, AliPID::ParticleName(abinitioPID)));
- /*
- TObjArray *cutlist = fCFM->GetParticleCutsList(icut + AliHFEcuts::kNcutStepsMCTrack);
- if(!cutlist){
+ /*
+ TObjArray *cutlist = fCFM->GetParticleCutsList(icut + AliHFEcuts::kNcutStepsMCTrack);
+ if(!cutlist){
AliDebug(2, Form("No cuts for step %d set", icut + AliHFEcuts::kNcutStepsMCTrack));
- } else {
+ } else {
AliDebug(2, Form("Cut Collection %s", cutlist->GetName()));
TIter cutiter(cutlist);
AliCFCutBase *cut;
while((cut = dynamic_cast<AliCFCutBase *>(cutiter()))){
- AliDebug(2, Form("Cut object %s, QA on? %s", cut->GetName(), cut->IsQAOn() ? "yes" : "no"));
+ AliDebug(2, Form("Cut object %s, QA on? %s", cut->GetName(), cut->IsQAOn() ? "yes" : "no"));
}
- }
- */
+ }
+ */
//if(!fCFM->CheckParticleCuts(icut + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)){
if(!fCuts->CheckParticleCuts(icut + AliHFEcuts::kNcutStepsMCTrack, (TObject *)track)){
AliDebug(2, Form("Track didn' survive cut %d", icut + AliHFEcuts::kNcutStepsMCTrack));
}
if(survived){
- AliDebug(2, "Use track in the PID");
- // Try a loose cut to reject pion contamination
- if(fClean) {
- if(abinitioPID == AliPID::kElectron){
- AliHFEpidTPC *pidTPC = (AliHFEpidTPC *) fPID->GetDetPID(AliHFEpid::kTPCpid);
- if(pidTPC) {
- Double_t numberOfSigmaTPC = pidTPC->NumberOfSigmas(track,AliPID::kElectron,AliHFEpidObject::kESDanalysis);
- if(numberOfSigmaTPC < -5) return;
- }
- }
- }
- // temporarily monitoring of the number of TPC clusters
- AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
- if(esdtrack && abinitioPID == AliPID::kElectron){
- fQAhistos->Fill("TPCclusters2_1", track->Pt(), esdtrack->GetTPCClusterInfo(2,1));
- fQAhistos->Fill("TPCclusters2_0", track->Pt(), esdtrack->GetTPCClusterInfo(2,0));
- }
- // Apply PID
- AliHFEpidObject hfetrack;
- hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
- hfetrack.SetRecTrack(track);
- hfetrack.SetAbInitioPID(abinitioPID);
- fPID->SetVarManager(fVarManager);
- fPID->IsSelected(&hfetrack, fContainer, "taggedTrackContainer", fPIDqa);
- }
+ AliDebug(2, "Use track in the PID");
+ // Apply PID
+ AliHFEpidObject hfetrack;
+ hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
+ hfetrack.SetRecTrack(track);
+ hfetrack.SetAbInitioPID(abinitioPID);
+ fPID->SetVarManager(fVarManager);
+ fPID->IsSelected(&hfetrack, fContainer, "taggedTrackContainer", fPIDqa);
+ }
}
//____________________________________________________________
return fCuts->GetQAhistograms();
}
+
void InitContainer();
void ProcessTrack(AliVParticle *track, Int_t abinitioPID);
-
+
AliHFEcontainer *GetContainer() const { return fContainer; }
AliHFEpidQAmanager *GetPIDqa() const { return fPIDqa; }
TList * GetPIDQA() const;
TList * GetCutQA() const;
AliHFEcollection * GetQAcollection() const { return fQAhistos; }
Bool_t GetClean() const { return fClean; };
+ Double_t GetMagneticField() const { return fMagneticField; };
void SetCuts(AliHFEcuts *cuts);
void SetPID(AliHFEpid *pid);
void SetClean(Bool_t clean) { fClean = clean; };
+ void SetMagneticField(Double_t magneticField) { fMagneticField = magneticField; };
+ void SetVariablesTRD(Bool_t variablesTRD) { fVariablesTRD = variablesTRD; };
private:
enum{
AliCFManager *fCFM; // CF Manager used for the track filtering
AliHFEcollection *fQAhistos; // QA histos
Bool_t fClean; // Clean
+ Double_t fMagneticField; // Magnetic field
+ Bool_t fVariablesTRD; // Use phi angle at the first plane of the TRD
ClassDef(AliHFEtaggedTrackAnalysis, 0)
};
THnSparseF *hSignal = dynamic_cast<THnSparseF *>(fHistos->Get("tofnSigma"));
if(!hSignal) return NULL;
hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
- if(species > 0 && species < AliPID::kSPECIES)
+ 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"),
Int_t pdg = 0;
if(!TString(track->IsA()->GetName()).CompareTo("AliMCParticle")){
const AliMCParticle *mctrack = dynamic_cast<const AliMCParticle *>(track);
- pdg = mctrack->Particle()->GetPdgCode();
+ pdg = mctrack ? mctrack->Particle()->GetPdgCode() : 0;
} else if(!TString(track->IsA()->GetName()).CompareTo("AliAODMCParticle")){
const AliAODMCParticle *aodmctrack = dynamic_cast<const AliAODMCParticle *>(track);
- pdg = aodmctrack->GetPdgCode();
+ pdg = aodmctrack ? aodmctrack->GetPdgCode() : 0;
}
return pdg;
}
goodTrack = kTRUE;
for(Int_t icut = 0; icut < nStep; icut++){
cutStep = dynamic_cast<AliHFEcutStep *>(fCutSteps->UncheckedAt(icut));
- if(!cutStep->IsSelected(track)){
+ if(cutStep && (!cutStep->IsSelected(track))){
// track cut away
goodTrack = kFALSE;
break;
fCutStep->AddCut(trackQuality);
AliHFEextraCuts *hfecuts = new AliHFEextraCuts("HFETPC","Extra cuts from the HFE group");
- hfecuts->SetClusterRatioTPC(0.6);
+ hfecuts->SetClusterRatioTPC(0.6, AliHFEextraCuts::kFoundOverCR);
fCutStep->AddCut(hfecuts);
AliCFTrackKineCuts *kineCuts = new AliCFTrackKineCuts((Char_t *)"RecKine", (Char_t *)"REC Kine Cuts");
fPionEfficiencies(NULL),
fProtonEfficiencies(NULL),
fKaonEfficiencies(NULL),
- fThresholds(NULL)
+ fThresholds(NULL),
+ fShowMessage(kFALSE)
{
//
// Default Constructor
fPionEfficiencies(NULL),
fProtonEfficiencies(NULL),
fKaonEfficiencies(NULL),
- fThresholds(NULL)
+ fThresholds(NULL),
+ fShowMessage(kFALSE)
{
//
// Main Constructor
fPionEfficiencies(NULL),
fProtonEfficiencies(NULL),
fKaonEfficiencies(NULL),
- fThresholds(NULL)
+ fThresholds(NULL),
+ fShowMessage(kFALSE)
{
//
// Copy constructor
}
for(Int_t itr = 4; itr <= 6; itr++){
- printf("========================================\n");
- printf("Analysing %d trackltes\n", itr);
- printf("========================================\n");
+ if(fShowMessage){
+ printf("========================================\n");
+ printf("Analysing %d trackltes\n", itr);
+ printf("========================================\n");
+ }
AnalyseNTracklets(itr);
}
}
return;
}
- printf("========================================\n");
- printf("Calculating threshold parameters\n");
- printf("========================================\n");
+ if(fShowMessage){
+ printf("========================================\n");
+ printf("Calculating threshold parameters\n");
+ printf("========================================\n");
+ }
TList *outlist = new TList;
outlist->SetName("thresholdTRD");
TList *lHistos = NULL, *lFormulas = NULL;
for(Int_t itracklet = 4; itracklet <= 6; itracklet++){
- printf("-------------------------------\n");
- printf("Processing %d tracklets\n", itracklet);
- printf("-------------------------------\n");
+ if(fShowMessage){
+ printf("-------------------------------\n");
+ printf("Processing %d tracklets\n", itracklet);
+ printf("-------------------------------\n");
+ }
lHistos = dynamic_cast<TList *>(fThresholds->FindObject(Form("%dTracklets", itracklet)));
if(!lHistos){
for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
- printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n");
- printf("Processing Electron Efficiency %f\n", fgkElectronEff[ieff]);
- printf(">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n");
+ if(fShowMessage){
+ printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n");
+ printf("Processing Electron Efficiency %f\n", fgkElectronEff[ieff]);
+ printf(">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n");
+ }
threshhist = dynamic_cast<TGraph *>(lHistos->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
if(!threshhist) continue;
Int_t threshbin = 0;
Double_t noElEff[2]; // value and error
for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
- printf("-----------------------------------------\n");
- printf("Doing Electron Efficiency %f\n", fgkElectronEff[ieff]);
- printf("-----------------------------------------\n");
+
+ if(fShowMessage){
+ printf("-----------------------------------------\n");
+ printf("Doing Electron Efficiency %f\n", fgkElectronEff[ieff]);
+ printf("-----------------------------------------\n");
+ }
effPi = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
effPi->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
effPr = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
}
//__________________________________________________________________
-void AliHFEtrdPIDqa::DrawTracklet(Int_t itracklet, Bool_t doFit){
+void AliHFEtrdPIDqa::DrawTracklet(Int_t itracklet, Double_t pmin, Double_t pmax, Bool_t doFit){
//
// Draw efficiencies and threshold as function of p
//
tr->GetXaxis()->SetTitle("p / GeV/c");
tr->GetYaxis()->SetTitle("Efficiency");
// Axis Range
- pi->GetYaxis()->SetRangeUser(0., 1.);
- pr->GetYaxis()->SetRangeUser(0., 1.);
- tr->GetYaxis()->SetRangeUser(0., 1.);
+ pi->GetYaxis()->SetRangeUser(1e-3, 1.);
+ pr->GetYaxis()->SetRangeUser(1e-3, 1.);
+ tr->GetYaxis()->SetRangeUser(1e-3, 1.);
+ if(pmin > 0 && pmax > 0.){
+ pi->GetXaxis()->SetRangeUser(pmin, pmax);
+ pr->GetXaxis()->SetRangeUser(pmin, pmax);
+ tr->GetXaxis()->SetRangeUser(pmin, pmax);
+ }
// Marker
pi->SetMarkerColor(kRed);
pi->SetMarkerStyle(20);
fThresholds = NULL;
}
}
+
+//__________________________________________________________________
+Double_t AliHFEtrdPIDqa::EvalPionEfficiency(Int_t ntls, Int_t eEff, Double_t p){
+ TList *graphs = dynamic_cast<TList *>(fPionEfficiencies->FindObject(Form("%dTracklets", ntls)));
+ if(!graphs) return -1.;
+ TGraph *measurement = dynamic_cast<TGraph *>(graphs->FindObject(Form("eff%d", eEff)));
+ if(!measurement) return -1.;
+ return measurement->Eval(p);
+}
+
+//__________________________________________________________________
+Double_t AliHFEtrdPIDqa::EvalProtonEfficiency(Int_t ntls, Int_t eEff, Double_t p){
+ TList *graphs = dynamic_cast<TList *>(fProtonEfficiencies->FindObject(Form("%dTracklets", ntls)));
+ if(!graphs) return -1.;
+ TGraph *measurement = dynamic_cast<TGraph *>(graphs->FindObject(Form("eff%d", eEff)));
+ if(!measurement) return -1.;
+ return measurement->Eval(p);
+}
+
+//__________________________________________________________________
+Double_t AliHFEtrdPIDqa::EvalThreshold(Int_t ntls, Int_t eEff, Double_t p){
+ TList *graphs = dynamic_cast<TList *>(fThresholds->FindObject(Form("%dTracklets", ntls)));
+ if(!graphs) return -1.;
+ TGraph *measurement = dynamic_cast<TGraph *>(graphs->FindObject(Form("eff%d", eEff)));
+ if(!measurement) return -1.;
+ return measurement->Eval(p);
+}
+
void Init();
void FinishAnalysis();
+ void ShowMessages() { fShowMessage = kTRUE; }
void StoreResults(const Char_t *filename = "HFEtrdPIDqa.root");
void SaveThresholdParameters(const Char_t * filename = "TRD.Thresholds.root");
- void DrawTracklet(Int_t tracklet, Bool_t doFit = kFALSE);
+ void DrawTracklet(Int_t tracklet, Double_t pmin = 0., Double_t pmax = 0., Bool_t doFit = kFALSE);
void ClearLists();
+ Double_t EvalPionEfficiency(Int_t ntls, Int_t eEff, Double_t p);
+ Double_t EvalProtonEfficiency(Int_t ntls, Int_t eEff, Double_t p);
+ Double_t EvalThreshold(Int_t ntls, Int_t eEff, Double_t p);
+
//---------------------------------------------------
// Getters for Histograms
THnSparseF *GetLikelihoodHistogram() const { return fHistos ? dynamic_cast<THnSparseF *>(fHistos->Get("fLikeTRD")) : NULL; }
TList *fKaonEfficiencies; //! List for Kaon efficiencies
TList *fThresholds; //! List for Threshold Graphs
+
+ Bool_t fShowMessage; // Display debug messages
ClassDef(AliHFEtrdPIDqa, 3) // QA class for TRD PID
};
//
#include <TAxis.h>
+#include <TBrowser.h>
#include <TH2.h>
#include <THnSparse.h>
#include <TString.h>
+#include "AliESDtrack.h"
#include "AliLog.h"
#include "AliPID.h"
return count + 1;
}
+//_________________________________________________________
+void AliHFEtrdPIDqaV1::Browse(TBrowser *b){
+ //
+ // Browse the PID QA
+ //
+ if(b){
+ if(fHistos){
+ b->Add(fHistos, fHistos->GetName());
+
+ // Make Projections of the dE/dx Spectra and add them to a new Folder
+ TString specnames[4] = {"All", "Electrons", "Pions", "Protons"};
+ Int_t specind[4] = {-1, AliPID::kElectron, AliPID::kPion, AliPID::kProton};
+ TList *listTM = new TList;
+ listTM->SetOwner();
+ TList *listLike = new TList;
+ listLike->SetOwner();
+ TList *listCharge = new TList;
+ listCharge->SetOwner();
+ TList *listTPCnsigma = new TList;
+ 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]);
+ hptr->SetName(Form("hTRDtm%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
+ listTM->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");
+ b->Add(listTPCnsigma, "Projections TPC spectra");
+ }
+ }
+}
+
//____________________________________________________________
void AliHFEtrdPIDqaV1::Initialize(){
//
Double_t maxTRDcharge[4] = {kMaxPID, kMaxP, 100000., 2.};
fHistos->CreateTHnSparse("hTRDcharge", "Total TRD charge; species; p [GeV/c]; TRD charge [a.u.]; selection step", 4, nBinsTRDcharge, minTRDcharge, maxTRDcharge);
fHistos->Sumw2("hTRDcharge");
+ // Monitoring of the TRD truncated mean according to version 1
+ const Int_t kTRDtmBins = 1000;
+ Int_t nBinsTRDtm[4] = {kPIDbins, kPbins, kTRDtmBins, kSteps};
+ Double_t minTRDtm[4] = {kMinPID, kMinP, 0., 0.};
+ Double_t maxTRDtm[4] = {kMaxPID, kMaxP, 20000., 2.};
+ fHistos->CreateTHnSparse("hTRDtruncatedMean", "TRD truncated Mean; species; p [GeV/c]; TRD signal [a.u.]; selection step", 4, nBinsTRDtm, minTRDtm, maxTRDtm);
+ fHistos->Sumw2("hTRDcharge");
}
//____________________________________________________________
//
AliDebug(1, Form("QA started for TRD PID for step %d", (Int_t)step));
Int_t species = track->GetAbInitioPID();
+ if(species >= AliPID::kSPECIES) species = -1;
AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
AliHFEpidTRD *trdpid = dynamic_cast<AliHFEpidTRD *>(fQAmanager->GetDetectorPID(AliHFEpid::kTRDpid));
container[2] = trdpid->GetElectronLikelihood(track->GetRecTrack(), anatype);
fHistos->Fill("hTRDlikelihood", container);
+ if(track->IsESDanalysis()){
+ container[2] = trdpid->GetTRDSignalV1(dynamic_cast<const AliESDtrack *>(track->GetRecTrack()));
+ fHistos->Fill("hTRDtruncatedMean", container);
+ }
for(UInt_t ily = 0; ily < 6; ily++){
container[2] = trdpid->GetChargeLayer(track->GetRecTrack(), ily, anatype);
fHistos->Fill("hTRDcharge", container);
}
- if(species >= AliPID::kSPECIES) species = -1;
}
//_________________________________________________________
//
// Get the TPC control histogram for the TRD selection step (either before or after PID)
//
- printf("histos :%p\n", fHistos);
THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTPCsigma"));
if(!histo){
AliError("QA histogram monitoring TPC nSigma not available");
return hSpec;
}
+//_________________________________________________________
+TH2 *AliHFEtrdPIDqaV1::MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species){
+ //
+ // Get the TPC control histogram for the TRD selection step (either before or after PID)
+ //
+ THnSparseF *histo = dynamic_cast<THnSparseF *>(fHistos->Get("hTRDtruncatedMean"));
+ if(!histo){
+ AliError("QA histogram monitoring TPC nSigma not available");
+ return NULL;
+ }
+ if(species > -1 && species < AliPID::kSPECIES){
+ // cut on species (if available)
+ histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
+ }
+ histo->GetAxis(3)->SetRange(step + 1, step + 1);
+
+ TH2 *hSpec = histo->Projection(2, 1);
+ // construct title and name
+ 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());
+ 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());
+ return hSpec;
+}
+
//_________________________________________________________
TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
//
class AliAODTrack;
class AliHFEcollection;
class AliHFEpidObject;
+class TBrowser;
class TCollection;
class AliHFEtrdPIDqaV1 : public AliHFEdetPIDqa{
AliHFEtrdPIDqaV1(const AliHFEtrdPIDqaV1 &c);
AliHFEtrdPIDqaV1 &operator=(const AliHFEtrdPIDqaV1 &o);
~AliHFEtrdPIDqaV1(){}
+ virtual Long64_t Merge(TCollection *coll);
+ virtual void Browse(TBrowser *b);
+ virtual Bool_t IsFolder() const { return kTRUE; };
virtual void Initialize();
- virtual Long64_t Merge(TCollection *coll);
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);
protected:
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff