From: wiechula Date: Thu, 12 Apr 2012 22:33:25 +0000 (+0000) Subject: o small update for signal extraction X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=b9d223bb5930e2127d31899ba6f78abff074556d;p=u%2Fmrichter%2FAliRoot.git o small update for signal extraction o fix for coverity 18194 19858 19859 --- diff --git a/PWGDQ/dielectron/AliDielectron.cxx b/PWGDQ/dielectron/AliDielectron.cxx index 3532e3c5640..73adbcc52b8 100644 --- a/PWGDQ/dielectron/AliDielectron.cxx +++ b/PWGDQ/dielectron/AliDielectron.cxx @@ -568,7 +568,6 @@ void AliDielectron::FillTrackArrays(AliVEvent * const ev, Int_t eventNr) AliESDtrack *track=static_cast(particle); track->SetESDEvent(static_cast(ev)); //only in trunk... } - //apply track cuts if (fTrackFilter.IsSelected(particle)!=selectedMask) continue; @@ -588,24 +587,34 @@ void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTra // 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 *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(qref->Clone()); - TVector2 *qssub1 = dynamic_cast(qrsub1->Clone()); - TVector2 *qssub2 = dynamic_cast(qrsub2->Clone()); + 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; + } - // subevents by charge separation sub1+ sub2- - if(fLikeSignSubEvents) { - qssub1 = dynamic_cast(qstd->Clone()); - qssub2 = dynamic_cast(qstd->Clone()); + // subevent separation + if(fLikeSignSubEvents || etagap) { + qcsub1 = new TVector2(*qcorr); + qcsub2 = new TVector2(*qcorr); Int_t ntracks=ev->GetNumberOfTracks(); @@ -618,39 +627,25 @@ void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTra 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; iCutGetEntries();++iCut) { - TString cutName=fEventPlanePreFilter.GetCuts()->At(iCut)->GetName(); - if(cutName.Contains("eta") || cutName.Contains("Eta")) { - etagap=kTRUE; - - qssub1 = dynamic_cast(qstd->Clone()); - qssub2 = dynamic_cast(qstd->Clone()); - - Int_t ntracks=ev->GetNumberOfTracks(); - - // track removals - for (Int_t itrack=0; itrackGetTrack(itrack); - AliVTrack *track= static_cast(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()) { @@ -678,19 +673,12 @@ void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTra 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(qstd->Clone()); - TVector2 *qcsub1 = dynamic_cast(qssub1->Clone()); - TVector2 *qcsub2 = dynamic_cast(qssub2->Clone()); - - // POI (particle of interest) rejection Int_t pairIndex=GetPairIndex(arr1,arr2); @@ -769,6 +757,7 @@ void AliDielectron::EventPlanePreFilter(Int_t arr1, Int_t arr2, TObjArray arrTra 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); diff --git a/PWGDQ/dielectron/AliDielectronSignalBase.h b/PWGDQ/dielectron/AliDielectronSignalBase.h index 29b700a07c8..edb802b5d94 100644 --- a/PWGDQ/dielectron/AliDielectronSignalBase.h +++ b/PWGDQ/dielectron/AliDielectronSignalBase.h @@ -36,6 +36,8 @@ public: kFitted, kLikeSign, kLikeSignArithm, + kLikeSignRcorr, + kLikeSignArithmRcorr, kEventMixing, kRotation }; @@ -49,6 +51,7 @@ public: void SetFitRange(Double_t min, Double_t max) {fFitMin=min; fFitMax=max;} void SetRebin(Int_t factor) {fRebin = factor;} void SetMethod(EBackgroundMethod method) {fMethod = method;} + Int_t GetMethod() const { return (Int_t)fMethod; } const TVectorD& GetValues() const {return fValues;} const TVectorD& GetErrors() const {return fErrors;} diff --git a/PWGDQ/dielectron/AliDielectronSignalExt.cxx b/PWGDQ/dielectron/AliDielectronSignalExt.cxx index c3195a2aa86..323cfe37121 100644 --- a/PWGDQ/dielectron/AliDielectronSignalExt.cxx +++ b/PWGDQ/dielectron/AliDielectronSignalExt.cxx @@ -78,6 +78,8 @@ void AliDielectronSignalExt::Process(TObjArray* const arrhist) switch ( fMethod ){ case kLikeSign : case kLikeSignArithm : + case kLikeSignRcorr: + case kLikeSignArithmRcorr: ProcessLS(arrhist); // process like-sign subtraction method break; @@ -139,7 +141,7 @@ void AliDielectronSignalExt::ProcessLS(TObjArray* const arrhist) Float_t background = 2*TMath::Sqrt(pp*mm); Float_t ebackground = TMath::Sqrt(mm+pp); - if (fMethod==kLikeSignArithm){ + if (fMethod==kLikeSignArithm || fMethod==kLikeSignArithmRcorr ){ //Arithmetic mean instead of geometric background=(pp+mm); ebackground=TMath::Sqrt(pp+mm); @@ -151,13 +153,23 @@ void AliDielectronSignalExt::ProcessLS(TObjArray* const arrhist) } //correct LS spectrum bin-by-bin with R factor obtained in mixed events - if(fMixingCorr) { + if(fMixingCorr || fMethod==kLikeSignRcorr || fMethod==kLikeSignArithmRcorr) { TH1* histMixPP = (TH1*)(arrhist->At(AliDielectron::kEv1PEv2P))->Clone("mixPP"); // ++ ME TH1* histMixMM = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2M))->Clone("mixMM"); // -- ME - TH1* histMixPM = (TH1*)(arrhist->At(AliDielectron::kEv1PEv2M))->Clone("mixPM"); // +- ME + + TH1* histMixPM = 0x0; + if (arrhist->At(AliDielectron::kEv1MEv2P)){ + histMixPM = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2P))->Clone("mixPM"); // -+ ME + } + + if (arrhist->At(AliDielectron::kEv1PEv2M)){ + TH1 *htmp=(TH1*)(arrhist->At(AliDielectron::kEv1PEv2M)); + if (!histMixPM) fHistDataME = (TH1*)htmp->Clone("mixPM"); // +- ME + else histMixPM->Add(htmp); + } + histMixPM->Sumw2(); - histMixPM->Add((TH1*)arrhist->At(AliDielectron::kEv1MEv2P)); // -+ ME // rebin the histograms if (fRebin>1) { @@ -236,59 +248,59 @@ void AliDielectronSignalExt::ProcessEM(TObjArray* const arrhist) // // event mixing of +- and -+ // - fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPMSE"); // +- SE - fHistDataME = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2P))->Clone("histMPME"); // -+ ME - fHistDataPM->Sumw2(); - fHistDataME->Sumw2(); - fHistDataPM->SetDirectory(0); - fHistDataME->SetDirectory(0); - - fHistDataME->Add((TH1*)arrhist->At(AliDielectron::kEv1PEv2M)); // +- ME - // rebin the histograms - if (fRebin>1) { - fHistDataPM->Rebin(fRebin); - fHistDataME->Rebin(fRebin); + if (!arrhist->At(AliDielectron::kEv1PM) || !(arrhist->At(AliDielectron::kEv1MEv2P) || arrhist->At(AliDielectron::kEv1PEv2M)) ){ + AliError("Either OS or mixed histogram missing"); + return; } - fHistSignal = new TH1D("HistSignal", "Mixed events background substracted signal", - fHistDataPM->GetXaxis()->GetNbins(), - fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax()); - fHistSignal->SetDirectory(0); - fHistBackground = new TH1D("HistBackground", "background contribution from mixed events", - fHistDataPM->GetXaxis()->GetNbins(), - fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax()); - fHistBackground->SetDirectory(0); + delete fHistDataPM; fHistDataPM=0x0; + delete fHistDataME; fHistDataME=0x0; + delete fHistBackground; fHistBackground=0x0; + + fHistDataPM = (TH1*)(arrhist->At(AliDielectron::kEv1PM))->Clone("histPM"); // +- SE + fHistDataPM->Sumw2(); + fHistDataPM->SetDirectory(0x0); - // fill out background and subtracted histogram - for(Int_t ibin=1; ibin<=fHistDataPM->GetXaxis()->GetNbins(); ibin++) { - Float_t pm = fHistDataPM->GetBinContent(ibin); - Float_t epm = fHistDataPM->GetBinError(ibin); - Float_t background = fHistDataME->GetBinContent(ibin); - Float_t ebackground = fHistDataME->GetBinError(ibin); + if (arrhist->At(AliDielectron::kEv1MEv2P)){ + fHistDataME = (TH1*)(arrhist->At(AliDielectron::kEv1MEv2P))->Clone("histMPME"); // -+ ME + } + + if (arrhist->At(AliDielectron::kEv1PEv2M)){ + TH1 *htmp=(TH1*)(arrhist->At(AliDielectron::kEv1PEv2M)); + if (!fHistDataME) fHistDataME = (TH1*)htmp->Clone("histMPME"); // -+ ME + else fHistDataME->Add(htmp); + } - fHistSignal->SetBinContent(ibin, pm); - fHistSignal->SetBinError(ibin, epm); - fHistBackground->SetBinContent(ibin, background); - fHistBackground->SetBinError(ibin, ebackground); + fHistBackground = (TH1*)fHistDataME->Clone("ME_Background"); + fHistBackground->SetDirectory(0x0); + fHistBackground->Sumw2(); + + // rebin the histograms + if (fRebin>1) { + fHistDataPM->Rebin(fRebin); + fHistDataME->Rebin(fRebin); + fHistBackground->Rebin(fRebin); } + //scale histograms to match integral between fScaleMin and fScaleMax + // or if fScaleMax < fScaleMin use fScaleMin as scale factor if (fScaleMax>fScaleMin) fScaleFactor=ScaleHistograms(fHistDataPM,fHistBackground,fScaleMin,fScaleMax); - else if (fScaleMin>0.){ + else if (fScaleMin>0.){ fScaleFactor=fScaleMin; fHistBackground->Scale(fScaleFactor); } - //subract background - fHistSignal->Add(fHistBackground,-1); + fHistSignal=(TH1*)fHistDataPM->Clone("Signal"); + fHistSignal->Add(fHistBackground,-1.); - // signal + // signal fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin), - fHistSignal->FindBin(fIntMax), fErrors(0)); + fHistSignal->FindBin(fIntMax), fErrors(0)); // background fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin), - fHistBackground->FindBin(fIntMax), - fErrors(1)); + fHistBackground->FindBin(fIntMax), + fErrors(1)); // S/B and significance SetSignificanceAndSOB(); diff --git a/PWGDQ/dielectron/AliDielectronVarCuts.cxx b/PWGDQ/dielectron/AliDielectronVarCuts.cxx index d1ebc865067..10f061f0c82 100644 --- a/PWGDQ/dielectron/AliDielectronVarCuts.cxx +++ b/PWGDQ/dielectron/AliDielectronVarCuts.cxx @@ -69,6 +69,7 @@ AliDielectronVarCuts::AliDielectronVarCuts(const char* name, const char* title) fActiveCuts[i]=0; fCutMin[i]=0; fCutMax[i]=0; + fCutExclude[i]=kFALSE; } }