#include <TMath.h>
#include <TObject.h>
-#include <AliESDEvent.h>
-#include <AliESDtrack.h>
#include <AliKFParticle.h>
#include <AliEventplane.h>
fPdgLeg2(11),
fSignalsMC(0x0),
fNoPairing(kFALSE),
+ fUseKF(kTRUE),
+ fHistoArray(0x0),
fHistos(0x0),
fPairCandidates(new TObjArray(11)),
fCfManagerPair(0x0),
fPreFilterAllSigns(kFALSE),
fHasMC(kFALSE),
fStoreRotatedPairs(kFALSE),
+ fDontClearArrays(kFALSE),
fEstimatorFilename(""),
- fTRDpidCorrectionFilename("")
+ fTRDpidCorrectionFilename(""),
+ fVZEROCalibrationFilename(""),
+ fVZERORecenteringFilename("")
{
//
// Default constructor
fPdgLeg2(11),
fSignalsMC(0x0),
fNoPairing(kFALSE),
+ fUseKF(kTRUE),
+ fHistoArray(0x0),
fHistos(0x0),
fPairCandidates(new TObjArray(11)),
fCfManagerPair(0x0),
fPreFilterAllSigns(kFALSE),
fHasMC(kFALSE),
fStoreRotatedPairs(kFALSE),
+ fDontClearArrays(kFALSE),
fEstimatorFilename(""),
- fTRDpidCorrectionFilename("")
+ fTRDpidCorrectionFilename(""),
+ fVZEROCalibrationFilename(""),
+ fVZERORecenteringFilename("")
{
//
// Named constructor
//
// Default destructor
//
+ if (fHistoArray) delete fHistoArray;
if (fHistos) delete fHistos;
if (fPairCandidates) delete fPairCandidates;
if (fDebugTree) delete fDebugTree;
if (fDebugTree) fDebugTree->SetDielectron(this);
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 (fMixing) fMixing->Init(this);
+ if (fHistoArray) {
+ fHistoArray->SetSignalsMC(fSignalsMC);
+ fHistoArray->Init();
+ }
}
//________________________________________________________________
return;
}
AliDielectronVarManager::SetEvent(ev1);
+ if (fMixing){
+ //set mixing bin to event data
+ Int_t bin=fMixing->FindBin(AliDielectronVarManager::GetData());
+ AliDielectronVarManager::SetValue(AliDielectronVarManager::kMixingBin,bin);
+ }
//in case we have MC load the MC event and process the MC particles
- if (AliDielectronMC::Instance()->HasMC()) {
- if (!AliDielectronMC::Instance()->ConnectMCEvent()){
- AliError("Could not properly connect the MC event, skipping this event!");
- return;
- }
+ if (AliDielectronMC::Instance()->ConnectMCEvent()){
ProcessMC(ev1);
}
if ((ev1&&fEventFilter.IsSelected(ev1)!=selectedMask) ||
(ev2&&fEventFilter.IsSelected(ev2)!=selectedMask)) return;
- AliDielectronVarManager::SetEvent(ev1);
+// AliDielectronVarManager::SetEvent(ev1); // why a second time???
//fill track arrays for the first event
if (ev1){
// TPC event plane correction
AliEventplane *cevplane = new AliEventplane();
- if (ev1 && fPreFilterEventPlane && ( fEventPlanePreFilter.GetCuts()->GetEntries()>0 || fEventPlanePOIPreFilter.GetCuts()->GetEntries()>0))
+ if (ev1 && cevplane && fPreFilterEventPlane && ( fEventPlanePreFilter.GetCuts()->GetEntries()>0 || fEventPlanePOIPreFilter.GetCuts()->GetEntries()>0))
EventPlanePreFilter(0, 1, fTracks[0], fTracks[1], ev1, cevplane);
-
+
if (!fNoPairing){
// create pairs and fill pair candidate arrays
for (Int_t itrackArr1=0; itrackArr1<4; ++itrackArr1){
}
//in case there is a histogram manager, fill the QA histograms
+ if (fHistos && fSignalsMC) FillMCHistograms(ev1);
if (fHistos) FillHistograms(ev1);
+
// clear arrays
- ClearArrays();
+ if (!fDontClearArrays) ClearArrays();
AliDielectronVarManager::SetTPCEventPlane(0x0);
delete cevplane;
}
if(!fSignalsMC) return;
//loop over all MC data and Fill the CF container if it exist
- if (!fCfManagerPair) return;
- fCfManagerPair->SetPdgMother(fPdgMother);
- if(!fCfManagerPair->GetStepForMCtruth()) return;
+ if (!fCfManagerPair && !fHistoArray) return;
+ if(fCfManagerPair) fCfManagerPair->SetPdgMother(fPdgMother);
+
+ Bool_t bFillCF = (fCfManagerPair ? fCfManagerPair->GetStepForMCtruth() : kFALSE);
+ Bool_t bFillHF = (fHistoArray ? fHistoArray->GetStepForMCGenerated() : kFALSE);
+ if(!bFillCF && !bFillHF) return;
// signals to be studied
Int_t nSignals = fSignalsMC->GetEntries();
// mix the particles which satisfy only one of the signal branches
for(Int_t i1=0;i1<indexes1[isig];++i1) {
for(Int_t i2=0;i2<indexes2[isig];++i2) {
- fCfManagerPair->FillMC(labels1[isig][i1], labels2[isig][i2], isig);
+ if(bFillCF) fCfManagerPair->FillMC(labels1[isig][i1], labels2[isig][i2], isig);
+ if(bFillHF) fHistoArray->Fill(labels1[isig][i1], labels2[isig][i2], isig);
}
}
// mix the particles which satisfy both branches
for(Int_t i1=0;i1<indexes12[isig];++i1) {
for(Int_t i2=0; i2<i1; ++i2) {
- fCfManagerPair->FillMC(labels12[isig][i1], labels12[isig][i2], isig);
+ if(bFillCF) fCfManagerPair->FillMC(labels12[isig][i1], labels12[isig][i2], isig);
+ if(bFillHF) fHistoArray->Fill(labels12[isig][i1], labels12[isig][i2], isig);
}
}
} // end loop over signals
TString className,className2;
Double_t values[AliDielectronVarManager::kNMaxValues]={0.};
//Fill event information
- if (ev){
- AliDielectronVarManager::Fill(ev, values);
+ if (ev){ //TODO: Why not use GetData() ??? See below event plane stuff!!!
+ AliDielectronVarManager::Fill(ev, values); //data should already be stored in AliDielectronVarManager from SetEvent, does EV plane correction rely on this???
+ if (fMixing){
+ //set mixing bin to event data
+ Int_t bin=fMixing->FindBin(values);
+ values[AliDielectronVarManager::kMixingBin]=bin;
+ }
+
if (fHistos->GetHistogramList()->FindObject("Event"))
- fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, values);
+// fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, AliDielectronVarManager::GetData());
+ fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, values); //check event plane stuff and replace with above...
}
//Fill track information, separately for the track array candidates
for (Int_t itrack=0; itrack<ntracks; ++itrack){
//get particle
AliVParticle *particle=ev->GetTrack(itrack);
- //TODO: temporary solution, perhaps think about a better implementation
- // This is needed to use AliESDpidCuts, which relies on the ESD event
- // is set as a AliESDtrack attribute... somehow ugly!
- if (ev->IsA()==AliESDEvent::Class()){
- AliESDtrack *track=static_cast<AliESDtrack*>(particle);
- track->SetESDEvent(static_cast<AliESDEvent*>(ev)); //only in trunk...
- }
-
+
//apply track cuts
if (fTrackFilter.IsSelected(particle)!=selectedMask) continue;
// remove contribution of all tracks to the Q-vector that are in invariant mass window
//
AliEventplane *evplane = const_cast<AliVEvent *>(ev)->GetEventplane();
+// AliEventplane *evplane = ev->GetEventplane();
if(!evplane) return;
// do not change these vectors qref
- TVector2 *qref = evplane->GetQVector();
+ TVector2 * const qref = evplane->GetQVector();
if(!qref) return;
// random subevents
TVector2 *qrsub1 = evplane->GetQsub1();
TVector2 *qrsub2 = evplane->GetQsub2();
// copy references
- TVector2 *qstd = dynamic_cast<TVector2 *>(qref->Clone());
- TVector2 *qssub1 = dynamic_cast<TVector2 *>(qrsub1->Clone());
- TVector2 *qssub2 = dynamic_cast<TVector2 *>(qrsub2->Clone());
+ TVector2 *qcorr = new TVector2(*qref);
+ TVector2 *qcsub1 = 0x0;
+ TVector2 *qcsub2 = 0x0;
+ // printf("qrsub1 %p %f \n",qrsub1,qrsub1->X());
- // subevents by charge separation sub1+ sub2-
- if(fLikeSignSubEvents) {
- qssub1 = dynamic_cast<TVector2 *>(qstd->Clone());
- qssub2 = dynamic_cast<TVector2 *>(qstd->Clone());
+
+ // 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;
+ }
+
+ // subevent separation
+ if(fLikeSignSubEvents || etagap) {
+ qcsub1 = new TVector2(*qcorr);
+ qcsub2 = new TVector2(*qcorr);
Int_t ntracks=ev->GetNumberOfTracks();
Double_t cQX = evplane->GetQContributionX(track);
Double_t cQY = evplane->GetQContributionY(track);
- Short_t charge=track->Charge();
- if (charge<0) qssub1->Set(qssub1->X()-cQX, qssub1->Y()-cQY);
- if (charge>0) qssub2->Set(qssub2->X()-cQX, qssub2->Y()-cQY);
- }
- }
-
- // subevents eta division sub1+ sub2-
- 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;
-
- qssub1 = dynamic_cast<TVector2 *>(qstd->Clone());
- qssub2 = dynamic_cast<TVector2 *>(qstd->Clone());
-
- Int_t ntracks=ev->GetNumberOfTracks();
-
- // track removals
- for (Int_t itrack=0; itrack<ntracks; ++itrack){
- AliVParticle *particle=ev->GetTrack(itrack);
- AliVTrack *track= static_cast<AliVTrack*>(particle);
- if (!track) continue;
-
- Double_t cQX = evplane->GetQContributionX(track);
- Double_t cQY = evplane->GetQContributionY(track);
-
+ // by charge sub1+ sub2-
+ if(fLikeSignSubEvents) {
+ Short_t charge=track->Charge();
+ if (charge<0) qcsub1->Set(qcsub1->X()-cQX, qcsub1->Y()-cQY);
+ if (charge>0) qcsub2->Set(qcsub2->X()-cQX, qcsub2->Y()-cQY);
+ }
+ // by eta sub1+ sub2-
+ if(etagap) {
Double_t eta=track->Eta();
- if (eta<0.0) qssub1->Set(qssub1->X()-cQX, qssub1->Y()-cQY);
- if (eta>0.0) qssub2->Set(qssub2->X()-cQX, qssub2->Y()-cQY);
+ if (eta<0.0) qcsub1->Set(qcsub1->X()-cQX, qcsub1->Y()-cQY);
+ if (eta>0.0) qcsub2->Set(qcsub2->X()-cQX, qcsub2->Y()-cQY);
}
}
}
+ else {
+ // by a random
+ qcsub1 = new TVector2(*qrsub1);
+ qcsub2 = new TVector2(*qrsub2);
+ }
// apply cuts, e.g. etagap
if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
// update Q vectors
- qstd->Set(qstd->X()-cQX, qstd->Y()-cQY);
- qssub1->Set(qssub1->X()-cQXsub1, qssub1->Y()-cQYsub1);
- qssub2->Set(qssub2->X()-cQXsub2, qssub2->Y()-cQYsub2);
+ qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
+ qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
+ qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
}
}
- if(!qstd || !qssub1 || !qssub2) return;
-
- TVector2 *qcorr = dynamic_cast<TVector2 *>(qstd->Clone());
- TVector2 *qcsub1 = dynamic_cast<TVector2 *>(qssub1->Clone());
- TVector2 *qcsub2 = dynamic_cast<TVector2 *>(qssub2->Clone());
-
-
// POI (particle of interest) rejection
Int_t pairIndex=GetPairIndex(arr1,arr2);
Int_t ntrack1=arrTracks1.GetEntriesFast();
Int_t ntrack2=arrTracks2.GetEntriesFast();
AliDielectronPair candidate;
+ candidate.SetKFUsage(fUseKF);
UInt_t selectedMask=(1<<fEventPlanePOIPreFilter.GetCuts()->GetEntries())-1;
for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1){
qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
}
+ // printf("qrsub1 %p %f \t qcsub1 %p %f \n",qrsub1,qrsub1->X(),qcsub1,qcsub1->X());
// set AliEventplane with corrected values
cevplane->SetQVector(qcorr);
cevplane->SetQsub(qcsub1, qcsub2);
Int_t ntrack1=arrTracks1.GetEntriesFast();
Int_t ntrack2=arrTracks2.GetEntriesFast();
AliDielectronPair candidate;
-
+ candidate.SetKFUsage(fUseKF);
// flag arrays for track removal
Bool_t *bTracks1 = new Bool_t[ntrack1];
for (Int_t itrack1=0; itrack1<ntrack1; ++itrack1) bTracks1[itrack1]=kFALSE;
Int_t ntrack2=arrTracks2.GetEntriesFast();
AliDielectronPair *candidate=new AliDielectronPair;
+ candidate->SetKFUsage(fUseKF);
UInt_t selectedMask=(1<<fPairFilter.GetCuts()->GetEntries())-1;
candidate->SetTracks(static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1,
static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2);
candidate->SetType(pairIndex);
- candidate->SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother));
+ Int_t label=AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother);
+ candidate->SetLabel(label);
+ if (label>-1) candidate->SetPdgCode(fPdgMother);
//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);
//apply cut
if (cutMask!=selectedMask) continue;
PairArray(pairIndex)->Add(candidate);
//get a new candidate
candidate=new AliDielectronPair;
+ candidate->SetKFUsage(fUseKF);
}
}
//delete the surplus candidate
while ( fTrackRotator->NextCombination() ){
AliDielectronPair candidate;
+ candidate.SetKFUsage(fUseKF);
candidate.SetTracks(&fTrackRotator->GetKFTrackP(), &fTrackRotator->GetKFTrackN(),
fTrackRotator->GetVTrackP(),fTrackRotator->GetVTrackN());
candidate.SetType(kEv1PMRot);
//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(10)->Add(new AliDielectronPair(candidate));
+ if(fStoreRotatedPairs) PairArray(kEv1PMRot)->Add(new AliDielectronPair(candidate));
}
}
}
}
fSignalsMC->Add(signal);
}
+//________________________________________________________________
+void AliDielectron::FillMCHistograms(const AliVEvent *ev) {
+ //
+ // fill QA MC histograms for pairs and legs of all added mc signals
+ //
+
+ if (!fSignalsMC) return;
+ TString className,className2;
+ 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++) {
+
+ className.Form("Pair_%s",fSignalsMC->At(isig)->GetName());
+ className2.Form("Track_Legs_%s",fSignalsMC->At(isig)->GetName());
+ Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
+ Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
+ if(!pairClass && !legClass) return;
+
+ 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));
+
+ 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
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff