]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGPP/MUON/lite/AliAnalysisTaskMuonCuts.cxx
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGPP / MUON / lite / AliAnalysisTaskMuonCuts.cxx
index e154109d3f6bf59d49ba5342aac3f1eba753f672..a85663022472c3b5f0bddad10b3b9e9073187f77 100644 (file)
@@ -43,6 +43,7 @@
 #include "TObjArray.h"
 #include "TF1.h"
 #include "TStyle.h"
+#include "TArrayI.h"
 //#include "TMCProcess.h"
 
 // STEER includes
@@ -60,7 +61,9 @@
 
 // PWG includes
 #include "AliMergeableCollection.h"
+#include "AliMuonEventCuts.h"
 #include "AliMuonTrackCuts.h"
+#include "AliAnalysisMuonUtility.h"
 
 
 /// \cond CLASSIMP
@@ -141,7 +144,7 @@ void AliAnalysisTaskMuonCuts::MyUserCreateOutputObjects()
       AddObjectToCollection(histo);
 
       histoName = GetHistoName(kDCAVsPCheck, itheta, isrc);
-      histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 800., 100, 0., 200.);
+      histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 800., 400, 0., 200.);
       histo->SetXTitle("p (GeV/c)");
       histo->SetYTitle("DCA (cm)");
       AddObjectToCollection(histo);
