X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGDQ%2Fdielectron%2FAliDielectron.cxx;h=65287b1173de3c00ecebdaa044d363db7af7d817;hb=8f8a1f7485bb3f8336ff87007742e83a909c9a5b;hp=37c9d8d097009cdd16c1ee9395fb8bd451cacff1;hpb=7ab4dc00cc4869b0e7e7bcf146c85bf4ea509bd6;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGDQ/dielectron/AliDielectron.cxx b/PWGDQ/dielectron/AliDielectron.cxx index 37c9d8d0970..65287b1173d 100644 --- a/PWGDQ/dielectron/AliDielectron.cxx +++ b/PWGDQ/dielectron/AliDielectron.cxx @@ -46,9 +46,13 @@ The names are available via the function PairClassName(Int_t i) #include #include #include +#include #include +#include +#include +#include #include #include #include @@ -62,6 +66,8 @@ The names are available via the function PairClassName(Int_t i) #include "AliDielectronDebugTree.h" #include "AliDielectronSignalMC.h" #include "AliDielectronMixingHandler.h" +#include "AliDielectronV0Cuts.h" +#include "AliDielectronPID.h" #include "AliDielectron.h" @@ -91,6 +97,10 @@ const char* AliDielectron::fgkPairClassNames[11] = { //________________________________________________________________ AliDielectron::AliDielectron() : TNamed("AliDielectron","AliDielectron"), + fCutQA(kFALSE), + fQAmonitor(0x0), + fPostPIDCntrdCorr(0x0), + fPostPIDWdthCorr(0x0), fEventFilter("EventFilter"), fTrackFilter("TrackFilter"), fPairPreFilter("PairPreFilter"), @@ -103,6 +113,7 @@ AliDielectron::AliDielectron() : fPdgLeg2(11), fSignalsMC(0x0), fNoPairing(kFALSE), + fProcessLS(kTRUE), fUseKF(kTRUE), fHistoArray(0x0), fHistos(0x0), @@ -121,7 +132,10 @@ AliDielectron::AliDielectron() : fEstimatorFilename(""), fTRDpidCorrectionFilename(""), fVZEROCalibrationFilename(""), - fVZERORecenteringFilename("") + fVZERORecenteringFilename(""), + fEffMapFilename(""), + fZDCRecenteringFilename("") + { // // Default constructor @@ -132,6 +146,10 @@ AliDielectron::AliDielectron() : //________________________________________________________________ AliDielectron::AliDielectron(const char* name, const char* title) : TNamed(name,title), + fCutQA(kFALSE), + fQAmonitor(0x0), + fPostPIDCntrdCorr(0x0), + fPostPIDWdthCorr(0x0), fEventFilter("EventFilter"), fTrackFilter("TrackFilter"), fPairPreFilter("PairPreFilter"), @@ -144,6 +162,7 @@ AliDielectron::AliDielectron(const char* name, const char* title) : fPdgLeg2(11), fSignalsMC(0x0), fNoPairing(kFALSE), + fProcessLS(kTRUE), fUseKF(kTRUE), fHistoArray(0x0), fHistos(0x0), @@ -162,7 +181,9 @@ AliDielectron::AliDielectron(const char* name, const char* title) : fEstimatorFilename(""), fTRDpidCorrectionFilename(""), fVZEROCalibrationFilename(""), - fVZERORecenteringFilename("") + fVZERORecenteringFilename(""), + fEffMapFilename(""), + fZDCRecenteringFilename("") { // // Named constructor @@ -176,13 +197,16 @@ AliDielectron::~AliDielectron() // // Default destructor // - if (fHistoArray) delete fHistoArray; + if (fQAmonitor) delete fQAmonitor; + 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; } //________________________________________________________________ @@ -205,17 +229,39 @@ void AliDielectron::Init() 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); + fQAmonitor->AddPairFilter(&fPairFilter); + fQAmonitor->AddEventFilter(&fEventFilter); + fQAmonitor->Init(); + } +} //________________________________________________________________ void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2) @@ -229,6 +275,16 @@ void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2) AliError("At least first event must be set!"); 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 @@ -237,10 +293,11 @@ void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2) } //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); } - + //if candidate array doesn't exist, create it if (!fPairCandidates->UncheckedAt(0)) { InitPairCandidateArrays(); @@ -252,10 +309,11 @@ void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2) UInt_t selectedMask=(1<GetEntries())-1; //apply event cuts - if ((ev1&&fEventFilter.IsSelected(ev1)!=selectedMask) || - (ev2&&fEventFilter.IsSelected(ev2)!=selectedMask)) return; - -// AliDielectronVarManager::SetEvent(ev1); // why a second time??? + 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){ @@ -271,15 +329,15 @@ void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2) } // 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); + 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 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); } } @@ -296,9 +354,15 @@ void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2) //process event mixing if (fMixing) { fMixing->Fill(ev1,this); -// FillHistograms(0x0,kTRUE); + // 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); @@ -306,8 +370,16 @@ void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2) // 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; iCutGetEntries();++iCut) { + if ( fTrackFilter.GetCuts()->At(iCut)->IsA() == AliDielectronV0Cuts::Class() ) + ((AliDielectronV0Cuts*)fTrackFilter.GetCuts()->At(iCut))->ResetUniqueEventNumbers(); + } + } + } //________________________________________________________________ @@ -321,17 +393,33 @@ void AliDielectron::ProcessMC(AliVEvent *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;isigAt(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 @@ -390,11 +478,14 @@ void AliDielectron::ProcessMC(AliVEvent *ev1) // Do the pairing and fill the CF container with pure MC info for(Int_t isig=0; isigFillMC(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 @@ -402,6 +493,7 @@ void AliDielectron::ProcessMC(AliVEvent *ev1) for(Int_t i2=0; i2FillMC(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 @@ -469,29 +561,29 @@ void AliDielectron::FillHistograms(const AliVEvent *ev, Bool_t pairInfoOnly) TString className,className2; Double_t values[AliDielectronVarManager::kNMaxValues]={0.}; + //Fill event information - 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 (ev){ + if (fHistos->GetHistogramList()->FindObject("Event")) { + fHistos->FillClass("Event", AliDielectronVarManager::kNMaxValues, AliDielectronVarManager::GetData()); + } } - if (fHistos->GetHistogramList()->FindObject("Event")) -// 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 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; itrackFillClass(className, AliDielectronVarManager::kNMaxValues, values); + if(trkClass) + fHistos->FillClass(className, AliDielectronVarManager::kNMaxValues, values); + if(mergedtrkClass && i<2) + fHistos->FillClass(className2, AliDielectronVarManager::kNMaxValues, values); //only ev1 } } } @@ -588,209 +680,312 @@ void AliDielectron::FillTrackArrays(AliVEvent * const ev, Int_t eventNr) Int_t ntracks=ev->GetNumberOfTracks(); - // fTrackFilter.Init(); UInt_t selectedMask=(1<GetEntries())-1; for (Int_t itrack=0; itrackGetTrack(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 // Needed for rejection in the Q-Vector of the event plane // remove contribution of all tracks to the Q-vector that are in invariant mass window // + AliEventplane *evplane = const_cast(ev)->GetEventplane(); -// AliEventplane *evplane = 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; iCutGetEntries();++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(man->GetTask("EventplaneSelection")); + if(!eptask) return; + + // track mapping + TMap mapRemovedTracks; + + + Double_t cQX=0., cQY=0.; + // apply cuts to the tracks, e.g. etagap + if(fEventPlanePreFilter.GetCuts()->GetEntries()) { + UInt_t selectedMask=(1<GetEntries())-1; + Int_t ntracks=ev->GetNumberOfTracks(); + for (Int_t itrack=0; itrackGetTrack(itrack); + AliVTrack *track= static_cast(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())); + } + } - // subevent separation - if(fLikeSignSubEvents || etagap) { - qcsub1 = new TVector2(*qcorr); - qcsub2 = new TVector2(*qcorr); + // POI (particle of interest) rejection + Int_t pairIndex=GetPairIndex(arr1,arr2); - Int_t ntracks=ev->GetNumberOfTracks(); - - // track removals - for (Int_t itrack=0; itrackGetTrack(itrack); - AliVTrack *track= static_cast(particle); - if (!track) continue; + Int_t ntrack1=arrTracks1.GetEntriesFast(); + Int_t ntrack2=arrTracks2.GetEntriesFast(); + AliDielectronPair candidate; + candidate.SetKFUsage(fUseKF); - 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); + UInt_t selectedMask=(1<GetEntries())-1; + for (Int_t itrack1=0; itrack1(track1), fPdgLeg1, + static_cast(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<GetEntries())-1; - Int_t ntracks=ev->GetNumberOfTracks(); - for (Int_t itrack=0; itrackGetTrack(itrack); - AliVTrack *track= static_cast(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(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(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<GetEntries())-1; - for (Int_t itrack1=0; itrack1(track1), fPdgLeg1, - static_cast(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; iCutGetEntries();++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; itrackGetTrack(itrack); + AliVTrack *track= static_cast(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<GetEntries())-1; + Int_t ntracks=ev->GetNumberOfTracks(); + for (Int_t itrack=0; itrackGetTrack(itrack); + AliVTrack *track= static_cast(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<GetEntries())-1; + for (Int_t itrack1=0; itrack1(track1), fPdgLeg1, + static_cast(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(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(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(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(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) { @@ -897,7 +1092,7 @@ void AliDielectron::PairPreFilter(Int_t arr1, Int_t arr2, TObjArray &arrTracks1, 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(); @@ -953,11 +1148,12 @@ void AliDielectron::FillPairArrays(Int_t arr1, Int_t arr2) for (Int_t itrack2=0; itrack2SetTracks(static_cast(arrTracks1.UncheckedAt(itrack1)), fPdgLeg1, - static_cast(arrTracks2.UncheckedAt(itrack2)), fPdgLeg2); + static_cast(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); @@ -969,15 +1165,22 @@ void AliDielectron::FillPairArrays(Int_t arr1, Int_t arr2) //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==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 @@ -1009,14 +1212,16 @@ void AliDielectron::FillPairArrayTR() //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)); + } } } @@ -1057,6 +1262,66 @@ void AliDielectron::AddSignalMC(AliDielectronSignalMC* signal) { } 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) { // @@ -1064,38 +1329,78 @@ 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; isigGetEntries(); 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(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(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; itrackGetLabel(); + 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; +}