#include <TString.h>
#include <TTree.h>
+#include "AliAODInputHandler.h"
+#include "AliAODMCParticle.h"
+#include "AliAODTrack.h"
#include "AliCFContainer.h"
#include "AliCFManager.h"
#include "AliCFEffGrid.h"
#include "AliStack.h"
#include "AliHFEpid.h"
+#include "AliHFEcollection.h"
#include "AliHFEcuts.h"
#include "AliHFEmcQA.h"
+#include "AliHFEpairs.h"
+#include "AliHFEsecVtxs.h"
#include "AliHFEsecVtx.h"
+#include "AliHFEelecbackground.h"
#include "AliAnalysisTaskHFE.h"
//____________________________________________________________
, fQAlevel(0)
, fPIDdetectors("")
, fPIDstrategy(0)
- , fESD(0x0)
- , fMC(0x0)
- , fCFM(0x0)
- , fCorrelation(0x0)
- , fPIDperformance(0x0)
- , fPID(0x0)
- , fCuts(0x0)
- , fSecVtx(0x0)
- , fMCQA(0x0)
- , fNEvents(0x0)
- , fNElectronTracksEvent(0x0)
- , fQA(0x0)
- , fOutput(0x0)
- , fHistMCQA(0x0)
- , fHistSECVTX(0x0)
+ , fPlugins(0)
+ , fRecEvent(NULL)
+ , fMC(NULL)
+ , fCFM(NULL)
+ , fCorrelation(NULL)
+ , fPIDperformance(NULL)
+ , fSignalToBackgroundMC(NULL)
+ , fPID(NULL)
+ , fCuts(NULL)
+ , fSecVtx(NULL)
+ , fElecBackGround(NULL)
+ , fMCQA(NULL)
+ , fNEvents(NULL)
+ , fNElectronTracksEvent(NULL)
+ , fQA(NULL)
+ , fOutput(NULL)
+ , fHistMCQA(NULL)
+ , fHistSECVTX(NULL)
+ , fHistELECBACKGROUND(NULL)
+// , fQAcoll(0x0)
{
//
// Dummy constructor
DefineOutput(0, TH1I::Class());
DefineOutput(1, TList::Class());
DefineOutput(2, TList::Class());
+// DefineOutput(3, TList::Class());
// Initialize cuts
fPID = new AliHFEpid;
-
- SetHasMCData();
}
//____________________________________________________________
, fQAlevel(0)
, fPIDdetectors("")
, fPIDstrategy(0)
- , fESD(0x0)
- , fMC(0x0)
- , fCFM(0x0)
- , fCorrelation(0x0)
- , fPIDperformance(0x0)
- , fPID(0x0)
- , fCuts(0x0)
- , fSecVtx(0x0)
- , fMCQA(0x0)
- , fNEvents(0x0)
- , fNElectronTracksEvent(0x0)
- , fQA(0x0)
- , fOutput(0x0)
- , fHistMCQA(0x0)
- , fHistSECVTX(0x0)
+ , fPlugins(0)
+ , fRecEvent(NULL)
+ , fMC(NULL)
+ , fCFM(NULL)
+ , fCorrelation(NULL)
+ , fPIDperformance(NULL)
+ , fSignalToBackgroundMC(NULL)
+ , fPID(NULL)
+ , fCuts(NULL)
+ , fSecVtx(NULL)
+ , fElecBackGround(NULL)
+ , fMCQA(NULL)
+ , fNEvents(NULL)
+ , fNElectronTracksEvent(NULL)
+ , fQA(NULL)
+ , fOutput(NULL)
+ , fHistMCQA(NULL)
+ , fHistSECVTX(NULL)
+ , fHistELECBACKGROUND(NULL)
+// , fQAcoll(0x0)
{
//
// Default constructor
- //
+ //
DefineInput(0, TChain::Class());
DefineOutput(0, TH1I::Class());
DefineOutput(1, TList::Class());
DefineOutput(2, TList::Class());
+// DefineOutput(3, TList::Class());
// Initialize cuts
fPID = new AliHFEpid;
-
- SetHasMCData();
}
//____________________________________________________________
, fQAlevel(ref.fQAlevel)
, fPIDdetectors(ref.fPIDdetectors)
, fPIDstrategy(ref.fPIDstrategy)
- , fESD(ref.fESD)
+ , fPlugins(ref.fPlugins)
+ , fRecEvent(ref.fRecEvent)
, fMC(ref.fMC)
, fCFM(ref.fCFM)
, fCorrelation(ref.fCorrelation)
, fPIDperformance(ref.fPIDperformance)
+ , fSignalToBackgroundMC(ref.fSignalToBackgroundMC)
, fPID(ref.fPID)
, fCuts(ref.fCuts)
, fSecVtx(ref.fSecVtx)
+ , fElecBackGround(ref.fElecBackGround)
, fMCQA(ref.fMCQA)
, fNEvents(ref.fNEvents)
, fNElectronTracksEvent(ref.fNElectronTracksEvent)
, fOutput(ref.fOutput)
, fHistMCQA(ref.fHistMCQA)
, fHistSECVTX(ref.fHistSECVTX)
+ , fHistELECBACKGROUND(ref.fHistELECBACKGROUND)
+// , fQAcoll(ref.fQAcoll)
{
//
// Copy Constructor
fQAlevel = ref.fQAlevel;
fPIDdetectors = ref.fPIDdetectors;
fPIDstrategy = ref.fPIDstrategy;
- fESD = ref.fESD;
+ fPlugins = ref.fPlugins;
+ fRecEvent = ref.fRecEvent;
fMC = ref.fMC;
fCFM = ref.fCFM;
fCorrelation = ref.fCorrelation;
fPIDperformance = ref.fPIDperformance;
+ fSignalToBackgroundMC = ref.fSignalToBackgroundMC;
fPID = ref.fPID;
fCuts = ref.fCuts;
fSecVtx = ref.fSecVtx;
+ fElecBackGround = ref.fElecBackGround;
fMCQA = ref.fMCQA;
fNEvents = ref.fNEvents;
fNElectronTracksEvent = ref.fNElectronTracksEvent;
fOutput = ref.fOutput;
fHistMCQA = ref.fHistMCQA;
fHistSECVTX = ref.fHistSECVTX;
+ fHistELECBACKGROUND = ref.fHistELECBACKGROUND;
+
+// fQAcoll = ref.fQAcoll;
return *this;
}
//
// Destructor
//
- if(fESD) delete fESD;
+ if(fRecEvent) delete fRecEvent;
if(fMC) delete fMC;
if(fPID) delete fPID;
if(fQA){
fHistSECVTX->Clear();
delete fHistSECVTX;
}
+ if(fHistELECBACKGROUND){
+ fHistELECBACKGROUND->Clear();
+ delete fHistELECBACKGROUND;
+ }
if(fSecVtx) delete fSecVtx;
+ if(fElecBackGround) delete fElecBackGround;
if(fMCQA) delete fMCQA;
if(fNEvents) delete fNEvents;
if(fCorrelation){
delete fCorrelation;
}
if(fPIDperformance) delete fPIDperformance;
+ if(fSignalToBackgroundMC) delete fSignalToBackgroundMC;
+// if(fQAcoll) delete fQAcoll;
}
//____________________________________________________________
esdchain->SetBranchStatus("Tracks", 1);
}
*/
- AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
- if(!esdH){
- AliError("No ESD input handler");
+ AliDebug(3, "Connecting input data");
+ AliVEventHandler *inputH = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+ if(!inputH){
+ AliError("Analysis Manager does not have any input event handler");
return;
- } else {
- fESD = esdH->GetEvent();
}
- AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
- if(!mcH){
- AliError("No MC truth handler");
- return;
+ if(IsAODanalysis()){
+ AliInfo("Analysis Mode: AOD Analysis");
+ AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler *>(inputH);
+ fRecEvent = aodH->GetEvent();
+ if(fRecEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName()))
+ fMC = aodH->MCEvent();
} else {
- fMC = mcH->MCEvent();
+ AliInfo("Analysis Mode: ESD Analysis");
+ AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler *>(inputH);
+ fRecEvent = esdH->GetEvent();
+ if(HasMCData()){
+ AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+ if(!mcH){
+ AliError("No MC truth handler");
+ return;
+ } else {
+ fMC = mcH->MCEvent();
+ }
+ }
}
}
// QA histograms are created if requested
// Called once per worker
//
+ AliDebug(3, "Creating Output Objects");
+ printf("Analysis Mode: %s Analysis\n", IsAODanalysis() ? "AOD" : "ESD");
+ printf("MC Data available %s\n", HasMCData() ? "Yes" : "No");
+
+ // example how to use the AliHFEcollection
+ //fQAcoll = new AliHFEcollection("fQAcoll", "QA");
+ //fQAcoll->CreateTH1F("fNevents", "Number of Events in the Analysis", 2, 0, 2);
+ //fQAcoll->CreateProfile("fNtrdclusters", "Number of TRD clusters as function of momentum; p[GeV/c]", 20, 0, 20);
+
fNEvents = new TH1I("nEvents", "Number of Events in the Analysis", 2, 0, 2); // Number of Events neccessary for the analysis and not a QA histogram
fNElectronTracksEvent = new TH1I("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100);
// First Step: TRD alone
// Initialize correction Framework and Cuts
fCFM = new AliCFManager;
MakeParticleContainer();
- // Temporary fix: Initialize particle cuts with 0x0
+ // Temporary fix: Initialize particle cuts with NULL
for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++)
- fCFM->SetParticleCutsList(istep, 0x0);
+ fCFM->SetParticleCutsList(istep, NULL);
if(!fCuts){
AliWarning("Cuts not available. Default cuts will be used");
fCuts = new AliHFEcuts;
fCuts->CreateStandardCuts();
}
+ if(IsAODanalysis()) fCuts->SetAOD();
fCuts->Initialize(fCFM);
if(fCuts->IsInDebugMode()) fQA->Add(fCuts->GetQAhistograms());
fOutput->AddAt(fCFM->GetParticleContainer(), 0);
fOutput->AddAt(fCorrelation, 1);
fOutput->AddAt(fPIDperformance, 2);
- fOutput->AddAt(fNElectronTracksEvent, 3);
+ fOutput->AddAt(fSignalToBackgroundMC, 3);
+ fOutput->AddAt(fNElectronTracksEvent, 4);
// Initialize PID
if(IsQAOn(kPIDqa)){
fPID->InitializePID(fPIDdetectors.Data()); // Only restrictions to TPC allowed
// mcQA----------------------------------
- if (IsQAOn(kMCqa)) {
+ if (HasMCData() && IsQAOn(kMCqa)) {
AliInfo("MC QA on");
if(!fMCQA) fMCQA = new AliHFEmcQA;
if(!fHistMCQA) fHistMCQA = new TList();
}
// secvtx----------------------------------
- if (IsSecVtxOn()) {
+ if (GetPlugin(kSecVtx)) {
AliInfo("Secondary Vertex Analysis on");
fSecVtx = new AliHFEsecVtx;
+ fSecVtx->SetHasMCData(HasMCData());
if(!fHistSECVTX) fHistSECVTX = new TList();
- fHistSECVTX->SetName("SecVtx");
- fSecVtx->CreateHistograms("secvtx_");
- TIter next(gDirectory->GetList());
- TObject *obj;
- int counter = 0;
- TString objname;
- while ((obj = next.Next())) {
- objname = obj->GetName();
- TObjArray *toks = objname.Tokenize("_");
- if (toks->GetEntriesFast()){
- TObjString *fpart = (TObjString *)(toks->UncheckedAt(0));
- if ((fpart->String()).CompareTo("secvtx") == 0) fHistSECVTX->AddAt(obj, counter++);
- }
- }
+ fSecVtx->CreateHistograms(fHistSECVTX);
fOutput->Add(fHistSECVTX);
- }
+ }
+
+ // background----------------------------------
+ if (GetPlugin(kIsElecBackGround)) {
+ AliInfo("Electron BackGround Analysis on");
+ fElecBackGround = new AliHFEelecbackground;
+ fElecBackGround->SetHasMCData(HasMCData());
+
+ if(!fHistELECBACKGROUND) fHistELECBACKGROUND = new TList();
+ fElecBackGround->CreateHistograms(fHistELECBACKGROUND);
+ fOutput->Add(fHistELECBACKGROUND);
+ }
}
//____________________________________________________________
//
// Run the analysis
//
- if(!fESD){
- AliError("No ESD Event");
+ AliDebug(3, "Starting Single Event Analysis");
+ if(!fRecEvent){
+ AliError("Reconstructed Event not available");
return;
}
- if(!fMC){
- AliError("No MC Event");
- return;
+ if(HasMCData()){
+ AliDebug(4, Form("MC Event: %p", fMC));
+ if(!fMC){
+ AliError("No MC Event, but MC Data required");
+ return;
+ }
}
if(!fCuts){
AliError("HFE cuts not available");
return;
}
- //fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC);
- Double_t container[6];
- // container for the output THnSparse
- Double_t dataE[5]; // [pT, eta, Phi, type, 'C' or 'B']
-
- // Loop over the Monte Carlo tracks to see whether we have overlooked any track
- AliMCParticle *mctrack = 0x0;
- Int_t nElectrons = 0;
+ if(HasMCData()) ProcessMC(); // Run the MC loop + MC QA in case MC Data are available
- if (IsSecVtxOn()) {
- fSecVtx->SetEvent(fESD);
- fSecVtx->SetStack(fMC->Stack());
- }
+ if(IsAODanalysis()) ProcessAOD();
+ else ProcessESD();
+ // Done!!!
+ PostData(0, fNEvents);
+ PostData(1, fOutput);
+ PostData(2, fQA);
+// PostData(3, fQAcoll->GetList());
+}
- // run mc QA ------------------------------------------------
- if (IsQAOn(kMCqa)) {
- AliDebug(2, "Running MC QA");
-
- fMCQA->SetStack(fMC->Stack());
- fMCQA->Init();
-
- Int_t nMCTracks = fMC->Stack()->GetNtrack();
-
- // loop over all tracks for decayed electrons
- for (Int_t igen = 0; igen < nMCTracks; igen++){
- TParticle* mcpart = fMC->Stack()->Particle(igen);
- fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kCharm);
- fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kBeauty);
- fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kCharm);
- fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kBeauty);
- fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 0); // no accept cut
- fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0); // no accept cut
- if (TMath::Abs(mcpart->Eta()) < 0.9) {
- fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
- fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
- }
- if (TMath::Abs(GetRapidity(mcpart)) < 0.5) {
- fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
- fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
- }
+//____________________________________________________________
+void AliAnalysisTaskHFE::Terminate(Option_t *){
+ //
+ // Terminate not implemented at the moment
+ //
+ if(GetPlugin(kPostProcess)){
+ fOutput = dynamic_cast<TList *>(GetOutputData(1));
+ if(!fOutput){
+ AliError("Results not available");
+ return;
}
- fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
- fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
-
- } // end of MC QA loop
- // -----------------------------------------------------------------
+ PostProcess();
+ }
+}
+//____________________________________________________________
+void AliAnalysisTaskHFE::Load(TString filename){
//
- // Loop MC
+ // Load Results into the task
//
- fCFM->SetMCEventInfo(fMC);
- for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){
- mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(imc));
+ fQA = NULL; fOutput = NULL; fNEvents = NULL;
+ TFile *input = TFile::Open(filename.Data());
+ if(!input || input->IsZombie()){
+ AliError("Cannot read file");
+ return;
+ }
+ TH1 *htmp = dynamic_cast<TH1I *>(input->Get("nEvents"));
+ if(htmp)
+ fNEvents = dynamic_cast<TH1I *>(htmp->Clone());
+ else
+ AliError("Event Counter histogram not found");
+ TList *ltmp = dynamic_cast<TList *>(input->Get("Results"));
+ if(ltmp)
+ fOutput = dynamic_cast<TList *>(ltmp->Clone());
+ else
+ AliError("Output Histograms not found");
+ ltmp = dynamic_cast<TList *>(input->Get("QA"));
+ if(ltmp)
+ fQA = dynamic_cast<TList *>(ltmp->Clone());
+ else
+ AliError("QA histograms not found");
+ input->Close();
+ delete input;
+ Int_t nObjects = 0;
+ if(fNEvents) nObjects++;
+ if(fOutput) nObjects++;
+ if(fQA) nObjects++;
+ AliInfo(Form("Loaded %d Objects into task", nObjects));
+}
- container[0] = mctrack->Pt();
- container[1] = mctrack->Eta();
- container[2] = mctrack->Phi();
+//____________________________________________________________
+void AliAnalysisTaskHFE::ProcessMC(){
+ //
+ // Runs the MC Loop (filling the container for the MC Cut Steps with the observables pt, eta and phi)
+ // In case MC QA is on also MC QA loop is done
+ //
+ AliDebug(3, "Processing MC Information");
+ Int_t nElectrons = 0;
+ if(IsESDanalysis()){
+ if (HasMCData() && IsQAOn(kMCqa)) {
+ AliDebug(2, "Running MC QA");
+
+ fMCQA->SetStack(fMC->Stack());
+ fMCQA->SetGenEventHeader(fMC->GenEventHeader());
+ fMCQA->Init();
+
+ Int_t nMCTracks = fMC->Stack()->GetNtrack();
+
+ // loop over all tracks for decayed electrons
+ for (Int_t igen = 0; igen < nMCTracks; igen++){
+ TParticle* mcpart = fMC->Stack()->Particle(igen);
+ fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kCharm);
+ fMCQA->GetQuarkKine(mcpart, igen, AliHFEmcQA::kBeauty);
+ fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kCharm);
+ fMCQA->GetHadronKine(mcpart, AliHFEmcQA::kBeauty);
+ fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 0); // no accept cut
+ fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 0); // no accept cut
+ if (TMath::Abs(mcpart->Eta()) < 0.9) {
+ fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
+ fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 1); // accept |eta|<0.9
+ }
+ if (TMath::Abs(GetRapidity(mcpart)) < 0.5) {
+ fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
+ fMCQA->GetDecayedKine(mcpart, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 2); // accept |y|<0.5
+ }
+ }
+ fMCQA->EndOfEventAna(AliHFEmcQA::kCharm);
+ fMCQA->EndOfEventAna(AliHFEmcQA::kBeauty);
- if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) continue;
- fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated);
- if(IsSignalElectron(mctrack)) fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCsignal);
- (dynamic_cast<TH1F *>(fQA->At(2)))->Fill(mctrack->Phi() - TMath::Pi());
- if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, mctrack)) continue;
- // find the label in the vector
- fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance);
- nElectrons++;
+ } // end of MC QA loop
+ // -----------------------------------------------------------------
+ fCFM->SetMCEventInfo(fMC);
+ // fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
+ } else {
+ fCFM->SetMCEventInfo(fRecEvent);
+ }
+ // Run MC loop
+ for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){
+ if(ProcessMCtrack(fMC->GetTrack(imc))) nElectrons++;
}
- (dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons);
// fCFM->CheckEventCuts(AliCFManager::kEvtRecCuts, fESD);
-
+ (dynamic_cast<TH1F *>(fQA->At(3)))->Fill(nElectrons);
+}
+//____________________________________________________________
+void AliAnalysisTaskHFE::ProcessESD(){
//
- // Loop ESD
+ // Run Analysis of reconstructed event in ESD Mode
+ // Loop over Tracks, filter according cut steps defined in AliHFEcuts
//
+ AliDebug(3, "Processing ESD Event");
+ AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fRecEvent);
+ if(!fESD){
+ AliError("ESD Event required for ESD Analysis")
+ return;
+ }
+ if (GetPlugin(kIsElecBackGround)) {
+ fElecBackGround->SetEvent(fESD);
+ }
+ if (GetPlugin(kSecVtx)) {
+ fSecVtx->SetEvent(fESD);
+ fSecVtx->GetPrimaryCondition();
+ }
+
+ if(HasMCData()){
+ if (GetPlugin(kSecVtx)) {
+ fSecVtx->SetStack(fMC->Stack());
+ }
+ if (GetPlugin(kIsElecBackGround)) {
+ fElecBackGround->SetMCEvent(fMC);
+ }
+ }
+
+
+ Double_t container[6];
+ memset(container, 0, sizeof(Double_t) * 6);
+ // container for the output THnSparse
+ Double_t dataE[5]; // [pT, eta, Phi, type, 'C' or 'B']
Int_t nElectronCandidates = 0;
- AliESDtrack *track = 0x0, *htrack = 0x0;
+ AliESDtrack *track = NULL, *htrack = NULL;
+ AliMCParticle *mctrack = NULL;
+ TParticle* mctrack4QA = NULL;
Int_t pid = 0;
// For double counted tracks
LabelContainer cont(fESD->GetNumberOfTracks());
Bool_t signal = kTRUE;
fCFM->SetRecEventInfo(fESD);
+ //
+ // Loop ESD
+ //
for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
track = fESD->GetTrack(itrack);
// RecKine: ITSTPC cuts
if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
- // Check if it is electrons near the vertex
- if(!(mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
- TParticle* mctrack4QA = fMC->Stack()->Particle(TMath::Abs(track->GetLabel()));
-
- container[3] = mctrack->Pt();
- container[4] = mctrack->Eta();
- container[5] = mctrack->Phi();
+ if(HasMCData()){
+ // Check if it is electrons near the vertex
+ if(!(mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
+ mctrack4QA = fMC->Stack()->Particle(TMath::Abs(track->GetLabel()));
+
+ container[3] = mctrack->Pt();
+ container[4] = mctrack->Eta();
+ container[5] = mctrack->Phi();
- if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
+ if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE;
+ }
if(signal) {
alreadyseen = cont.Find(TMath::Abs(track->GetLabel()));
(dynamic_cast<TH1F *>(fQA->At(1)))->Fill(track->GetAlpha()); // Check the acceptance without tight cuts
(dynamic_cast<TProfile *>(fQA->At(4)))->Fill(container[0], track->GetTRDpidQuality());
(dynamic_cast<TProfile *>(fQA->At(5)))->Fill(container[0], track->GetTRDncls());
+ //fQAcoll->Fill("fNtrdclusters", container[0], track->GetTRDncls());
}
// RecPrim
- if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) continue;
- if(signal) {
- fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepRecPrim + 2*(AliHFEcuts::kNcutESDSteps + 1));
- fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim);
- if(alreadyseen) {
- fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutESDSteps + 1);
- }
- }
-
+ if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track, container, signal, alreadyseen)) continue;
// HFEcuts: ITS layers cuts
- if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) continue;
- if(signal) {
- fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsITS + 2*(AliHFEcuts::kNcutESDSteps + 1));
- fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsITS);
- if(alreadyseen) {
- fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutESDSteps + 1);
- }
- }
+ if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track, container, signal, alreadyseen)) continue;
// HFEcuts: Nb of tracklets TRD0
- if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD, track)) continue;
+ if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTRD, track, container, signal, alreadyseen)) continue;
if(signal) {
- fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD + 2*(AliHFEcuts::kNcutESDSteps + 1));
- fCFM->GetParticleContainer()->Fill(&container[3], AliHFEcuts::kStepHFEcutsTRD);
- if(alreadyseen) {
- fCFM->GetParticleContainer()->Fill(&container[3], (AliHFEcuts::kStepHFEcutsTRD + (AliHFEcuts::kNcutESDSteps + 1)));
- }
// dimensions 3&4&5 : pt,eta,phi (MC)
((THnSparseF *)fCorrelation->At(0))->Fill(container);
}
- if (IsQAOn(kMCqa)) {
- // mc qa for after the reconstruction cuts
- AliDebug(2, "Running MC QA");
- fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 3); // charm
- fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 3); // beauty
+ if(HasMCData() && IsQAOn(kMCqa)) {
+ // mc qa for after the reconstruction cuts
+ AliDebug(2, "Running MC QA");
+ fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 3); // charm
+ fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 3); // beauty
}
// track accepted, do PID
if(!fPID->IsSelected(&hfetrack)) continue;
nElectronCandidates++;
- if (IsQAOn(kMCqa)) {
- // mc qa for after the reconstruction and pid cuts
- AliDebug(2, "Running MC QA");
- fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 4); // charm
- fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 4); // beauty
+ if (HasMCData() && IsQAOn(kMCqa)) {
+ // mc qa for after the reconstruction and pid cuts
+ AliDebug(2, "Running MC QA");
+ fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 4); // charm
+ fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 4); // beauty
}
// Fill Containers
((THnSparseF *)fCorrelation->At(1))->Fill(container);
}
- // Track selected: distinguish between true and fake
- AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack->Particle()->GetPdgCode()));
- if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){
- Int_t type = IsSignalElectron(track);
- AliDebug(1, Form("Type: %d\n", type));
- if(type){
- dataE[4] = type; // beauty[1] or charm[2]
- dataE[3] = 2; // signal electron
- }
- else{
- dataE[3] = 1; // not a signal electron
- dataE[4] = 0;
- }
- // pair analysis [mj]
- if (IsSecVtxOn()) {
- AliDebug(2, "Running Secondary Vertex Analysis");
- fSecVtx->InitAnaPair();
+ if(GetPlugin(kSecVtx)) {
+ AliDebug(2, "Running Secondary Vertex Analysis");
+ if(track->Pt()>0.5){
+ fSecVtx->InitHFEpairs();
+ fSecVtx->InitHFEsecvtxs();
+ AliESDtrack *htrack = 0x0;
for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){
htrack = fESD->GetTrack(jtrack);
- if ( itrack == jtrack ) continue;
- //if( fPID->IsSelected(htrack) && (itrack < jtrack)) continue;
- if( abs(fSecVtx->GetMCPID(track)) == 11 && (itrack < jtrack)) continue;
- fSecVtx->AnaPair(track, htrack, jtrack);
+ if ( itrack == jtrack ) continue; // since it is for tagging single electron, don't need additional condition
+ if (htrack->Pt()<0.5) continue;
+ if (!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC, htrack)) continue;
+ if (!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim, htrack)) continue;
+ fSecVtx->PairAnalysis(track, htrack, jtrack); // e-h pairing
}
- // based on the partner of e info, you run secandary vertexing function
- fSecVtx->RunSECVTX(track);
- } // end of pair analysis
- }
- else {
- // Fill THnSparse with the information for Fake Electrons
- dataE[3] = 0;
- dataE[4] = 0;
+ /*for(int ip=0; ip<fSecVtx->HFEpairs()->GetEntriesFast(); ip++){
+ AliHFEpairs *pair = (AliHFEpairs*) (fSecVtx->HFEpairs()->UncheckedAt(ip));
+ if(HasMCData()){
+ if(!(pair->GetPairCode()>1. && pair->GetPairCode()<4.)) // apply various cuts
+ fSecVtx->HFEpairs()->RemoveAt(ip);
+ }
+ }*/
+ fSecVtx->HFEpairs()->Compress();
+ fSecVtx->RunSECVTX(track); // secondary vertexing with e,h1,h2,.. tracks
+ for(int ip=0; ip<fSecVtx->HFEsecvtxs()->GetEntriesFast(); ip++){
+ AliHFEsecVtxs *secvtx=0x0;
+ secvtx = (AliHFEsecVtxs*) (fSecVtx->HFEsecvtxs()->UncheckedAt(ip));
+ // here you apply cuts, then if it doesn't pass the cut, remove it from the fSecVtx->HFEsecvtxs()
+ }
+ fSecVtx->DeleteHFEpairs();
+ fSecVtx->DeleteHFEsecvtxs();
+ }
}
- // fill the performance THnSparse, if the mc origin could be defined
- if(dataE[3] > -1){
- AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4]));
- fPIDperformance->Fill(dataE);
+
+ if(HasMCData()){
+ // Track selected: distinguish between true and fake
+ AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack->Particle()->GetPdgCode()));
+ if((pid = TMath::Abs(mctrack->Particle()->GetPdgCode())) == 11){
+ Int_t type = IsSignalElectron(track);
+ AliDebug(1, Form("Type: %d\n", type));
+ if(type){
+ dataE[4] = type; // beauty[1] or charm[2]
+ dataE[3] = 2; // signal electron
+ }
+ else{
+ dataE[3] = 1; // not a signal electron
+ dataE[4] = 0;
+ }
+ }
+ else {
+ // Fill THnSparse with the information for Fake Electrons
+ dataE[3] = 0;
+ dataE[4] = 0;
+ }
+ // fill the performance THnSparse, if the mc origin could be defined
+ if(dataE[3] > -1){
+ AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4]));
+ fPIDperformance->Fill(dataE);
+ }
}
+ // Electron background analysis
+ if (GetPlugin(kIsElecBackGround)) {
+
+ AliDebug(2, "Running BackGround Analysis");
+
+ for(Int_t jtrack = 0; jtrack < fESD->GetNumberOfTracks(); jtrack++){
+ htrack = fESD->GetTrack(jtrack);
+ if ( itrack == jtrack ) continue;
+ fElecBackGround->PairAnalysis(track, htrack);
+ }
+ } // end of electron background analysis
+
}
-
-
fNEvents->Fill(1);
+ //fQAcoll->Fill("fNevents", 1);
fNElectronTracksEvent->Fill(nElectronCandidates);
-
- // Done!!!
- PostData(0, fNEvents);
- PostData(1, fOutput);
- PostData(2, fQA);
}
//____________________________________________________________
-void AliAnalysisTaskHFE::Terminate(Option_t *){
+void AliAnalysisTaskHFE::ProcessAOD(){
//
- // Terminate not implemented at the moment
+ // Run Analysis in AOD Mode
+ // Function is still in development
//
- if(IsRunningPostProcess()){
- fOutput = dynamic_cast<TList *>(GetOutputData(1));
- if(!fOutput){
- AliError("Results not available");
- return;
+ AliDebug(3, "Processing AOD Event");
+ AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fRecEvent);
+ if(!fAOD){
+ AliError("AOD Event required for AOD Analysis")
+ return;
+ }
+
+ AliAODTrack *track = NULL;
+ AliAODMCParticle *mctrack = NULL;
+ Double_t container[6]; memset(container, 0, sizeof(Double_t) * 6);
+ Double_t dataE[5]; // [pT, eta, Phi, type, 'C' or 'B']
+ Int_t nElectronCandidates = 0;
+ Int_t pid;
+ for(Int_t itrack = 0; itrack < fAOD->GetNumberOfTracks(); itrack++){
+ track = fAOD->GetTrack(itrack);
+ if(!track) continue;
+ if(track->GetFlags() != 1<<4) continue; // Only process AOD tracks where the HFE is set
+
+ container[0] = track->Pt();
+ container[1] = track->Eta();
+ container[2] = track->Phi();
+
+ dataE[0] = track->Pt();
+ dataE[1] = track->Eta();
+ dataE[2] = track->Phi();
+ dataE[3] = -1;
+ dataE[4] = -1;
+
+ if(HasMCData()){
+ Int_t label = TMath::Abs(track->GetLabel());
+ if(label){
+ mctrack = dynamic_cast<AliAODMCParticle *>(fMC->GetTrack(label));
+ container[3] = mctrack->Pt();
+ container[4] = mctrack->Eta();
+ container[5] = mctrack->Phi();
+ }
+ }
+ // track accepted, do PID
+ AliHFEpidObject hfetrack;
+ hfetrack.fAnalysisType = AliHFEpidObject::kAODanalysis;
+ hfetrack.fRecTrack = track;
+ if(HasMCData()) hfetrack.fMCtrack = mctrack;
+ //if(!fPID->IsSelected(&hfetrack)) continue; // we will do PID here as soon as possible
+ // Particle identified - Fill CF Container
+ fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepHFEcutsTRD +1 + 2*(AliHFEcuts::kNcutESDSteps + 1));
+ nElectronCandidates++;
+ if(HasMCData()){
+ // Track selected: distinguish between true and fake
+ AliDebug(1, Form("Candidate Selected, filling THnSparse, PID: %d\n", mctrack->GetPdgCode()));
+ if((pid = TMath::Abs(mctrack->GetPdgCode())) == 11){
+ Int_t type = IsSignalElectron(track);
+ AliDebug(1, Form("Type: %d\n", type));
+ if(type){
+ dataE[4] = type; // beauty[1] or charm[2]
+ dataE[3] = 2; // signal electron
+ }
+ else{
+ dataE[3] = 1; // not a signal electron
+ dataE[4] = 0;
+ }
+ }
+ else {
+ // Fill THnSparse with the information for Fake Electrons
+ dataE[3] = 0;
+ dataE[4] = 0;
+ }
+ // fill the performance THnSparse, if the mc origin could be defined
+ if(dataE[3] > -1){
+ AliDebug(1, Form("Entries: [%.3f|%.3f|%.3f|%f|%f]\n", dataE[0],dataE[1],dataE[2],dataE[3],dataE[4]));
+ fPIDperformance->Fill(dataE);
+ }
}
- PostProcess();
}
+ fNEvents->Fill(1);
+ fNElectronTracksEvent->Fill(nElectronCandidates);
}
//____________________________________________________________
-void AliAnalysisTaskHFE::Load(TString filename){
+Bool_t AliAnalysisTaskHFE::ProcessMCtrack(AliVParticle *track){
//
- // Load Results into the task
+ // Filter the Monte Carlo Track
+ // Additionally Fill a THnSparse for Signal To Background Studies
+ // Works for AOD and MC analysis Type
//
- fQA = NULL; fOutput = NULL; fNEvents = NULL;
- TFile *input = TFile::Open(filename.Data());
- if(!input || input->IsZombie()){
- AliError("Cannot read file");
- return;
+ Double_t container[3], signalContainer[5];
+ Double_t vertex[3]; // Production vertex cut to mask gammas which are NOT supposed to have hits in the first ITS layer(s)
+ if(IsESDanalysis()){
+ AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(track);
+ container[0] = mctrack->Pt();
+ container[1] = mctrack->Eta();
+ container[2] = mctrack->Phi();
+
+ signalContainer[0] = mctrack->Pt();
+ signalContainer[1] = mctrack->Eta();
+ signalContainer[2] = mctrack->Phi();
+
+ vertex[0] = mctrack->Particle()->Vx();
+ vertex[1] = mctrack->Particle()->Vy();
+ } else {
+ AliAODMCParticle *aodmctrack = dynamic_cast<AliAODMCParticle *>(track);
+ container[0] = aodmctrack->Pt();
+ container[1] = aodmctrack->Eta();
+ container[2] = aodmctrack->Phi();
+
+ signalContainer[0] = aodmctrack->Pt();
+ signalContainer[1] = aodmctrack->Eta();
+ signalContainer[2] = aodmctrack->Phi();
+
+ aodmctrack->XvYvZv(vertex);
}
- TH1 *htmp = dynamic_cast<TH1I *>(input->Get("nEvents"));
- if(htmp)
- fNEvents = dynamic_cast<TH1I *>(htmp->Clone());
- else
- AliError("Event Counter histogram not found");
- TList *ltmp = dynamic_cast<TList *>(input->Get("Results"));
- if(ltmp)
- fOutput = dynamic_cast<TList *>(ltmp->Clone());
- else
- AliError("Output Histograms not found");
- ltmp = dynamic_cast<TList *>(input->Get("QA"));
- if(ltmp)
- fQA = dynamic_cast<TList *>(ltmp->Clone());
- else
- AliError("QA histograms not found");
- input->Close();
- delete input;
- Int_t nObjects = 0;
- if(fNEvents) nObjects++;
- if(fOutput) nObjects++;
- if(fQA) nObjects++;
- AliInfo(Form("Loaded %d Objects into task", nObjects));
+
+ if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, track)) return kFALSE;
+ fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCGenerated);
+ if((signalContainer[3] = static_cast<Double_t >(IsSignalElectron(track))) > 1e-3) fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCsignal);
+ signalContainer[4] = 0;
+ // apply cut on the sqrt of the production vertex
+ Double_t radVertex = TMath::Sqrt(vertex[0]*vertex[0] + vertex[1] * vertex[1]);
+ if(radVertex < 3.5){
+ // Within first ITS layer(2) -> Background we cannot reject by ITS cut, let it pass
+ signalContainer[4] = 1;
+ } else if (radVertex < 7.5){
+ signalContainer[4] = 2;
+ }
+ fSignalToBackgroundMC->Fill(signalContainer);
+ (dynamic_cast<TH1F *>(fQA->At(2)))->Fill(container[2] - TMath::Pi());
+ //if(IsESDanalysis()){
+ if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCInAcceptance, track)) return kFALSE;
+ fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepMCInAcceptance);
+ //}
+ return kTRUE;
}
//____________________________________________________________
//delete electronPurity; delete signalPurity; delete fakeContamination;
}
+//____________________________________________________________
+void AliAnalysisTaskHFE::DrawMCSignal2Background(){
+ //
+ // Draw the MC signal/background plots
+ //
+ fSignalToBackgroundMC = dynamic_cast<THnSparseF *>(fOutput->FindObject("SignalToBackgroundMC"));
+ if(!fSignalToBackgroundMC) return;
+
+ // First Select everything within the first ITS Layer
+ fSignalToBackgroundMC->GetAxis(4)->SetRange(2,2);
+ // For Signal electrons we project axis 3 to everything > 0
+ fSignalToBackgroundMC->GetAxis(3)->SetRange(2,3);
+ TH1 *hSignal = fSignalToBackgroundMC->Projection(0);
+ hSignal->SetName("hSignal");
+ hSignal->SetTitle("Signal Electrons");
+ // For Background studies project axis 3 to bin 0
+ fSignalToBackgroundMC->GetAxis(3)->SetRange(1,1);
+ TH1 *hBackground = fSignalToBackgroundMC->Projection(0);
+ hBackground->SetName("hBackground");
+ hBackground->SetTitle("Background Electrons");
+ // For All electrons we undo the projection of axis 3
+ fSignalToBackgroundMC->GetAxis(3)->SetRange(0, fSignalToBackgroundMC->GetAxis(3)->GetNbins());
+ TH1 *hAll = fSignalToBackgroundMC->Projection(0);
+ hAll->SetName("hAll");
+ hAll->SetTitle("All Electrons");
+ // Prepare Canvas
+ TCanvas *cEff = new TCanvas("cEff", "MC Sig/Backg studies", 800, 400);
+ cEff->Divide(2);
+ // Plot Efficiency
+ TH1 *hEff = (TH1 *)hSignal->Clone();
+ hEff->Divide(hAll);
+ hEff->SetName("sigEff");
+ hEff->SetTitle("Signal/(Signal+Background)");
+ hEff->GetXaxis()->SetTitle("p_{T} / GeV/c");
+ hEff->GetYaxis()->SetTitle("Efficiency");
+ hEff->SetStats(kFALSE);
+ hEff->SetMarkerStyle(22);
+ hEff->SetMarkerColor(kBlue);
+ cEff->cd(1);
+ hEff->Draw("p");
+ // Plot Signal/Background
+ TH1 *hSB = (TH1 *)hSignal->Clone();
+ hSB->Divide(hBackground);
+ hSB->SetName("sigEff");
+ hSB->SetTitle("Signal/Background");
+ hSB->GetXaxis()->SetTitle("p_{T} / GeV/c");
+ hSB->GetYaxis()->SetTitle("Signal/Background");
+ hSB->SetStats(kFALSE);
+ hSB->SetMarkerStyle(22);
+ hSB->SetMarkerColor(kBlue);
+ cEff->cd(2);
+ hSB->Draw("p");
+
+ // Undo projections
+ fSignalToBackgroundMC->GetAxis(4)->SetRange(0, fSignalToBackgroundMC->GetAxis(4)->GetNbins());
+}
+
//____________________________________________________________
void AliAnalysisTaskHFE::MakeParticleContainer(){
//
for(Int_t i=0; i<=nBin[3]; i++) binEdges2[3][i] = i;
for(Int_t i=0; i<=nBin[4]; i++) binEdges2[4][i] = i;
- fPIDperformance = new THnSparseF("PIDperformance", "PID performance; pT [GeV/c]; theta [rad]; phi [rad] type (0 - not el, 1 - other el, 2 - HF el; flavor (0 - no, 1 - charm, 2 - bottom)", nDim, nBin);
- for(Int_t idim = 0; idim < nDim; idim++)
+ fPIDperformance = new THnSparseF("PIDperformance", "PID performance; pT [GeV/c]; theta [rad]; phi [rad]; type (0 - not el, 1 - other el, 2 - HF el; flavor (0 - no, 1 - charm, 2 - bottom)", nDim, nBin);
+ fSignalToBackgroundMC = new THnSparseF("SignalToBackgroundMC", "PID performance; pT [GeV/c]; theta [rad]; phi [rad]; flavor (0 - no, 1 - charm, 2 - bottom); ITS Cluster (0 - no, 1 - first (and maybe second), 2 - second)", nDim, nBin);
+ for(Int_t idim = 0; idim < nDim; idim++){
fPIDperformance->SetBinEdges(idim, binEdges2[idim]);
+ fSignalToBackgroundMC->SetBinEdges(idim, binEdges2[idim]);
+ }
}
//____________________________________________________________
// Print Analysis status
//
printf("\n\tAnalysis Settings\n\t========================================\n\n");
- printf("\tSecondary Vertex finding: %s\n", IsSecVtxOn() ? "YES" : "NO");
- printf("\tPrimary Vertex resolution: %s\n", IsPriVtxOn() ? "YES" : "NO");
+ printf("\tSecondary Vertex finding: %s\n", GetPlugin(kSecVtx) ? "YES" : "NO");
+ printf("\tPrimary Vertex resolution: %s\n", GetPlugin(kPriVtx) ? "YES" : "NO");
printf("\n");
printf("\tParticle Identification Detectors:\n");
TObjArray *detectors = fPIDdetectors.Tokenize(":");
kCharm = 1,
kBeauty = 2
};
- AliMCParticle *mctrack = NULL;
TString objname = fTrack->IsA()->GetName();
- if(!objname.CompareTo("AliESDtrack")){
- AliDebug(2, "Checking signal for ESD track");
- AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(fTrack);
- mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(esdtrack->GetLabel())));
- }
- else if(!objname.CompareTo("AliAODtrack")){
- AliError("AOD Analysis not implemented yet");
- return kNoSignal;
- }
- else if(!objname.CompareTo("AliMCParticle")){
- AliDebug(2, "Checking signal for MC track");
- mctrack = dynamic_cast<AliMCParticle *>(fTrack);
- }
- else{
- AliError("Input object not supported");
- return kNoSignal;
+ Int_t pid = 0;
+ if(IsESDanalysis()){
+ // ESD Analysis
+ AliMCParticle *mctrack = NULL;
+ if(!objname.CompareTo("AliESDtrack")){
+ AliDebug(2, "Checking signal for ESD track");
+ AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(fTrack);
+ mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(esdtrack->GetLabel())));
+ }
+ else if(!objname.CompareTo("AliMCParticle")){
+ AliDebug(2, "Checking signal for MC track");
+ mctrack = dynamic_cast<AliMCParticle *>(fTrack);
+ }
+ else{
+ AliError("Input object not supported");
+ return kNoSignal;
+ }
+ if(!mctrack) return kNoSignal;
+ TParticle *ecand = mctrack->Particle();
+ if(TMath::Abs(ecand->GetPdgCode()) != 11) return kNoSignal; // electron candidate not true electron
+ Int_t motherLabel = TMath::Abs(ecand->GetFirstMother());
+ AliDebug(3, Form("mother label: %d\n", motherLabel));
+ if(!motherLabel) return kNoSignal; // mother track unknown
+ AliMCParticle *motherTrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(motherLabel));
+ if(!motherTrack) return kNoSignal;
+ TParticle *mparticle = motherTrack->Particle();
+ pid = TMath::Abs(mparticle->GetPdgCode());
+ } else {
+ // AOD Analysis - Different Data handling
+ AliAODMCParticle *aodmc = NULL;
+ if(!objname.CompareTo("AliAODTrack")){
+ AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(fTrack);
+ Int_t aodlabel = TMath::Abs(aodtrack->GetLabel());
+ if(aodlabel >= fMC->GetNumberOfTracks()) return kNoSignal;
+ aodmc = dynamic_cast<AliAODMCParticle *>(fMC->GetTrack(aodlabel));
+ } else if(!objname.CompareTo("AliAODMCParticle")){
+ aodmc = dynamic_cast<AliAODMCParticle *>(fTrack);
+ } else{
+ AliError("Input object not supported");
+ return kNoSignal;
+ }
+ if(!aodmc) return kNoSignal;
+ Int_t motherLabel = TMath::Abs(aodmc->GetMother());
+ AliDebug(3, Form("mother label: %d\n", motherLabel));
+ if(!motherLabel || motherLabel >= fMC->GetNumberOfTracks()) return kNoSignal;
+ AliAODMCParticle *aodmother = dynamic_cast<AliAODMCParticle *>(fMC->GetTrack(motherLabel));
+ pid = aodmother->GetPdgCode();
}
- if(!mctrack) return kNoSignal;
- TParticle *ecand = mctrack->Particle();
- if(TMath::Abs(ecand->GetPdgCode()) != 11) return kNoSignal; // electron candidate not true electron
- Int_t motherLabel = TMath::Abs(ecand->GetFirstMother());
- AliDebug(3, Form("mother label: %d\n", motherLabel));
- if(!motherLabel) return kNoSignal; // mother track unknown
- AliMCParticle *motherTrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(motherLabel));
- if(!motherTrack) return kNoSignal;
- TParticle *mparticle = motherTrack->Particle();
- Int_t pid = TMath::Abs(mparticle->GetPdgCode());
+ // From here the two analysis modes go together
AliDebug(3, Form("PDG code: %d\n", pid));
// identify signal according to Pdg Code
return rapidity;
}
+//__________________________________________
+void AliAnalysisTaskHFE::SwitchOnPlugin(Int_t plug){
+ //
+ // Switch on Plugin
+ // Available:
+ // - Primary vertex studies
+ // - Secondary vertex Studies
+ // - Post Processing
+ //
+ switch(plug){
+ case kPriVtx: SETBIT(fPlugins, plug); break;
+ case kSecVtx: SETBIT(fPlugins, plug); break;
+ case kIsElecBackGround: SETBIT(fPlugins, plug); break;
+ case kPostProcess: SETBIT(fPlugins, plug); break;
+ default: AliError("Unknown Plugin");
+ };
+}
+
+//__________________________________________
+Bool_t AliAnalysisTaskHFE::ProcessCutStep(Int_t cutStep, AliVParticle *track, Double_t *container, Bool_t signal, Bool_t alreadyseen){
+ //
+ // Check single track cuts for a given cut step
+ // Fill the particle container
+ //
+ if(!fCFM->CheckParticleCuts(cutStep, track)) return kFALSE;
+ if(signal) {
+ fCFM->GetParticleContainer()->Fill(container, cutStep + 2*(AliHFEcuts::kNcutESDSteps + 1));
+ fCFM->GetParticleContainer()->Fill(&container[3], cutStep);
+ if(alreadyseen) {
+ fCFM->GetParticleContainer()->Fill(&container[3], cutStep + AliHFEcuts::kNcutESDSteps + 1);
+ }
+ }
+ return kTRUE;
+}
class AliHFEcuts;
class AliHFEmcQA;
class AliHFEsecVtx;
+class AliHFEelecbackground;
+class AliHFEcollection;
class AliCFManager;
-class AliESDEvent;
-class AliESDtrackCuts;
+class AliVEvent;
class AliMCEvent;
class AliVParticle;
class TH1I;
class AliAnalysisTaskHFE : public AliAnalysisTask{
public:
- enum{
- kPIDqa = 0,
- kMCqa =1
- };
+ enum{
+ kPIDqa = 0,
+ kMCqa =1
+ };
+ enum{
+ kPriVtx = 0,
+ kSecVtx = 1,
+ kIsElecBackGround = 2,
+ kPostProcess = 3
+ };
AliAnalysisTaskHFE();
AliAnalysisTaskHFE(const char * name);
AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref);
virtual void Terminate(Option_t *);
Bool_t IsQAOn(Int_t qaLevel) const { return TESTBIT(fQAlevel, qaLevel); };
- Bool_t IsSecVtxOn() const { return TestBit(kIsSecVtxOn); };
- Bool_t IsPriVtxOn() const { return TestBit(kIsPriVtxOn); };
- Bool_t IsRunningPostProcess() const { return TestBit(kIsRunningPostProcess); };
+ Bool_t IsAODanalysis() const { return TestBit(kAODanalysis); };
+ Bool_t IsESDanalysis() const { return !TestBit(kAODanalysis); };
Bool_t HasMCData() const { return TestBit(kHasMCdata); }
+ Bool_t GetPlugin(Int_t plug) const { return TESTBIT(fPlugins, plug); };
Int_t IsSignalElectron(AliVParticle *fTrack) const;
void Load(TString filename = "HFEtask.root");
void PostProcess();
void SetHFECuts(AliHFEcuts * const cuts) { fCuts = cuts; };
void SetQAOn(Int_t qaLevel) { SETBIT(fQAlevel, qaLevel); };
+ void SwitchOnPlugin(Int_t plug);
void SetHasMCData(Bool_t hasMC = kTRUE) { SetBit(kHasMCdata, hasMC); };
- void SetPriVtxOn(Bool_t option = kTRUE) { SetBit(kIsPriVtxOn, option); };
- void SetSecVtxOn(Bool_t option = kTRUE) { SetBit(kIsSecVtxOn, option); };
- void SetRunPostProcess(Bool_t option = kTRUE) { SetBit(kIsRunningPostProcess, option); };
void SetPIDdetectors(Char_t * const detectors){ fPIDdetectors = detectors; }
void SetPIDStrategy(UInt_t strategy) { fPIDstrategy = strategy; }
void AddPIDdetector(Char_t *detector);
+ void SetAODAnalysis() { SetBit(kAODanalysis, kTRUE); };
+ void SetESDAnalysis() { SetBit(kAODanalysis, kFALSE); };
void PrintStatus() const;
Float_t GetRapidity(TParticle *part) const;
private:
enum{
- kIsSecVtxOn = BIT(19),
- kIsPriVtxOn = BIT(20),
- kIsRunningPostProcess = BIT(21),
- kHasMCdata = BIT(22)
+ kHasMCdata = BIT(19),
+ kAODanalysis = BIT(20)
};
class LabelContainer{
public:
Int_t *fCurrent; // Current entry to mimic an iterator
};
void MakeParticleContainer();
+ void ProcessMC();
+ void ProcessESD();
+ void ProcessAOD();
+ Bool_t ProcessMCtrack(AliVParticle *track);
+ Bool_t ProcessCutStep(Int_t cutStep, AliVParticle *track, Double_t *container, Bool_t signal, Bool_t alreadyseen);
+ void DrawMCSignal2Background();
ULong_t fQAlevel; // QA level
TString fPIDdetectors; // Detectors for Particle Identification
UInt_t fPIDstrategy; // PID Strategy
- AliESDEvent *fESD; //! The ESD Event
+ UShort_t fPlugins; // Enabled Plugins
+ AliVEvent *fRecEvent; //! The ESD Event
AliMCEvent *fMC; //! The MC Event
AliCFManager *fCFM; //! Correction Framework Manager
TList *fCorrelation; //! response matrix for unfolding
THnSparseF *fPIDperformance; //! info on contamination and yield of electron spectra
+ THnSparseF *fSignalToBackgroundMC; //! Signal To Background Studies on pure MC information
AliHFEpid *fPID; //! PID
AliHFEcuts *fCuts; // Cut Collection
AliHFEsecVtx *fSecVtx; //! Secondary Vertex Analysis
+ AliHFEelecbackground *fElecBackGround;//! Background analysis
AliHFEmcQA *fMCQA; //! MC QA
TH1I *fNEvents; //! counter for the number of Events
TH1I *fNElectronTracksEvent; //! Number of Electron candidates after PID decision per Event
TList *fOutput; //! Container for Task Output
TList *fHistMCQA; //! Output container for MC QA histograms
TList *fHistSECVTX; //! Output container for sec. vertexing results
+ TList *fHistELECBACKGROUND; //! Output container for electron background analysis
+// AliHFEcollection *fQAcoll; //! collection class for basic QA histograms
ClassDef(AliAnalysisTaskHFE, 1) // The electron Analysis Task
};
#include <TH1F.h>
#include <TH2F.h>
+#include <THnSparse.h>
+#include <TProfile.h>
#include <TList.h>
#include <TString.h>
#include <TBrowser.h>
-#include <TIterator.h>
+#include <TMath.h>
#include "AliLog.h"
#include "AliHFEcollection.h"
//___________________________________________________________________
AliHFEcollection::AliHFEcollection():
TNamed()
- , fListE(0x0)
+ , fList(0x0)
{
+
//
// default constructor
//
- fListE = new TList();
- if(!fListE){
+ fList = new TList();
+ if(!fList){
AliError("Initialization of the list failed");
}
else{
// list is owner of the objects. Once list is deleted, the objects
// it contains will be deleted too
- fListE->SetOwner(kTRUE);
+ fList->SetOwner(kTRUE);
}
//Printf("%s:%d,%p",(char*)__FILE__,__LINE__,fInstance);
//___________________________________________________________________
AliHFEcollection::AliHFEcollection(char* name, char* title):
TNamed(name, title)
- , fListE(0x0)
+ , fList(0x0)
{
//
// constructor
//
- fListE = new TList();
- if(!fListE){
+ fList = new TList();
+ if(!fList){
AliError("Initialization of the list failed");
}
else{
// list is owner of the objects. Once list is deleted, the objects
// it contains will be deleted too
- fListE->SetOwner(kTRUE);
+ fList->SetOwner(kTRUE);
}
}
//___________________________________________________________________
AliHFEcollection::AliHFEcollection(const AliHFEcollection &c) :
TNamed(c)
- , fListE(0x0)
+ , fList(0x0)
{
//
}
//___________________________________________________________________
void AliHFEcollection::Copy(TObject &ref) const {
+
//
// Performs the copying of the object
//
+
AliHFEcollection &target = dynamic_cast<AliHFEcollection &>(ref);
- target.fListE = fListE;
+ target.fList = fList;
}
//___________________________________________________________________
AliHFEcollection::~AliHFEcollection(){
+
//
// Destructor
//
+
AliInfo("DESTRUCTOR");
}
//___________________________________________________________________
Bool_t AliHFEcollection::CreateTH1F(const char* name, const char* title, Int_t nBin, Float_t nMin, Float_t nMax){
+
//
// Creates a TH1F histogram for the collection
//
- if(!fListE){
+
+ if(!fList){
AliError("No TList pointer ! ");
return kFALSE;
}
else{
- fListE->Add(new TH1F(name, title, nBin, nMin, nMax));
+ fList->Add(new TH1F(name, title, nBin, nMin, nMax));
return CheckObject(name);
}
}
//___________________________________________________________________
Bool_t AliHFEcollection::CreateTH2F(const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY){
+
//
// Creates a TH2F histogram for the collection
//
- if(!fListE){
+
+ if(!fList){
AliError("No TList pointer ! ");
return kFALSE;
}
- fListE->Add(new TH2F(name, title, nBinX, nMinX, nMaxX, nBinY, nMinY, nMaxY));
+ fList->Add(new TH2F(name, title, nBinX, nMinX, nMaxX, nBinY, nMinY, nMaxY));
return CheckObject(name);
}
//___________________________________________________________________
Bool_t AliHFEcollection::CreateTH1Fvector1(Int_t X, const char* name, const char* title, Int_t nBin, Float_t nMin, Float_t nMax){
+
//
// create a 1 dimensional array of size [X]
//
- if(!fListE){
+
+ if(!fList){
AliError("No TList pointer ! ");
return kFALSE;
}
}
//___________________________________________________________________
Bool_t AliHFEcollection::CreateTH2Fvector1(Int_t X, const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY){
+
//
// create a 1 dimensinal array of TH2F histograms with size [X]
//
- if(!fListE){
+
+ if(!fList){
AliError("No TList pointer !");
return kFALSE;
}
}
//___________________________________________________________________
Bool_t AliHFEcollection::CreateTH1Fvector2(Int_t X, Int_t Y, const char* name, const char* title, Int_t nBin, Float_t nMin, Float_t nMax){
+
//
// create a 2 dimensional array of histograms of size [X, Y]
//
- if(!fListE){
+
+ if(!fList){
AliError("No TList pointer ! ");
return kFALSE;
}
}
}
}
- return kTRUE;
+ return kTRUE;
+}
+//___________________________________________________________________
+Bool_t AliHFEcollection::CreateProfile(const char* name, const char* title, Int_t nbins, Double_t xmin, Double_t xmax){
+
+ //
+ // create a simple TProfile
+ //
+
+ if(!fList){
+ AliError("No TList pointer ! ");
+ return kFALSE;
+ }
+ fList->Add(new TProfile(name, title, nbins, xmin, xmax));
+ return CheckObject(name);
+
+}
+//___________________________________________________________________
+Bool_t AliHFEcollection::CreateTHnSparse(const char* name, const char* title, Int_t dim, Int_t* nbins, Double_t* xmin, Double_t* xmax){
+
+ //
+ // create 'dim' dimensional THnSparse
+ //
+
+ if(!fList){
+ AliError("No TList pointer ! ");
+ return kFALSE;
+ }
+ fList->Add(new THnSparseF(name, title, dim, nbins, xmin, xmax));
+ return CheckObject(name);
+
+}
+//___________________________________________________________________
+TObject* AliHFEcollection::Get(const char* name){
+
+ //
+ // Get histogram with the required name
+ //
+
+ if(!CheckObject(name)){
+ AliError(Form("Not possible to return pointer to the object '%s'\n", name));
+ return 0;
+ }
+
+ return fList->FindObject(name);
}
//___________________________________________________________________
-TObject* AliHFEcollection::Get(const char* name, Int_t X){
+Bool_t AliHFEcollection::Fill(const char* name, Double_t v){
+
//
- // Get the histogram from the vector
- // Paramter:
- // name - vector name
- // X - Number of the desired histogram
+ // fill function for one TH1 histograms
//
- TString hname = name;
- hname.Append(Form("_[%d]", X));
- if(!CheckObject(hname.Data())){
- AliError("No such object found in the list");
- return 0x0;
+
+ if(!CheckObject(name)){
+ AliError(Form("Not possible to return pointer to the object '%s'\n", name));
+ return kFALSE;
}
- else{
- return Get(hname.Data());
+
+ // chack the possible object types
+ if(fList->FindObject(name)->InheritsFrom("TH1")){
+ (dynamic_cast<TH1F*>(fList->FindObject(name)))->Fill(v);
+ return kTRUE;
+ }
+
+ return kFALSE;
+
+}
+//___________________________________________________________________
+Bool_t AliHFEcollection::Fill(const char* name, Int_t X, Double_t v){
+
+ //
+ // fill function for one dimension arrays of TH1
+ //
+
+ const char* n = Form("%s_[%d]", name, X);
+ TObject *o = Get(n);
+ if(!o){
+ return kFALSE;
}
+ Fill(o->GetName(), v);
+ return kTRUE;
}
//___________________________________________________________________
-TObject* AliHFEcollection::Get(const char* name, Int_t X, Int_t Y){
+Bool_t AliHFEcollection::Fill(const char* name, Int_t X, Int_t Y, Double_t v){
+
+ const char* n = Form("%s_[%d][%d]", name, X, Y);
+ TObject *o = Get(n);
+ if(!o){
+ return kFALSE;
+ }
+ Fill(o->GetName(), v);
+ return kTRUE;
+}
+//___________________________________________________________________
+Bool_t AliHFEcollection::Fill(const char* name, Int_t X, Double_t v1, Double_t v2){
+
//
- // Get histogram from the 2D vector
- // Parameters:
- // name - Name of the vector
- // X,Y - Indices of the histogram
+ // fill function for one dimension array of TH2
//
- TString hname = name;
- hname.Append(Form("_[%d][%d]", X, Y));
- if(!CheckObject(hname.Data())){
- AliError("No such object found in the list");
- AliError(Form("name: %s", hname.Data()));
- return 0x0;
+
+ const char* n = Form("%s_[%d]", name, X);
+ TObject *o = Get(n);
+ if(!o){
+ return kFALSE;
}
- else{
- return Get(hname.Data());
+ Fill(o->GetName(), v1, v2);
+
+ return kTRUE;
+}
+//___________________________________________________________________
+Bool_t AliHFEcollection::Fill(const char* name, Double_t v1, Double_t v2){
+
+ //
+ // fill function for TH2 objects
+ //
+
+ if(!CheckObject(name)){
+ AliError(Form("Not possible to return pointer to the object '%s'\n", name));
+ return kFALSE;
}
+
+ // chack the possible object types
+ if(fList->FindObject(name)->InheritsFrom("TH2")){
+ (dynamic_cast<TH2F*>(fList->FindObject(name)))->Fill(v1, v2);
+ return kTRUE;
+ }
+ if(fList->FindObject(name)->InheritsFrom("TProfile")){
+ (dynamic_cast<TProfile*>(fList->FindObject(name)))->Fill(v1, v2);
+ return kTRUE;
+ }
+
+ return kFALSE;
+
}
+
//___________________________________________________________________
Bool_t AliHFEcollection::CheckObject(const char* name){
+
//
// check wheter the creation of the histogram was succesfull
//
- if(!fListE){
+ if(!fList){
AliError("No TList pointer ! ");
return kFALSE;
}
- if(!fListE->FindObject(name)){
- AliError("Creating or Finding the object failed");
+ if(!fList->FindObject(name)){
+ AliError(Form("Creating or Finding the object '%s' failed\n", name));
return kFALSE;
}
return kTRUE;
}
//___________________________________________________________________
-TObject* AliHFEcollection::Get(const char* name){
- //
- // Get histogram with the required name
+Bool_t AliHFEcollection::BinLogAxis(const char* name, Int_t dim){
+
//
- return fListE->FindObject(name);
+ // converts the axis (defined by the dimension) of THx or THnSparse
+ // object to Log scale. Number of bins and bin min and bin max are preserved
+ //
+
+
+ if(!CheckObject(name)){
+ return kFALSE;
+ }
+
+ TObject *o = Get(name);
+ TAxis *axis = 0x0;
+ if(o->InheritsFrom("TH1")){
+ axis = (dynamic_cast<TH1F*>(o))->GetXaxis();
+ }
+ if(o->InheritsFrom("TH2")){
+ if(0 == dim){
+ axis = (dynamic_cast<TH2F*>(o))->GetXaxis();
+ }
+ else if(1 == dim){
+ axis = (dynamic_cast<TH2F*>(o))->GetYaxis();
+ }
+ else{
+ AliError("Only dim = 0 or 1 possible for TH2F");
+ }
+ }
+ if(o->InheritsFrom("THnSparse")){
+ axis = (dynamic_cast<THnSparse*>(o))->GetAxis(dim);
+ }
+
+ if(!axis){
+ AliError(Form("Axis '%d' could not be identified in the object '%s'\n", dim, name));
+ return kFALSE;
+ }
+
+ Int_t bins = axis->GetNbins();
+
+ Double_t from = axis->GetXmin();
+ Double_t to = axis->GetXmax();
+ Double_t *newBins = new Double_t[bins+1];
+ newBins[0] = from;
+ Double_t factor = TMath::Power(to/from, 1./bins);
+ for(Int_t i=1; i<=bins; ++i){
+ newBins[i] = factor * newBins[i-1];
+ }
+ axis->Set(bins, newBins);
+ delete newBins;
+
+ return kTRUE;
+
+
}
//___________________________________________________________________
Long64_t AliHFEcollection::Merge(TCollection *list){
+
//
// Merge the collections
//
- if(!fListE){
+
+ if(!fList){
AliError("AliHFEcollection::Merge : No TList pointer ! ");
return 0;
}
- return fListE->Merge(list);
+ return fList->Merge(list);
}
//____________________________________________________________________
void AliHFEcollection::Browse(TBrowser *b)
{
+
//
// Browse the content of the directory.
//
if (b) {
TObject *obj = 0;
- TIter nextin(fListE);
+ TIter nextin(fList);
//Add objects that are only in memory
while ((obj = nextin())) {
Bool_t CreateTH2F(const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY);
Bool_t CreateTH1Fvector1(Int_t X, const char* name, const char* title, Int_t nBin, Float_t nMin, Float_t nMax);
+ Bool_t CreateTH1Fvector2(Int_t X, Int_t Y, const char* name, const char* title, Int_t nBin, Float_t nMin, Float_t nMax);
Bool_t CreateTH2Fvector1(Int_t X, const char* name, const char* title, Int_t nBinX, Float_t nMinX, Float_t nMaxX, Int_t nBinY, Float_t nMinY, Float_t nMaxY);
+ Bool_t CreateProfile(const char* name, const char* title, Int_t nbins, Double_t xmin, Double_t xmax);
+ Bool_t CreateTHnSparse(const char* name, const char* title, Int_t dim, Int_t* nbins, Double_t* xmin, Double_t* xmax);
- Bool_t CreateTH1Fvector2(Int_t X, Int_t Y, const char* name, const char* title, Int_t nBin, Float_t nMin, Float_t nMax);
-
+ Bool_t BinLogAxis(const char* name, Int_t dim);
+
Long64_t Merge(TCollection *list);
// Get functions
- TList* GetList() const { return fListE; }
+ TList* GetList() const { return fList; }
TObject* Get(const char* name);
- TObject* Get(const char* name, Int_t X);
- TObject* Get(const char* name, Int_t X, Int_t Y);
+ // Fill functions
+ Bool_t Fill(const char* name, Double_t v);
+ Bool_t Fill(const char* name, Int_t X, Double_t v);
+ Bool_t Fill(const char* name, Int_t X, Int_t Y, Double_t v);
+ Bool_t Fill(const char* name, Double_t v1, Double_t v2);
+ Bool_t Fill(const char* name, Int_t X, Double_t v1, Double_t v2);
private:
Bool_t CheckObject(const char* name);
void Copy(TObject &ref) const;
private:
- TList* fListE; //! Object container
+ TList* fList; //! Object container
ClassDef(AliHFEcollection, 1)
SetHFElectronITSCuts();
SetHFElectronTRDCuts();
+ // Publish to the cuts which analysis type they are (ESD Analysis by default)
+ if(IsAOD()){
+ AliInfo("Setting AOD Analysis");
+ TObjArray *genCuts = dynamic_cast<TObjArray *>(fCutList->FindObject("fPartGenCuts"));
+ if(genCuts){
+ AliCFParticleGenCuts *myGenCut = dynamic_cast<AliCFParticleGenCuts *>(genCuts->FindObject("fCutsGenMC"));
+ if(myGenCut) myGenCut->SetAODMC(kTRUE);
+ }
+ }
+
// Connect the event cuts
- SetEventCutList(kEventStepGenerated);
+ /*SetEventCutList(kEventStepGenerated);
SetEventCutList(kEventStepReconstructed);
cfm->SetEventCutsList(kEventStepGenerated, dynamic_cast<TObjArray *>(fCutList->FindObject("fEvGenCuts")));
- cfm->SetEventCutsList(kEventStepReconstructed, dynamic_cast<TObjArray *>(fCutList->FindObject("fEvRecCuts")));
+ cfm->SetEventCutsList(kEventStepReconstructed, dynamic_cast<TObjArray *>(fCutList->FindObject("fEvRecCuts")));*/
// Connect the particle cuts
cfm->SetParticleCutsList(kStepMCGenerated, dynamic_cast<TObjArray *>(fCutList->FindObject("fPartGenCuts")));
TList *GetQAhistograms() const { return fHistQA; }
void SetDebugMode();
+ void SetAOD() { SetBit(kAOD, kTRUE); }
+ void SetESD() { SetBit(kAOD, kFALSE); }
void UnsetDebugMode() { SetBit(kDebugMode, kFALSE); }
Bool_t IsInDebugMode() const { return TestBit(kDebugMode); };
+ Bool_t IsAOD() const { return TestBit(kAOD); }
+ Bool_t IsESD() const { return !TestBit(kAOD); }
// Getters
Bool_t IsRequireITSpixel() const { return TESTBIT(fRequirements, kITSPixel); };
private:
enum{
- kDebugMode = BIT(14)
+ kDebugMode = BIT(14),
+ kAOD = BIT(15)
};
typedef enum{
kPrimary = 0,
#include <AliLog.h>
#include <AliStack.h>
+#include <AliGenEventHeader.h>
#include <AliAODMCParticle.h>
#include "AliHFEmcQA.h"
//_______________________________________________________________________________________________
AliHFEmcQA::AliHFEmcQA() :
fStack(0x0)
+ ,fMCHeader(0x0)
,fMCArray(0x0)
,fNparents(0)
{
AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p):
TObject(p)
,fStack(0x0)
+ ,fMCHeader(0x0)
,fMCArray(0x0)
,fNparents(p.fNparents)
{
Int_t maPdgcode = partMother->GetPdgCode();
Int_t maPdgcodeCopy = maPdgcode;
+ // get mc primary vertex
+ TArrayF mcPrimVtx;
+ if(fMCHeader) fMCHeader->PrimaryVertex(mcPrimVtx);
+
// get electron production vertex
TLorentzVector ePoint;
mcpart->ProductionVertex(ePoint);
+ // calculated production vertex to primary vertex (in xy plane)
+ Float_t decayLxy = TMath::Sqrt((mcPrimVtx[0]-ePoint[0])*(mcPrimVtx[0]-ePoint[0])+(mcPrimVtx[1]-ePoint[1])*(mcPrimVtx[1]-ePoint[1]));
+
// if the mother is charmed hadron
Bool_t isMotherDirectCharm = kFALSE;
if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
// ratio between pT of electron and pT of mother B hadron
if(grandMa->Pt()) fHistComm[iq][icut].fDePtRatio->Fill(grandMa->Pt(),mcpart->Pt()/grandMa->Pt());
- // distance between electron production point and mother hadron production point
- TLorentzVector debPoint;
- grandMa->ProductionVertex(debPoint);
- TLorentzVector dedistance = ePoint - debPoint;
- fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),dedistance.P());
+ // distance between electron production point and primary vertex
+ fHistComm[iq][icut].fDeDistance->Fill(grandMa->Pt(),decayLxy);
return;
}
}
// ratio between pT of electron and pT of mother B hadron
if(partMotherCopy->Pt()) fHistComm[iq][icut].fePtRatio->Fill(partMotherCopy->Pt(),mcpart->Pt()/partMotherCopy->Pt());
- // distance between electron production point and mother hadron production point
- TLorentzVector ebPoint;
- partMotherCopy->ProductionVertex(ebPoint);
- TLorentzVector edistance = ePoint - ebPoint;
- fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),edistance.P());
+ // distance between electron production point and primary vertex
+ fHistComm[iq][icut].feDistance->Fill(partMotherCopy->Pt(),decayLxy);
}
}
} // end of if
}
//__________________________________________
-Int_t AliHFEmcQA::GetElectronSource(AliAODMCParticle *mcpart)
+Int_t AliHFEmcQA::GetSource(AliAODMCParticle * const mcpart)
{
- // decay electron's origin
+ // decay particle's origin
- if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
+ //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
Int_t origin = -1;
Bool_t isFinalOpenCharm = kFALSE;
}
//__________________________________________
-Int_t AliHFEmcQA::GetElectronSource(TParticle* mcpart)
+Int_t AliHFEmcQA::GetSource(TParticle * const mcpart)
{
- // decay electron's origin
+ // decay particle's origin
- if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
+ //if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
Int_t origin = -1;
Bool_t isFinalOpenCharm = kFALSE;
class TParticle;
class TString;
class AliStack;
+class AliGenEventHeader;
class AliAODMCParticle;
//________________________________________________________________
void PostAnalyze() const;
void CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt=""); // create histograms for mc qa analysis
void SetStack(AliStack* const stack){fStack=stack;} // set stack pointer
+ void SetGenEventHeader(AliGenEventHeader* const mcHeader){fMCHeader=mcHeader;} // set stack pointer
void SetMCArray(TClonesArray* const mcarry){fMCArray=mcarry;} // set mcarray pointer
void Init();
void GetDecayedKine(TParticle *part, const Int_t kquark, const Int_t kdecayed, Int_t icut); // get decay electron kinematics distribution
void GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut); // get decay electron kinematics for AOD
void EndOfEventAna(const Int_t kquark); // run analysis which should be done at the end of the event loop
- Int_t GetElectronSource(TParticle* mcpart); // return electron source id
- Int_t GetElectronSource(AliAODMCParticle *mcpart); // return electron source id for AOD
+ Int_t GetSource(TParticle * const mcpart); // return electron source id
+ Int_t GetSource(AliAODMCParticle * const mcpart); // return electron source id for AOD
protected:
void IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel); //
Float_t GetRapidity(AliAODMCParticle *part) const; // return rapidity
AliStack* fStack; // stack pointer
- TClonesArray *fMCArray; // mc array pointer
+ AliGenEventHeader* fMCHeader; // mcheader pointer
+ TClonesArray *fMCArray; // mc array pointer
static const Int_t fgkGluon; // gluon pdg code
static const Int_t fgkMaxGener; // ancester level wanted to be checked
//
// PID Steering Class
// Interface to the user task
+// Several strategies for Electron identification implemented.
+// In addition users can combine different detectors or use
+// single detector PID
//
// Authors:
// Markus Fasel <M.Fasel@gsi.de>
case 3: InitStrategy3(); break;
case 4: InitStrategy4(); break;
case 5: InitStrategy5(); break;
+ case 6: InitStrategy6(); break;
default: strategyStatus = kFALSE;
}
return strategyStatus;
fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRD PID");
SETBIT(fEnabledDetectors, kTRDpid);
} else if(det->String().CompareTo("TOF") == 0){
+ AliInfo("Doing TOF PID");
fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOF PID");
SETBIT(fEnabledDetectors, kTOFpid);
}
// MC Event
return (TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track)) == 11);
}
- if(fPIDstrategy > 0 && fPIDstrategy < 6){
+ if(fPIDstrategy > 0 && fPIDstrategy < 7){
Int_t pid = 0;
switch(fPIDstrategy){
case 1: pid = IdentifyStrategy1(track); break;
case 3: pid = IdentifyStrategy3(track); break;
case 4: pid = IdentifyStrategy4(track); break;
case 5: pid = IdentifyStrategy5(track); break;
+ case 6: pid = IdentifyStrategy6(track); break;
default: break;
}
return pid;
//
// Combines TPC and TOF PID decision
//
- if(fDetectorPID[kTOFpid]->IsSelected(track)) return fDetectorPID[kTPCpid]->IsSelected(track);
- return kFALSE;
+ if(track->fAnalysisType != AliHFEpidObject::kESDanalysis) return kFALSE;
+ AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(track->fRecTrack);
+ if(!esdTrack) return kFALSE;
+
+ AliHFEpidTOF *tofPID = dynamic_cast<AliHFEpidTOF*>(fDetectorPID[kTOFpid]);
+ if(!tofPID){
+ AliWarning("TOF pid object is NULL");
+ return kFALSE;
+ }
+
+
+ AliHFEpidTPC *tpcPID = dynamic_cast<AliHFEpidTPC*>(fDetectorPID[kTPCpid]);
+ if(!tpcPID){
+ AliWarning("TPC pid object is NULL");
+ return kFALSE;
+ }
+
+ // charge of the particle
+ Int_t charge = 0;
+ if(esdTrack->GetOuterParam()) charge = esdTrack->GetOuterParam()->Charge();
+ else charge = esdTrack->GetInnerParam()->Charge();
+
+ // setup the TPC PID
+ const Int_t cTPCsigma = 2;
+ // set the number of sigmas in the TPC
+ tpcPID->SetTPCnSigma(cTPCsigma);
+ // turn on the the dEdx line crossing for kaons and protons
+ tpcPID->AddTPCdEdxLineCrossing(3, cTPCsigma); // for kaons
+ tpcPID->AddTPCdEdxLineCrossing(4, cTPCsigma); // for protons
+ // request the tpc PID information
+ Int_t pidTPC = tpcPID->IsSelected(track);
+
+ // without TPC pid it makes no sense to look at the TOF pid with this detector combination
+ //if(0 == pidTPC){
+ // return kFALSE;
+ //}
+
+ // is the TPC in the crossing line region? 0 - no crossing, otherwise the AliPID of the crossing
+ // particle returned - 3 for Kaon and 4 for Proton
+ Int_t crossing = tpcPID->GetCrossingType();
+
+ // setup the TOF pid
+ const Int_t cTOFsigma = 3;
+ tofPID->SetTOFnSigma(cTOFsigma);
+ // request TOF PID information
+ Int_t pidTOF = tofPID->IsSelected(track);
+
+ // case that the TOF for some reason does not deliver a PID information
+ // or the TPC is not in the crossing point region, only the TPC will be used
+ if(0 != pidTOF && 0 != crossing){
+ if(3 == crossing){
+ return (321 == pidTOF) ? kFALSE : kTRUE;
+ }
+ else if(4 == crossing){
+ return (2212 == pidTOF) ? kFALSE : kTRUE;
+ }
+ else{
+ // something went wrong
+ AliError("Wrong crossing type returned from AliHFEpidTPC - check your code!!");
+ return kFALSE;
+ }
+ }
+ else{
+ // tpc ONLY
+ return (11 == pidTPC) ? kTRUE : kFALSE;
+ }
}
//____________________________________________________________
//
AliHFEpidTPC *pid = new AliHFEpidTPC("strat3TPCpid");
pid->SetTPCnSigma(3);
- pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 3.);
+ pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 1.);
Bool_t status = pid->InitializePID();
if(IsQAOn() && status) pid->SetQAOn(fQAlist);
if(HasMCData() && status) pid->SetHasMCData();
fDetectorPID[kTRDpid] = trdpid;
}
+//____________________________________________________________
+void AliHFEpid::InitStrategy6(){
+ //
+ // Combined TPC-TOF PID, combination is discribed in the funtion MakePidTpcTof
+ //
+ AliHFEpidTPC *tpcpid = new AliHFEpidTPC("strat6TPCpid");
+ AliHFEpidTOF *tofpid = new AliHFEpidTOF("strat6TOFpid");
+ Bool_t status = tpcpid->InitializePID();
+ if(!status)
+ AliError("Initialization of TPC PID failed");
+ Bool_t status1 = tofpid->InitializePID();
+ if(!status1)
+ AliError("Initialization of TOF PID failed");
+ status &= status1;
+ if(IsQAOn() && status){
+ tpcpid->SetQAOn(fQAlist);
+ tofpid->SetQAOn(fQAlist);
+ }
+ if(HasMCData() && status){
+ tpcpid->SetHasMCData();
+ tofpid->SetHasMCData();
+ }
+ fDetectorPID[kTPCpid] = tpcpid;
+ fDetectorPID[kTOFpid] = tofpid;
+}
+
//____________________________________________________________
Bool_t AliHFEpid::IdentifyStrategy1(AliHFEpidObject *track){
return TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) == 11;
Bool_t AliHFEpid::IdentifyStrategy5(AliHFEpidObject *track){
return MakePidTpcTrd(track);
}
+
+//____________________________________________________________
+Bool_t AliHFEpid::IdentifyStrategy6(AliHFEpidObject *track){
+ return MakePidTpcTof(track);
+}
+
void InitStrategy3();
void InitStrategy4();
void InitStrategy5();
+ void InitStrategy6();
Bool_t IdentifyStrategy1(AliHFEpidObject *track);
Bool_t IdentifyStrategy2(AliHFEpidObject *track);
Bool_t IdentifyStrategy3(AliHFEpidObject *track);
Bool_t IdentifyStrategy4(AliHFEpidObject *track);
Bool_t IdentifyStrategy5(AliHFEpidObject *track);
+ Bool_t IdentifyStrategy6(AliHFEpidObject *track);
private:
enum{
kIsQAOn = BIT(14),
#include "AliLog.h"
#include "AliMCParticle.h"
#include "AliPID.h"
+#include "AliTOFpidESD.h"
#include "AliHFEpidTOF.h"
#include "AliHFEpidBase.h"
AliHFEpidBase(name)
, fPID(0x0)
, fQAList(0x0)
+ , fPIDtofESD(0x0)
+ , fNsigmaTOF(3)
{
//
// Constructor
//
+
+ fPIDtofESD = new AliTOFpidESD();
+
}
//___________________________________________________________________
AliHFEpidTOF::AliHFEpidTOF(const AliHFEpidTOF &c):
AliHFEpidBase("")
, fPID(0x0)
, fQAList(0x0)
+ , fPIDtofESD(0x0)
+ , fNsigmaTOF(3)
{
//
// Copy operator
// Destructor
//
if(fPID) delete fPID;
+ if(fPIDtofESD) delete fPIDtofESD;
if(fQAList){
fQAList->Delete();
delete fQAList;
target.fPID = fPID;
target.fQAList = fQAList;
+ target.fPIDtofESD = new AliTOFpidESD(*fPIDtofESD);
AliHFEpidBase::Copy(ref);
}
// the ESD PID will be checked and if necessary improved
// in the mean time. Best of luck
//
- // returns 10 (== kUnknown)if PID can not be assigned
+ // returns 10 (== kUnknown)if PID can not be assigned for whatever reason
//
+
if(vtrack->fAnalysisType == AliHFEpidObject::kESDanalysis){
AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(vtrack->fRecTrack);
if(!esdTrack) return 0;
Long_t status = 0;
status = track->GetStatus();
- if(!(status & AliESDtrack::kTOFout)) return AliPID::kUnknown;
+ if(!(status & AliESDtrack::kTOFout)) return 0;
- (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(0.);
+ if(IsQAon())(dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(0.);
Double_t tItrackL = track->GetIntegratedLength();
Double_t tTOFsignal = track->GetTOFsignal();
- if(tItrackL > 0)
- (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(1.);
+ if(IsQAon()){
+ if(tItrackL > 0)
+ (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(1.);
- if(tTOFsignal > 0)
- (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(2.);
+ if(tTOFsignal > 0)
+ (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(2.);
+ }
- if(tItrackL <=0 || tTOFsignal <=0) return AliPID::kUnknown;
+ if(tItrackL <=0 || tTOFsignal <=0) return 0;
- (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(3.);
- (dynamic_cast<TH1F *>(fQAList->At(kHistTOFsignal)))->Fill(tTOFsignal/1000.);
- (dynamic_cast<TH1F *>(fQAList->At(kHistTOFlength)))->Fill(tItrackL);
+ if(IsQAon()){
+ (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(3.);
+ (dynamic_cast<TH1F *>(fQAList->At(kHistTOFsignal)))->Fill(tTOFsignal/1000.);
+ (dynamic_cast<TH1F *>(fQAList->At(kHistTOFlength)))->Fill(tItrackL);
+ }
// get the TOF pid probabilities
Double_t tESDpid[5] = {0., 0., 0., 0., 0.};
Float_t tTOFpidSum = 0.;
tMAXindex = i;
}
}
+
+ Int_t pdg = 0;
+
+ switch(tMAXindex){
+ case 0: pdg = 11; break;
+ case 1: pdg = 13; break;
+ case 2: pdg = 211; break;
+ case 3: pdg = 321; break;
+ case 4: pdg = 2212; break;
+ default: pdg = 0;
+ };
+
Double_t p = track->GetOuterParam()->P();
Double_t beta = (tItrackL/100.)/(TMath::C()*(tTOFsignal/1e12));
- if(TMath::Abs(tTOFpidSum - 1) > 0.01) return AliPID::kUnknown;
- else{
- // should be the same as AliPID flags
-
+ // sanity check, should not be necessary
+ if(TMath::Abs(tTOFpidSum - 1) > 0.01) return 0;
+
+ Double_t nSigmas = fPIDtofESD->GetNumberOfSigmas(track, (AliPID::EParticleType)tMAXindex);
+ if(TMath::Abs(nSigmas) > fNsigmaTOF) return 0;
+
+
+ // should be the same as AliPID flags
+
+ if(IsQAon()){
(dynamic_cast<TH2F *>(fQAList->At(kHistTOFpid0+tMAXindex)))->Fill(beta, p);
(dynamic_cast<TH2F *>(fQAList->At(kHistTOFpidBetavP)))->Fill(beta, p);
- return tMAXindex;
}
+ //return tMAXindex;
+ return pdg;
+
+}
+//___________________________________________________________________
+Double_t AliHFEpidTOF::Likelihood(const AliESDtrack *track, Int_t species, Float_t rsig){
+
+ //gives probability for track to come from a certain particle species;
+ //no priors considered -> probabilities for equal abundances of all species!
+ //optional restriction to r-sigma bands of different particle species; default: rsig = 2. (see below)
+
+ //IMPORTANT: Tracks which are judged to be outliers get negative likelihoods -> unusable for combination with further detectors!
+
+ if(!track) return -1.;
+ Bool_t outlier = kTRUE;
+ // Check whether distance from the respective particle line is smaller than r sigma
+ for(Int_t hypo = 0; hypo < AliPID::kSPECIES; hypo++){
+ if(TMath::Abs(fPIDtofESD->GetNumberOfSigmas(track, (AliPID::EParticleType)hypo)) > rsig)
+ outlier = kTRUE;
+ else {
+ outlier = kFALSE;
+ break;
+ }
+ }
+ if(outlier)
+ return -2.;
+
+ Double_t tofProb[5];
+
+ track->GetTOFpid(tofProb);
+
+ return tofProb[species];
}
-
//___________________________________________________________________
Int_t AliHFEpidTOF::MakePIDaod(AliAODTrack * /*aodTrack*/, AliAODMCParticle * /*mcTrack*/){
AliError("AOD PID not yet implemented");
class AliAODMCParticle;
class AliESDtrack;
class AliMCParticle;
+class AliTOFpidESD;
class AliHFEpidTOF : public AliHFEpidBase{
public:
virtual Int_t IsSelected(AliHFEpidObject *track);
virtual Bool_t HasQAhistos() const { return kTRUE; };
-
+ void SetTOFnSigma(Short_t nSigma) { fNsigmaTOF = nSigma; };
+
+ Double_t Likelihood(const AliESDtrack *track, Int_t species, Float_t rsig = 2.);
+
protected:
void Copy(TObject &ref) const;
void AddQAhistograms(TList *qaHist);
kHistTOFpid2 = 6,
kHistTOFpid3 = 7,
kHistTOFpid4 = 8
-
} QAHist_t;
+
+ AliPID *fPID; //! PID Object
+ TList *fQAList; //! QA histograms
+ AliTOFpidESD *fPIDtofESD; //! TOF pid object
+
+ Short_t fNsigmaTOF; // TOF sigma band
- AliPID *fPID; //! PID Object
- TList *fQAList; //! QA histograms
- ClassDef(AliHFEpidTOF, 1) // TOF PID class
+ ClassDef(AliHFEpidTOF, 1)
};
#endif
//___________________________________________________________________
AliHFEpidTPC::AliHFEpidTPC(const char* name) :
// add a list here
- AliHFEpidBase(name)
+ AliHFEpidBase(name)
+ , fLineCrossingType(0)
, fLineCrossingsEnabled(0)
, fNsigmaTPC(3)
, fRejectionEnabled(0)
//___________________________________________________________________
AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
- AliHFEpidBase("")
+ AliHFEpidBase("")
+ , fLineCrossingType(0)
, fLineCrossingsEnabled(0)
, fNsigmaTPC(2)
, fRejectionEnabled(0)
return MakePIDaod(aodtrack, aodmctrack);
}
}
-
//___________________________________________________________________
Int_t AliHFEpidTPC::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle *mctrack){
//
// Doing TPC PID as explained in IsSelected for ESD tracks
//
if(IsQAon()) FillTPChistograms(esdTrack, mctrack);
+ Float_t nsigma = fPIDtpcESD->GetNumberOfSigmas(esdTrack, AliPID::kElectron);
// exclude crossing points:
// Determine the bethe values for each particle species
Bool_t isLineCrossing = kFALSE;
+ fLineCrossingType = 0; // default value
for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
if(ispecies == AliPID::kElectron) continue;
if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
- if(fPIDtpcESD->GetNumberOfSigmas(esdTrack, (AliPID::EParticleType)ispecies) < fLineCrossingSigma[ispecies]){
- // Point in a line crossing region, no PID possible
- isLineCrossing = kTRUE;
+ if(TMath::Abs(fPIDtpcESD->GetNumberOfSigmas(esdTrack, (AliPID::EParticleType)ispecies)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
+ // Point in a line crossing region, no PID possible, but !PID still possible ;-)
+ isLineCrossing = kTRUE;
+ fLineCrossingType = ispecies;
break;
}
}
// Check particle rejection
if(HasParticleRejection()){
Int_t reject = Reject(esdTrack);
- AliDebug(1, Form("PID code from Rejection: %d", reject));
if(reject != 0) return reject;
}
// Check whether distance from the electron line is smaller than n-sigma
// Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
- Float_t nsigma = fPIDtpcESD->GetNumberOfSigmas(esdTrack, AliPID::kElectron);
Float_t p = 0.;
Int_t pdg = 0;
if(HasAsymmetricSigmaCut() && (p = esdTrack->P()) >= fPAsigCut[0] && p <= fPAsigCut[1]){
if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
}
if(IsQAon() && pdg != 0) (dynamic_cast<TH2I *>(fQAList->At(kHistTPCselected)))->Fill(esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P(), esdTrack->GetTPCsignal());
+
return pdg;
}
Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
if(!TESTBIT(fRejectionEnabled, ispec)) continue;
- AliDebug(1, Form("Particle Rejection enabled for species %d", ispec));
// Particle rejection enabled
if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
Double_t sigma = fPIDtpcESD->GetNumberOfSigmas(track, static_cast<AliPID::EParticleType>(ispec));
- AliDebug(1, Form("Sigma %f, min %f, max %f", sigma, fRejection[4*ispec + 1], fRejection[4*ispec+3]));
if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
}
return 0;
fQAList = new TList;
fQAList->SetName("fTPCqaHistos");
- fQAList->AddAt(new TH2I("fHistTPCelectron","TPC signal for Electrons", 200, 0, 20, 200, 0, 200), kHistTPCelectron);
- fQAList->AddAt(new TH2I("fHistTPCmuon","TPC signal for Muons", 200, 0, 20, 200, 0, 200), kHistTPCmuon);
- fQAList->AddAt(new TH2I("fHistTPCpion","TPC signal for Pions", 200, 0, 20, 200, 0, 200), kHistTPCpion);
- fQAList->AddAt(new TH2I("fHistTPCkaon","TPC signal for Kaons", 200, 0, 20, 200, 0, 200), kHistTPCkaon);
- fQAList->AddAt(new TH2I("fHistTPCproton","TPC signal for Protons", 200, 0, 20, 200, 0, 200), kHistTPCproton);
- fQAList->AddAt(new TH2I("fHistTPCothers","TPC signal for other species", 200, 0, 20, 200, 0, 200), kHistTPCothers);
- fQAList->AddAt(new TH2I("fHistTPCall","TPC signal for all species", 200, 0, 20, 200, 0, 200), kHistTPCall);
- fQAList->AddAt(new TH2I("fHistTPCselected","TPC signal for all selected particles", 200, 0, 20, 200, 0, 200), kHistTPCselected);
+ fQAList->AddAt(new TH2I("fHistTPCelectron","TPC signal for Electrons", 200, 0, 20, 60, 0, 600), kHistTPCelectron);
+ fQAList->AddAt(new TH2I("fHistTPCmuon","TPC signal for Muons", 200, 0, 20, 60, 0, 600), kHistTPCmuon);
+ fQAList->AddAt(new TH2I("fHistTPCpion","TPC signal for Pions", 200, 0, 20, 60, 0, 600), kHistTPCpion);
+ fQAList->AddAt(new TH2I("fHistTPCkaon","TPC signal for Kaons", 200, 0, 20, 60, 0, 600), kHistTPCkaon);
+ fQAList->AddAt(new TH2I("fHistTPCproton","TPC signal for Protons", 200, 0, 20, 60, 0, 600), kHistTPCproton);
+ fQAList->AddAt(new TH2I("fHistTPCothers","TPC signal for other species", 200, 0, 20, 60, 0, 600), kHistTPCothers);
+ fQAList->AddAt(new TH2I("fHistTPCall","TPC signal for all species", 200, 0, 20, 60, 0, 600), kHistTPCall);
+ fQAList->AddAt(new TH2I("fHistTPCselected","TPC signal for all selected particles", 200, 0, 20, 60, 0, 600), kHistTPCselected);
fQAList->AddAt(new TH2F("fHistTPCprobEl","TPC likelihood for electrons to be an electron vs. p", 200, 0.,20.,200,0.,1.), kHistTPCprobEl);
fQAList->AddAt(new TH2F("fHistTPCprobPi","TPC likelihood for pions to be an electron vs. p", 200, 0.,20.,200, 0.,1.), kHistTPCprobPi);
virtual Int_t IsSelected(AliHFEpidObject *track);
virtual Bool_t HasQAhistos() const { return kTRUE; };
+ Int_t GetCrossingType() const {return fLineCrossingType; }
+
void AddTPCdEdxLineCrossing(Int_t species, Double_t sigma);
Bool_t HasAsymmetricSigmaCut() const { return TestBit(kAsymmetricSigmaCut);}
Bool_t HasParticleRejection() const { return TestBit(kRejection); }
kRejection = BIT(21)
};
Double_t fLineCrossingSigma[AliPID::kSPECIES]; // with of the exclusion point
+ Int_t fLineCrossingType; // 0 for no line crossing, otherwise AliPID of the particle crossing the electron dEdx band
UChar_t fLineCrossingsEnabled; // Bitmap showing which line crossing is set
Float_t fPAsigCut[2]; // Momentum region where to perform asymmetric sigma cut
Float_t fNAsigmaTPC[2]; // Asymmetric TPC Sigma band
* provided "as is" without express or implied warranty. *
**************************************************************************/
//
-// QA class of primary vertex study for Heavy Flavor electrons
-//
+// QA class of primary vertex study for Heavy Flavor electrons
+// this has functionality to reject electrons from primary vertex
+// and check primary vertex characteristics
+//
// Authors:
// MinJung Kweon <minjung@physi.uni-heidelberg.de>
//
**************************************************************************/
//
// QA class of primary vertex study for Heavy Flavor electrons
+// this has functionality to reject electrons from primary vertex
+// and check primary vertex characteristics
//
#ifndef ALIHFEPRIVTX_H
-/**************************************************************************
+ /**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
* *
* Author: The ALICE Off-line Project. *
#include <TParticle.h>
#include <AliESDEvent.h>
+#include <AliAODEvent.h>
+#include <AliVTrack.h>
#include <AliESDtrack.h>
+#include <AliAODTrack.h>
#include "AliHFEsecVtx.h"
#include <AliKFParticle.h>
#include <AliKFVertex.h>
#include <AliLog.h>
#include <AliStack.h>
+#include <AliAODMCParticle.h>
+#include "AliHFEpairs.h"
+#include "AliHFEsecVtxs.h"
ClassImp(AliHFEsecVtx)
//_______________________________________________________________________________________________
AliHFEsecVtx::AliHFEsecVtx():
- fESD1(0x0)
- ,fStack(0x0)
- ,fNparents(0)
- ,fHistTagged()
- ,fPairTagged(0)
- ,fdistTwoSecVtx(-1)
- ,fcosPhi(-1)
- ,fsignedLxy(-1)
- ,finvmass(-1)
- ,finvmassSigma(-1)
- ,fBTagged(0)
- ,fBElec(0)
+ fESD1(0x0)
+ ,fAOD1(0x0)
+ ,fStack(0x0)
+ ,fUseMCPID(kFALSE)
+ ,fkSourceLabel()
+ ,fNparents(0)
+ ,fParentSelect()
+ ,fNoOfHFEpairs(0)
+ ,fNoOfHFEsecvtxs(0)
+ ,fHFEpairs(0x0)
+ ,fHFEsecvtxs(0x0)
+ ,fMCArray(0x0)
+ ,fPVx(0)
+ ,fPVy(0)
+ ,fCosPhi(-1)
+ ,fSignedLxy(-1)
+ ,fKFchi2(-1)
+ ,fInvmass(-1)
+ ,fInvmassSigma(-1)
+ ,fKFip(0)
+ ,fPairQA(0x0)
+ ,fSecvtxQA(0x0)
+ ,fSecVtxList(0x0)
{
- // Default constructor
-
- Init();
+ //
+ // Default constructor
+ //
+ Init();
}
//_______________________________________________________________________________________________
AliHFEsecVtx::AliHFEsecVtx(const AliHFEsecVtx &p):
- TObject(p)
- ,fESD1(0x0)
- ,fStack(0x0)
- ,fNparents(p.fNparents)
- ,fHistTagged()
- ,fPairTagged(p.fPairTagged)
- ,fdistTwoSecVtx(p.fdistTwoSecVtx)
- ,fcosPhi(p.fcosPhi)
- ,fsignedLxy(p.fsignedLxy)
- ,finvmass(p.finvmass)
- ,finvmassSigma(p.finvmassSigma)
- ,fBTagged(p.fBTagged)
- ,fBElec(p.fBElec)
+ TObject(p)
+ ,fESD1(0x0)
+ ,fAOD1(0x0)
+ ,fStack(0x0)
+ ,fUseMCPID(p.fUseMCPID)
+ ,fkSourceLabel()
+ ,fNparents(p.fNparents)
+ ,fParentSelect()
+ ,fNoOfHFEpairs(p.fNoOfHFEpairs)
+ ,fNoOfHFEsecvtxs(p.fNoOfHFEsecvtxs)
+ ,fHFEpairs(0x0)
+ ,fHFEsecvtxs(0x0)
+ ,fMCArray(0x0)
+ ,fPVx(p.fPVx)
+ ,fPVy(p.fPVy)
+ ,fCosPhi(p.fCosPhi)
+ ,fSignedLxy(p.fSignedLxy)
+ ,fKFchi2(p.fKFchi2)
+ ,fInvmass(p.fInvmass)
+ ,fInvmassSigma(p.fInvmassSigma)
+ ,fKFip(p.fKFip)
+ ,fPairQA(0x0)
+ ,fSecvtxQA(0x0)
+ ,fSecVtxList(0x0)
{
- // Copy constructor
+ //
+ // Copy constructor
+ //
}
//_______________________________________________________________________________________________
AliHFEsecVtx&
AliHFEsecVtx::operator=(const AliHFEsecVtx &)
{
+ //
// Assignment operator
+ //
AliInfo("Not yet implemented.");
return *this;
//_______________________________________________________________________________________________
AliHFEsecVtx::~AliHFEsecVtx()
{
- // Destructor
+ //
+ // Destructor
+ //
- //cout << "Analysis Done." << endl;
+ //cout << "Analysis Done." << endl;
}
//__________________________________________
void AliHFEsecVtx::Init()
{
-
- // set pdg code and index
-
- fNparents = 7;
-
- fParentSelect[0][0] = 411;
- fParentSelect[0][1] = 421;
- fParentSelect[0][2] = 431;
- fParentSelect[0][3] = 4122;
- fParentSelect[0][4] = 4132;
- fParentSelect[0][5] = 4232;
- fParentSelect[0][6] = 4332;
-
- fParentSelect[1][0] = 511;
- fParentSelect[1][1] = 521;
- fParentSelect[1][2] = 531;
- fParentSelect[1][3] = 5122;
- fParentSelect[1][4] = 5132;
- fParentSelect[1][5] = 5232;
- fParentSelect[1][6] = 5332;
-
-/*
- fid[0][0] = 0;
- fid[0][1] = 1;
- fid[0][2] = 2;
-
- fid[1][0] = 0;
- fid[1][1] = 1;
- fid[1][2] = 3;
-
- fid[2][0] = 0;
- fid[2][1] = 2;
- fid[2][2] = 3;
-
- fid[3][0] = 1;
- fid[3][1] = 2;
- fid[3][2] = 3;
-
- fia[0][0][0] = 0;
- fia[0][0][1] = 1;
- fia[0][1][0] = 0;
- fia[0][1][1] = 2;
- fia[0][2][0] = 1;
- fia[0][2][1] = 2;
-
- fia[1][0][0] = 0;
- fia[1][0][1] = 1;
- fia[1][1][0] = 0;
- fia[1][1][1] = 3;
- fia[1][2][0] = 1;
- fia[1][2][1] = 3;
-
- fia[2][0][0] = 0;
- fia[2][0][1] = 2;
- fia[2][1][0] = 0;
- fia[2][1][1] = 3;
- fia[2][2][0] = 2;
- fia[2][2][1] = 3;
-
- fia[3][0][0] = 1;
- fia[3][0][1] = 2;
- fia[3][1][0] = 1;
- fia[3][1][1] = 3;
- fia[3][2][0] = 2;
- fia[3][2][1] = 3;
-*/
-
+ //
+ // set pdg code and index
+ //
+
+ fNparents = 7;
+
+ fParentSelect[0][0] = 411;
+ fParentSelect[0][1] = 421;
+ fParentSelect[0][2] = 431;
+ fParentSelect[0][3] = 4122;
+ fParentSelect[0][4] = 4132;
+ fParentSelect[0][5] = 4232;
+ fParentSelect[0][6] = 4332;
+
+ fParentSelect[1][0] = 511;
+ fParentSelect[1][1] = 521;
+ fParentSelect[1][2] = 531;
+ fParentSelect[1][3] = 5122;
+ fParentSelect[1][4] = 5132;
+ fParentSelect[1][5] = 5232;
+ fParentSelect[1][6] = 5332;
}
-//__________________________________________
-void AliHFEsecVtx::ResetTagVar()
-{
- // reset tag variables
-
- fdistTwoSecVtx = -1;
- fcosPhi = -1;
- fsignedLxy = -1;
- finvmass = -1;
- finvmassSigma = -1;
- fBTagged = kFALSE;
- fBElec = kFALSE;
+//_______________________________________________________________________________________________
+void AliHFEsecVtx::CreateHistograms(TList *qaList)
+{
+ //
+ // create histograms
+ //
+
+ fSecVtxList = new TList;
+ fSecVtxList->SetName("SecVtx");
+
+ MakeContainer();
+ /*
+ fkSourceLabel[kAll]="all";
+ fkSourceLabel[kDirectCharm]="directCharm";
+ fkSourceLabel[kDirectBeauty]="directBeauty";
+ fkSourceLabel[kBeautyCharm]="beauty2charm";
+ fkSourceLabel[kGamma]="gamma";
+ fkSourceLabel[kPi0]="pi0";
+ fkSourceLabel[kElse]="others";
+ fkSourceLabel[kBeautyGamma]="beauty22gamma";
+ fkSourceLabel[kBeautyPi0]="beauty22pi0";
+ fkSourceLabel[kBeautyElse]="beauty22others";
+
+ TString hname;
+ TString hnopt="secvtx_";
+ for (Int_t isource = 0; isource < 10; isource++ ){
+ }
+ */
+
+ qaList->AddLast(fSecVtxList);
}
-//__________________________________________
-void AliHFEsecVtx::InitAnaPair()
+//_______________________________________________________________________________________________
+void AliHFEsecVtx::GetPrimaryCondition()
{
- // initialize pair tagging variables
-
- fPairTagged = 0;
- for (Int_t i=0; i<20; i++){
- fpairedTrackID[i] = -1;
- fpairedChi2[i] = -1;
- fpairedInvMass[i] = -1;
- fpairedSignedLxy[i] = -1;
- }
-
- fid[0][0] = 0;
- fid[0][1] = 1;
- fid[0][2] = 2;
-
- fid[1][0] = 0;
- fid[1][1] = 1;
- fid[1][2] = 3;
-
- fid[2][0] = 0;
- fid[2][1] = 2;
- fid[2][2] = 3;
-
- fid[3][0] = 1;
- fid[3][1] = 2;
- fid[3][2] = 3;
-
- fia[0][0][0] = 0;
- fia[0][0][1] = 1;
- fia[0][1][0] = 0;
- fia[0][1][1] = 2;
- fia[0][2][0] = 1;
- fia[0][2][1] = 2;
-
- fia[1][0][0] = 0;
- fia[1][0][1] = 1;
- fia[1][1][0] = 0;
- fia[1][1][1] = 3;
- fia[1][2][0] = 1;
- fia[1][2][1] = 3;
-
- fia[2][0][0] = 0;
- fia[2][0][1] = 2;
- fia[2][1][0] = 0;
- fia[2][1][1] = 3;
- fia[2][2][0] = 2;
- fia[2][2][1] = 3;
-
- fia[3][0][0] = 1;
- fia[3][0][1] = 2;
- fia[3][1][0] = 1;
- fia[3][1][1] = 3;
- fia[3][2][0] = 2;
- fia[3][2][1] = 3;
-
+ //
+ // get primary characteristics and set
+ //
+
+ if (fESD1) {
+ AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
+ if( primVtxCopy.GetNDF() <1 ) return;
+ fPVx = primVtxCopy.GetX();
+ fPVx = primVtxCopy.GetY();
+ }
+ else if(fAOD1) {
+ AliKFVertex primVtxCopy(*(fAOD1->GetPrimaryVertex()));
+ if( primVtxCopy.GetNDF() <1 ) return;
+ fPVx = primVtxCopy.GetX();
+ fPVx = primVtxCopy.GetY();
+ }
}
//_______________________________________________________________________________________________
-void AliHFEsecVtx::CreateHistograms(TString hnopt)
-{
- // create histograms
-
- fkSourceLabel[kAll]="all";
- fkSourceLabel[kDirectCharm]="directCharm";
- fkSourceLabel[kDirectBeauty]="directBeauty";
- fkSourceLabel[kBeautyCharm]="beauty2charm";
- fkSourceLabel[kGamma]="gamma";
- fkSourceLabel[kPi0]="pi0";
- fkSourceLabel[kElse]="others";
- fkSourceLabel[kBeautyGamma]="beauty22gamma";
- fkSourceLabel[kBeautyPi0]="beauty22pi0";
- fkSourceLabel[kBeautyElse]="beauty22others";
-
-
- TString hname;
- for (Int_t isource = 0; isource < 10; isource++ ){
-
- hname=hnopt+"InvMass_"+fkSourceLabel[isource];
- fHistPair[isource].fInvMass = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
- hname=hnopt+"InvMassCut1_"+fkSourceLabel[isource];
- fHistPair[isource].fInvMassCut1 = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
- hname=hnopt+"InvMassCut2_"+fkSourceLabel[isource];
- fHistPair[isource].fInvMassCut2 = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
- hname=hnopt+"KFChi2_"+fkSourceLabel[isource];
- fHistPair[isource].fKFChi2 = new TH1F(hname,hname,200,0,20);
- hname=hnopt+"OpenAngle_"+fkSourceLabel[isource];
- fHistPair[isource].fOpenAngle = new TH1F(hname,hname,100,0,3.14);
- hname=hnopt+"CosOpenAngle_"+fkSourceLabel[isource];
- fHistPair[isource].fCosOpenAngle = new TH1F(hname,hname,100,-1.1,1.1);
- hname=hnopt+"SignedLxy_"+fkSourceLabel[isource];
- fHistPair[isource].fSignedLxy = new TH2F(hname,hname,1000,-5,5,120,-2,10);
- hname=hnopt+"KFIP_"+fkSourceLabel[isource];
- fHistPair[isource].fKFIP = new TH1F(hname,hname,1000,-5,5);
- hname=hnopt+"IPMax_"+fkSourceLabel[isource];
- fHistPair[isource].fIPMax= new TH1F(hname,hname,500,0,5);
-
- }
+void AliHFEsecVtx::PairAnalysis(AliVTrack* track1, AliVTrack* track2, Int_t index2)
+{
+ //
+ // calculate e-h pair characteristics and tag pair
+ //
- hname=hnopt+"pt_beautyelec";
- fHistTagged.fPtBeautyElec= new TH1F(hname,hname,250,0,50);
- hname=hnopt+"pt_taggedelec";
- fHistTagged.fPtTaggedElec= new TH1F(hname,hname,250,0,50);
- hname=hnopt+"pt_righttaggedelec";
- fHistTagged.fPtRightTaggedElec = new TH1F(hname,hname,250,0,50);
- hname=hnopt+"pt_wrongtaggedelec";
- fHistTagged.fPtWrongTaggedElec = new TH1F(hname,hname,250,0,50);
-
- hname=hnopt+"InvmassBeautyElecSecVtx";
- fHistTagged.fInvmassBeautyElecSecVtx= new TH1F(hname,hname,120,-2,10);
- hname=hnopt+"InvmassNotBeautyElecSecVtx";
- fHistTagged.fInvmassNotBeautyElecSecVtx= new TH1F(hname,hname,120,-2,10);
-
- hname=hnopt+"SignedLxyBeautyElecSecVtx";
- fHistTagged.fSignedLxyBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
- hname=hnopt+"SignedLxyNotBeautyElecSecVtx";
- fHistTagged.fSignedLxyNotBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
-
- hname=hnopt+"DistTwoVtxBeautyElecSecVtx";
- fHistTagged.fDistTwoVtxBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
- hname=hnopt+"DistTwoVtxNotBeautyElecSecVtx";
- fHistTagged.fDistTwoVtxNotBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
-
- hname=hnopt+"CosPhiBeautyElecSecVtx";
- fHistTagged.fCosPhiBeautyElecSecVtx= new TH1F(hname,hname,100,-1.1,1.1);
- hname=hnopt+"CosPhiNotBeautyElecSecVtx";
- fHistTagged.fCosPhiNotBeautyElecSecVtx= new TH1F(hname,hname,100,-1.1,1.1);
-
- hname=hnopt+"Chi2BeautyElecSecVtx";
- fHistTagged.fChi2BeautyElecSecVtx= new TH1F(hname,hname,200,0,20);
- hname=hnopt+"Chi2NotBeautyElecSecVtx";
- fHistTagged.fChi2NotBeautyElecSecVtx= new TH1F(hname,hname,200,0,20);
-
- hname=hnopt+"InvmassBeautyElec2trkVtx";
- fHistTagged.fInvmassBeautyElec2trkVtx= new TH1F(hname,hname,120,-2,10);
- hname=hnopt+"InvmassNotBeautyElec2trkVtx";
- fHistTagged.fInvmassNotBeautyElec2trkVtx= new TH1F(hname,hname,120,-2,10);
-
- hname=hnopt+"SignedLxyBeautyElec2trkVtx";
- fHistTagged.fSignedLxyBeautyElec2trkVtx= new TH1F(hname,hname,1000,-5,5);
- hname=hnopt+"SignedLxyNotBeautyElec2trkVtx";
- fHistTagged.fSignedLxyNotBeautyElec2trkVtx= new TH1F(hname,hname,1000,-5,5);
-
- hname=hnopt+"Chi2BeautyElec2trkVtx";
- fHistTagged.fChi2BeautyElec2trkVtx= new TH1F(hname,hname,200,0,20);
- hname=hnopt+"Chi2NotBeautyElec2trkVtx";
- fHistTagged.fChi2NotBeautyElec2trkVtx= new TH1F(hname,hname,200,0,20);
+ // get KF particle input pid
+ Int_t pdg1 = GetPDG(track1);
+ Int_t pdg2 = GetPDG(track2);
+
+ if(pdg1==-1 || pdg2==-1) {
+ //printf("out if considered pid range \n");
+ return;
+ }
+
+ // create KF particle of pair
+ if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
+ else AliKFParticle::SetField(fESD1->GetMagneticField());
+ AliKFParticle kfTrack1(*track1, pdg1);
+ AliKFParticle kfTrack2(*track2, pdg2);
+
+ AliKFParticle kfSecondary(kfTrack1,kfTrack2);
+
+ //secondary vertex point from kf particle
+ Double_t kfx = kfSecondary.GetX();
+ Double_t kfy = kfSecondary.GetY();
+ //Double_t kfz = kfSecondary.GetZ();
+
+ //momentum at the decay point from kf particle
+ Double_t kfpx = kfSecondary.GetPx();
+ Double_t kfpy = kfSecondary.GetPy();
+ //Double_t kfpz = kfSecondary.GetPz();
+
+ Double_t dx = kfx-fPVx;
+ Double_t dy = kfy-fPVy;
+
+ // discriminating variables
+ // invariant mass of the KF particle
+ Double_t invmass = -1;
+ Double_t invmassSigma = -1;
+ kfSecondary.GetMass(invmass,invmassSigma);
+
+ // chi2 of the KF particle
+ Double_t kfchi2 = -1;
+ if(kfSecondary.GetNDF()>0) kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
+
+ // opening angle between two particles in XY plane
+ Double_t phi = kfTrack1.GetAngleXY(kfTrack2);
+ Double_t cosphi = TMath::Cos(phi);
+
+ // projection of kf vertex vector to the kf momentum direction
+ Double_t signedLxy=-999.;
+ Double_t psqr = kfpx*kfpx+kfpy*kfpy;
+ /*
+ //[the other way to think about]
+ if(psqr>0) {
+ if((dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr)>0) signedLxy = TMath::Sqrt(dx*dx+dy*dy);
+ if((dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr)<0) signedLxy = -1*TMath::Sqrt(dx*dx+dy*dy);
+ }*/
+ if(psqr>0) signedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
+
+ // DCA from primary to e-h KF particle (impact parameter of KF particle)
+ Double_t vtx[2]={fPVx, fPVy};
+ Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx);
+
+ /* //later consider the below
+ Float_t dcaR1=-999., dcaR2=-999.;
+ Float_t dcaZ1=-999., dcaZ2=-999.;
+
+ if(IsAODanalysis()){
+ Double_t trkIpPar1[2];
+ Double_t trkIpCov1[3];
+ Double_t trkIpPar2[2];
+ Double_t trkIpCov2[3];
+
+ AliESDtrack esdTrk1(track1);
+ AliESDtrack esdTrk2(track2);
+ esdTrk1.PropagateToDCA(fAOD1->GetPrimaryVertex(),0,10000,trkIpPar1,trkIpCov1);
+ esdTrk2.PropagateToDCA(fAOD1->GetPrimaryVertex(),0,10000,trkIpPar2,trkIpCov2);
+
+ dcaR1=trkIpPar1[0];
+ dcaZ1=trkIpPar1[1];
+ dcaR2=trkIpPar2[0];
+ dcaZ2=trkIpPar2[1];
+ }
+ else {
+ ((AliESDtrack*)track1)->GetImpactParameters(dcaR1,dcaZ1);
+ ((AliESDtrack*)track2)->GetImpactParameters(dcaR2,dcaZ2);
+ }*/
+
+ Int_t paircode = -1;
+ if (HasMCData()) paircode = GetPairCode(track1,track2);
+
+ AliHFEpairs hfepair;
+ hfepair.SetTrkLabel(index2);
+ hfepair.SetInvmass(invmass);
+ hfepair.SetKFChi2(kfchi2);
+ hfepair.SetOpenangle(phi);
+ hfepair.SetCosOpenangle(cosphi);
+ hfepair.SetSignedLxy(signedLxy);
+ hfepair.SetKFIP(kfip);
+ hfepair.SetPairCode(paircode);
+ AddHFEpairToArray(&hfepair);
+ fNoOfHFEpairs++;
+
+ // fill into container for later QA
+ Double_t dataE[6];
+ dataE[0]=invmass;
+ dataE[1]=kfchi2;
+ dataE[2]=phi;
+ dataE[3]=signedLxy;
+ dataE[4]=kfip;
+ dataE[5]=paircode;
+ fPairQA->Fill(dataE);
}
//_______________________________________________________________________________________________
-void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t index2)
+void AliHFEsecVtx::RunSECVTX(AliVTrack *track)
{
- // calculate e-h pair characteristics and tag pair
-
- Int_t sourcePart = PairCode(track1,track2);
-
- // get KF particle input pid
- Int_t pdg1 = GetMCPID(track1);
- Int_t pdg2 = GetMCPID(track2);
-
-
- // create KF particle of pair
- AliKFParticle::SetField(fESD1->GetMagneticField());
- AliKFParticle kfTrack1(*track1, pdg1);
- AliKFParticle kfTrack2(*track2, pdg2);
-
- AliKFParticle kfSecondary(kfTrack1,kfTrack2);
+ //
+ // run secondary vertexing algorithm and do tagging
+ //
- // copy primary vertex from ESD
- AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
- if( primVtxCopy.GetNDF() <1 ) return;
-
- //primary vertex point
- Double_t pvx = primVtxCopy.GetX();
- Double_t pvy = primVtxCopy.GetY();
- //Double_t pvz = primVtxCopy.GetZ();
+ //printf("number of considered pairs= %d\n",HFEpairs()->GetEntriesFast());
+ FindSECVTXCandid(track);
+}
- //secondary vertex point from kf particle
- Double_t kfx = kfSecondary.GetX();
- Double_t kfy = kfSecondary.GetY();
- //Double_t kfz = kfSecondary.GetZ();
-
- //momentum at the decay point from kf particle
- Double_t kfpx = kfSecondary.GetPx();
- Double_t kfpy = kfSecondary.GetPy();
- //Double_t kfpz = kfSecondary.GetPz();
+//_______________________________________________________________________________________________
+void AliHFEsecVtx::FindSECVTXCandid(AliVTrack *track)
+{
+ //
+ // find secondary vertex candidate and store them into container
+ //
+
+ AliVTrack *htrack[20];
+ Int_t htracklabel[20];
+ Double_t vtxchi2cut=3.; // testing cut
+ Double_t dataE[6]={-999.,-999.,-999.,-999.,-1.,0};
+ if (HFEpairs()->GetEntriesFast()>20){
+ AliDebug(3, "number of paired hadron is over maximum(20)");
+ return;
+ }
-
- Double_t dx = kfx-pvx;
- Double_t dy = kfy-pvy;
-
-
-
- // discriminating variables ----------------------------------------------------------
-
- // invariant mass of the KF particle
- Double_t invmass = -1;
- Double_t invmassSigma = -1;
- kfSecondary.GetMass(invmass,invmassSigma);
-
- // chi2 of the KF particle
- Double_t kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
-
- // opening angle between two particles in XY plane
- Double_t phi = kfTrack1.GetAngleXY(kfTrack2);
- Double_t cosphi = TMath::Cos(phi);
-
- // projection of kf vertex vector to the kf momentum direction
- Double_t costheta = ( dx*kfpx + dy*kfpy)/TMath::Sqrt(dx*dx+dy*dy)*TMath::Sqrt(kfpx*kfpx + kfpy*kfpy);
- Double_t signedLxy = TMath::Sqrt(dx*dx+dy*dy)*costheta;
-
- // DCA from primary to e-h KF particle (impact parameter of KF particle)
- Double_t vtx[2]={pvx, pvy};
- Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx);
-
-
- Float_t dcaR=-1;
- Float_t dcaR1=-1, dcaR2=-1;
- Float_t dcaZ1=-1, dcaZ2=-1;
- track1->GetImpactParameters(dcaR1,dcaZ1);
- track2->GetImpactParameters(dcaR2,dcaZ2);
-
- if (TMath::Abs(dcaR1) >= TMath::Abs(dcaR2)) dcaR=dcaR1;
- else dcaR=dcaR2;
-
- // fill histograms
- fHistPair[sourcePart].fInvMass->Fill(invmass,invmassSigma);
- fHistPair[sourcePart].fKFChi2->Fill(kfchi2);
- fHistPair[sourcePart].fOpenAngle->Fill(phi);
- fHistPair[sourcePart].fCosOpenAngle->Fill(cosphi);
- fHistPair[sourcePart].fSignedLxy->Fill(signedLxy,invmass);
- fHistPair[sourcePart].fKFIP->Fill(kfip);
- fHistPair[sourcePart].fIPMax->Fill(TMath::Abs(dcaR));
-
- // pair cuts
- if( kfchi2 >2. ) return;
-
- if ( signedLxy > 0.05 && cosphi > 0.5) fHistPair[sourcePart].fInvMassCut1->Fill(invmass,invmassSigma);
-
- // pair tagging if it passed the above cuts
- if(signedLxy > 0. && cosphi > 0.) fHistPair[sourcePart].fInvMassCut2->Fill(invmass,invmassSigma);
-
-
- // pair tagging condition
- if ( signedLxy > 0.0 && cosphi > 0) { // testing loose cut
- //if ( signedLxy > 0.06 && cosphi > 0) {
- fpairedTrackID[fPairTagged] = index2;
- fpairedChi2[fPairTagged] = kfchi2;
- fpairedInvMass[fPairTagged] = invmass;
- fpairedSignedLxy[fPairTagged] = signedLxy;
- fPairTagged++;
+ // get paired track objects
+ AliHFEpairs *pair=0x0;
+ if (IsAODanalysis()){
+ for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
+ pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
+ htracklabel[ip] = pair->GetTrkLabel();
+ htrack[ip] = fAOD1->GetTrack(pair->GetTrkLabel());
+ }
+ }
+ else{
+ for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
+ pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
+ htracklabel[ip] = pair->GetTrkLabel();
+ htrack[ip] = fESD1->GetTrack(pair->GetTrkLabel());
+ }
+ }
+ // in case there is only one paired track with the electron, put pair characteristics into secvtx container
+ // for the moment, I only apply pair vertex chi2 cut
+ if (HFEpairs()->GetEntriesFast() == 1){
+ if (pair->GetKFChi2()<vtxchi2cut) { // you can also put single track cut
+ AliHFEsecVtxs hfesecvtx;
+ hfesecvtx.SetTrkLabel1(pair->GetTrkLabel());
+ hfesecvtx.SetTrkLabel2(-999);
+ hfesecvtx.SetInvmass(pair->GetInvmass());
+ hfesecvtx.SetKFChi2(pair->GetKFChi2());
+ hfesecvtx.SetSignedLxy(pair->GetSignedLxy());
+ hfesecvtx.SetKFIP(pair->GetKFIP());
+ AddHFEsecvtxToArray(&hfesecvtx);
+ fNoOfHFEsecvtxs++;
+
+ dataE[0]=pair->GetInvmass();
+ dataE[1]=pair->GetKFChi2();
+ dataE[2]=pair->GetSignedLxy();
+ dataE[3]=pair->GetKFIP();
+ if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
+ dataE[5]=2;
+ fSecvtxQA->Fill(dataE);
+ }
+ return;
+ }
+
+ // in case there are multiple paired track with the electron, calculate secvtx characteristics
+ // put the secvtx characteristics into container if it passes cuts
+ for (int i=0; i<HFEpairs()->GetEntriesFast()-1; i++){
+ for (int j=i+1; j<HFEpairs()->GetEntriesFast(); j++){
+ CalcSECVTXProperty(track, htrack[i], htrack[j]);
+ if (fKFchi2<vtxchi2cut) {
+ AliHFEsecVtxs hfesecvtx;
+ hfesecvtx.SetTrkLabel1(htracklabel[i]);
+ hfesecvtx.SetTrkLabel2(htracklabel[j]);
+ hfesecvtx.SetKFChi2(fKFchi2);
+ hfesecvtx.SetInvmass(fInvmass);
+ hfesecvtx.SetSignedLxy(fSignedLxy);
+ hfesecvtx.SetKFIP(fKFip);
+ AddHFEsecvtxToArray(&hfesecvtx);
+ fNoOfHFEsecvtxs++;
+
+ dataE[0]=pair->GetInvmass();
+ dataE[1]=pair->GetKFChi2();
+ dataE[2]=pair->GetSignedLxy();
+ dataE[3]=pair->GetKFIP();
+ if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
+ dataE[5]=3;
+ fSecvtxQA->Fill(dataE);
}
-
+ }
+ }
}
//_______________________________________________________________________________________________
-void AliHFEsecVtx::RunSECVTX(AliESDtrack *track)
+void AliHFEsecVtx::CalcSECVTXProperty(AliVTrack* track1, AliVTrack* track2, AliVTrack* track3)
{
- // run secondary vertexing algorithm and do tagging
+ //
+ // calculate secondary vertex properties
+ //
+
+ // get KF particle input pid
+ Int_t pdg1 = GetPDG(track1);
+ Int_t pdg2 = GetPDG(track2);
+ Int_t pdg3 = GetPDG(track3);
+
+ if(pdg1==-1 || pdg2==-1 || pdg3==-1) {
+ //printf("out if considered pid range \n");
+ return;
+ }
+
+ // create KF particle of pair
+ if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
+ else AliKFParticle::SetField(fESD1->GetMagneticField());
+ AliKFParticle kfTrack1(*track1, pdg1);
+ AliKFParticle kfTrack2(*track2, pdg2);
+ AliKFParticle kfTrack3(*track3, pdg3);
+
+ AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
+
+ //secondary vertex point from kf particle
+ Double_t kfx = kfSecondary.GetX();
+ Double_t kfy = kfSecondary.GetY();
+ //Double_t kfz = kfSecondary.GetZ();
- ResetTagVar();
+ //momentum at the decay point from kf particle
+ Double_t kfpx = kfSecondary.GetPx();
+ Double_t kfpy = kfSecondary.GetPy();
+ //Double_t kfpz = kfSecondary.GetPz();
- Int_t imclabel = TMath::Abs(track->GetLabel());
- if(imclabel<0) return;
- TParticle* mcpart = fStack->Particle(imclabel);
- Int_t esource = GetElectronSource(imclabel);
- if (esource == kDirectBeauty || esource == kBeautyCharm || esource == kBeautyGamma || esource == kBeautyPi0 || esource == kBeautyElse){
- fHistTagged.fPtBeautyElec->Fill(mcpart->Pt());
- fBElec = kTRUE;
- }
+ Double_t dx = kfx-fPVx;
+ Double_t dy = kfy-fPVy;
+ // discriminating variables ----------------------------------------------------------
- if (fPairTagged >= 4) {
- FindSECVTXCandid4Tracks(track);
- }
- else if (fPairTagged == 3) {
- FindSECVTXCandid3Tracks(track);
- }
- else if (fPairTagged == 2) {
- FindSECVTXCandid2Tracks(track);
- }
- else if (fPairTagged == 1) {
- ApplyPairTagCut();
- }
+ if(kfSecondary.GetNDF()>0) fKFchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
+ // invariant mass of the KF particle
+ kfSecondary.GetMass(fInvmass,fInvmassSigma);
- if (fBTagged) {
- fHistTagged.fPtTaggedElec->Fill(mcpart->Pt());
- if (esource == kDirectBeauty || esource == kBeautyCharm || esource == kBeautyGamma || esource == kBeautyPi0 || esource == kBeautyElse){
- fHistTagged.fPtRightTaggedElec->Fill(mcpart->Pt());
- }
- else fHistTagged.fPtWrongTaggedElec->Fill(mcpart->Pt());
- }
+ // DCA from primary to e-h KF particle (impact parameter of KF particle)
+ Double_t vtx[2]={fPVx, fPVy};
+ fKFip = kfSecondary.GetDistanceFromVertexXY(vtx);
+ // projection of kf vertex vector to the kf momentum direction
+ Double_t psqr = kfpx*kfpx+kfpy*kfpy;
+ if(psqr>0) fSignedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
}
//_______________________________________________________________________________________________
-void AliHFEsecVtx::ApplyPairTagCut()
+Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track)
{
- // apply tagging cut for e-h pair
+ //
+ // return mc pid
+ //
- if (fBElec){
- fHistTagged.fInvmassBeautyElec2trkVtx->Fill(fpairedInvMass[0]);
- fHistTagged.fSignedLxyBeautyElec2trkVtx->Fill(fpairedSignedLxy[0]);
- fHistTagged.fChi2BeautyElec2trkVtx->Fill(fpairedChi2[0]);
- }
- else {
- fHistTagged.fInvmassNotBeautyElec2trkVtx->Fill(fpairedInvMass[0]);
- fHistTagged.fSignedLxyNotBeautyElec2trkVtx->Fill(fpairedSignedLxy[0]);
- fHistTagged.fChi2NotBeautyElec2trkVtx->Fill(fpairedChi2[0]);
- }
-
- if (fpairedChi2[0] > 2.0) return;
- if (fpairedInvMass[0] < 1.5) return;
- if (fpairedSignedLxy[0] < 0.05) return;
+ Int_t label = TMath::Abs(track->GetLabel());
+ TParticle* mcpart = fStack->Particle(label);
+ if ( !mcpart ) return 0;
+ Int_t pdgCode = mcpart->GetPdgCode();
- fBTagged = kTRUE;
+ return pdgCode;
}
//_______________________________________________________________________________________________
-Bool_t AliHFEsecVtx::ApplyPairTagCut(Int_t id)
+Int_t AliHFEsecVtx::GetPairOriginESD(AliESDtrack* trk1, AliESDtrack* trk2)
{
- // apply tagging cut for e-h pair of given indexed hadron
-
- if (fpairedChi2[id] > 2.0) return kFALSE;
- if (fpairedInvMass[id] < 1.5) return kFALSE;
- if (fpairedSignedLxy[id] < 0.05) return kFALSE;
-
- fBTagged = kTRUE;
- return kTRUE;
-
+ //
+ // return pdg code of the origin(source) of the pair
+ //
+ //
+ // ---*---*---*-----ancester A----- track1
+ // |____*______
+ // |______ track2
+ // => if they originated from same ancester,
+ // then return "the absolute value of pdg code of ancester A"
+ //
+ // ---*---*---B hadron-----ancester A----- track1
+ // |____*______
+ // |______ track2
+ // => if they originated from same ancester, and this ancester originally comes from B hadrons
+ // then return -1*"the absolute value of pdg code of ancester A"
+ //
+ // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
+ //
+
+ if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
+ TParticle* part1 = fStack->Particle(trk1->GetLabel());
+ TParticle* part2 = fStack->Particle(trk2->GetLabel());
+ TParticle* part2cp = part2;
+ if (!(part1) || !(part2)) return 0;
+
+ Int_t srcpdg = 0;
+
+ //if the two tracks' mother's label is same, get the mother info
+ //in case of charm, check if it originated from beauty
+ for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
+ Int_t label1 = part1->GetFirstMother();
+ if (label1 < 0) return 0;
+
+ for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
+ Int_t label2 = part2->GetFirstMother();
+ if (label2 < 0) break;
+
+ if (label1 == label2){ //check if two tracks are originated from same mother
+ TParticle* commonmom = fStack->Particle(label2);
+ srcpdg = abs(commonmom->GetPdgCode());
+
+ //check ancester to see if it is originally from beauty
+ for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
+ Int_t ancesterlabel = commonmom->GetFirstMother();
+ if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg
+
+ TParticle* commonancester = fStack->Particle(ancesterlabel);
+ Int_t ancesterpdg = abs(commonancester->GetPdgCode());
+
+ for (Int_t l=0; l<fNparents; l++){
+ if (abs(ancesterpdg)==fParentSelect[1][l]){
+ srcpdg = -1*srcpdg; //multiply -1 for hadron from bottom
+ return srcpdg;
+ }
+ }
+ commonmom = commonancester;
+ }
+ }
+ part2 = fStack->Particle(label2); //if their mother is different, go to earlier generation of 2nd particle
+ if (!(part2)) break;
+ }
+ part1 = fStack->Particle(label1); //if their mother is different, go to earlier generation of 1st particle
+ part2 = part2cp;
+ if (!(part1)) return 0;
+ }
+
+ return srcpdg;
}
//_______________________________________________________________________________________________
-Bool_t AliHFEsecVtx::ApplyTagCut(Double_t chi2)
+Int_t AliHFEsecVtx::GetPairOriginAOD(AliAODTrack* trk1, AliAODTrack* trk2)
{
- // apply tagging cut for secondary vertex
- if (chi2 > 2.0) return kFALSE;
- if (finvmass < 1.5) return kFALSE;
- if (fsignedLxy < 0.05) return kFALSE;
- if (fcosPhi < 0.90) return kFALSE;
- if (fdistTwoSecVtx > 0.1) return kFALSE;
-
- fBTagged = kTRUE;
- return kTRUE;
+ //
+ // return pdg code of the origin(source) of the pair
+ //
+ //
+ // ---*---*---*-----ancester A----- track1
+ // |____*______
+ // |______ track2
+ // => if they originated from same ancester,
+ // then return "the absolute value of pdg code of ancester A"
+ //
+ // ---*---*---B hadron-----ancester A----- track1
+ // |____*______
+ // |______ track2
+ // => if they originated from same ancester, and this ancester originally comes from B hadrons
+ // then return -1*"the absolute value of pdg code of ancester A"
+ //
+ // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
+ //
+
+ if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
+ AliAODMCParticle *part1 = (AliAODMCParticle*)fMCArray->At(trk1->GetLabel());
+ AliAODMCParticle *part2 = (AliAODMCParticle*)fMCArray->At(trk2->GetLabel());
+ AliAODMCParticle *part2cp = part2;
+ if (!(part1) || !(part2)) return 0;
+
+ Int_t srcpdg = 0;
+
+ //if the two tracks' mother's label is same, get the mother info
+ //in case of charm, check if it originated from beauty
+ for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
+ Int_t label1 = part1->GetMother();
+ if (label1 < 0) return 0;
+
+ for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
+ Int_t label2 = part2->GetMother();
+ if (label2 < 0) return 0;
+
+ if (label1 == label2){ //check if two tracks are originated from same mother
+ AliAODMCParticle *commonmom = (AliAODMCParticle*)fMCArray->At(label1);
+ srcpdg = abs(commonmom->GetPdgCode());
+
+ //check ancester to see if it is originally from beauty
+ for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
+ Int_t ancesterlabel = commonmom->GetMother();
+ if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg
+
+ AliAODMCParticle *commonancester = (AliAODMCParticle*)fMCArray->At(ancesterlabel);
+ Int_t ancesterpdg = abs(commonancester->GetPdgCode());
+
+ for (Int_t l=0; l<fNparents; l++){
+ if (abs(ancesterpdg)==fParentSelect[1][l]){
+ srcpdg = -1*srcpdg; //multiply -1 for charm from bottom
+ return srcpdg;
+ }
+ }
+ commonmom = commonancester;
+ }
+ }
+ part2 = (AliAODMCParticle*)fMCArray->At(label2); //if their mother is different, go to earlier generation of 2nd particle
+ if (!(part2)) break;
+ }
+ part1 = (AliAODMCParticle*)fMCArray->At(label1); //if their mother is different, go to earlier generation of 2nd particle
+ part2 = part2cp;
+ if (!(part1)) return 0;
+ }
+
+ return srcpdg;
}
//_______________________________________________________________________________________________
-void AliHFEsecVtx::ApplyVtxTagCut(Double_t chi2, Int_t id1, Int_t id2)
+Int_t AliHFEsecVtx::GetPairCode(const AliVTrack* const trk1, const AliVTrack* const trk2)
{
- // apply tagging cut for e-h pair of given indexed hadron
+ //
+ // return pair code which is predefinded as:
+ // kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0,
+ // kElse, kBeautyGamma, kBeautyPi0, kBeautyElse
+ //
+
+ Int_t srcpdg = -1;
+ Int_t srccode = kElse;
+
+ if(IsAODanalysis()) srcpdg = GetPairOriginAOD((AliAODTrack*)trk1,(AliAODTrack*)trk2);
+ else srcpdg = GetPairOriginESD((AliESDtrack*)trk1,(AliESDtrack*)trk2);
+
+ if (srcpdg < 0) srccode = kBeautyElse;
+ for (Int_t i=0; i<fNparents; i++){
+ if (abs(srcpdg)==fParentSelect[0][i]) {
+ if (srcpdg>0) srccode = kDirectCharm;
+ if (srcpdg<0) srccode = kBeautyCharm;
+ }
+ if (abs(srcpdg)==fParentSelect[1][i]) {
+ if (srcpdg>0) srccode = kDirectBeauty;
+ if (srcpdg<0) return kElse;
+ }
+ }
+ if (srcpdg == 22) srccode = kGamma;
+ if (srcpdg == -22) srccode = kBeautyGamma;
+ if (srcpdg == 111) srccode = kPi0;
+ if (srcpdg == -111) srccode = kBeautyPi0;
- if(!ApplyTagCut(chi2)){
- if(!ApplyPairTagCut(id1)) ApplyPairTagCut(id2);
- }
-
+ return srccode;
}
//_______________________________________________________________________________________________
-void AliHFEsecVtx::FindSECVTXCandid4Tracks(AliESDtrack *track)
+Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack)
{
- // find secondary vertex for >= 4 e-h pairs
+ //
+ // return decay electron's origin
+ //
- // sort pair in increasing order (kFALSE-increasing order)
- Int_t index[20];
- Int_t indexA[4];
- Int_t indexB[3];
- Double_t sevchi2[4];
- AliESDtrack *htrack[4];
+ if (iTrack < 0) {
+ AliDebug(1, "Stack label is negative, return\n");
+ return -1;
+ }
- if(fPairTagged>20) return; // protection
+ TParticle* mcpart = fStack->Particle(iTrack);
- TMath::Sort(fPairTagged,fpairedChi2,index,kFALSE);
+// if ( abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !
- // select 4 partner tracks retruning smallest pair chi2
+ Int_t iLabel = mcpart->GetFirstMother();
+ if (iLabel<0){
+ AliDebug(1, "Stack label is negative, return\n");
+ return -1;
+ }
- for (Int_t i=0; i<4; i++){
- htrack[i] = fESD1->GetTrack(fpairedTrackID[index[i]]);
- }
-
- // calculated chi2 with 1 electron and 3 partner tracks
- for (Int_t i=0; i<4; i++){
- sevchi2[i] = GetSecVtxChi2(track, htrack[fid[i][0]], htrack[fid[i][1]], htrack[fid[i][2]]);
- }
+ Int_t origin = -1;
+ Bool_t isFinalOpenCharm = kFALSE;
- // select 3 partner tracks retruning smallest pair chi2
- // [think] if two smallest chi2 are similar, have to think about better handling of selection
- TMath::Sort(4,sevchi2,indexA,kFALSE);
+ TParticle *partMother = fStack->Particle(iLabel);
+ Int_t maPdgcode = partMother->GetPdgCode();
- // calculated chi2 with 1 electron and 2 partner tracks
- for (Int_t i=0; i<3; i++){
- sevchi2[i] = GetSecVtxChi2(track, htrack[fia[indexA[0]][i][0]], htrack[fia[indexA[0]][i][1]]);
- }
+ // if the mother is charmed hadron
+ if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
- // select 2 partner tracks retruning smallest pair chi2
- TMath::Sort(3,sevchi2,indexB,kFALSE);
+ for (Int_t i=0; i<fNparents; i++){
+ if (abs(maPdgcode)==fParentSelect[0][i]){
+ isFinalOpenCharm = kTRUE;
+ }
+ }
+ if (!isFinalOpenCharm) return -1;
- // calculate secondary vertex quality variables with 1 electron and 2 hadrons
- CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
+ // iterate until you find B hadron as a mother or become top ancester
+ for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
- if (fBElec){
- fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
- }
- else {
- fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
+ Int_t jLabel = partMother->GetFirstMother();
+ if (jLabel == -1){
+ origin = kDirectCharm;
+ return origin;
+ }
+ if (jLabel < 0){ // safety protection
+ AliDebug(1, "Stack label is negative, return\n");
+ return -1;
+ }
+
+ // if there is an ancester
+ TParticle* grandMa = fStack->Particle(jLabel);
+ Int_t grandMaPDG = grandMa->GetPdgCode();
+
+ for (Int_t j=0; j<fNparents; j++){
+ if (abs(grandMaPDG)==fParentSelect[1][j]){
+ origin = kBeautyCharm;
+ return origin;
}
+ }
+
+ partMother = grandMa;
+ } // end of iteration
+ } // end of if
+ else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
+ for (Int_t i=0; i<fNparents; i++){
+ if (abs(maPdgcode)==fParentSelect[1][i]){
+ origin = kDirectBeauty;
+ return origin;
+ }
+ }
+ } // end of if
- ApplyVtxTagCut(sevchi2[indexB[0]],index[fia[indexA[0]][indexB[0]][0]],index[fia[indexA[0]][indexB[0]][1]]);
+ //============ gamma ================
+ else if ( abs(maPdgcode) == 22 ) {
+ origin = kGamma;
-}
-
-//_______________________________________________________________________________________________
-void AliHFEsecVtx::FindSECVTXCandid3Tracks(AliESDtrack *track)
-{
- // find secondary vertex for 3 e-h pairs
+ // iterate until you find B hadron as a mother or become top ancester
+ for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
- // sort pair in increasing order (kFALSE-increasing order)
- Int_t indexA[1] = { 0 };
- Int_t indexB[3];
- Double_t sevchi2[3];
- AliESDtrack *htrack[3];
+ Int_t jLabel = partMother->GetFirstMother();
+ if (jLabel == -1){
+ origin = kGamma;
+ return origin;
+ }
+ if (jLabel < 0){ // safety protection
+ AliDebug(1, "Stack label is negative, return\n");
+ return -1;
+ }
+
+ // if there is an ancester
+ TParticle* grandMa = fStack->Particle(jLabel);
+ Int_t grandMaPDG = grandMa->GetPdgCode();
+
+ for (Int_t j=0; j<fNparents; j++){
+ if (abs(grandMaPDG)==fParentSelect[1][j]){
+ origin = kBeautyGamma;
+ return origin;
+ }
+ }
- // select 4 partner tracks retruning smallest pair chi2
+ partMother = grandMa;
+ } // end of iteration
- for (Int_t i=0; i<3; i++){
- htrack[i] = fESD1->GetTrack(fpairedTrackID[i]);
- }
-
- // calculated chi2 with 1 electron and 2 partner tracks
- for (Int_t i=0; i<3; i++){
- sevchi2[i] = GetSecVtxChi2(track, htrack[fia[indexA[0]][i][0]], htrack[fia[indexA[0]][i][1]]);
- }
+ return origin;
+ } // end of if
- // select 2 partner tracks retruning smallest pair chi2
- TMath::Sort(3,sevchi2,indexB,kFALSE);
+ //============ pi0 ================
+ else if ( abs(maPdgcode) == 111 ) {
+ origin = kPi0;
- // calculate secondary vertex quality variables with 1 electron and 2 hadrons
- CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
+ // iterate until you find B hadron as a mother or become top ancester
+ for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
- if (fBElec){
- fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
- }
- else {
- fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
+ Int_t jLabel = partMother->GetFirstMother();
+ if (jLabel == -1){
+ origin = kPi0;
+ return origin;
+ }
+ if (jLabel < 0){ // safety protection
+ AliDebug(1, "Stack label is negative, return\n");
+ return -1;
+ }
+
+ // if there is an ancester
+ TParticle* grandMa = fStack->Particle(jLabel);
+ Int_t grandMaPDG = grandMa->GetPdgCode();
+
+ for (Int_t j=0; j<fNparents; j++){
+ if (abs(grandMaPDG)==fParentSelect[1][j]){
+ origin = kBeautyPi0;
+ return origin;
}
+ }
- ApplyVtxTagCut(sevchi2[indexB[0]],fia[indexA[0]][indexB[0]][0],fia[indexA[0]][indexB[0]][1]);
-}
+ partMother = grandMa;
+ } // end of iteration
-//_______________________________________________________________________________________________
-void AliHFEsecVtx::FindSECVTXCandid2Tracks(AliESDtrack *track)
-{
- // find secondary vertex for 2 e-h pairs
+ return origin;
+ } // end of if
- Double_t sevchi2[1];
- AliESDtrack *htrack[2];
+ else {
+ origin = kElse;
+ // iterate until you find B hadron as a mother or become top ancester
+ for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
- for (Int_t i=0; i<2; i++){
- htrack[i] = fESD1->GetTrack(fpairedTrackID[i]);
+ Int_t jLabel = partMother->GetFirstMother();
+ if (jLabel == -1){
+ origin = kElse;
+ return origin;
+ }
+ if (jLabel < 0){ // safety protection
+ AliDebug(1, "Stack label is negative, return\n");
+ return -1;
+ }
+
+ // if there is an ancester
+ TParticle* grandMa = fStack->Particle(jLabel);
+ Int_t grandMaPDG = grandMa->GetPdgCode();
+
+ for (Int_t j=0; j<fNparents; j++){
+ if (abs(grandMaPDG)==fParentSelect[1][j]){
+ origin = kBeautyElse;
+ return origin;
}
-
- sevchi2[0] = GetSecVtxChi2(track, htrack[0], htrack[1]);
-
- // calculate secondary vertex quality variables with 1 electron and 2 hadrons
- CalcSECVTXProperty(track,htrack[0],htrack[1]);
+ }
- if (fBElec){
- fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[0]);
- }
- else {
- fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[0]);
- }
+ partMother = grandMa;
+ } // end of iteration
+ }
- ApplyVtxTagCut(sevchi2[0],0,1);
+ return origin;
}
//_______________________________________________________________________________________________
-void AliHFEsecVtx::CalcSECVTXProperty(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
+Int_t AliHFEsecVtx::GetPDG(AliVTrack *track)
{
- // calculate secondary vertex properties
-
- Int_t pdg1 = GetMCPID(track1);
- Int_t pdg2 = GetMCPID(track2);
- Int_t pdg3 = GetMCPID(track3);
-
- AliKFParticle::SetField(fESD1->GetMagneticField());
- AliKFParticle kfTrack1(*track1, pdg1);
- AliKFParticle kfTrack2(*track2, pdg2);
- AliKFParticle kfTrack3(*track3, pdg3);
-
- AliKFParticle kfSecondary12(kfTrack1,kfTrack2);
- AliKFParticle kfSecondary13(kfTrack1,kfTrack3);
- AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
-
- // copy primary vertex from ESD
- AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
- //printf("primary ndf= %d\n",primVtxCopy.GetNDF());
- if( primVtxCopy.GetNDF() <1 ) return;
-
- Double_t kdx12 = kfSecondary12.GetX()-primVtxCopy.GetX();
- Double_t kdy12 = kfSecondary12.GetY()-primVtxCopy.GetY();
- //Double_t kdz12 = kfSecondary12.GetZ()-primVtxCopy.GetZ();
-
- Double_t kdx13 = kfSecondary13.GetX()-primVtxCopy.GetX();
- Double_t kdy13 = kfSecondary13.GetY()-primVtxCopy.GetY();
- //Double_t kdz13 = kfSecondary13.GetZ()-primVtxCopy.GetZ();
-
- Double_t kdx = kfSecondary.GetX()-primVtxCopy.GetX();
- Double_t kdy = kfSecondary.GetY()-primVtxCopy.GetY();
- //Double_t kdz = kfSecondary.GetZ()-primVtxCopy.GetZ();
-
- // calculate distance and angle between two secvtxes
- fdistTwoSecVtx = TMath::Sqrt((kdx12-kdx13)*(kdx12-kdx13) + (kdy12-kdy13)*(kdy12-kdy13));
- fcosPhi = ( kdx12*kdx13 + kdy12*kdy13 ) / ( TMath::Sqrt(kdx12*kdx12+kdy12*kdy12)*TMath::Sqrt(kdx13*kdx13+kdy13*kdy13) );
- //Double_t lengthdiff = TMath::Abs(TMath::Sqrt(kdx12*kdx12+kdy12*kdy12) - TMath::Sqrt(kdx13*kdx13+kdy13*kdy13));
-
- // calculate angle between secondary vertex vector and secondary particle momentum vector in transverse plane
- Double_t cosTheta = ( kdx*kfSecondary.GetPx() + kdy*kfSecondary.GetPy()) / TMath::Sqrt(kdx*kdx+kdy*kdy)*TMath::Sqrt(kfSecondary.GetPx()*kfSecondary.GetPx()+kfSecondary.GetPy()*kfSecondary.GetPy());
- // calculate signed Lxy
- fsignedLxy = TMath::Sqrt(kdx*kdx+kdy*kdy)*cosTheta;
-
- // calculate invariant mass of the kf particle
- kfSecondary.GetMass(finvmass,finvmassSigma);
-
- if (fBElec){
- fHistTagged.fInvmassBeautyElecSecVtx->Fill(finvmass);
- fHistTagged.fSignedLxyBeautyElecSecVtx->Fill(fsignedLxy);
- fHistTagged.fDistTwoVtxBeautyElecSecVtx->Fill(fdistTwoSecVtx);
- fHistTagged.fCosPhiBeautyElecSecVtx->Fill(fcosPhi);
- }
- else {
- fHistTagged.fInvmassNotBeautyElecSecVtx->Fill(finvmass);
- fHistTagged.fSignedLxyNotBeautyElecSecVtx->Fill(fsignedLxy);
- fHistTagged.fDistTwoVtxNotBeautyElecSecVtx->Fill(fdistTwoSecVtx);
- fHistTagged.fCosPhiNotBeautyElecSecVtx->Fill(fcosPhi);
- }
+ //
+ // get KF particle input pdg for mass hypothesis
+ //
+
+ Int_t pdgCode=-1;
+
+ if (fUseMCPID && HasMCData()){
+ pdgCode = GetMCPDG(track);
+ }
+ else if(fESD1){
+ Int_t pid=0;
+ Double_t prob;
+ GetESDPID((AliESDtrack*)track, pid, prob);
+ switch(pid){
+ case 0: pdgCode = 11; break;
+ case 1: pdgCode = 13; break;
+ case 2: pdgCode = 211; break;
+ case 3: pdgCode = 321; break;
+ case 4: pdgCode = 2212; break;
+ default: pdgCode = -1;
+ }
+ }
+ else if(fAOD1){
+ Int_t pid = ((AliAODTrack*)track)->GetMostProbablePID();
+ switch(pid){
+ case 0: pdgCode = 11; break;
+ case 1: pdgCode = 13; break;
+ case 2: pdgCode = 211; break;
+ case 3: pdgCode = 321; break;
+ case 4: pdgCode = 2212; break;
+ default: pdgCode = -1;
+ }
+ }
+ return pdgCode;
}
//_______________________________________________________________________________________________
-Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track)
+Int_t AliHFEsecVtx::GetMCPDG(AliVTrack *track)
{
- // return mc pid
-
- Int_t label = TMath::Abs(track->GetLabel());
- TParticle* mcpart = fStack->Particle(label);
- if ( !mcpart ) return 0;
- Int_t pdgCode = mcpart->GetPdgCode();
-
- return pdgCode;
+ //
+ // return mc pdg code
+ //
+
+ Int_t label = TMath::Abs(track->GetLabel());
+ Int_t pdgCode;
+
+ if (IsAODanalysis()) {
+ AliAODMCParticle *mcpart = (AliAODMCParticle*)fMCArray->At(label);
+ if ( !mcpart ) return 0;
+ pdgCode = mcpart->GetPdgCode();
+ }
+ else {
+ TParticle* mcpart = fStack->Particle(label);
+ if ( !mcpart ) return 0;
+ pdgCode = mcpart->GetPdgCode();
+ }
+
+ return pdgCode;
}
-//_______________________________________________________________________________________________
-Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3, AliESDtrack* track4)
+//______________________________________________________________________________
+void AliHFEsecVtx::GetESDPID(AliESDtrack *track, Int_t &recpid, Double_t &recprob)
{
- // return 4 track secondary vertex chi2
-
- Int_t pdg1 = GetMCPID(track1);
- Int_t pdg2 = GetMCPID(track2);
- Int_t pdg3 = GetMCPID(track3);
- Int_t pdg4 = GetMCPID(track4);
+ //
+ // calculate likehood for esd pid
+ //
- AliKFParticle::SetField(fESD1->GetMagneticField());
- AliKFParticle kfTrack1(*track1, pdg1);
- AliKFParticle kfTrack2(*track2, pdg2);
- AliKFParticle kfTrack3(*track3, pdg3);
- AliKFParticle kfTrack4(*track4, pdg4);
- AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3,kfTrack4);
+ recpid = -1;
+ recprob = -1;
- return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
+ Int_t ipid=-1;
+ Double_t max=0.;
-}
-
-//_______________________________________________________________________________________________
-Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
-{
- // return 3 track secondary vertex chi2
+ Double_t probability[5];
- Int_t pdg1 = GetMCPID(track1);
- Int_t pdg2 = GetMCPID(track2);
- Int_t pdg3 = GetMCPID(track3);
+ // get probability of the diffenrent particle types
+ track->GetESDpid(probability);
- AliKFParticle::SetField(fESD1->GetMagneticField());
- AliKFParticle kfTrack1(*track1, pdg1);
- AliKFParticle kfTrack2(*track2, pdg2);
- AliKFParticle kfTrack3(*track3, pdg3);
- AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
-
- return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
+ // find most probable particle in ESD pid
+ // 0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
+ ipid = TMath::LocMax(5,probability);
+ max = TMath::MaxElement(5,probability);
+ recpid = ipid;
+ recprob = max;
}
-//_______________________________________________________________________________________________
-Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2)
+//_____________________________________________________________________________
+void AliHFEsecVtx::AddHFEpairToArray(const AliHFEpairs* const pair)
{
- // return 2 track secondary vertex chi2
-
- Int_t pdg1 = GetMCPID(track1);
- Int_t pdg2 = GetMCPID(track2);
-
- AliKFParticle::SetField(fESD1->GetMagneticField());
- AliKFParticle kfTrack1(*track1, pdg1);
- AliKFParticle kfTrack2(*track2, pdg2);
- AliKFParticle kfSecondary(kfTrack1,kfTrack2);
-
- return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
+ //
+ // Add a HFE pair to the array
+ //
+ Int_t n = HFEpairs()->GetEntriesFast();
+ if(n!=fNoOfHFEpairs)AliError(Form("fNoOfHFEpairs != HFEpairs()->GetEntriesFast %i != %i \n", fNoOfHFEpairs, n));
+ new((*HFEpairs())[n]) AliHFEpairs(*pair);
}
-//_______________________________________________________________________________________________
-Int_t AliHFEsecVtx::PairOrigin(AliESDtrack* track1, AliESDtrack* track2)
+//_____________________________________________________________________________
+TClonesArray *AliHFEsecVtx::HFEpairs()
{
-
- //
- // return pdg code of the origin(source) of the pair
- //
- //
- // ---*---*---*-----ancester A----- track1
- // |____*______
- // |______ track2
- // => if they originated from same ancester,
- // then return "the absolute value of pdg code of ancester A"
- //
- // ---*---*---B hadron-----ancester A----- track1
- // |____*______
- // |______ track2
- // => if they originated from same ancester, and this ancester originally comes from B hadrons
- // then return -1*"the absolute value of pdg code of ancester A"
- //
- // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
- //
-
- if (track1->GetLabel()<0 || track2->GetLabel()<0) return 0;
- TParticle* part1 = fStack->Particle(TMath::Abs(track1->GetLabel()));
- TParticle* part2 = fStack->Particle(TMath::Abs(track2->GetLabel()));
- if (!(part1)) return 0;
- if (!(part2)) return 0;
-
- TParticle* part1Crtgen = part1; // copy track into current generation particle
- TParticle* part2Crtgen = part2; // copy track into current generation particle
-
-
- Int_t sourcePDG = 0;
-
- // if the two tracks' mother's label is same, get the mother info
- // in case of charm, check if it originated from beauty
- for (Int_t i=0; i<100; i++){ // iterate 100
- Int_t iLabel = part1Crtgen->GetFirstMother(); //label of mother of current generation for 1st partilce
- if (iLabel < 0) return 0;
-
- for (Int_t j=0; j<100; j++){ // iterate 100
- Int_t jLabel = part2Crtgen->GetFirstMother(); //label of mother of current generation for 2nd partilce
- if (jLabel < 0) return 0; // if jLabel == -1
-
- if (iLabel == jLabel){ // check if two tracks are originated from same mother
- TParticle* thismother = fStack->Particle(jLabel); // if yes, get "thismother" info
- sourcePDG = abs(thismother->GetPdgCode()); // get the pdg code of "this mother"
-
- // check ancester to see if it is originally from beauty
- for (Int_t k=0; k<10; k++){ // check up to 10 ancesters
- Int_t ancesterLabel = thismother->GetFirstMother();
- if (ancesterLabel < 0) return sourcePDG; // if "thismoter" doesn't have mother anymore, return thismother's pdg
-
- TParticle* thisancester = fStack->Particle(ancesterLabel);
- Int_t ancesterPDG = abs(thisancester->GetPdgCode());
-
- for (Int_t l=0; l<fNparents; l++){
- if (abs(ancesterPDG)==fParentSelect[1][l]){
- sourcePDG = -1*sourcePDG; // multiply -1 for charm from bottom
- return sourcePDG;
- }
- }
- thismother = thisancester;
- }
-
- }
- part2Crtgen = fStack->Particle(jLabel); // if their mother is different, go up to previous generation of 2nd particle
- }
- part1Crtgen = fStack->Particle(iLabel); // if their mother is different, go up to previous generation of 2nd particle
-
- // if you don't find additionional(2nd particle) track originated from a given beauty hadron, break to save time
- Int_t motherPDGtmp = abs(part1Crtgen->GetPdgCode());
- for (Int_t l=0; l<fNparents; l++){
- if (abs(motherPDGtmp)==fParentSelect[1][l]) return sourcePDG;
- }
-
- }
-
- return sourcePDG;
-
+ //
+ // Returns the list of HFE pairs
+ //
+
+ if (!fHFEpairs) {
+ fHFEpairs = new TClonesArray("AliHFEpairs", 200);
+ }
+ return fHFEpairs;
}
-//_______________________________________________________________________________________________
-Int_t AliHFEsecVtx::PairCode(AliESDtrack* track1,AliESDtrack* track2)
+//_____________________________________________________________________________
+void AliHFEsecVtx::DeleteHFEpairs()
{
-
- //
- // return pair code which is predefinded as:
- // kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0, kElse, kBeautyGamma, kBeautyPi0, kBeautyElse
- //
-
- Int_t pairOriginsPDG = PairOrigin(track1,track2);
-
- Int_t sourcePart = kElse;
-
- if (pairOriginsPDG < 0) {
- sourcePart = kBeautyElse;
- }
- for (Int_t i=0; i<fNparents; i++){
- if (abs(pairOriginsPDG)==fParentSelect[0][i]) {
- if (pairOriginsPDG>0) sourcePart = kDirectCharm;
- if (pairOriginsPDG<0) {
- sourcePart = kBeautyCharm;
- }
- }
- if (abs(pairOriginsPDG)==fParentSelect[1][i]) {
- if (pairOriginsPDG>0) {
- sourcePart = kDirectBeauty;
- }
- if (pairOriginsPDG<0) return kElse;
- }
- }
- if (pairOriginsPDG == 22) sourcePart = kGamma;
- if (pairOriginsPDG == -22) {
- sourcePart = kBeautyGamma;
- }
- if (pairOriginsPDG == 111) sourcePart = kPi0;
- if (pairOriginsPDG == -111) {
- sourcePart = kBeautyPi0;
- }
-
- return sourcePart;
-
+ //
+ // Delete the list of HFE pairs
+ //
+
+ if (fHFEpairs) {
+ fHFEpairs->Delete();
+ //delete fHFEpairs;
+ }
}
-//_______________________________________________________________________________________________
-Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack)
+//_____________________________________________________________________________
+void AliHFEsecVtx::InitHFEpairs()
{
+ //
+ // Initialization should be done before make all possible pairs for a given electron candidate
+ //
- // return decay electron's origin
-
- if (iTrack < 0) {
- AliDebug(1, "Stack label is negative, return\n");
- return -1;
- }
-
- TParticle* mcpart = fStack->Particle(iTrack);
-
- if ( abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !!!!!!!!!!!!!!!!!
-
- Int_t iLabel = mcpart->GetFirstMother();
- if (iLabel<0){
- AliDebug(1, "Stack label is negative, return\n");
- return -1;
- }
-
- Int_t origin = -1;
- Bool_t isFinalOpenCharm = kFALSE;
-
- TParticle *partMother = fStack->Particle(iLabel);
- Int_t maPdgcode = partMother->GetPdgCode();
-
- // if the mother is charmed hadron
- if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
-
- for (Int_t i=0; i<fNparents; i++){
- if (abs(maPdgcode)==fParentSelect[0][i]){
- isFinalOpenCharm = kTRUE;
- }
- }
- if (!isFinalOpenCharm) return -1;
-
-
- // iterate until you find B hadron as a mother or become top ancester
- for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
-
- Int_t jLabel = partMother->GetFirstMother();
- if (jLabel == -1){
- origin = kDirectCharm;
- return origin;
- }
- if (jLabel < 0){ // safety protection
- AliDebug(1, "Stack label is negative, return\n");
- return -1;
- }
-
- // if there is an ancester
- TParticle* grandMa = fStack->Particle(jLabel);
- Int_t grandMaPDG = grandMa->GetPdgCode();
-
- for (Int_t j=0; j<fNparents; j++){
- if (abs(grandMaPDG)==fParentSelect[1][j]){
-
- origin = kBeautyCharm;
- return origin;
- }
- }
-
- partMother = grandMa;
- } // end of iteration
- } // end of if
- else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
- for (Int_t i=0; i<fNparents; i++){
- if (abs(maPdgcode)==fParentSelect[1][i]){
- origin = kDirectBeauty;
- return origin;
- }
- }
- } // end of if
-
-
- //============ gamma ================
- else if ( abs(maPdgcode) == 22 ) {
- origin = kGamma;
-
- // iterate until you find B hadron as a mother or become top ancester
- for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
-
- Int_t jLabel = partMother->GetFirstMother();
- if (jLabel == -1){
- origin = kGamma;
- return origin;
- }
- if (jLabel < 0){ // safety protection
- AliDebug(1, "Stack label is negative, return\n");
- return -1;
- }
-
- // if there is an ancester
- TParticle* grandMa = fStack->Particle(jLabel);
- Int_t grandMaPDG = grandMa->GetPdgCode();
-
- for (Int_t j=0; j<fNparents; j++){
- if (abs(grandMaPDG)==fParentSelect[1][j]){
- origin = kBeautyGamma;
- return origin;
- }
- }
-
- partMother = grandMa;
- } // end of iteration
-
- return origin;
- } // end of if
-
- //============ pi0 ================
- else if ( abs(maPdgcode) == 111 ) {
- origin = kPi0;
-
- // iterate until you find B hadron as a mother or become top ancester
- for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
-
- Int_t jLabel = partMother->GetFirstMother();
- if (jLabel == -1){
- origin = kPi0;
- return origin;
- }
- if (jLabel < 0){ // safety protection
- AliDebug(1, "Stack label is negative, return\n");
- return -1;
- }
-
- // if there is an ancester
- TParticle* grandMa = fStack->Particle(jLabel);
- Int_t grandMaPDG = grandMa->GetPdgCode();
-
- for (Int_t j=0; j<fNparents; j++){
- if (abs(grandMaPDG)==fParentSelect[1][j]){
- origin = kBeautyPi0;
- return origin;
- }
- }
-
- partMother = grandMa;
- } // end of iteration
-
- return origin;
- } // end of if
-
-
- else {
- origin = kElse;
-
- // iterate until you find B hadron as a mother or become top ancester
- for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
-
- Int_t jLabel = partMother->GetFirstMother();
- if (jLabel == -1){
- origin = kElse;
- return origin;
- }
- if (jLabel < 0){ // safety protection
- AliDebug(1, "Stack label is negative, return\n");
- return -1;
- }
-
- // if there is an ancester
- TParticle* grandMa = fStack->Particle(jLabel);
- Int_t grandMaPDG = grandMa->GetPdgCode();
-
- for (Int_t j=0; j<fNparents; j++){
- if (abs(grandMaPDG)==fParentSelect[1][j]){
- origin = kBeautyElse;
- return origin;
- }
- }
-
- partMother = grandMa;
- } // end of iteration
-
- }
-
- return origin;
-
+ fNoOfHFEpairs = 0;
}
-/*
-//_______________________________________________________________________________________________
-Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
+//_____________________________________________________________________________
+void AliHFEsecVtx::AddHFEsecvtxToArray(const AliHFEsecVtxs* const secvtx)
{
+ //
+ // Add a HFE secondary vertex to the array
+ //
- //
- // decay electron's origin
- //
-
- if (iTrackLabel < 0) {
- AliDebug(1, "Stack label is negative, return\n");
- return -1;
- }
-
- TParticle* mcpart = fStack->Particle(iTrackLabel);
- Int_t iLabel = mcpart->GetFirstMother();
- if (iLabel<0){
- AliDebug(1, "Stack label is negative, return\n");
- return -1;
- }
-
- Int_t origin = -1;
- Bool_t isFinalOpenCharm = kFALSE;
-
- TParticle *partMother = fStack->Particle(iLabel);
- Int_t maPdgcode = partMother->GetPdgCode(); // get mother's pdg code
-
- //beauty --------------------------
- if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
- for (Int_t i=0; i<fNparents; i++){
- if (abs(maPdgcode)==fParentSelect[1][i]){
- origin = kDirectBeauty;
- return origin;
- }
- else return -1; // this track is originated beauties not in the final B hadron list => excited beauties
- }
- } // end of if
+ Int_t n = HFEsecvtxs()->GetEntriesFast();
+ if(n!=fNoOfHFEsecvtxs)AliError(Form("fNoOfHFEsecvtxs != HFEsecvtxs()->GetEntriesFast %i != %i \n", fNoOfHFEsecvtxs, n));
+ new((*HFEsecvtxs())[n]) AliHFEsecVtxs(*secvtx);
+}
- //charm --------------------------
- else if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
-
- for (Int_t i=0; i<fNparents; i++){
- if (abs(maPdgcode)==fParentSelect[0][i])
- isFinalOpenCharm = kTRUE;
- }
- if (!isFinalOpenCharm) return -1; // this track is originated charms not in the final D hadron list => excited charms
- // to prevent any possible double counting
-
- for (Int_t i=0; i<100; i++){ // iterate 100 until you find B hadron as a mother or become top ancester
-
- Int_t jLabel = partMother->GetFirstMother();
- if (jLabel == -1){
- origin = kDirectCharm;
- return origin;
- }
- if (jLabel < 0){ // safety protection even though not really necessary here
- AliDebug(1, "Stack label is negative, return\n");
- return -1;
- }
-
- // if there is an ancester, check if it in the final B hadron list
- TParticle* grandMa = fStack->Particle(jLabel);
- Int_t grandMaPDG = grandMa->GetPdgCode();
-
- for (Int_t j=0; j<fNparents; j++){
- if (abs(grandMaPDG)==fParentSelect[1][j]){
- origin = kBeautyCharm;
- return origin;
- }
- }
-
- partMother = grandMa;
- } // end of iteration
- } // end of if
-
- //gamma --------------------------
- else if ( abs(maPdgcode) == 22 ) {
- origin = kGamma;
-
- // iterate until you find B hadron as a mother or become top ancester
- for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
-
- Int_t jLabel = partMother->GetFirstMother();
- if (jLabel == -1){
- origin = kGamma;
- return origin;
- }
- if (jLabel < 0){ // safety protection
- AliDebug(1, "Stack label is negative, return\n");
- return -1;
- }
-
- // if there is an ancester
- TParticle* grandMa = fStack->Particle(jLabel);
- Int_t grandMaPDG = grandMa->GetPdgCode();
-
- for (Int_t j=0; j<fNparents; j++){
- if (abs(grandMaPDG)==fParentSelect[1][j]){
- origin = kBeautyGamma;
- return origin;
- }
- }
-
- partMother = grandMa;
- } // end of iteration
-
- return origin;
- } // end of if
-
- //pi0 --------------------------
- else if ( abs(maPdgcode) == 111 ) {
- origin = kPi0;
-
- // iterate until you find B hadron as a mother or become top ancester
- for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
-
- Int_t jLabel = partMother->GetFirstMother();
- if (jLabel == -1){
- origin = kPi0;
- return origin;
- }
- if (jLabel < 0){ // safety protection
- AliDebug(1, "Stack label is negative, return\n");
- return -1;
- }
-
- // if there is an ancester
- TParticle* grandMa = fStack->Particle(jLabel);
- Int_t grandMaPDG = grandMa->GetPdgCode();
-
- for (Int_t j=0; j<fNparents; j++){
- if (abs(grandMaPDG)==fParentSelect[1][j]){
- origin = kBeautyPi0;
- return origin;
- }
- }
-
- partMother = grandMa;
- } // end of iteration
-
- return origin;
- } // end of if
-
-
- //else --------------------------
- else {
- origin = kElse;
-
- // iterate until you find B hadron as a mother or become top ancester
- for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
-
- Int_t jLabel = partMother->GetFirstMother();
- if (jLabel == -1){
- origin = kElse;
- return origin;
- }
- if (jLabel < 0){ // safety protection
- AliDebug(1, "Stack label is negative, return\n");
- return -1;
- }
-
- // if there is an ancester
- TParticle* grandMa = fStack->Particle(jLabel);
- Int_t grandMaPDG = grandMa->GetPdgCode();
-
- for (Int_t j=0; j<fNparents; j++){
- if (abs(grandMaPDG)==fParentSelect[1][j]){
- origin = kBeautyElse;
- return origin;
- }
- }
-
- partMother = grandMa;
- } // end of iteration
+//_____________________________________________________________________________
+TClonesArray *AliHFEsecVtx::HFEsecvtxs()
+{
+ //
+ // Returns the list of HFE secvtx
+ //
+
+ if (!fHFEsecvtxs) {
+ fHFEsecvtxs = new TClonesArray("AliHFEsecVtxs", 200);
+ }
+ return fHFEsecvtxs;
+}
- }
+//_____________________________________________________________________________
+void AliHFEsecVtx::DeleteHFEsecvtxs()
+{
+ //
+ // Delete the list of HFE pairs
+ //
+
+ if (fHFEsecvtxs) {
+ fHFEsecvtxs->Delete();
+ //delete fHFEsecvtx;
+ }
+}
- return origin;
+//_____________________________________________________________________________
+void AliHFEsecVtx::InitHFEsecvtxs()
+{
+ //
+ // Initialization should be done
+ //
+ fNoOfHFEsecvtxs = 0;
}
-*/
//_______________________________________________________________________________________________
Bool_t AliHFEsecVtx::SingleTrackCut(AliESDtrack* track) const
{
- // test cuts
-
- //if (track->Pt() < 1.0) return kFALSE;
- //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE;
- //if (!(track->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
- //if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
- if (!(TESTBIT(track->GetITSClusterMap(),0))) return kFALSE; // ask hit on the first pixel layer
- //if (!(TESTBIT(track->GetITSClusterMap(),0) | TESTBIT(track->GetITSClusterMap(),1))) return kFALSE;
-
+ //if (track->Pt() < 1.0) return kFALSE;
+ //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE;
+ //if (!(track->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
+ //if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
+ if (!(TESTBIT(track->GetITSClusterMap(),0))) return kFALSE; // ask hit on the first pixel layer
+ //if (!(TESTBIT(track->GetITSClusterMap(),0) | TESTBIT(track->GetITSClusterMap(),1))) return kFALSE;
/*
- Float_t dcaR=-1;
- Float_t dcaZ=-1;
- track->GetImpactParameters(dcaR,dcaZ);
- if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) < 0.005) return kFALSE;
- if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) > 0.3) return kFALSE;
+ Float_t dcaR=-1;
+ Float_t dcaZ=-1;
+ track->GetImpactParameters(dcaR,dcaZ);
+ if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) < 0.005) return kFALSE;
+ if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) > 0.3) return kFALSE;
*/
- return kTRUE;
+ return kTRUE;
+}
+//____________________________________________________________
+void AliHFEsecVtx::MakeContainer(){
+
+ //
+ // make container
+ //
+
+ const Int_t nDimPair=6;
+ Int_t nBinPair[nDimPair] = {200, 500, 314, 2000, 2000, 11};
+ const Double_t kInvmassmin = 0., kInvmassmax = 20.;
+ const Double_t kKFChi2min = 0, kKFChi2max= 50;
+ const Double_t kOpenanglemin = 0, kOpenanglemax = 3.14;
+ const Double_t kSignedLxymin = -10, kSignedLxymax= 10;
+ const Double_t kKFIPmin = -10, kKFIPmax= 10;
+ const Double_t kPairCodemin = -1, kPairCodemax= 10;
+
+ Double_t* binEdgesPair[nDimPair];
+ for(Int_t ivar = 0; ivar < nDimPair; ivar++)
+ binEdgesPair[ivar] = new Double_t[nBinPair[ivar] + 1];
+
+ for(Int_t i=0; i<=nBinPair[0]; i++) binEdgesPair[0][i]=(Double_t)kInvmassmin + (kInvmassmax - kInvmassmin)/nBinPair[0]*(Double_t)i;
+ for(Int_t i=0; i<=nBinPair[1]; i++) binEdgesPair[1][i]=(Double_t)kKFChi2min + (kKFChi2max - kKFChi2min)/nBinPair[1]*(Double_t)i;
+ for(Int_t i=0; i<=nBinPair[2]; i++) binEdgesPair[2][i]=(Double_t)kOpenanglemin + (kOpenanglemax - kOpenanglemin)/nBinPair[2]*(Double_t)i;
+ for(Int_t i=0; i<=nBinPair[3]; i++) binEdgesPair[3][i]=(Double_t)kSignedLxymin + (kSignedLxymax - kSignedLxymin)/nBinPair[3]*(Double_t)i;
+ for(Int_t i=0; i<=nBinPair[4]; i++) binEdgesPair[4][i]=(Double_t)kKFIPmin + (kKFIPmax - kKFIPmin)/nBinPair[4]*(Double_t)i;
+ for(Int_t i=0; i<=nBinPair[5]; i++) binEdgesPair[5][i]=(Double_t)kPairCodemin + (kPairCodemax - kPairCodemin)/nBinPair[5]*(Double_t)i;
+
+ fPairQA = new THnSparseF("pairQA", "QA for Pair; invmass[GeV/c^2]; KF chi2; opening angle; signed Lxy; KF ip; pair code ", nDimPair, nBinPair);
+ for(Int_t idim = 0; idim < nDimPair; idim++){
+ fPairQA->SetBinEdges(idim, binEdgesPair[idim]);
+ }
+
+ fSecVtxList->AddAt(fPairQA,0);
+
+ const Int_t nDimSecvtx=6;
+ Double_t* binEdgesSecvtx[nDimSecvtx];
+ Int_t nBinSecvtx[nDimSecvtx] = {200, 500, 2000, 2000, 11, 3};
+ const Double_t kNtrksmin = 0, kNtrksmax= 3;
+ for(Int_t ivar = 0; ivar < nDimSecvtx; ivar++)
+ binEdgesSecvtx[ivar] = new Double_t[nBinSecvtx[ivar] + 1];
+
+ for(Int_t i=0; i<=nBinSecvtx[0]; i++) binEdgesSecvtx[0][i]=binEdgesPair[0][i];
+ for(Int_t i=0; i<=nBinSecvtx[1]; i++) binEdgesSecvtx[1][i]=binEdgesPair[1][i];
+ for(Int_t i=0; i<=nBinSecvtx[2]; i++) binEdgesSecvtx[2][i]=binEdgesPair[3][i];
+ for(Int_t i=0; i<=nBinSecvtx[3]; i++) binEdgesSecvtx[3][i]=binEdgesPair[4][i];
+ for(Int_t i=0; i<=nBinSecvtx[4]; i++) binEdgesSecvtx[4][i]=binEdgesPair[5][i];
+ for(Int_t i=0; i<=nBinSecvtx[5]; i++) binEdgesSecvtx[5][i]=(Double_t)kNtrksmin + (kNtrksmax - kNtrksmin)/nBinSecvtx[5]*(Double_t)i;
+
+ fSecvtxQA = new THnSparseF("secvtxQA", "QA for Secvtx; invmass[GeV/c^2]; KF chi2; signed Lxy; KF ip; pair code; n tracks ", nDimSecvtx, nBinSecvtx);
+ for(Int_t idim = 0; idim < nDimSecvtx; idim++){
+ fSecvtxQA->SetBinEdges(idim, binEdgesSecvtx[idim]);
+ }
+
+ fSecVtxList->AddAt(fSecvtxQA,1);
}
//#include <TObject.h>
#endif
+#ifndef ROOT_THnSparse
+#include <THnSparse.h>
+#endif
+
class TH1F;
class TH2F;
class TString;
class AliESDEvent;
+class AliAODEvent;
+class AliVTrack;
class AliESDtrack;
+class AliAODTrack;
class AliStack;
+class AliHFEpairs;
+class AliHFEsecVtxs;
//________________________________________________________________
class AliHFEsecVtx : public TObject {
AliHFEsecVtx &operator=(const AliHFEsecVtx &); // assignment operator
virtual ~AliHFEsecVtx();
- void CreateHistograms(TString hnopt="");
+ void CreateHistograms(TList *qaList);
+
+ Bool_t HasMCData() const { return TestBit(kHasMCData); };
+ Bool_t IsAODanalysis() const { return TestBit(kAODanalysis); };
+ Bool_t IsESDanalysis() const { return !TestBit(kAODanalysis); };
+
+ void SetHasMCData(Bool_t hasMCdata = kTRUE) { SetBit(kHasMCData,hasMCdata); };
+ void SetAODAnalysis() { SetBit(kAODanalysis, kTRUE); };
+ void SetESDAnalysis() { SetBit(kAODanalysis, kFALSE); };
+ void SetEvent(AliESDEvent* const ESD){fESD1=ESD;}; // set ESD pointer
+ void SetEventAOD(AliAODEvent* const AOD){fAOD1=AOD;}; // set ESD pointer
+ void SetStack(AliStack* const stack){fStack=stack;}; // set stack pointer
+ void SetMCArray(TClonesArray* const mcarry){fMCArray=mcarry;} // set mcarray pointer
+ void SetUseMCPID(Bool_t usemcpid){fUseMCPID=usemcpid;};
- void SetEvent(AliESDEvent* const ESD){fESD1=ESD;}; // set ESD pointer
- void SetStack(AliStack* const stack){fStack=stack;} // set stack pointer
- void InitAnaPair(); // initialize default parameters
- void AnaPair(AliESDtrack* ESDtrack1, AliESDtrack* ESDtrack2, Int_t index2); // do e-h analysis
- void RunSECVTX(AliESDtrack *track); // run secondary vertexing algorithm
Int_t GetMCPID(AliESDtrack *track); // return MC pid
- Int_t PairOrigin(AliESDtrack* track1, AliESDtrack* track2); // return pair origin as a pdg code
- Int_t PairCode(AliESDtrack* track1,AliESDtrack* track2); // return corresponding pair code to pdg code
+ Int_t GetMCPDG(AliVTrack *track); // return MC pid
+ Int_t GetPairOriginESD(AliESDtrack* track1, AliESDtrack* track2); // return pair origin as a pdg code
+ Int_t GetPairOriginAOD(AliAODTrack* track1, AliAODTrack* track2); // return pair origin as a pdg code
+ Int_t GetPairCode(const AliVTrack* const track1, const AliVTrack* const track2); // return corresponding pair code to pdg code
Int_t GetElectronSource(Int_t mclabel); // return origin of the electron
+ Int_t GetPDG(AliVTrack *track); // return pdg
+ void GetESDPID(AliESDtrack *track, Int_t &recpid, Double_t &recprob); //return esd pid likelihood
+ void GetPrimaryCondition();
- Bool_t SingleTrackCut(AliESDtrack* track1) const; // single track cut
+ TClonesArray *HFEpairs();
+ TClonesArray *HFEsecvtxs();
- protected:
+ void AddHFEpairToArray(const AliHFEpairs* const pair);
+ void AddHFEsecvtxToArray(const AliHFEsecVtxs* const secvtx);
- void Init();
- void FindSECVTXCandid4Tracks(AliESDtrack* track1); // find secondary vertex for 4 tracks
- void FindSECVTXCandid3Tracks(AliESDtrack* track1); // find secondary vertex for 3 tracks
- void FindSECVTXCandid2Tracks(AliESDtrack* track1); // find secondary vertex for 2 tracks
- void CalcSECVTXProperty(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3); // calculated distinctive variables
- void ApplyPairTagCut(); // apply strict tagging condition for 1 pair secondary vertex
- void ApplyVtxTagCut(Double_t chi2, Int_t id1, Int_t id2); // apply tagging condition for 3 track secondary vertex
- void ResetTagVar(); // reset tagging variables
+ void InitHFEpairs();
+ void InitHFEsecvtxs();
- Bool_t ApplyPairTagCut(Int_t id); // apply strict tagging condition for 1 pair secondary vertex
- Bool_t ApplyTagCut(Double_t chi2); // apply tagging condition
+ void DeleteHFEpairs();
+ void DeleteHFEsecvtxs();
+
+ Bool_t SingleTrackCut(AliESDtrack* track1) const; // single track cut
+ void PairAnalysis(AliVTrack* ESDtrack1, AliVTrack* ESDtrack2, Int_t index2); // do e-h analysis
+ void RunSECVTX(AliVTrack *track); // run secondary vertexing algorithm
+ void MakeContainer(); // make containers
- Double_t GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3, AliESDtrack* track4); // get secondary vertex chi2 for 4 tracks
- Double_t GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3); // get secondary vertex chi2 for 3 tracks
- Double_t GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2); // get secondary vertex chi2 for 2 tracks
+ protected:
+ void Init();
+ void FindSECVTXCandid(AliVTrack *track);
+ void CalcSECVTXProperty(AliVTrack* track1, AliVTrack* track2, AliVTrack* track3); // calculated distinctive variables
private:
+ enum{
+ kHasMCData = BIT(15), // bitset for mc data usage
+ kAODanalysis = BIT(16) // bitset for aod analysis
+ };
+ enum {kAll, kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0, kElse, kBeautyGamma, kBeautyPi0, kBeautyElse}; // electron origin
+ enum {kCharm=4, kBeauty=5}; // quark flavor
AliESDEvent* fESD1; // ESD pointer
- AliStack* fStack; // stack pointer
+ AliAODEvent* fAOD1; // AOD pointer
+ AliStack* fStack; // stack pointer
- Int_t fParentSelect[2][7]; // heavy hadron species
- Int_t fNparents; // number of heavy hadrons to be considered
+ Bool_t fUseMCPID; // if use MC pid
TString fkSourceLabel[10]; // electron source label
- enum {kAll, kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0, kElse, kBeautyGamma, kBeautyPi0, kBeautyElse};
- enum {kCharm=4, kBeauty=5};
-
- struct AliHistspair{
- TH2F *fInvMass; // histogram to store pair invariant mass
- TH2F *fInvMassCut1; // histogram to store pair invariant mass after cut1
- TH2F *fInvMassCut2; // histogram to store pair invariant mass after cut2
- TH1F *fKFChi2; // histogram to store pair vertex chi2
- TH1F *fOpenAngle; // histogram to store pair opening angle
- TH1F *fCosOpenAngle; // histogram to store cosine of the pair opening angle
- TH2F *fSignedLxy; // histogram to store signed Lxy
- TH1F *fKFIP; // histogram to store pair impact parameter
- TH1F *fIPMax; // histogram to store larger impact parameter among pair tracks
-
- AliHistspair()
- : fInvMass()
- , fInvMassCut1()
- , fInvMassCut2()
- , fKFChi2()
- , fOpenAngle()
- , fCosOpenAngle()
- , fSignedLxy()
- , fKFIP()
- , fIPMax()
- {
- // define constructor
- }
- AliHistspair(const AliHistspair & p)
- : fInvMass(p.fInvMass)
- , fInvMassCut1(p.fInvMassCut1)
- , fInvMassCut2(p.fInvMassCut2)
- , fKFChi2(p.fKFChi2)
- , fOpenAngle(p.fOpenAngle)
- , fCosOpenAngle(p.fCosOpenAngle)
- , fSignedLxy(p.fSignedLxy)
- , fKFIP(p.fKFIP)
- , fIPMax(p.fIPMax)
- {
- // copy constructor
- }
- AliHistspair &operator=(const AliHistspair &)
- {
- // assignment operator, not yet implemented
- return *this;
- }
- };
-
- struct AliHiststag{
- TH1F *fPtBeautyElec; // histogram for electrons of single track cut passed
- TH1F *fPtTaggedElec; // histogram for total b tagged electrons
- TH1F *fPtRightTaggedElec; // histogram for right b tagged electrons
- TH1F *fPtWrongTaggedElec; // histogram for wrong b tagged electrons
- TH1F *fInvmassBeautyElecSecVtx; // histogram for right-tagged b invariant mass
- TH1F *fInvmassNotBeautyElecSecVtx; // histogram for mis-tagged b invariant mass
- TH1F *fSignedLxyBeautyElecSecVtx; // histogram for right-tagged b signed Lxy
- TH1F *fSignedLxyNotBeautyElecSecVtx; // histogram for mis-tagged b signed Lxy
- TH1F *fDistTwoVtxBeautyElecSecVtx; // histogram for right-tagged b two vertex distance
- TH1F *fDistTwoVtxNotBeautyElecSecVtx; // histogram for mis-tagged b two vertex distance
- TH1F *fCosPhiBeautyElecSecVtx; // histogram for right-tagged b cos of opening angle
- TH1F *fCosPhiNotBeautyElecSecVtx; // histogram for mis-tagged b cos of opening angle
- TH1F *fChi2BeautyElecSecVtx; // histogram for right-tagged b chi2 of secondary vertex
- TH1F *fChi2NotBeautyElecSecVtx; // histogram for mis-tagged b chi2 of secondary vertex
- TH1F *fInvmassBeautyElec2trkVtx; // histogram for right-tagged b invariant mass of two track vertex
- TH1F *fInvmassNotBeautyElec2trkVtx; // histogram for mis-tagged b invariant mass of two track vertex
- TH1F *fSignedLxyBeautyElec2trkVtx; // histogram for right-tagged b signed Lxy of two track vertex
- TH1F *fSignedLxyNotBeautyElec2trkVtx; // histogram for mis-tagged b signed Lxy of two track vertex
- TH1F *fChi2BeautyElec2trkVtx; // histogram for right-tagged b chi2 of two track vertex
- TH1F *fChi2NotBeautyElec2trkVtx; // histogram for mis-tagged b chi2 of two track vertex
-
- AliHiststag()
- : fPtBeautyElec()
- , fPtTaggedElec()
- , fPtRightTaggedElec()
- , fPtWrongTaggedElec()
- , fInvmassBeautyElecSecVtx()
- , fInvmassNotBeautyElecSecVtx()
- , fSignedLxyBeautyElecSecVtx()
- , fSignedLxyNotBeautyElecSecVtx()
- , fDistTwoVtxBeautyElecSecVtx()
- , fDistTwoVtxNotBeautyElecSecVtx()
- , fCosPhiBeautyElecSecVtx()
- , fCosPhiNotBeautyElecSecVtx()
- , fChi2BeautyElecSecVtx()
- , fChi2NotBeautyElecSecVtx()
- , fInvmassBeautyElec2trkVtx()
- , fInvmassNotBeautyElec2trkVtx()
- , fSignedLxyBeautyElec2trkVtx()
- , fSignedLxyNotBeautyElec2trkVtx()
- , fChi2BeautyElec2trkVtx()
- , fChi2NotBeautyElec2trkVtx()
- {
- // define constructor
- }
- AliHiststag(const AliHiststag & p)
- : fPtBeautyElec(p.fPtBeautyElec)
- , fPtTaggedElec(p.fPtTaggedElec)
- , fPtRightTaggedElec(p.fPtRightTaggedElec)
- , fPtWrongTaggedElec(p.fPtWrongTaggedElec)
- , fInvmassBeautyElecSecVtx(p.fInvmassBeautyElecSecVtx)
- , fInvmassNotBeautyElecSecVtx(p.fInvmassNotBeautyElecSecVtx)
- , fSignedLxyBeautyElecSecVtx(p.fSignedLxyBeautyElecSecVtx)
- , fSignedLxyNotBeautyElecSecVtx(p.fSignedLxyNotBeautyElecSecVtx)
- , fDistTwoVtxBeautyElecSecVtx(p.fDistTwoVtxBeautyElecSecVtx)
- , fDistTwoVtxNotBeautyElecSecVtx(p.fDistTwoVtxNotBeautyElecSecVtx)
- , fCosPhiBeautyElecSecVtx(p.fCosPhiBeautyElecSecVtx)
- , fCosPhiNotBeautyElecSecVtx(p.fCosPhiNotBeautyElecSecVtx)
- , fChi2BeautyElecSecVtx(p.fChi2BeautyElecSecVtx)
- , fChi2NotBeautyElecSecVtx(p.fChi2NotBeautyElecSecVtx)
- , fInvmassBeautyElec2trkVtx(p.fInvmassBeautyElec2trkVtx)
- , fInvmassNotBeautyElec2trkVtx(p.fInvmassNotBeautyElec2trkVtx)
- , fSignedLxyBeautyElec2trkVtx(p.fSignedLxyBeautyElec2trkVtx)
- , fSignedLxyNotBeautyElec2trkVtx(p.fSignedLxyNotBeautyElec2trkVtx)
- , fChi2BeautyElec2trkVtx(p.fChi2BeautyElec2trkVtx)
- , fChi2NotBeautyElec2trkVtx(p.fChi2NotBeautyElec2trkVtx)
- {
- // copy constructor
- }
- AliHiststag &operator=(const AliHiststag &)
- {
- // assignment operator, not yet implemented
- return *this;
- }
- };
-
- AliHistspair fHistPair[10]; // struct of above given histogram for different electron sources
- AliHiststag fHistTagged; // struct of above given histogram
+ Int_t fNparents; // number of heavy hadrons to be considered
+ Int_t fParentSelect[2][7]; // heavy hadron species
- Int_t fPairTagged; // number of tagged e-h pairs
- Int_t fpairedTrackID[20]; // paird hadron track track
- Double_t fpairedChi2[20]; // pair chi2
- Double_t fpairedInvMass[20]; // pair invariant mass
- Double_t fpairedSignedLxy[20]; // pair signed Lxy
+ Int_t fNoOfHFEpairs; // number of e-h pairs
+ Int_t fNoOfHFEsecvtxs; // number of secondary vertexes
- Int_t fid[4][3]; // index to store sorting result
- Int_t fia[4][3][2]; // index to store sorting result
+ TClonesArray *fHFEpairs; //! Array of pair
+ TClonesArray *fHFEsecvtxs; //! Array of secondary vertexes
+ TClonesArray *fMCArray; //! mc array pointer
- Double_t fdistTwoSecVtx; // distance between two pair vertex
- Double_t fcosPhi; // cos of opening angle of two pair vertex
- Double_t fsignedLxy; // signed Lxy of secondary vertex
- Double_t finvmass; // invariant mass of secondary vertex
- Double_t finvmassSigma; // invariant mass sigma of secondary vertex
+ Double_t fPVx; // primary vertex copy x
+ Double_t fPVy; // primary vertex copy y
+ Double_t fCosPhi; // cos of opening angle of two pair vertex
+ Double_t fSignedLxy; // signed Lxy of secondary vertex
+ Double_t fKFchi2; // chi2 of secondary vertex
+ Double_t fInvmass; // invariant mass of secondary vertex
+ Double_t fInvmassSigma; // invariant mass sigma of secondary vertex
+ Double_t fKFip; // impact parameter of secondary vertex track
- Bool_t fBTagged; // if b tagged, set true
- Bool_t fBElec; // if set true for b electron, set true
+ THnSparseF *fPairQA; // qa histos for pair analysis
+ THnSparseF *fSecvtxQA; // qa histos for secvtx
+ TList *fSecVtxList; // list for secondary vertexing outputs
- ClassDef(AliHFEsecVtx,0);
+ ClassDef(AliHFEsecVtx,0);
};
#endif