@@ -192,18 +195,15 @@ void AliAnalysisTaskMuonCuts::ProcessEvent(TString physSel, const TObjArray& sel
   /// Fill histogram
   //
 
-  if ( GetVertexSPD()->GetNContributors() < fMinNvtxContirbutors ) return;
-
   TString histoName = "";
   AliVParticle* track = 0x0;
-  Int_t nTracks = GetNTracks();
+  Int_t nTracks = AliAnalysisMuonUtility::GetNTracks(InputEvent());
   for (Int_t itrack = 0; itrack < nTracks; itrack++) {
-    track = GetTrack(itrack);
-    fMuonTrackCuts->SetNSigmaPdca(1.e10);
+    track = AliAnalysisMuonUtility::GetTrack(itrack, InputEvent());
+    fMuonTrackCuts->CustomParam()->SetNSigmaPdca(1.e10);
     if ( ! fMuonTrackCuts->IsSelected(track) ) continue;
 
-    Double_t rAbsEnd =  ( fAODEvent ) ? ((AliAODTrack*)track)->GetRAtAbsorberEnd(): ((AliESDMuonTrack*)track)->GetRAtAbsorberEnd();
-    Double_t thetaAbsEndDeg = TMath::ATan( rAbsEnd / 505. ) * TMath::RadToDeg();
+    Double_t thetaAbsEndDeg = AliAnalysisMuonUtility::GetThetaAbsDeg(track);
     Int_t thetaAbsBin = ( thetaAbsEndDeg < 3. ) ? kThetaAbs23 : kThetaAbs310;
 
     Int_t trackSrc = GetParticleType(track);
@@ -213,7 +213,8 @@ void AliAnalysisTaskMuonCuts::ProcessEvent(TString physSel, const TObjArray& sel
     Double_t dca = dcaAtVz.Mag();
     Double_t pDca = pTotMean * dca;
 
-    Double_t chi2 = pDca / fMuonTrackCuts->GetSigmaPdca(rAbsEnd) ;
+    Double_t sigmaPdca = fMuonTrackCuts->IsThetaAbs23(track) ? fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca23() : fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca310();
+    Double_t chi2 = pDca / sigmaPdca;
     chi2 *= chi2;
     Double_t chiProb = TMath::Prob(chi2, 2);
 
@@ -243,7 +244,7 @@ void AliAnalysisTaskMuonCuts::ProcessEvent(TString physSel, const TObjArray& sel
     } // loop on selected trigger classes
 
     for ( Int_t isigma=0; isigma<fSigmaCuts.GetSize(); ++isigma) {
-      fMuonTrackCuts->SetNSigmaPdca(fSigmaCuts[isigma]);
+      fMuonTrackCuts->CustomParam()->SetNSigmaPdca(fSigmaCuts[isigma]);
       if ( ! fMuonTrackCuts->IsSelected(track) ) continue;
       for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
         TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
@@ -306,8 +307,10 @@ void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
     }
   }
   delete optArray;
+  
+  fMuonTrackCuts->Print("param");
 
-  Int_t srcColors[kNtrackSources] = {kBlack, kRed, kGreen, kBlue, kViolet, 7, kOrange};
+  Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange, kGray};
 
   TCanvas* can = 0x0;
   Int_t xshift = 100;
@@ -335,7 +338,7 @@ void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
     for ( Int_t ihisto=0; ihisto<2; ++ihisto ) {
       TH2* histo = 0x0;
       histoPattern = "";
-      histoPattern = Form("%s & %s", fHistoTypeKeys->At(recoDcaHisto[ihisto])->GetName(), fThetaAbsKeys->At(itheta)->GetName());
+      histoPattern = Form("%s%s*", fHistoTypeKeys->At(recoDcaHisto[ihisto])->GetName(), fThetaAbsKeys->At(itheta)->GetName());
       histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
       if ( ! histo ) continue;
 
@@ -421,10 +424,10 @@ void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
   //////////////
   // Fit pDCA //
   //////////////
-  Double_t nSigmaCut = fMuonTrackCuts->GetNSigmaPdca(); //6.;
-  Double_t sigmaMeasCut[2] = { fMuonTrackCuts->GetSigmaPdca(AliMuonTrackCuts::kThetaAbs23), fMuonTrackCuts->GetSigmaPdca(AliMuonTrackCuts::kThetaAbs23)}; //{99., 54.}; //{120., 63.};
-  Double_t relPResolution = fMuonTrackCuts->GetRelPResolution(); //4.5e-4; //6.e-3;//8.e-4;
-  Double_t angleResolution = 535.*fMuonTrackCuts->GetSlopeResolution(); //6.e-4;
+  Double_t nSigmaCut = fMuonTrackCuts->GetMuonTrackCutsParam().GetNSigmaPdca(); //6.;
+  Double_t sigmaMeasCut[2] = { fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca23(), fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca310()}; //{99., 54.}; //{120., 63.};
+  Double_t relPResolution = fMuonTrackCuts->GetMuonTrackCutsParam().GetRelPResolution(); //4.5e-4; //6.e-3;//8.e-4;
+  Double_t angleResolution = 535.*fMuonTrackCuts->GetMuonTrackCutsParam().GetSlopeResolution(); //6.e-4;
   Double_t pMinCut = 0.1;
   Double_t pMaxCut =  800.;
   const Int_t kNCutFuncs = 2;
@@ -441,7 +444,7 @@ void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
   fitFunc->SetParNames("Norm", "Mean", "Sigma");
   gStyle->SetOptFit(1111);
   Double_t xMinFit[2] = {0., 0.};
-  Double_t xMaxFit[2] = {320., 150.}; // {360., 180.};
+  Double_t xMaxFit[2] = {260., 150.}; //{320., 150.}; // {360., 180.};
   printf("\nSigma p x DCA:\n");
   Double_t averageSigmaPdca[kNtrackSources*kNthetaAbs] = {0.};
   for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
@@ -473,7 +476,7 @@ void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
       Int_t ican = 0;
 
       for ( Int_t ibin=2; ibin<=nPbins; ++ibin ) {
-        currName = Form("hPDCA_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data());
+        currName = Form("hPDCA_%s_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data(), fThetaAbsKeys->At(itheta)->GetName());
         Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
         Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
         if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
@@ -521,28 +524,6 @@ void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
         printf("  %s %s: Sigma = %g +- %g  (chi2/%i = %g)\n",  fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName(), trendFit->GetParameter(0), trendFit->GetParError(0), trendFit->GetNDF(), trendFit->GetChisquare()/((Double_t)trendFit->GetNDF()));
       }
       leg->AddEntry(sigmaVsP, fSrcKeys->At(isrc)->GetName(), "lp");
-
-      // Plot 2D function for check!
-      histoPattern = GetHistoName(kPDCAVsPCheck, itheta, isrc);
-      TH2* histoCheck = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
-      if ( ! histoCheck ) histoCheck = histo; // needed for old data
-      currName = histoCheck->GetName();
-      currName.Append("_plotCut");
-      TCanvas* pdcaVsPcan = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,(igroup2+1)*yshift,600,600);
-      pdcaVsPcan->SetLogz();
-      pdcaVsPcan->SetRightMargin(0.12);
-      histoCheck->Draw("COLZ");
-
-      for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
-        currName = Form("%s_cutFunc%i", histo->GetName(), icut);
-        TF1* cutFunction = new TF1(currName.Data(),cutFormula.Data(), pMinCut, pMaxCut);
-        cutParam[icut][0] = sigmaMeasCut[itheta];
-        cutParam[icut][1] = nSigmaCut;
-        cutFunction->SetParameters(cutParam[icut]);
-        cutFunction->SetLineWidth(2);
-        cutFunction->SetLineColor(cutColor[icut]);
-        cutFunction->Draw("same");
-      } // loop on cut func
     } // loop on src
     can->cd(itheta+1);
     leg->Draw("same");
@@ -565,33 +546,37 @@ void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
   printf("\n");
 
   igroup2++;
-  for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
-    for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
-      histoPattern = GetHistoName(kDCAVsPCheck, itheta, isrc);
-      TH2* histoCheck = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
-      if ( ! histoCheck ) continue;
-      currName = histoCheck->GetName();
-      currName.Append("_plotCut");
-      TCanvas* pdcaVsPcan = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,(igroup2+1)*yshift,600,600);
-      pdcaVsPcan->SetRightMargin(0.12);
-      pdcaVsPcan->SetLogz();
-      histoCheck->Draw("COLZ");
-
-      for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
+  for ( Int_t icheck=0; icheck<2; icheck++ ) {
+    Int_t hDraw = ( icheck == 0 ) ? kPDCAVsPCheck : kDCAVsPCheck;
+    for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
+      for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
+        histoPattern = GetHistoName(hDraw, itheta, isrc);
+        TH2* histoCheck = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
+        if ( ! histoCheck ) continue;
         currName = histoCheck->GetName();
-        currName.Append(Form("_cutFunc%i",icut));
-        TString currFormula = cutFormula;
-        currFormula.Append("/x");
-        TF1* cutFunction = new TF1(currName.Data(),currFormula.Data(), pMinCut, pMaxCut);
-        cutParam[icut][0] = sigmaMeasCut[itheta];
-        cutParam[icut][1] = nSigmaCut;
-        cutFunction->SetParameters(cutParam[icut]);
-        cutFunction->SetLineWidth(2);
-        cutFunction->SetLineColor(cutColor[icut]);
-        cutFunction->Draw("same");
-      } // loop on cut functions
-    } // loop on src
-  } //loop on theta
+        currName.Append("_plotCut");
+        TCanvas* pdcaVsPcan = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,(igroup2+icheck)*yshift,600,600);
+        pdcaVsPcan->SetRightMargin(0.12);
+        pdcaVsPcan->SetLogy();
+        pdcaVsPcan->SetLogz();
+        histoCheck->Draw("COLZ");
+
+        for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
+          currName = histoCheck->GetName();
+          currName.Append(Form("_cutFunc%i",icut));
+          TString currFormula = cutFormula;
+          if ( icheck == 1 ) currFormula.Append("/x");
+          TF1* cutFunction = new TF1(currName.Data(),currFormula.Data(), pMinCut, pMaxCut);
+          cutParam[icut][0] = sigmaMeasCut[itheta];
+          cutParam[icut][1] = nSigmaCut;
+          cutFunction->SetParameters(cutParam[icut]);
+          cutFunction->SetLineWidth(2);
+          cutFunction->SetLineColor(cutColor[icut]);
+          cutFunction->Draw("same");
+        } // loop on cut functions
+      } // loop on src
+    } //loop on theta
+  } // loop on check
 
 
   ///////////////////////////
