#include <AliKFParticle.h>
+#include <AliESDInputHandler.h>
+#include <AliAnalysisManager.h>
+#include <AliEPSelectionTask.h>
#include <AliEventplane.h>
#include <AliVEvent.h>
#include <AliVParticle.h>
#include "AliDielectronDebugTree.h"
#include "AliDielectronSignalMC.h"
#include "AliDielectronMixingHandler.h"
+#include "AliDielectronV0Cuts.h"
#include "AliDielectron.h"
//________________________________________________________________
AliDielectron::AliDielectron() :
TNamed("AliDielectron","AliDielectron"),
+ fCutQA(kFALSE),
+ fQAmonitor(0x0),
fEventFilter("EventFilter"),
fTrackFilter("TrackFilter"),
fPairPreFilter("PairPreFilter"),
//________________________________________________________________
AliDielectron::AliDielectron(const char* name, const char* title) :
TNamed(name,title),
+ fCutQA(kFALSE),
+ fQAmonitor(0x0),
fEventFilter("EventFilter"),
fTrackFilter("TrackFilter"),
fPairPreFilter("PairPreFilter"),
//
// Default destructor
//
+ if (fQAmonitor) delete fQAmonitor;
if (fHistoArray) delete fHistoArray;
if (fHistos) delete fHistos;
if (fPairCandidates) delete fPairCandidates;
fHistoArray->SetSignalsMC(fSignalsMC);
fHistoArray->Init();
}
-}
+
+ if (fCutQA) {
+ fQAmonitor = new AliDielectronCutQA(Form("QAcuts_%s",GetName()),"QAcuts");
+ fQAmonitor->AddTrackFilter(&fTrackFilter);
+ fQAmonitor->AddPairFilter(&fPairFilter);
+ fQAmonitor->AddEventFilter(&fEventFilter);
+ fQAmonitor->Init();
+ }
+}
//________________________________________________________________
void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2)
return;
}
+ // modify event numbers in MC so that we can identify new events
+ // in AliDielectronV0Cuts (not neeeded for collision data)
+ if(GetHasMC()) {
+ ev1->SetBunchCrossNumber(1);
+ ev1->SetOrbitNumber(1);
+ ev1->SetPeriodNumber(1);
+ }
+
+ // set event
+ 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
+ // why do not apply the event cuts first ????
if (AliDielectronMC::Instance()->ConnectMCEvent()){
ProcessMC(ev1);
}
UInt_t selectedMask=(1<<fEventFilter.GetCuts()->GetEntries())-1;
//apply event cuts
- if ((ev1&&fEventFilter.IsSelected(ev1)!=selectedMask) ||
- (ev2&&fEventFilter.IsSelected(ev2)!=selectedMask)) return;
+ UInt_t cutmask = fEventFilter.IsSelected(ev1);
+ if(fCutQA) fQAmonitor->FillAll(ev1);
+ if(fCutQA) fQAmonitor->Fill(cutmask,ev1);
+ if ((ev1&&cutmask!=selectedMask) ||
+ (ev2&&fEventFilter.IsSelected(ev2)!=selectedMask)) return;
//fill track arrays for the first event
if (ev1){
}
// TPC event plane correction
- AliEventplane *cevplane = new AliEventplane();
- if (ev1 && cevplane && fPreFilterEventPlane && ( fEventPlanePreFilter.GetCuts()->GetEntries()>0 || fEventPlanePOIPreFilter.GetCuts()->GetEntries()>0))
- EventPlanePreFilter(0, 1, fTracks[0], fTracks[1], ev1, cevplane);
-
- // set event
- AliDielectronVarManager::SetEvent(ev1);
- if (fMixing){
- //set mixing bin to event data
- Int_t bin=fMixing->FindBin(AliDielectronVarManager::GetData());
- AliDielectronVarManager::SetValue(AliDielectronVarManager::kMixingBin,bin);
- }
+ if (ev1 && fPreFilterEventPlane && ( fEventPlanePreFilter.GetCuts()->GetEntries()>0 || fEventPlanePOIPreFilter.GetCuts()->GetEntries()>0))
+ EventPlanePreFilter(0, 1, fTracks[0], fTracks[1], ev1);
if (!fNoPairing){
// create pairs and fill pair candidate arrays
// clear arrays
if (!fDontClearArrays) ClearArrays();
+
+ // reset TPC EP and unique identifiers for v0 cut class
AliDielectronVarManager::SetTPCEventPlane(0x0);
- delete cevplane;
+ if(GetHasMC()) { // only for MC needed
+ for (Int_t iCut=0; iCut<fTrackFilter.GetCuts()->GetEntries();++iCut) {
+ if ( fTrackFilter.GetCuts()->At(iCut)->IsA() == AliDielectronV0Cuts::Class() )
+ ((AliDielectronV0Cuts*)fTrackFilter.GetCuts()->At(iCut))->ResetUniqueEventNumbers();
+ }
+ }
+
}
//________________________________________________________________
if (fHistos) FillHistogramsMC(dieMC->GetMCEvent(), ev1);
if(!fSignalsMC) return;
- //loop over all MC data and Fill the CF container if it exist
- if (!fCfManagerPair && !fHistoArray) 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;
-
// signals to be studied
Int_t nSignals = fSignalsMC->GetEntries();
+ Bool_t bFillCF = (fCfManagerPair ? fCfManagerPair->GetStepForMCtruth() : kFALSE);
+ Bool_t bFillHF = (fHistoArray ? fHistoArray->GetStepForMCGenerated() : kFALSE);
+ Bool_t bFillHist = kFALSE;
+ for(Int_t isig=0;isig<nSignals;isig++) {
+ TString sigName = fSignalsMC->At(isig)->GetName();
+ const THashList *histlist = fHistos->GetHistogramList();
+ 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;
+
+
// initialize 2D arrays of labels for particles from each MC signal
Int_t** labels1; // labels for particles satisfying branch 1
Int_t** labels2; // labels for particles satisfying branch 2
// 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
//Fill event information
if (ev){
if (fHistos->GetHistogramList()->FindObject("Event")) {
- if(fMixing) AliDielectronVarManager::SetEvent(ev); // data can be overwritten by mixed event
fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, AliDielectronVarManager::GetData());
}
}
//Fill track information, separately for the track array candidates
if (!pairInfoOnly){
+ className2.Form("Track_%s",fgkPairClassNames[1]); // unlike sign, SE only
for (Int_t i=0; i<4; ++i){
className.Form("Track_%s",fgkTrackClassNames[i]);
- if (!fHistos->GetHistogramList()->FindObject(className.Data())) continue;
+ Bool_t mergedtrkClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0;
+ Bool_t trkClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0;
+ if (!trkClass && !mergedtrkClass) continue;
Int_t ntracks=fTracks[i].GetEntriesFast();
for (Int_t itrack=0; itrack<ntracks; ++itrack){
AliDielectronVarManager::Fill(fTracks[i].UncheckedAt(itrack), values);
- fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
+ if(trkClass)
+ fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values);
+ if(mergedtrkClass && i<2)
+ fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values); //only ev1
}
}
}
AliVParticle *particle=ev->GetTrack(itrack);
//apply track cuts
- if (fTrackFilter.IsSelected(particle)!=selectedMask) continue;
+ UInt_t cutmask=fTrackFilter.IsSelected(particle);
+ //fill cut QA
+ if(fCutQA) fQAmonitor->FillAll(particle);
+ if(fCutQA) fQAmonitor->Fill(cutmask,particle);
+
+ if (cutmask!=selectedMask) continue;
//fill selected particle into the corresponding track arrays
Short_t charge=particle->Charge();
if (charge>0) fTracks[eventNr*2].Add(particle);
else if (charge<0) fTracks[eventNr*2+1].Add(particle);
+
+
}
}
//________________________________________________________________
-void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTracks1, TObjArray arrTracks2, const AliVEvent *ev, AliEventplane *cevplane)
+void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTracks1, TObjArray arrTracks2, const AliVEvent *ev)
{
//
// Prefilter tracks and tracks from pairs
// remove contribution of all tracks to the Q-vector that are in invariant mass window
//
- // TODO: how should we implement the filtering for the nanoADOs
-
AliEventplane *evplane = const_cast<AliVEvent *>(ev)->GetEventplane();
- if(!evplane) return;
-
- // do not change these vectors qref
- TVector2 * const qref = evplane->GetQVector();
- if(!qref) return;
- // random subevents
- TVector2 *qrsub1 = evplane->GetQsub1();
- TVector2 *qrsub2 = evplane->GetQsub2();
-
- // copy references
- TVector2 *qcorr = new TVector2(*qref);
- TVector2 *qcsub1 = 0x0;
- TVector2 *qcsub2 = 0x0;
- // printf("qrsub1 %p %f \n",qrsub1,qrsub1->X());
-
-
- // 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;
- }
+ if(!evplane) { // nanoAODs , here we do NOT have sub event reaction planes
+ // if(1) {
+ // get the EPselectionTask for recalculation of weighting factors
+ AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+ AliEPSelectionTask *eptask = dynamic_cast<AliEPSelectionTask *>(man->GetTask("EventplaneSelection"));
+ if(!eptask) return;
+
+ // 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;
+ }
- // subevent separation
- if(fLikeSignSubEvents || etagap) {
- qcsub1 = new TVector2(*qcorr);
- qcsub2 = new TVector2(*qcorr);
+ Double_t cQX=0., cQY=0.;
+ // apply cuts to the tracks, e.g. etagap
+ if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
+ UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
+ Int_t ntracks=ev->GetNumberOfTracks();
+ for (Int_t itrack=0; itrack<ntracks; ++itrack){
+ AliVParticle *particle=ev->GetTrack(itrack);
+ AliVTrack *track= static_cast<AliVTrack*>(particle);
+ if (!track) continue;
+ //event plane cuts
+ UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
+ //apply cut
+ if (cutMask==selectedMask) continue;
+
+ mapRemovedTracks.Add(track,track);
+ cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
+ cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
+ }
+ }
- 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;
+ // POI (particle of interest) rejection
+ Int_t pairIndex=GetPairIndex(arr1,arr2);
- 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) qcsub1->Set(qcsub1->X()-cQX, qcsub1->Y()-cQY);
- if (eta>0.0) qcsub2->Set(qcsub2->X()-cQX, qcsub2->Y()-cQY);
+ 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){
+ Int_t end=ntrack2;
+ if (arr1==arr2) end=itrack1;
+ Bool_t accepted=kFALSE;
+ for (Int_t itrack2=0; itrack2<end; ++itrack2){
+ TObject *track1=arrTracks1.UncheckedAt(itrack1);
+ TObject *track2=arrTracks2.UncheckedAt(itrack2);
+ if (!track1 || !track2) continue;
+ //create the pair
+ candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
+ static_cast<AliVTrack*>(track2), fPdgLeg2);
+ candidate.SetType(pairIndex);
+ candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
+
+ //event plane pair cuts
+ UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
+ //apply cut
+ if (cutMask==selectedMask) continue;
+
+ accepted=kTRUE;
+ //remove the tracks from the Track arrays
+ arrTracks2.AddAt(0x0,itrack2);
}
+ if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
}
- }
- else {
- // by a random
- qcsub1 = new TVector2(*qrsub1);
- qcsub2 = new TVector2(*qrsub2);
- }
-
- // apply cuts, e.g. etagap
- if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
- UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
- Int_t ntracks=ev->GetNumberOfTracks();
- for (Int_t itrack=0; itrack<ntracks; ++itrack){
- AliVParticle *particle=ev->GetTrack(itrack);
- AliVTrack *track= static_cast<AliVTrack*>(particle);
+ //compress the track arrays
+ arrTracks1.Compress();
+ arrTracks2.Compress();
+
+ //Modify the components: subtract the tracks
+ ntrack1=arrTracks1.GetEntriesFast();
+ ntrack2=arrTracks2.GetEntriesFast();
+ // remove leg1 contribution
+ for (Int_t itrack=0; itrack<ntrack1; ++itrack){
+ AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
if (!track) continue;
-
- //event plane cuts
- UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
- //apply cut
- if (cutMask==selectedMask) continue;
+ // track contribution was already removed
+ if (mapRemovedTracks.FindObject(track)) continue;
+ else mapRemovedTracks.Add(track,track);
- Double_t cQX = 0.0;
- Double_t cQY = 0.0;
- if(!etagap) {
- cQX = evplane->GetQContributionX(track);
- cQY = evplane->GetQContributionY(track);
- }
- Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
- Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
- Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
- Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
-
- // update Q vectors
- qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
- qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
- qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
+ cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
+ cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
}
- }
+ // remove leg2 contribution
+ for (Int_t itrack=0; itrack<ntrack2; ++itrack){
+ AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
+ if (!track) continue;
+ // track contribution was already removed
+ if (mapRemovedTracks.FindObject(track)) continue;
+ else mapRemovedTracks.Add(track,track);
- // 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){
- Int_t end=ntrack2;
- if (arr1==arr2) end=itrack1;
- Bool_t accepted=kFALSE;
- for (Int_t itrack2=0; itrack2<end; ++itrack2){
- TObject *track1=arrTracks1.UncheckedAt(itrack1);
- TObject *track2=arrTracks2.UncheckedAt(itrack2);
- if (!track1 || !track2) continue;
- //create the pair
- candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
- static_cast<AliVTrack*>(track2), fPdgLeg2);
-
- candidate.SetType(pairIndex);
- candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
-
- //event plane cuts
- UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
- //apply cut
- if (cutMask==selectedMask) continue;
+ cQX += (eptask->GetWeight(track) * TMath::Cos(2*track->Phi()));
+ cQY += (eptask->GetWeight(track) * TMath::Sin(2*track->Phi()));
+ }
- accepted=kTRUE;
- //remove the tracks from the Track arrays
- arrTracks2.AddAt(0x0,itrack2);
+ // build a corrected alieventplane using the values from the var manager
+ // these uncorrected values are filled using the stored magnitude and angle in the header
+ TVector2 qcorr;
+ qcorr.Set(AliDielectronVarManager::GetValue(AliDielectronVarManager::kTPCxH2uc)-cQX,
+ AliDielectronVarManager::GetValue(AliDielectronVarManager::kTPCyH2uc)-cQY);
+ // fill alieventplane
+ AliEventplane cevplane;
+ cevplane.SetQVector(&qcorr);
+ AliDielectronVarManager::SetTPCEventPlane(&cevplane);
+ cevplane.SetQVector(0);
+ return;
+ } //end: nanoAODs
+ else
+ {
+ // this is done in case of ESDs or AODs
+ Bool_t isESD=(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class());
+ // copy event plane object
+ AliEventplane cevplane(*evplane);
+ // Int_t nMaxID=cevplane->GetQContributionXArray()->GetSize();
+
+ TVector2 *qcorr = cevplane.GetQVector();
+ if(!qcorr) return;
+ TVector2 *qcsub1 = 0x0;
+ TVector2 *qcsub2 = 0x0;
+
+ // 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 configuration for eta gap or LS (default is rndm)
+ if(fLikeSignSubEvents && etagap) {
+ // start with the full Qvector/event in both sub events
+ qcsub1 = new TVector2(*qcorr);
+ qcsub2 = new TVector2(*qcorr);
+ cevplane.SetQsub(qcsub1,qcsub2);
+
+ 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;
+ if (track->GetID()>=0 && !isESD) continue;
+ Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
+
+ // set contributions to zero
+ // charge sub1+ sub2-
+ if(fLikeSignSubEvents) {
+ Short_t charge=track->Charge();
+ if (charge<0) {
+ cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
+ }
+ if (charge>0) {
+ cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
+ }
+ }
+ // eta sub1+ sub2-
+ if(etagap) {
+ Double_t eta=track->Eta();
+ if (eta<0.0) {
+ cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
+ }
+ if (eta>0.0) {
+ cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
+ }
+ }
+ } // end: loop over tracks
+ } // end: sub event configuration
+
+ // apply cuts, e.g. etagap
+ if(fEventPlanePreFilter.GetCuts()->GetEntries()) {
+ UInt_t selectedMask=(1<<fEventPlanePreFilter.GetCuts()->GetEntries())-1;
+ Int_t ntracks=ev->GetNumberOfTracks();
+ for (Int_t itrack=0; itrack<ntracks; ++itrack){
+ AliVParticle *particle=ev->GetTrack(itrack);
+ AliVTrack *track= static_cast<AliVTrack*>(particle);
+ if (!track) continue;
+ if (track->GetID()>=0 && !isESD) continue;
+ Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
+
+ //event plane cuts
+ UInt_t cutMask=fEventPlanePreFilter.IsSelected(track);
+ //apply cut
+ if (cutMask==selectedMask) continue;
+
+ // set contributions to zero
+ cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
+ }
+ } // end: track cuts
+
+ // 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){
+ Int_t end=ntrack2;
+ if (arr1==arr2) end=itrack1;
+ Bool_t accepted=kFALSE;
+ for (Int_t itrack2=0; itrack2<end; ++itrack2){
+ TObject *track1=arrTracks1.UncheckedAt(itrack1);
+ TObject *track2=arrTracks2.UncheckedAt(itrack2);
+ if (!track1 || !track2) continue;
+ //create the pair
+ candidate.SetTracks(static_cast<AliVTrack*>(track1), fPdgLeg1,
+ static_cast<AliVTrack*>(track2), fPdgLeg2);
+
+ candidate.SetType(pairIndex);
+ candidate.SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(&candidate,fPdgMother));
+
+ //event plane cuts
+ UInt_t cutMask=fEventPlanePOIPreFilter.IsSelected(&candidate);
+ //apply cut
+ if (cutMask==selectedMask) continue;
+
+ accepted=kTRUE;
+ //remove the tracks from the Track arrays
+ arrTracks2.AddAt(0x0,itrack2);
+ }
if ( accepted ) arrTracks1.AddAt(0x0,itrack1);
- }
- //compress the track arrays
- arrTracks1.Compress();
- arrTracks2.Compress();
-
-
- //Modify the components: subtract the tracks
- ntrack1=arrTracks1.GetEntriesFast();
- ntrack2=arrTracks2.GetEntriesFast();
-
- // remove leg1 contribution
- for (Int_t itrack=0; itrack<ntrack1; ++itrack){
- AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
- if (!track) continue;
-
- Double_t cQX = evplane->GetQContributionX(track);
- Double_t cQY = evplane->GetQContributionY(track);
- Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
- Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
- Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
- Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
-
- // update Q vectors
- qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
- qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
- qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
- }
- // remove leg2 contribution
- for (Int_t itrack=0; itrack<ntrack2; ++itrack){
- AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
- if (!track) continue;
-
- Double_t cQX = evplane->GetQContributionX(track);
- Double_t cQY = evplane->GetQContributionY(track);
- Double_t cQXsub1 = evplane->GetQContributionXsub1(track);
- Double_t cQYsub1 = evplane->GetQContributionYsub1(track);
- Double_t cQXsub2 = evplane->GetQContributionXsub2(track);
- Double_t cQYsub2 = evplane->GetQContributionYsub2(track);
-
- // update Q vectors
- qcorr->Set(qcorr->X()-cQX, qcorr->Y()-cQY);
- qcsub1->Set(qcsub1->X()-cQXsub1, qcsub1->Y()-cQYsub1);
- qcsub2->Set(qcsub2->X()-cQXsub2, qcsub2->Y()-cQYsub2);
- }
+ }
+ //compress the track arrays
+ arrTracks1.Compress();
+ arrTracks2.Compress();
+
+ //Modify the components: subtract the tracks
+ ntrack1=arrTracks1.GetEntriesFast();
+ ntrack2=arrTracks2.GetEntriesFast();
+ // remove leg1 contribution
+ for (Int_t itrack=0; itrack<ntrack1; ++itrack){
+ AliVTrack *track= static_cast<AliVTrack*>(arrTracks1.UncheckedAt(itrack));
+ if (!track) continue;
+ if (track->GetID()>=0 && !isESD) continue;
+ Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
+ // set contributions to zero
+ cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
+ }
+ // remove leg2 contribution
+ for (Int_t itrack=0; itrack<ntrack2; ++itrack){
+ AliVTrack *track= static_cast<AliVTrack*>(arrTracks2.UncheckedAt(itrack));
+ if (!track) continue;
+ if (track->GetID()>=0 && !isESD) continue;
+ Int_t tmpID = isESD ? track->GetID() : track->GetID()*-1 - 1;
+ // set contributions to zero
+ cevplane.GetQContributionXArray()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionYArray()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionXArraysub1()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionYArraysub1()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionXArraysub2()->SetAt(0.0, tmpID);
+ cevplane.GetQContributionYArraysub2()->SetAt(0.0, tmpID);
+ }
- // 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);
- AliDielectronVarManager::SetTPCEventPlane(cevplane);
-}
+ // set corrected AliEventplane and fill variables with corrected values
+ AliDielectronVarManager::SetTPCEventPlane(&cevplane);
+ delete qcsub1;
+ delete qcsub2;
+ } // end: ESD or AOD case
+}
//________________________________________________________________
void AliDielectron::PairPreFilter(Int_t arr1, Int_t arr2, TObjArray &arrTracks1, TObjArray &arrTracks2)
{
if (fCfManagerPair) fCfManagerPair->Fill(cutMask,candidate);
//histogram array for the pair
if (fHistoArray) fHistoArray->Fill(pairIndex,candidate);
+ // cut qa
+ if(pairIndex==kEv1PM && fCutQA) {
+ fQAmonitor->FillAll(candidate);
+ fQAmonitor->Fill(cutMask,candidate);
+ }
//apply cut
if (cutMask!=selectedMask) continue;
}
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
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff