o small update for signal extraction
authorwiechula <wiechula@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Apr 2012 22:33:25 +0000 (22:33 +0000)
committerwiechula <wiechula@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Apr 2012 22:33:25 +0000 (22:33 +0000)
o fix for coverity 18194 19858 19859

PWGDQ/dielectron/AliDielectron.cxx
PWGDQ/dielectron/AliDielectronSignalBase.h
PWGDQ/dielectron/AliDielectronSignalExt.cxx
PWGDQ/dielectron/AliDielectronVarCuts.cxx

index 3532e3c..73adbcc 100644 (file)
@@ -568,7 +568,6 @@ void AliDielectron::FillTrackArrays(AliVEvent * const ev, Int_t eventNr)
       AliESDtrack *track=static_cast<AliESDtrack*>(particle);
       track->SetESDEvent(static_cast<AliESDEvent*>(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<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());
+
+
+  // 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;
+  }
 
-  // subevents by charge separation sub1+ sub2-
-  if(fLikeSignSubEvents) {
-    qssub1 = dynamic_cast<TVector2 *>(qstd->Clone());
-    qssub2 = dynamic_cast<TVector2 *>(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; 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()) {
@@ -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<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);
   
@@ -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);
index 29b700a..edb802b 100644 (file)
@@ -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;}
index c3195a2..323cfe3 100644 (file)
@@ -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();
 
index d1ebc86..10f061f 100644 (file)
@@ -69,6 +69,7 @@ AliDielectronVarCuts::AliDielectronVarCuts(const char* name, const char* title)
     fActiveCuts[i]=0;
     fCutMin[i]=0;
     fCutMax[i]=0;
+    fCutExclude[i]=kFALSE;
   }
 }