#include <TList.h>
#include <TMath.h>
#include <TObject.h>
+#include <TGrid.h>
#include <AliKFParticle.h>
#include "AliDielectronSignalMC.h"
#include "AliDielectronMixingHandler.h"
#include "AliDielectronV0Cuts.h"
+#include "AliDielectronPID.h"
#include "AliDielectron.h"
TNamed("AliDielectron","AliDielectron"),
fCutQA(kFALSE),
fQAmonitor(0x0),
+ fPostPIDCntrdCorr(0x0),
+ fPostPIDWdthCorr(0x0),
fEventFilter("EventFilter"),
fTrackFilter("TrackFilter"),
fPairPreFilter("PairPreFilter"),
fPdgLeg2(11),
fSignalsMC(0x0),
fNoPairing(kFALSE),
+ fProcessLS(kTRUE),
fUseKF(kTRUE),
fHistoArray(0x0),
fHistos(0x0),
fEstimatorFilename(""),
fTRDpidCorrectionFilename(""),
fVZEROCalibrationFilename(""),
- fVZERORecenteringFilename("")
+ fVZERORecenteringFilename(""),
+ fEffMapFilename(""),
+ fZDCRecenteringFilename("")
+
{
//
// Default constructor
TNamed(name,title),
fCutQA(kFALSE),
fQAmonitor(0x0),
+ fPostPIDCntrdCorr(0x0),
+ fPostPIDWdthCorr(0x0),
fEventFilter("EventFilter"),
fTrackFilter("TrackFilter"),
fPairPreFilter("PairPreFilter"),
fPdgLeg2(11),
fSignalsMC(0x0),
fNoPairing(kFALSE),
+ fProcessLS(kTRUE),
fUseKF(kTRUE),
fHistoArray(0x0),
fHistos(0x0),
fEstimatorFilename(""),
fTRDpidCorrectionFilename(""),
fVZEROCalibrationFilename(""),
- fVZERORecenteringFilename("")
+ fVZERORecenteringFilename(""),
+ fEffMapFilename(""),
+ fZDCRecenteringFilename("")
{
//
// Named constructor
// Default destructor
//
if (fQAmonitor) delete fQAmonitor;
- if (fHistoArray) delete fHistoArray;
+ if (fPostPIDCntrdCorr) delete fPostPIDCntrdCorr;
+ if (fPostPIDWdthCorr) delete fPostPIDWdthCorr;
if (fHistos) delete fHistos;
if (fPairCandidates) delete fPairCandidates;
if (fDebugTree) delete fDebugTree;
if (fMixing) delete fMixing;
if (fSignalsMC) delete fSignalsMC;
if (fCfManagerPair) delete fCfManagerPair;
+ if (fHistoArray) delete fHistoArray;
}
//________________________________________________________________
fTrackRotator->SetPdgLegs(fPdgLeg1,fPdgLeg2);
}
if (fDebugTree) fDebugTree->SetDielectron(this);
- if(fEstimatorFilename.Contains(".root")) AliDielectronVarManager::InitEstimatorAvg(fEstimatorFilename.Data());
+
+ TString allfiles = fEstimatorFilename;
+ allfiles+=fTRDpidCorrectionFilename;
+ allfiles+=fVZEROCalibrationFilename;
+ allfiles+=fVZERORecenteringFilename;
+ allfiles+=fEffMapFilename;
+ allfiles+=fZDCRecenteringFilename;
+ if(allfiles.Contains("alien://")) TGrid::Connect("alien://",0,0,"t");
+
+ if(fEstimatorFilename.Contains(".root")) AliDielectronVarManager::InitEstimatorAvg(fEstimatorFilename.Data());
if(fTRDpidCorrectionFilename.Contains(".root")) AliDielectronVarManager::InitTRDpidEffHistograms(fTRDpidCorrectionFilename.Data());
if(fVZEROCalibrationFilename.Contains(".root")) AliDielectronVarManager::SetVZEROCalibrationFile(fVZEROCalibrationFilename.Data());
if(fVZERORecenteringFilename.Contains(".root")) AliDielectronVarManager::SetVZERORecenteringFile(fVZERORecenteringFilename.Data());
-
+ if(fEffMapFilename.Contains(".root")) AliDielectronVarManager::InitEffMap(fEffMapFilename.Data());
+ if(fZDCRecenteringFilename.Contains(".root")) AliDielectronVarManager::SetZDCRecenteringFile(fZDCRecenteringFilename.Data());
+
if (fMixing) fMixing->Init(this);
if (fHistoArray) {
fHistoArray->SetSignalsMC(fSignalsMC);
fHistoArray->Init();
}
+ if(fPostPIDCntrdCorr) AliDielectronPID::SetCentroidCorrFunction(fPostPIDCntrdCorr);
+ if(fPostPIDWdthCorr) AliDielectronPID::SetWidthCorrFunction(fPostPIDWdthCorr);
+
if (fCutQA) {
fQAmonitor = new AliDielectronCutQA(Form("QAcuts_%s",GetName()),"QAcuts");
fQAmonitor->AddTrackFilter(&fTrackFilter);
// create pairs and fill pair candidate arrays
for (Int_t itrackArr1=0; itrackArr1<4; ++itrackArr1){
for (Int_t itrackArr2=itrackArr1; itrackArr2<4; ++itrackArr2){
- FillPairArrays(itrackArr1, itrackArr2);
+ if(!fProcessLS && GetPairIndex(itrackArr1,itrackArr2)!=kEv1PM) continue;
+ FillPairArrays(itrackArr1, itrackArr2);
}
}
// FillHistograms(0x0,kTRUE);
}
+ // fill candidate variables
+ Double_t ntracks = fTracks[0].GetEntriesFast() + fTracks[1].GetEntriesFast();
+ Double_t npairs = PairArray(AliDielectron::kEv1PM)->GetEntriesFast();
+ AliDielectronVarManager::SetValue(AliDielectronVarManager::kTracks, ntracks);
+ AliDielectronVarManager::SetValue(AliDielectronVarManager::kPairs, npairs);
+
//in case there is a histogram manager, fill the QA histograms
if (fHistos && fSignalsMC) FillMCHistograms(ev1);
if (fHistos) FillHistograms(ev1);
if (fHistos) FillHistogramsMC(dieMC->GetMCEvent(), ev1);
+ // mc tracks
+ if(!dieMC->GetNMCTracks()) return;
+
+ // signals to be studied
if(!fSignalsMC) return;
- //loop over all MC data and Fill the CF container if it exist
- if (!fCfManagerPair && !fHistoArray) return;
+ Int_t nSignals = fSignalsMC->GetEntries();
+ if(!nSignals) return;
+
+ //loop over all MC data and Fill the HF, CF containers and histograms if they exist
if(fCfManagerPair) fCfManagerPair->SetPdgMother(fPdgMother);
- Bool_t bFillCF = (fCfManagerPair ? fCfManagerPair->GetStepForMCtruth() : kFALSE);
- Bool_t bFillHF = (fHistoArray ? fHistoArray->GetStepForMCGenerated() : kFALSE);
- if(!bFillCF && !bFillHF) return;
+ Bool_t bFillCF = (fCfManagerPair ? fCfManagerPair->GetStepForMCtruth() : kFALSE);
+ Bool_t bFillHF = (fHistoArray ? fHistoArray->GetStepForMCGenerated() : kFALSE);
+ Bool_t bFillHist = kFALSE;
+ if(fHistos) {
+ const THashList *histlist = fHistos->GetHistogramList();
+ for(Int_t isig=0;isig<nSignals;isig++) {
+ TString sigName = fSignalsMC->At(isig)->GetName();
+ bFillHist |= histlist->FindObject(Form("Pair_%s_MCtruth",sigName.Data()))!=0x0;
+ bFillHist |= histlist->FindObject(Form("Track_Leg_%s_MCtruth",sigName.Data()))!=0x0;
+ bFillHist |= histlist->FindObject(Form("Track_%s_%s_MCtruth",fgkPairClassNames[1],sigName.Data()))!=0x0;
+ if(bFillHist) break;
+ }
+ }
+ // check if there is anything to fill
+ if(!bFillCF && !bFillHF && !bFillHist) return;
- // signals to be studied
- Int_t nSignals = fSignalsMC->GetEntries();
// initialize 2D arrays of labels for particles from each MC signal
Int_t** labels1; // labels for particles satisfying branch 1
// Do the pairing and fill the CF container with pure MC info
for(Int_t isig=0; isig<nSignals; ++isig) {
+ // printf("INDEXES: %d-%d both%d\n",indexes1[isig],indexes2[isig],indexes12[isig]);
// mix the particles which satisfy only one of the signal branches
for(Int_t i1=0;i1<indexes1[isig];++i1) {
+ if(!indexes2[isig]) FillMCHistograms(labels1[isig][i1], -1, isig); // (e.g. single electrons only, no pairs)
for(Int_t i2=0;i2<indexes2[isig];++i2) {
if(bFillCF) fCfManagerPair->FillMC(labels1[isig][i1], labels2[isig][i2], isig);
if(bFillHF) fHistoArray->Fill(labels1[isig][i1], labels2[isig][i2], isig);
+ FillMCHistograms(labels1[isig][i1], labels2[isig][i2], isig);
}
}
// mix the particles which satisfy both branches
for(Int_t i2=0; i2<i1; ++i2) {
if(bFillCF) fCfManagerPair->FillMC(labels12[isig][i1], labels12[isig][i2], isig);
if(bFillHF) fHistoArray->Fill(labels12[isig][i1], labels12[isig][i2], isig);
+ FillMCHistograms(labels12[isig][i1], labels12[isig][i2], isig);
}
}
} // end loop over signals
// track mapping
TMap mapRemovedTracks;
- // eta gap ?
- Bool_t etagap = kFALSE;
- for (Int_t iCut=0; iCut<fEventPlanePreFilter.GetCuts()->GetEntries();++iCut) {
- TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName();
- if(cutName.Contains("eta") || cutName.Contains("Eta")) etagap=kTRUE;
- }
Double_t cQX=0., cQY=0.;
// apply cuts to the tracks, e.g. etagap
UInt_t cutMask=fPairPreFilterLegs.IsSelected(arrTracks1.UncheckedAt(itrack));
//apply cut
- if (cutMask!=selectedMask) arrTracks1.AddAt(0x0,itrack);;
+ if (cutMask!=selectedMask) arrTracks1.AddAt(0x0,itrack);
}
arrTracks1.Compress();
for (Int_t itrack2=0; itrack2<end; ++itrack2){
//create the pair
candidate->SetTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
- static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
+ static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
candidate->SetType(pairIndex);
Int_t label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother);
candidate->SetLabel(label);
if (label>-1) candidate->SetPdgCode(fPdgMother);
+ else candidate->SetPdgCode(0);
// check for gamma kf particle
label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,22);
//pair cuts
UInt_t cutMask=fPairFilter.IsSelected(candidate);
-
+
//CF manager for the pair
if (fCfManagerPair) fCfManagerPair->Fill(cutMask,candidate);
- //histogram array for the pair
- if (fHistoArray) fHistoArray->Fill(pairIndex,candidate);
+
// cut qa
- if(pairIndex==AliDielectron::kEv1PM && fCutQA) {
+ if(pairIndex==kEv1PM && fCutQA) {
fQAmonitor->FillAll(candidate);
fQAmonitor->Fill(cutMask,candidate);
}
//apply cut
if (cutMask!=selectedMask) continue;
+ //histogram array for the pair
+ if (fHistoArray) fHistoArray->Fill(pairIndex,candidate);
+
//add the candidate to the candidate array
PairArray(pairIndex)->Add(candidate);
//get a new candidate
//CF manager for the pair
if (fCfManagerPair) fCfManagerPair->Fill(cutMask,&candidate);
- //histogram array for the pair
- if (fHistoArray) fHistoArray->Fill((Int_t)kEv1PMRot,&candidate);
//apply cut
if (cutMask==selectedMask) {
- if(fHistos) FillHistogramsPair(&candidate);
- if(fStoreRotatedPairs) PairArray(kEv1PMRot)->Add(new AliDielectronPair(candidate));
- }
+
+ //histogram array for the pair
+ if (fHistoArray) fHistoArray->Fill((Int_t)kEv1PMRot,&candidate);
+
+ if(fHistos) FillHistogramsPair(&candidate);
+ if(fStoreRotatedPairs) PairArray(kEv1PMRot)->Add(new AliDielectronPair(candidate));
+ }
}
}
}
fSignalsMC->Add(signal);
}
+
+//________________________________________________________________
+void AliDielectron::FillMCHistograms(Int_t label1, Int_t label2, Int_t nSignal) {
+ //
+ // fill QA MC TRUTH histograms for pairs and legs of all added mc signals
+ //
+
+ TString className,className2,className3;
+ className.Form("Pair_%s_MCtruth",fSignalsMC->At(nSignal)->GetName());
+ className2.Form("Track_Legs_%s_MCtruth",fSignalsMC->At(nSignal)->GetName());
+ className3.Form("Track_%s_%s_MCtruth",fgkPairClassNames[1],fSignalsMC->At(nSignal)->GetName());
+ Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
+ Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
+ Bool_t trkClass=fHistos->GetHistogramList()->FindObject(className3.Data())!=0x0;
+ // printf("fill signal %d: pair %d legs %d trk %d \n",nSignal,pairClass,legClass,trkClass);
+ if(!pairClass && !legClass && !trkClass) return;
+
+ // printf("leg labels: %d-%d \n",label1,label2);
+ AliVParticle* part1 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label1);
+ AliVParticle* part2 = AliDielectronMC::Instance()->GetMCTrackFromMCEvent(label2);
+ if(!part1 && !part2) return;
+ if(part1&&part2) {
+ // fill only unlike sign (and only SE)
+ if(part1->Charge()*part2->Charge()>=0) return;
+ }
+
+
+ AliDielectronMC* dieMC = AliDielectronMC::Instance();
+
+ Int_t mLabel1 = dieMC->GetMothersLabel(label1); // should work for both ESD and AOD
+ Int_t mLabel2 = dieMC->GetMothersLabel(label2);
+
+ // check the same mother option
+ AliDielectronSignalMC* sigMC = (AliDielectronSignalMC*)fSignalsMC->At(nSignal);
+ if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kSame && mLabel1!=mLabel2) return;
+ if(sigMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent && mLabel1==mLabel2) return;
+
+ // fill event values
+ Double_t values[AliDielectronVarManager::kNMaxValues];
+ AliDielectronVarManager::Fill(dieMC->GetMCEvent(), values); // get event informations
+
+ // fill the leg variables
+ // printf("leg:%d trk:%d part1:%p part2:%p \n",legClass,trkClass,part1,part2);
+ if (legClass || trkClass) {
+ if(part1) AliDielectronVarManager::Fill(part1,values);
+ if(part1 && trkClass) fHistos->FillClass(className3, AliDielectronVarManager::kNMaxValues, values);
+ if(part1 && part2 && legClass) fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
+ if(part2) AliDielectronVarManager::Fill(part2,values);
+ if(part2 && trkClass) fHistos->FillClass(className3, AliDielectronVarManager::kNMaxValues, values);
+ if(part1 && part2 && legClass) fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
+ }
+
+ //fill pair information
+ if (pairClass && part1 && part2) {
+ AliDielectronVarManager::FillVarMCParticle2(part1,part2,values);
+ fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
+ }
+
+}
+
//________________________________________________________________
void AliDielectron::FillMCHistograms(const AliVEvent *ev) {
//
//
if (!fSignalsMC) return;
- TString className,className2;
+ TString className,className2,className3;
Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
AliDielectronVarManager::Fill(ev, values); // get event informations
//loop over all added mc signals
for(Int_t isig=0; isig<fSignalsMC->GetEntries(); isig++) {
+ //check if and what to fill
className.Form("Pair_%s",fSignalsMC->At(isig)->GetName());
className2.Form("Track_Legs_%s",fSignalsMC->At(isig)->GetName());
+ className3.Form("Track_%s_%s",fgkPairClassNames[1],fSignalsMC->At(isig)->GetName()); // unlike sign, SE only
Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
- if(!pairClass && !legClass) return;
+ Bool_t mergedtrkClass=fHistos->GetHistogramList()->FindObject(className3.Data())!=0x0;
+ if(!pairClass && !legClass && !mergedtrkClass) continue;
+
+ // fill pair and/or their leg variables
+ if(pairClass || legClass) {
+ Int_t npairs=PairArray(AliDielectron::kEv1PM)->GetEntriesFast(); // only SE +-
+ for (Int_t ipair=0; ipair<npairs; ++ipair){
+ AliDielectronPair *pair=static_cast<AliDielectronPair*>(PairArray(AliDielectron::kEv1PM)->UncheckedAt(ipair));
+
+ Bool_t isMCtruth = AliDielectronMC::Instance()->IsMCTruth(pair, (AliDielectronSignalMC*)fSignalsMC->At(isig));
+ if(isMCtruth) {
+ //fill pair information
+ if (pairClass){
+ AliDielectronVarManager::Fill(pair, values);
+ fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
+ }
+ //fill leg information, both + and - in the same histo
+ if (legClass){
+ AliDielectronVarManager::Fill(pair->GetFirstDaughter(),values);
+ fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
+ AliDielectronVarManager::Fill(pair->GetSecondDaughter(),values);
+ fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
+ }
+ } //is signal
+ } //loop: pairs
+ }
- Int_t ntracks=PairArray(AliDielectron::kEv1PM)->GetEntriesFast(); // only SE +-
- for (Int_t ipair=0; ipair<ntracks; ++ipair){
- AliDielectronPair *pair=static_cast<AliDielectronPair*>(PairArray(AliDielectron::kEv1PM)->UncheckedAt(ipair));
+ // fill single tracks of signals
+ if(!mergedtrkClass) continue;
+ // loop over SE track arrays
+ for (Int_t i=0; i<2; ++i){
+ Int_t ntracks=fTracks[i].GetEntriesFast();
+ for (Int_t itrack=0; itrack<ntracks; ++itrack){
+ Int_t label=((AliVParticle*)fTracks[i].UncheckedAt(itrack))->GetLabel();
+ Bool_t isMCtruth1 = AliDielectronMC::Instance()->IsMCTruth(label, (AliDielectronSignalMC*)fSignalsMC->At(isig), 1);
+ Bool_t isMCtruth2 = AliDielectronMC::Instance()->IsMCTruth(label, (AliDielectronSignalMC*)fSignalsMC->At(isig), 2);
+ // skip if track does not correspond to the signal
+ if(!isMCtruth1 && !isMCtruth2) continue;
+ AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
+ fHistos->FillClass(className3, AliDielectronVarManager::kNMaxValues, values);
+ } //loop: tracks
+ } //loop: arrays
- Bool_t isMCtruth = AliDielectronMC::Instance()->IsMCTruth(pair, (AliDielectronSignalMC*)fSignalsMC->At(isig));
- if(isMCtruth) {
- //fill pair information
- if (pairClass){
- AliDielectronVarManager::Fill(pair, values);
- fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
- }
- //fill leg information, both + and - in the same histo
- if (legClass){
- AliDielectronVarManager::Fill(pair->GetFirstDaughter(),values);
- fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
- AliDielectronVarManager::Fill(pair->GetSecondDaughter(),values);
- fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values);
- }
- } //is signal
- } //loop: pairs
} //loop: MCsignals
}
+
+//______________________________________________
+void AliDielectron::SetCentroidCorrFunction(TF1 *fun, UInt_t varx, UInt_t vary, UInt_t varz)
+{
+ fun->GetHistogram()->GetXaxis()->SetUniqueID(varx);
+ fun->GetHistogram()->GetYaxis()->SetUniqueID(vary);
+ fun->GetHistogram()->GetZaxis()->SetUniqueID(varz);
+ fPostPIDCntrdCorr=fun;
+}
+//______________________________________________
+void AliDielectron::SetWidthCorrFunction(TF1 *fun, UInt_t varx, UInt_t vary, UInt_t varz)
+{
+ fun->GetHistogram()->GetXaxis()->SetUniqueID(varx);
+ fun->GetHistogram()->GetYaxis()->SetUniqueID(vary);
+ fun->GetHistogram()->GetZaxis()->SetUniqueID(varz);
+ fPostPIDWdthCorr=fun;
+}