@@ -616,7 +601,7 @@ void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
       Int_t ican = 0;
 
       for ( Int_t ibin=0; ibin<=nPbins; ++ibin ) {
-        currName = Form("hChiProb_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data());
+        currName = Form("hChiProb_%s_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data(), fThetaAbsKeys->At(itheta)->GetName());
         Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
         Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
         if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
@@ -646,8 +631,7 @@ void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
   Float_t fracOfHeight = 0.35;
   Float_t rightMargin = 0.03;
   Int_t cutColors[14] = {kBlack, kRed, kBlue, kGreen, kCyan, kMagenta, kOrange, kViolet, kSpring, kGray, kSpring, kAzure, kPink, kYellow};
-  Int_t* orderCuts = 0x0;
-  Int_t nSigmaCuts = 0;
+  TArrayI orderCuts;
   Int_t checkHistos[2] = {kSigmaVsPt, kSigmaVsEta};
   Bool_t useCustomSigma = furtherOpt.Contains("CUSTOMSIGMA");
 
@@ -665,11 +649,11 @@ void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
           delete histo;
           continue;
         }
-        if ( ! orderCuts ) {
+        if ( orderCuts.GetSize() == 0 ) {
           // Re-order axis
           TAxis* axis = histo->GetYaxis();
-          nSigmaCuts = ( useCustomSigma ) ? fSigmaCuts.GetSize() : axis->GetNbins();
-          orderCuts = new Int_t[nSigmaCuts];
+          Int_t nSigmaCuts = ( useCustomSigma ) ? fSigmaCuts.GetSize() : axis->GetNbins();
+          orderCuts.Set(nSigmaCuts);
           Int_t countBin = 0;
           for ( Int_t isigma=0; isigma<axis->GetNbins(); ++isigma ) {
             TString currLabel = axis->GetBinLabel(isigma+1);
@@ -708,7 +692,7 @@ void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
 
         TH1* refCutHisto = 0x0;
         TString legTitle = "";
-        for ( Int_t isigma=0; isigma<nSigmaCuts; ++isigma ) {
+        for ( Int_t isigma=0; isigma<orderCuts.GetSize(); ++isigma ) {
           currName = histo->GetName();
           currName.Append(Form("_sigma%i", isigma));
           Int_t currBin = orderCuts[isigma];
@@ -733,7 +717,24 @@ void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
           currName.Append("_ratio");
           TH1* histoRatio = (TH1*)projectHisto->Clone(currName.Data());
           histoRatio->Sumw2();
-          histoRatio->Divide(refCutHisto);
+          TH1* auxNum = projectHisto;
+          TH1* auxDen = refCutHisto;
+          Bool_t mustInvert = kFALSE;
+          if ( auxNum->Integral()>auxDen->Integral() ) {
+            auxNum = refCutHisto;
+            auxDen = projectHisto;
+            mustInvert = kTRUE;
+          }
+          histoRatio->Divide(auxNum,auxDen,1.,1.,"B");
+          if ( mustInvert ) {
+            TH1* auxHistoOne = static_cast<TH1*>(histoRatio->Clone("auxHistoOne"));
+            for ( Int_t ibin=1; ibin<=auxHistoOne->GetXaxis()->GetNbins(); ibin++ ) {
+              auxHistoOne->SetBinContent(ibin,1.);
+              auxHistoOne->SetBinError(ibin,0.);
+            }
+            histoRatio->Divide(auxHistoOne,histoRatio);
+            delete auxHistoOne;
+          }
           histoRatio->SetTitle("");
           histoRatio->GetYaxis()->SetTitle(legTitle.Data());
           histoRatio->GetXaxis()->SetLabelSize(0.04/fracOfHeight);
@@ -757,7 +758,7 @@ void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
 
   igroup1++;
   igroup2=0;
-  Double_t ptMin[] = {0., 2., 4., 15.};
+  Double_t ptMin[] = {0., 2., 4., 15., 60.};
   Int_t nPtMins = sizeof(ptMin)/sizeof(ptMin[0]);
   for ( Int_t iptmin=0; iptmin<nPtMins; ++iptmin) {
     for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
@@ -765,13 +766,13 @@ void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
       for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
         TLegend* leg = 0x0;
         histoName = GetHistoName(kSigmaVsPt, itheta, isrc);
-        for ( Int_t icent=1; icent<=fCentralityClasses->GetNbins(); ++icent ) {
-          TH2* histo = (TH2*)GetSum(physSel, trigClassName, fCentralityClasses->GetBinLabel(icent), histoName);
+        for ( Int_t icent=1; icent<=GetCentralityClasses()->GetNbins(); ++icent ) {
+          TH2* histo = (TH2*)GetSum(physSel, trigClassName, GetCentralityClasses()->GetBinLabel(icent), histoName);
           if ( ! histo ) continue;
           Int_t ptMinBin = histo->GetXaxis()->FindBin(ptMin[iptmin]);
           Int_t ptMaxBin = histo->GetXaxis()->GetNbins()+1;
           currName = histo->GetName();
-          currName += Form("_%s_ptMin%i", fCentralityClasses->GetBinLabel(icent), TMath::Nint(ptMin[iptmin]));
+          currName += Form("_%s_ptMin%i", GetCentralityClasses()->GetBinLabel(icent), TMath::Nint(ptMin[iptmin]));
           TH1* projectHisto = histo->ProjectionY(currName.Data(), ptMinBin, ptMaxBin);
           if ( projectHisto->GetMaximum() < 50. ) {
             delete projectHisto;
@@ -795,17 +796,16 @@ void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
           drawOpt = "e";
           if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
           projectHisto->Draw(drawOpt.Data());
-          leg->AddEntry(projectHisto, fCentralityClasses->GetBinLabel(icent), "lp");
+          leg->AddEntry(projectHisto, GetCentralityClasses()->GetBinLabel(icent), "lp");
           Double_t keptEvents = projectHisto->GetBinContent(orderCuts[0]);
-          Double_t totalEvents = projectHisto->GetBinContent(orderCuts[nSigmaCuts-1]);
+          Double_t totalEvents = projectHisto->GetBinContent(orderCuts[orderCuts.GetSize()-1]);
           Double_t accepted = ( totalEvents > 0. ) ? keptEvents / totalEvents : 1.;
           Double_t acceptedErr = ( totalEvents > 0. ) ? TMath::Sqrt(accepted*(1.-accepted)/totalEvents) : 1.;
-          printf("%12s %11s %6s (pt>%g) rejected evts : %6.2f +- %6.3f %%\n", fSrcKeys->At(isrc)->GetName(), fThetaAbsKeys->At(itheta)->GetName(), fCentralityClasses->GetBinLabel(icent), ptMin[iptmin], (1.-accepted)*100., acceptedErr*100.);
-          //printf("  rejected %g  total %g   (%s vs %s)\n",totalEvents-keptEvents,totalEvents,projectHisto->GetXaxis()->GetBinLabel(orderCuts[0]),projectHisto->GetXaxis()->GetBinLabel(orderCuts[nSigmaCuts-1]));
+          printf("%12s %11s %6s (pt>%g) rejected evts : %6.2f +- %6.3f %%\n", fSrcKeys->At(isrc)->GetName(), fThetaAbsKeys->At(itheta)->GetName(), GetCentralityClasses()->GetBinLabel(icent), ptMin[iptmin], (1.-accepted)*100., acceptedErr*100.);
+          //printf("  rejected %g  total %g   (%s vs %s)\n",totalEvents-keptEvents,totalEvents,projectHisto->GetXaxis()->GetBinLabel(orderCuts[0]),projectHisto->GetXaxis()->GetBinLabel(orderCuts[orderCuts.GetSize()-1]));
         } // loop on centrality
         if ( leg ) leg->Draw("same");
       } // loop on theta abs
     } // loop on sources
   } // loop on pt min
-  delete [] orderCuts;
 }