1) Adding class AliAnalysisMuonUtility which contains static methods allowing to...
authorpcrochet <pcrochet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 21 Sep 2012 14:10:15 +0000 (14:10 +0000)
committerpcrochet <pcrochet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 21 Sep 2012 14:10:15 +0000 (14:10 +0000)
15 files changed:
PWG/CMakelibPWGmuon.pkg
PWG/PWGmuonLinkDef.h
PWG/muon/AddTaskFlowSingleMu.C
PWG/muon/AliAnalysisTaskFlowSingleMu.cxx
PWG/muon/AliAnalysisTaskFlowSingleMu.h
PWG/muon/AliAnalysisTaskSingleMu.cxx
PWG/muon/AliMuonPairCuts.cxx
PWG/muon/AliMuonPairCuts.h
PWG/muon/AliMuonTrackCuts.cxx
PWG/muon/AliMuonTrackCuts.h
PWG/muon/AliVAnalysisMuon.cxx
PWG/muon/AliVAnalysisMuon.h
PWGPP/MUON/lite/AliAnalysisTaskMuonCuts.cxx
PWGPP/MUON/lite/AliAnalysisTaskTrigChEff.cxx
PWGPP/macros/AddTaskMTRchamberEfficiency.C

index 73e2c79..44baae1 100644 (file)
@@ -57,11 +57,13 @@ set ( SRCS
     muon/AliAODMuonReplicator.cxx 
     muon/AliAnalysisTaskMuonCollisionMultiplicity.cxx 
     muon/AliCFMuonResUpsilon.cxx 
+    muon/AliMuonEventCuts.cxx 
     muon/AliMuonTrackCuts.cxx 
     muon/AliMuonPairCuts.cxx 
     muon/AliMergeableCollection.cxx 
     muon/AliVAnalysisMuon.cxx 
-    muon/AliAnalysisTaskFlowSingleMu.cxx
+    muon/AliAnalysisTaskFlowSingleMu.cxx 
+    muon/AliAnalysisMuonUtility.cxx
     )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
index 66e18c1..f5195b1 100644 (file)
 #pragma link C++ class AliAnalysisMuMu+;
 #pragma link C++ class AliAnalysisMuMuFromAOD+;
 #pragma link C++ class AliAnalysisMuMuFromESD+;
+#pragma link C++ class AliMuonEventCuts+;
 #pragma link C++ class AliMuonTrackCuts+;
 #pragma link C++ class AliMuonPairCuts+;
 #pragma link C++ class AliMergeableCollection+;
 #pragma link C++ class AliMergeableCollectionIterator+;
 #pragma link C++ class AliVAnalysisMuon+;
 #pragma link C++ class AliAnalysisTaskFlowSingleMu+;
+#pragma link C++ class AliAnalysisMuonUtility+;
 #endif
 
index ab95a5e..88a061a 100644 (file)
@@ -41,6 +41,7 @@ AliAnalysisTaskFlowSingleMu* AddTaskFlowSingleMu(Bool_t isMC = kFALSE)
   // Create task
   AliAnalysisTaskFlowSingleMu *flowSingleMuTask = new AliAnalysisTaskFlowSingleMu("FlowSingleMuTask", *muonTrackCuts);
   if ( isMC ) flowSingleMuTask->SetTrigClassPatterns("ANY");
+  flowSingleMuTask->GetMuonEventCuts()->SetVertexVzLimits(-10., 10.);
   mgr->AddTask(flowSingleMuTask);
 
   // Connect containers
index 0a325a8..1cdc932 100644 (file)
@@ -76,7 +76,9 @@
 #include "AliVAnalysisMuon.h"
 #include "AliMergeableCollection.h"
 #include "AliCounterCollection.h"
+#include "AliMuonEventCuts.h"
 #include "AliMuonTrackCuts.h"
+#include "AliAnalysisMuonUtility.h"
 
 
 /// \cond CLASSIMP
@@ -226,10 +228,8 @@ void AliAnalysisTaskFlowSingleMu::ProcessEvent(TString physSel, const TObjArray&
   //
   // Base event cuts
   //
-  if ( GetVertexSPD()->GetNContributors() < fMinNvtxContirbutors ) return;
   //Bool_t isPileupFromSPD = ( fAODEvent && ! fAODEvent->GetTracklets() ) ? InputEvent()->IsPileupFromSPD(3, 0.8, 3., 2., 5.) : InputEvent()->IsPileupFromSPDInMultBins(); // Avoid break when reading Muon AODs (tracklet info is not present and IsPileupFromSPDInMultBins crashes
   //if ( isPileupFromSPD ) return;
-  if ( TMath::Abs(GetVertexSPD()->GetZ()) > 10. ) return;
 
   TArrayD psiPlane(fEPKeys->GetEntriesFast());
   psiPlane.Reset(-999.);
@@ -275,6 +275,7 @@ void AliAnalysisTaskFlowSingleMu::ProcessEvent(TString physSel, const TObjArray&
         if ( qsub1 && qsub2 )
           ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hResoSub2%s",currEP.Data())))->Fill(TMath::Cos(qsub1->Phi()/2.-qsub2->Phi()/2.));
       }
+      //else ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, Form("hResoSub2%s",currEP.Data())))->Fill(TMath::Cos(2.*psiPlane[iep])); // REMEMBER TO CUT
     } // loop on trigger
   } // loop on EP
 
@@ -284,9 +285,9 @@ void AliAnalysisTaskFlowSingleMu::ProcessEvent(TString physSel, const TObjArray&
   Double_t containerInput[kNvars];
   AliVParticle* track = 0x0;
   for ( Int_t istep = 0; istep<2; ++istep ) {
-    Int_t nTracks = ( istep == kStepReconstructed ) ? GetNTracks() : GetNMCTracks();
+    Int_t nTracks = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetNTracks(InputEvent()) : AliAnalysisMuonUtility::GetNMCTracks(InputEvent(),MCEvent());
     for (Int_t itrack = 0; itrack < nTracks; itrack++) {
-      track = ( istep == kStepReconstructed ) ? GetTrack(itrack) : GetMCTrack(itrack);
+      track = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetTrack(itrack,InputEvent()) : AliAnalysisMuonUtility::GetMCTrack(itrack,InputEvent(),MCEvent());
       
       Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 );
       if ( ! isSelected ) continue;
@@ -316,7 +317,7 @@ void AliAnalysisTaskFlowSingleMu::ProcessEvent(TString physSel, const TObjArray&
 
         for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
           TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
-          if ( istep == kStepReconstructed && ! TrackPtCutMatchTrigClass(track, trigClassName) ) continue;
+          if ( istep == kStepReconstructed && ! fMuonTrackCuts->TrackPtCutMatchTrigClass(track, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) ) continue;
           ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, Form("FlowSingleMuContainer%s",currEP.Data())))->Fill(containerInput,istep);
         } // loop on selected trigger classes
       } // loop on event plane
@@ -396,8 +397,9 @@ void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
   TAxis* customBins = 0x0;
   
   //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 5., 6., 8., 10.};
-  Double_t myBins[] = {4., 5., 6., 8., 10.}; // REMEMBER TO CHOOSE
-  //Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10.};
+  //Double_t myBins[] = {4., 5., 6., 8., 10.}; // REMEMBER TO CHOOSE
+  //Double_t myBins[] = {6., 8., 10.};
+  Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10.};
   //Double_t myBins[] = {0., 0.5, 1., 1.5, 2., 3., 4., 10.};
   //Double_t myBins[] = {2., 3., 4., 5., 6., 8., 10., 80.};
   Int_t nMyBins = sizeof(myBins)/sizeof(myBins[0])-1;
@@ -408,6 +410,8 @@ void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
   Int_t yshift = 20;
   Int_t igroup1 = 0;
   Int_t igroup2 = 0;
+  
+  TList outList;
 
   //
   /// Read plane resolution file (if present)
@@ -458,15 +462,28 @@ void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
       resoCentrality.Add(new TObjString(centralityClasses->At(icent)->GetName()));
       Int_t icos = 0;
       Double_t sumRelErrSquare = 0.;
+      TString hResoTitle = "";
       for ( Int_t isub=0; isub<resoEParray.GetEntries(); isub++ ) {
-        for ( Int_t jsub=0; jsub<resoEParray.GetEntries(); jsub++ ) {
-          TH1* histo = static_cast<TH1*>(GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),Form("hResoSub3%s%s", resoEParray.At(isub)->GetName(),resoEParray.At(jsub)->GetName())));
+        TString resoDet1 = resoEParray.At(isub)->GetName();
+        for ( Int_t jsub=isub+1; jsub<resoEParray.GetEntries(); jsub++ ) {
+          TString resoDet2 = resoEParray.At(jsub)->GetName();
+          TH1* histo = 0x0;
+          // Try both combinations: we want the first two values to contain the selectedEP!
+          TString resoDet = "";
+          for ( Int_t icomb=0; icomb<2; icomb++ ) {
+            TString hName = "hResoSub3";
+            resoDet = ( icomb == 0 ) ? resoDet1 + resoDet2 : resoDet2 + resoDet1;
+            hName += resoDet;
+            histo = static_cast<TH1*>(GetSum(physSel,trigClassName,resoCentrality.At(icent)->GetName(),hName));
+            if ( histo ) break;
+          }
           if ( ! histo || histo->Integral() == 0 ) continue;
           currCos[icos] = histo->GetMean();
           if ( currCos[icos] <= 0. ) continue;
           Double_t relErr = histo->GetMeanError() / currCos[icos];
           sumRelErrSquare += relErr * relErr;
           icos++;
+          hResoTitle += Form("%s ", resoDet.Data());
         } // loop on sub events
       } // loop on sub events
       if ( icos != 3  ) {
@@ -476,7 +493,7 @@ void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
       }
       resolution[0][icent] = TMath::Sqrt((currCos[0]*currCos[1]/currCos[2]));
       resolution[1][icent] = resolution[0][icent]/2.*TMath::Sqrt(sumRelErrSquare);
-      printf("Resolution %s  %g +- %g\n", centralityClasses->At(icent)->GetName(), resolution[0][icent], resolution[1][icent]);
+      printf("Resolution (%s) %s  %g +- %g\n", hResoTitle.Data(), centralityClasses->At(icent)->GetName(), resolution[0][icent], resolution[1][icent]);
     }
   }
 
@@ -490,22 +507,7 @@ void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
     if ( icent == 0 ) resoCentrArray.AddAt(range[0],ib++);
     resoCentrArray.AddAt(range[1], ib++);
   }
-
-  if ( hasResolution ) {
-    currName = externalReso ? Form("externalReso") : Form("reso_%s_%s_%s", resoEParray.At(0)->GetName(), resoEParray.At(1)->GetName(), resoEParray.At(2)->GetName());
-    can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
-    TGraphErrors* resoGraph = new TGraphErrors(centralityClasses->GetEntries());
-    for ( Int_t icent=0; icent<resoCentrality.GetEntries(); icent++ ) {
-
-      resoGraph->SetPoint(icent, 0.5*(resoCentrArray[icent+1] + resoCentrArray[icent]), resolution[0][icent]);
-      resoGraph->SetPointError(icent, 0.5*(resoCentrArray[icent+1] - resoCentrArray[icent]), resolution[1][icent]);
-    }
-    resoGraph->SetTitle(Form("Resolution %s", selectEP.Data()));
-    resoGraph->GetXaxis()->SetTitle("Centrality");
-    resoGraph->GetYaxis()->SetTitle("Resolution");
-    resoGraph->Draw("apz");
-  }
-
+  
   const Int_t kNresoCentr = resoCentrality.GetEntries();
 
   //
@@ -523,87 +525,166 @@ void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
     for ( Int_t istep=0; istep<kNsteps; istep++ ) {
       stepName[istep] = cfContainer->GetStepTitle(istep);
       gridSparseArray[istep*kNresoCentr+icent] = cfContainer->GetGrid(istep);
+      if ( ! customBins ) customBins = gridSparseArray[istep*kNresoCentr+icent]->GetAxis(kHvarPt);
     } // loop on step
   } // loop on centrality
-
-  Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange};
+  
   
   Bool_t isMC = furtherOpt.Contains("MC");
   Int_t firstSrc = ( isMC ) ? 0 : kUnidentified;
   Int_t lastSrc  = ( isMC ) ? kNtrackSources - 1 : kUnidentified;
+  Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange};
   if ( ! isMC ) srcColors[kUnidentified] = 1;
-
-  TList outList;
-
-  TString graphTypeName[3] = {"raw","correctStat","correctSyst"};
-  Int_t drawOrder[3] = {0, 2, 1};
-  TString drawOpt[3] = {"apz","pz","a2"};
-  Int_t nTypes = ( hasResolution ) ? 3 : 1;
-
-  TString histoName = "", histoPattern = "";
-  ///////////
-  // Flow  //
-  ///////////
-  gStyle->SetOptFit(1111);
-  TString funcFormula = ( v2only ) ? "[0] * (1. + 2.*[1]*TMath::Cos(2*x))" : "[0] * (1. + 2.*([1]*TMath::Cos(x) + [2]*TMath::Cos(2*x) + [3]*TMath::Cos(3*x)))";
+  
+  
   //
-  TF1* func = new TF1("funcFlow",funcFormula.Data(), 0., TMath::Pi());
-  if  ( v2only ) func->SetParNames("scale", "v2");
-  else func->SetParNames("scale","v1", "v2", "v3");
-  Int_t v2par = ( v2only ) ? 1 : 2;
-
+  // Get histograms
+  //
+  const Int_t kNcentralities = centralityClasses->GetEntries();
+  const Int_t kNhistos = kNsteps*kNtrackSources*kNcentralities;
+  TObjArray ptVsDeltaPhiList(kNhistos);
+  TObjArray ptVsPhiList(kNhistos);
+  TArrayD wgtReso[3];
+  for ( Int_t ires=0; ires<3; ires++ ) {
+    wgtReso[ires].Set(kNhistos);
+    wgtReso[ires].Reset(0.);
+  }
   for ( Int_t istep=0; istep<kNsteps; istep++ ) {
     for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
       TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
-      TGraphErrors* flowVsCentrality[3];
-      for ( Int_t itype=0; itype<nTypes; itype++ ) {
-        histoName = Form("v2VsCentrality_%s_%s", baseName.Data(), graphTypeName[itype].Data());
-        flowVsCentrality[itype] = new TGraphErrors();
-        flowVsCentrality[itype]->SetName(histoName.Data());
-        flowVsCentrality[itype]->SetTitle(histoName.Data());
-      }
       for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
-        TH2* hDeltaPhi2D = 0x0;
+        Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
+        TH1* hPtVsDeltaPhi = 0x0;
+        TH1* hPtVsPhi = 0x0;
         TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
-        TArrayD weightedReso(3);
-        weightedReso.Reset(0.);
         Double_t nMuons = 0.;
         TString debugString = Form("\n%s", centralityClasses->At(icent)->GetName());
         for ( Int_t rcent=0; rcent<kNresoCentr; rcent++ ) {
           if ( resoCentrArray[rcent] < centralityArray[icent] ||
-               resoCentrArray[rcent+1] > centralityArray[icent+1] ) continue;
+              resoCentrArray[rcent+1] > centralityArray[icent+1] ) continue;
           AliCFGridSparse* gridSparse = gridSparseArray[istep*kNresoCentr+rcent];
           if ( ! gridSparse ||  gridSparse->GetEntries() == 0. ) continue;
-
+          
           SetSparseRange(gridSparse, kHvarEta, "", -3.999, -2.501);
           SetSparseRange(gridSparse, kHvarMotherType, "", isrc+1, isrc+1, "USEBIN");
           TH1* auxHisto = gridSparse->Project(kHvarDeltaPhi,kHvarPt);
-
+          
           Double_t currYield = auxHisto->Integral();
           if ( currYield == 0. ) continue;
           nMuons += currYield;
           debugString += Form("\nAdding %s  yield %g  reso", resoCentrality.At(rcent)->GetName(), currYield);
           for ( Int_t ires=0; ires<3; ires++ ) {
-            weightedReso[ires] += resolution[ires][rcent] * currYield;
+            wgtReso[ires][idx] += resolution[ires][rcent] * currYield;
             debugString += Form("  %g", resolution[ires][rcent]);
           }
-          histoName = Form("ptVsDeltaPhi_%s", baseNameCent.Data());
-          if ( ! hDeltaPhi2D ) hDeltaPhi2D = static_cast<TH2*>(auxHisto->Clone(histoName.Data()));
-          else hDeltaPhi2D->Add(auxHisto);
+          if ( ! hPtVsDeltaPhi ) hPtVsDeltaPhi = static_cast<TH1*>(auxHisto->Clone(Form("ptVsDeltaPhi_%s", baseNameCent.Data())));
+          else hPtVsDeltaPhi->Add(auxHisto);
           delete auxHisto;
-        }
-
-        if ( ! hDeltaPhi2D ) continue;
-
+          
+          auxHisto = gridSparse->Project(kHvarPhi,kHvarPt);
+          if ( ! hPtVsPhi ) hPtVsPhi = static_cast<TH1*>(auxHisto->Clone(Form("ptVsPhi_%s", baseNameCent.Data())));
+          else hPtVsPhi->Add(auxHisto);
+          delete auxHisto;
+          
+        } // loop on resolution centralities
+        
+        if ( nMuons == 0. ) continue;
         debugString += Form("\nWeighted reso ");
         for ( Int_t ires=0; ires<3; ires++ ) {
-          weightedReso[ires] = weightedReso[ires]/nMuons;
-          debugString += Form("  %g", weightedReso[ires]);
+          wgtReso[ires][idx] = wgtReso[ires][idx]/nMuons;
+          debugString += Form("  %g", wgtReso[ires][idx]);
         }
         AliDebug(1,debugString.Data());
+        
+        ptVsDeltaPhiList.AddAt(hPtVsDeltaPhi,idx);
+        ptVsPhiList.AddAt(hPtVsPhi,idx);
+      } // loop on centralities
+    } // loop on sources
+  } // loop on steps
+  
+  
+  
+  /////////////////////
+  // Plot resolution //
+  /////////////////////
+  if ( hasResolution ) {
+    currName = Form("reso_%s", selectEP.Data());
+    currName += externalReso ? Form("_external") : Form("_sub_%s_%s_%s", resoEParray.At(0)->GetName(), resoEParray.At(1)->GetName(), resoEParray.At(2)->GetName());
+    TGraphErrors* resoGraph = new TGraphErrors(kNresoCentr);
+    resoGraph->SetName(currName.Data());
+    currName.Append("_can");
+    can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
+    
+    for ( Int_t icent=0; icent<kNresoCentr; icent++ ) {
+      resoGraph->SetPoint(icent, 0.5*(resoCentrArray[icent+1] + resoCentrArray[icent]), resolution[0][icent]);
+      resoGraph->SetPointError(icent, 0.5*(resoCentrArray[icent+1] - resoCentrArray[icent]), resolution[1][icent]);
+    }
+    resoGraph->SetTitle(Form("Resolution %s", selectEP.Data()));
+    resoGraph->GetXaxis()->SetTitle("Centrality");
+    resoGraph->GetYaxis()->SetTitle("Resolution");
+    resoGraph->Draw("apz");
+    outList.Add(resoGraph);
+    
+    for ( Int_t istep=0; istep<kNsteps; istep++ ) {
+      for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
+        currName = Form("wgtReso_%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
+        resoGraph = new TGraphErrors();
+        resoGraph->SetName(currName.Data());
+        for ( Int_t icent=0; icent<kNcentralities; icent++ ) {
+          Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
+          if ( wgtReso[0][idx] == 0. ) continue;
+          resoGraph->SetPoint(icent, 0.5*(centralityArray[icent+1] + centralityArray[icent]), wgtReso[0][idx]);
+          resoGraph->SetPointError(icent, 0.5*(centralityArray[icent+1] - centralityArray[icent]), wgtReso[1][idx]);
+        } // loop on centralities
+        if ( resoGraph->GetN() == 0 ) continue;
+        printf("New reso %s\n", resoGraph->GetName());
+        resoGraph->SetTitle(Form("Resolution %s", selectEP.Data()));
+        resoGraph->GetXaxis()->SetTitle("Centrality");
+        resoGraph->GetYaxis()->SetTitle("Resolution");
+        resoGraph->Draw("pz");
+        outList.Add(resoGraph);
+      } // loop on sources
+    } // loop on steps
+  }
+  
+  
+  
+  TString graphTypeName[3] = {"raw","correctStat","correctSyst"};
+  Int_t drawOrder[3] = {0, 2, 1};
+  TString drawOpt[3] = {"apz","pz","a2"};
+  Int_t nTypes = ( hasResolution ) ? 3 : 1;
 
-        TAxis* ptAxis = hDeltaPhi2D->GetYaxis();
+  TString histoName = "", histoPattern = "";
+  ///////////
+  // Flow  //
+  ///////////
+  
+  gStyle->SetOptFit(1111);
+  TString funcFormula = ( v2only ) ? "[0] * (1. + 2.*[1]*TMath::Cos(2*x))" : "[0] * (1. + 2.*([1]*TMath::Cos(x) + [2]*TMath::Cos(2*x) + [3]*TMath::Cos(3*x)))";
+  //
+  TF1* func = new TF1("funcFlow",funcFormula.Data(), 0., TMath::Pi());
+  if  ( v2only ) func->SetParNames("scale", "v2");
+  else func->SetParNames("scale","v1", "v2", "v3");
+  Int_t v2par = ( v2only ) ? 1 : 2;
+  
+  for ( Int_t istep=0; istep<kNsteps; istep++ ) {
+    for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
+      TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
+      TGraphErrors* flowVsCentrality[3];
+      for ( Int_t itype=0; itype<nTypes; itype++ ) {
+        histoName = Form("v2VsCentrality_%s_%s", baseName.Data(), graphTypeName[itype].Data());
+        flowVsCentrality[itype] = new TGraphErrors();
+        flowVsCentrality[itype]->SetName(histoName.Data());
+        flowVsCentrality[itype]->SetTitle(histoName.Data());
+      }
+      for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
+        Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
+        TH2* histo2D = static_cast<TH2*> (ptVsDeltaPhiList.At(idx));
+        if ( ! histo2D ) continue;
+        TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
 
+        TAxis* ptAxis = histo2D->GetYaxis();
+        
         Int_t ipad = 0;
         TGraphErrors* flowVsPt[3];
         for ( Int_t itype=0; itype<nTypes; itype++ ) {
@@ -612,10 +693,7 @@ void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
           flowVsPt[itype]->SetName(histoName.Data());
           flowVsPt[itype]->SetTitle(histoName.Data());
         }
-
-        if ( ! customBins ) customBins = static_cast<TAxis*>(ptAxis->Clone());
-
-
+        
         for ( Int_t ipt=0; ipt<=customBins->GetNbins(); ipt++ ) {
           Int_t minCustomBin = ( ipt == 0 ) ? 1 : ipt;
           Int_t maxCustomBin = ( ipt == 0 ) ? customBins->GetNbins() : ipt;
@@ -623,7 +701,7 @@ void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
           Int_t ptMaxBin = ptAxis->FindBin(customBins->GetBinUpEdge(maxCustomBin)-1.e-4);
           Double_t ptMin = ptAxis->GetBinLowEdge(ptMinBin);
           Double_t ptMax = ptAxis->GetBinUpEdge(ptMaxBin);
-          TH1* projHisto = hDeltaPhi2D->ProjectionX(Form("checkHisto_%s_ptBin%i", baseNameCent.Data(), ipt), ptMinBin, ptMaxBin);
+          TH1* projHisto = histo2D->ProjectionX(Form("%s_%s_ptBin%i", baseName.Data(), histo2D->GetName(), ipt), ptMinBin, ptMaxBin);
           projHisto->SetTitle(Form("%g < p_{t} (GeV/c) < %g", ptMin, ptMax));
           if ( projHisto->Integral() < 50 ) break;
           if ( ipad%4 == 0 ) {
@@ -642,8 +720,8 @@ void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
           projHisto->Fit(func,"R");
           for ( Int_t itype=0; itype<nTypes; itype++ ) {
             TGraphErrors* currGraph = ( ipt == 0 ) ? flowVsCentrality[itype] : flowVsPt[itype];
-            Double_t resoVal = ( itype == 0 ) ? 1. : weightedReso[0];
-            Double_t resoErr = ( itype == 0 ) ? 0. : weightedReso[itype];
+            Double_t resoVal = ( itype == 0 ) ? 1. : wgtReso[0][idx];
+            Double_t resoErr = ( itype == 0 ) ? 0. : wgtReso[itype][idx];
             Double_t rawVal = func->GetParameter(v2par);
             Double_t rawErr = (itype == 2 ) ? 0. : func->GetParError(v2par);
             Double_t yVal = rawVal/resoVal;
@@ -658,7 +736,7 @@ void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
               resoErr = 0.;
               TString graphTitle = currGraph->GetTitle();
               if ( ! graphTitle.Contains("|") ) {
-                graphTitle += Form(" | Rel. error on EP resolution: %.2g (stat.) %.2g (syst.)", weightedReso[1]/weightedReso[0], weightedReso[2]/weightedReso[0]);
+                graphTitle += Form(" | Rel. error on EP resolution: %.2g (stat.) %.2g (syst.)", wgtReso[1][idx]/wgtReso[0][idx], wgtReso[2][idx]/wgtReso[0][idx]);
                 currGraph->SetTitle(graphTitle.Data());
               }
             }
@@ -704,12 +782,9 @@ void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
       } // loop on type
       igroup1++;
       igroup2 = 0;
-    } // loop on track sources      
+    } // loop on track sources
   } // loop on steps
 
-  delete customBins;
-
-
   ///////////////////////
   // Event plane check //
   ///////////////////////
@@ -721,20 +796,47 @@ void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
   can->SetGridy();
   igroup2++;
 
+  const Int_t kNfits = 2;
+  TString fitTypeName[kNfits] = {"Cos","Sin"};
+
   Bool_t addOffsetToCentrality = kFALSE;
+  Double_t offsetStep = 0.1;
 
   TLegend* leg = new TLegend(0.7,0.7,0.92,0.92);
   leg->SetBorderSize(1);
+  
+  TGraphErrors* epCheck[kNfits];
+  for ( Int_t itype=0; itype<kNfits; itype++ ) {
+    histoName = Form("checkFourier_%s",fitTypeName[itype].Data());
+    epCheck[itype] = new TGraphErrors();
+    epCheck[itype]->SetName(histoName.Data());
+//    epCheck[itype]->SetTitle(histoName.Data());
+  }
   for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
     TH1* histo = static_cast<TH1*> ( GetSum(physSel,trigClassName,centralityClasses->At(icent)->GetName(),Form("hEventPlane%s", selectEP.Data())) );
     if ( ! histo ) continue;
     histo->SetName(Form("%s_%s", histo->GetName(), centralityClasses->At(icent)->GetName()));
     if ( histo->Integral() < 50. ) continue;
     if ( histo->GetSumw2N() == 0 ) histo->Sumw2();
+    
+    for ( Int_t itype=0; itype<kNfits; itype++ ) {
+      TString checkFormula = Form("[0]*(1.+2.*[1]*TMath::%s(2.*x))",fitTypeName[itype].Data());
+      TF1* checkFunc = new TF1(Form("checkFunc_%s", fitTypeName[itype].Data()), checkFormula.Data(), histo->GetXaxis()->GetBinLowEdge(1), histo->GetXaxis()->GetBinUpEdge(histo->GetXaxis()->GetNbins()));
+      checkFunc->SetParameters(histo->Integral(), 0.);
+      histo->Fit(checkFunc,"R0");
+      TGraphErrors* currGraph = epCheck[itype];
+      Double_t yVal = checkFunc->GetParameter(1);
+      Double_t yErr = checkFunc->GetParError(1);
+      Double_t xVal = 0.5*(centralityArray[icent+1] + centralityArray[icent]);
+      Double_t xErr = 0.5*(centralityArray[icent+1] - centralityArray[icent]);
+      Int_t ipoint = currGraph->GetN();
+      currGraph->SetPoint(ipoint,xVal,yVal);
+      currGraph->SetPointError(ipoint,xErr,yErr);
+    } // loop on type
     histo->Fit("pol0","R0");
     TF1* pol0 = histo->GetFunction("pol0");
     histo->Scale(1./pol0->GetParameter(0));
-    Int_t offset = ( addOffsetToCentrality ) ? icolor : 0;
+    Double_t offset = ( addOffsetToCentrality ) ? offsetStep*(Double_t)icolor : 0.;
     TGraphErrors* graph = new TGraphErrors();
     graph->SetTitle(Form("%s %s", histo->GetTitle(), trigClassName.Data()));
     for ( Int_t ipoint=0; ipoint<histo->GetXaxis()->GetNbins(); ipoint++ ) {
@@ -746,9 +848,9 @@ void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
     graph->GetYaxis()->SetTitle("Data/Fit (pol0)");
     TString legTitle = Form("%s %%", centralityClasses->At(icent)->GetName());
     if ( addOffsetToCentrality ) {
-      graph->GetYaxis()->SetRangeUser(0, 1.5*(Double_t)centralityClasses->GetEntries());
+      graph->GetYaxis()->SetRangeUser(1.-offsetStep, 1.+1.5*offsetStep*(Double_t)centralityClasses->GetEntries());
       graph->GetYaxis()->SetNdivisions(2*centralityClasses->GetEntries(),0,0);
-      legTitle += Form(" ( + %i)", offset);
+      legTitle += Form(" ( + %g)", offset);
     }
     Int_t currColor = ( icolor < kNtrackSources ) ? srcColors[icolor] : 20+icolor;
     graph->SetLineColor(currColor);
@@ -768,6 +870,94 @@ void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
     icolor++;
   }
   leg->Draw("same");
+  
+  currName = "epCheckCan";
+  can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
+  can->SetGridy();
+  igroup2++;
+  leg = new TLegend(0.7,0.7,0.95,0.95);
+  leg->SetBorderSize(1);
+  for ( Int_t itype=0; itype<kNfits; itype++ ) {
+    epCheck[itype]->SetLineColor(srcColors[itype]);
+    epCheck[itype]->SetMarkerColor(srcColors[itype]);
+    epCheck[itype]->SetMarkerStyle(20+itype);
+    epCheck[itype]->GetXaxis()->SetTitle("Centrality (%)");
+    TString currDrawOpt = "pz";
+    if ( can->GetListOfPrimitives()->GetEntries() == 0 ) currDrawOpt.Prepend("a");
+    epCheck[itype]->Draw(currDrawOpt.Data());
+    leg->AddEntry(epCheck[itype],Form("<%s(2 #psi_{%s})>", fitTypeName[itype].Data(), selectEP.Data()),"lp");
+  }
+  leg->Draw("same");
+  
+  for ( Int_t istep=0; istep<kNsteps; istep++ ) {
+    for ( Int_t isrc = firstSrc; isrc <= lastSrc; ++isrc ) {
+      TString baseName = Form("%s_%s", stepName[istep].Data(), fSrcKeys->At(isrc)->GetName());
+      for ( Int_t icent=0; icent<centralityClasses->GetEntries(); icent++ ) {
+        Int_t idx = kNcentralities * kNtrackSources * istep + kNcentralities * isrc + icent;
+        TH2* histo2D = static_cast<TH2*> (ptVsPhiList.At(idx));
+        if ( ! histo2D ) continue;
+        TString baseNameCent = Form("%s_%s", baseName.Data(), centralityClasses->At(icent)->GetName());
+        
+        TAxis* ptAxis = histo2D->GetYaxis();
+        
+        Int_t ipad = 0;
+        TGraphErrors* detEffect[kNfits];
+        for ( Int_t itype=0; itype<kNfits; itype++ ) {
+          histoName = Form("checkResoVsPt_%s_%s", baseNameCent.Data(), fitTypeName[itype].Data());
+          detEffect[itype] = new TGraphErrors();
+          detEffect[itype]->SetName(histoName.Data());
+          detEffect[itype]->SetTitle(histoName.Data());
+        }
+        
+        for ( Int_t ipt=0; ipt<=customBins->GetNbins(); ipt++ ) {
+          Int_t minCustomBin = ( ipt == 0 ) ? 1 : ipt;
+          Int_t maxCustomBin = ( ipt == 0 ) ? customBins->GetNbins() : ipt;
+          Int_t ptMinBin = ptAxis->FindBin(customBins->GetBinLowEdge(minCustomBin)+1.e-4);
+          Int_t ptMaxBin = ptAxis->FindBin(customBins->GetBinUpEdge(maxCustomBin)-1.e-4);
+          Double_t ptMin = ptAxis->GetBinLowEdge(ptMinBin);
+          Double_t ptMax = ptAxis->GetBinUpEdge(ptMaxBin);
+          TH1* projHisto = histo2D->ProjectionX(Form("checkReso_%s_%s_ptBin%i", baseName.Data(), histo2D->GetName(), ipt), ptMinBin, ptMaxBin);
+          projHisto->SetTitle(Form("%g < p_{t} (GeV/c) < %g", ptMin, ptMax));
+          if ( projHisto->Integral() < 50 ) break;
+          if ( ipad%4 == 0 ) {
+            currName = projHisto->GetName();
+            currName.Append("_can");
+            showCan = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
+            showCan->Divide(2,2);
+            ipad = 0;
+          }
+          showCan->cd(++ipad);
+          for ( Int_t itype=0; itype<kNfits; itype++ ) {
+            TString checkFormula = Form("[0]*(1.+2.*[1]*TMath::%s(2.*x))",fitTypeName[itype].Data());
+            TF1* checkFunc = new TF1(Form("checkFunc_%s", fitTypeName[itype].Data()), checkFormula.Data(), projHisto->GetXaxis()->GetBinLowEdge(1), projHisto->GetXaxis()->GetBinUpEdge(projHisto->GetXaxis()->GetNbins()));
+            checkFunc->SetParameters(projHisto->Integral(), 0.);
+            projHisto->Fit(checkFunc,"R");
+            TGraphErrors* currGraph = detEffect[itype];
+            Double_t yVal = checkFunc->GetParameter(1);
+            Double_t yErr = checkFunc->GetParError(1);
+            Double_t xVal = 0.5*(ptMax + ptMin);
+            Double_t xErr = 0.5*(ptMax - ptMin);
+            Int_t ipoint = currGraph->GetN();
+            currGraph->SetPoint(ipoint,xVal,yVal);
+            currGraph->SetPointError(ipoint,xErr,yErr);
+          } // loop on type
+        } // loop on pt bins
+        for ( Int_t itype=0; itype<kNfits; itype++ ) {
+          currName = detEffect[itype]->GetName();
+          currName.Append("Can");
+          can = new TCanvas(currName.Data(),currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
+          igroup2++;
+          can->cd();
+          detEffect[itype]->GetXaxis()->SetTitle(ptAxis->GetTitle());
+          detEffect[itype]->GetYaxis()->SetTitle("v2");
+//          detEffect[itype]->GetYaxis()->SetRangeUser(0., 0.25);
+          detEffect[itype]->SetFillStyle(0);
+          detEffect[itype]->Draw("ap");
+//          outList.Add(detEffect[itype]);
+        } // loop on types
+      } // loop on centrality
+    } // loop on track sources
+  } // loop on steps
 
 
   //////////////////////
@@ -794,4 +984,6 @@ void AliAnalysisTaskFlowSingleMu::Terminate(Option_t *) {
     outList.Write();
     file->Close();
   }
+  
+  delete customBins;
 }
index 8a2dfaf..fcdabda 100644 (file)
@@ -30,30 +30,12 @@ class AliAnalysisTaskFlowSingleMu : public AliVAnalysisMuon {
   void MyUserCreateOutputObjects();
   void ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality);
 
-  /*
-  enum {
-    kEPV0A, ///< EP form V0A
-    kEPTPC, ///< EP form TPC
-    kEPrandom ///< Random EP
-  };
-
-  void SetEPtype ( Int_t epType = kEPV0A ) { fEPtype = epType; }
-  */
-
  private:
 
   AliAnalysisTaskFlowSingleMu(const AliAnalysisTaskFlowSingleMu&);
   AliAnalysisTaskFlowSingleMu& operator=(const AliAnalysisTaskFlowSingleMu&);
 
   TArrayD GetCentralityRange(TString sRange);
-
-  /*
-  enum {
-    kTrackContainer, ///< CF container for tracks
-    kHistoEP,      ///< Event plane distribution
-    kNobjectTypes    ///< Number of objects
-  };
-  */
   
   enum {
     kStepReconstructed,  ///< Reconstructed tracks
index f59ef09..687d216 100644 (file)
@@ -70,7 +70,9 @@
 #include "AliVAnalysisMuon.h"
 #include "AliMergeableCollection.h"
 #include "AliCounterCollection.h"
+#include "AliMuonEventCuts.h"
 #include "AliMuonTrackCuts.h"
+#include "AliAnalysisMuonUtility.h"
 
 
 /// \cond CLASSIMP
@@ -191,17 +193,8 @@ void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& sel
   /// Fill output objects
   //
 
-  if ( GetVertexSPD()->GetNContributors() < fMinNvtxContirbutors ) return;
-
-  Double_t ipVz = GetVertexSPD()->GetZ();
-  Double_t ipVzMC = 0;
-  if ( IsMC() ) {
-    if ( fMCEvent ) ipVzMC = fMCEvent->GetPrimaryVertex()->GetZ();
-    else if ( fAODEvent ) {
-      AliAODMCHeader* aodMCHeader = (AliAODMCHeader *)fAODEvent->FindListObject(AliAODMCHeader::StdBranchName());
-      if ( aodMCHeader ) ipVzMC = aodMCHeader->GetVtxZ();
-    }
-  }
+  Double_t ipVz = AliAnalysisMuonUtility::GetVertexSPD(InputEvent())->GetZ();
+  Double_t ipVzMC = AliAnalysisMuonUtility::IsMCEvent(InputEvent(),MCEvent()) ? AliAnalysisMuonUtility::GetMCVertexZ(InputEvent(),MCEvent()) : 0.;
   
   for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
     TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
@@ -215,9 +208,9 @@ void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& sel
   AliVParticle* track = 0x0;
 
   for ( Int_t istep = 0; istep<2; ++istep ) {
-    Int_t nTracks = ( istep == kStepReconstructed ) ? GetNTracks() : GetNMCTracks();
+    Int_t nTracks = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetNTracks(InputEvent()) : AliAnalysisMuonUtility::GetNMCTracks(InputEvent(),MCEvent());
     for (Int_t itrack = 0; itrack < nTracks; itrack++) {
-      track = ( istep == kStepReconstructed ) ? GetTrack(itrack) : GetMCTrack(itrack);
+      track = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetTrack(itrack,InputEvent()) : AliAnalysisMuonUtility::GetMCTrack(itrack,InputEvent(),MCEvent());
       
       Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 );
       if ( ! isSelected ) continue;
@@ -225,12 +218,12 @@ void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& sel
       // In W simulations with Pythia, sometimes muon is stored twice.
       // Remove muon in case it has another muon as daugther
       if ( istep == kStepGeneratedMC ) {
-        Int_t firstDaughter = GetDaughterIndex(track, 0);
+        Int_t firstDaughter = AliAnalysisMuonUtility::GetDaughterIndex(track, 0);
         if ( firstDaughter >= 0 ) {
           Bool_t hasMuonDaughter = kFALSE;
-          Int_t lastDaughter = GetDaughterIndex(track, 1);
+          Int_t lastDaughter = AliAnalysisMuonUtility::GetDaughterIndex(track, 1);
           for ( Int_t idaugh=firstDaughter; idaugh<=lastDaughter; idaugh++ ) {
-            AliVParticle* currTrack = GetMCTrack(idaugh);
+            AliVParticle* currTrack = AliAnalysisMuonUtility::GetMCTrack(idaugh,InputEvent(),MCEvent());
             if ( currTrack->PdgCode() == track->PdgCode() ) {
               hasMuonDaughter = kTRUE;
               break;
@@ -247,8 +240,7 @@ void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& sel
       
       Double_t thetaAbsEndDeg = 0;
       if ( istep == kStepReconstructed ) {
-        Double_t rAbsEnd =  ( fAODEvent ) ? ((AliAODTrack*)track)->GetRAtAbsorberEnd(): ((AliESDMuonTrack*)track)->GetRAtAbsorberEnd();
-        thetaAbsEndDeg = TMath::ATan( rAbsEnd / 505. ) * TMath::RadToDeg();
+        thetaAbsEndDeg = AliAnalysisMuonUtility::GetThetaAbsDeg(track);
       }
       else {
         thetaAbsEndDeg = ( TMath::Pi()-track->Theta() ) * TMath::RadToDeg();
@@ -265,7 +257,7 @@ void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& sel
       
       for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
         TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
-        if ( istep == kStepReconstructed && ! TrackPtCutMatchTrigClass(track, trigClassName) ) continue;
+        if ( istep == kStepReconstructed && ! fMuonTrackCuts->TrackPtCutMatchTrigClass(track, fMuonEventCuts->GetTrigClassPtCutLevel(trigClassName)) ) continue;
         ((AliCFContainer*)GetMergeableObject(physSel, trigClassName, centrality, "SingleMuContainer"))->Fill(containerInput,istep);
       } // loop on selected trigger classes
     } // loop on tracks
index 7e01ac4..3340d8f 100644 (file)
@@ -18,6 +18,7 @@
 #include "TMath.h"
 #include "TList.h"
 #include "TLorentzVector.h"
+#include "TArrayI.h"
 
 #include "AliLog.h"
 #include "AliVParticle.h"
@@ -190,6 +191,28 @@ void AliMuonPairCuts::SetIsMC ( Bool_t isMC )
 }
 
 //________________________________________________________________________
+Bool_t AliMuonPairCuts::TrackPtCutMatchTrigClass ( const AliVParticle* track1, const AliVParticle* track2, const TArrayI ptCutFromClass ) const
+{
+  /// Check if track pair passes the trigger pt cut level used in the trigger class
+  Bool_t matchTrig1 = fMuonTrackCuts.TrackPtCutMatchTrigClass(track1, ptCutFromClass);
+  Bool_t matchTrig2 = fMuonTrackCuts.TrackPtCutMatchTrigClass(track2, ptCutFromClass);
+  
+  Bool_t matchTrackerPt1 = kTRUE, matchTrackerPt2 = kTRUE;
+  if ( IsApplySharpPtCutInMatching() ) {
+    matchTrackerPt1 = ( track1->Pt() >= fMuonTrackCuts.GetSharpPtCut(ptCutFromClass[0]-1,kFALSE) );
+    matchTrackerPt2 = ( track2->Pt() >= fMuonTrackCuts.GetSharpPtCut(ptCutFromClass[0]-1,kFALSE) );
+  }
+  
+  matchTrig1 = ( matchTrig1 && matchTrackerPt1 );
+  matchTrig2 = ( matchTrig2 && matchTrackerPt2 );
+  
+  Bool_t passCut = ( ptCutFromClass[1]>0 ) ? ( matchTrig1 && matchTrig2 ) : ( matchTrig1 || matchTrig2 );
+
+  AliDebug(1,Form("Class matchTrig %i %i  trackMatchTrig %i %i trackPt %g %g (required %i)  passCut %i", ptCutFromClass[0], ptCutFromClass[1], matchTrig1, matchTrig2, track1->Pt(), track2->Pt(), IsApplySharpPtCutInMatching(),passCut));
+  return passCut;
+}
+
+//________________________________________________________________________
 void AliMuonPairCuts::Print(Option_t* option) const
 {
   //
index 88bdd0f..fa77482 100644 (file)
@@ -3,10 +3,10 @@
 
 #include "AliAnalysisCuts.h"
 #include "AliMuonTrackCuts.h"
-#include "TArrayD.h"
 
 class TList;
 class TVector3;
+class TArrayI;
 
 class AliMuonPairCuts : public AliAnalysisCuts
 {
@@ -55,6 +55,7 @@ class AliMuonPairCuts : public AliAnalysisCuts
   /// Get flag to apply the sharp pt cut when matching with trigger
   Bool_t IsApplySharpPtCutInMatching () const { return fMuonTrackCuts.IsApplySharpPtCutInMatching(); }
   
+  Bool_t TrackPtCutMatchTrigClass ( const AliVParticle* track1, const AliVParticle* track2, const TArrayI ptCutFromClass ) const;
   
   /// Return the MuonTrackCuts (normally the standard user do not need to modify this data member)
   AliMuonTrackCuts& GetMuonTrackCuts () { return fMuonTrackCuts; }
index 47600ba..6a022d8 100644 (file)
@@ -18,6 +18,7 @@
 #include "TMath.h"
 #include "TList.h"
 #include "TArrayD.h"
+#include "TArrayI.h"
 #include "TObjArray.h"
 #include "TObjString.h"
 #include "TFile.h"
@@ -31,6 +32,7 @@
 #include "AliESDMuonTrack.h"
 #include "AliAODTrack.h"
 #include "AliAnalysisManager.h"
+#include "AliAnalysisMuonUtility.h"
 
 /// \cond CLASSIMP
 ClassImp(AliMuonTrackCuts) // Class implementation in ROOT context
@@ -40,7 +42,6 @@ ClassImp(AliMuonTrackCuts) // Class implementation in ROOT context
 //________________________________________________________________________
 AliMuonTrackCuts::AliMuonTrackCuts() :
   AliAnalysisCuts(),
-  fIsESD(kFALSE),
   fIsMC(kFALSE),
   fUseCustomParam(kFALSE),
   fSharpPtCut(kFALSE),
@@ -53,29 +54,12 @@ AliMuonTrackCuts::AliMuonTrackCuts() :
 //________________________________________________________________________
 AliMuonTrackCuts::AliMuonTrackCuts(const char* name, const char* title ) :
 AliAnalysisCuts(name, title),
-fIsESD(kFALSE),
-fIsMC(kFALSE),
-fUseCustomParam(kFALSE),
-fSharpPtCut(kFALSE),
-fParameters(TArrayD(kNParameters))
-{
-  /// Constructor
-  fParameters.Reset();
-  SetDefaultFilterMask();
-}
-
-
-//________________________________________________________________________
-AliMuonTrackCuts::AliMuonTrackCuts(const char* name, const char* title, Bool_t isESD ) :
-  AliAnalysisCuts(name, title),
-  fIsESD(isESD),
   fIsMC(kFALSE),
   fUseCustomParam(kFALSE),
   fSharpPtCut(kFALSE),
   fParameters(TArrayD(kNParameters))
 {
-  /// Obsolete Constructor
-  AliWarning(Form("\n\n *****  This constructor is obsolete and will be removed!\nPlease use AliMuonTrackCuts(name, title) instead!\n"));
+  /// Constructor
   fParameters.Reset();
   SetDefaultFilterMask();
 }
@@ -84,7 +68,6 @@ AliMuonTrackCuts::AliMuonTrackCuts(const char* name, const char* title, Bool_t i
 //________________________________________________________________________
 AliMuonTrackCuts::AliMuonTrackCuts(const AliMuonTrackCuts& obj) :
   AliAnalysisCuts(obj),
-  fIsESD(obj.fIsESD),
   fIsMC(obj.fIsMC),
   fUseCustomParam(obj.fUseCustomParam),
   fSharpPtCut(obj.fSharpPtCut),
@@ -100,7 +83,6 @@ AliMuonTrackCuts& AliMuonTrackCuts::operator=(const AliMuonTrackCuts& obj)
   /// Assignment operator
   if ( this != &obj ) { 
     AliAnalysisCuts::operator=(obj);
-    fIsESD = obj.fIsESD;
     fIsMC = obj.fIsMC;
     fUseCustomParam = obj.fUseCustomParam;
     fSharpPtCut = obj.fSharpPtCut;
@@ -117,13 +99,6 @@ AliMuonTrackCuts::~AliMuonTrackCuts()
 }
 
 //________________________________________________________________________
-Bool_t AliMuonTrackCuts::IsESDTrack( const AliVParticle* track ) const
-{
-  /// Check if track is from ESD or AOD
-  return ( track->IsA() != AliAODTrack::Class() );
-}
-
-//________________________________________________________________________
 Bool_t AliMuonTrackCuts::RunMatchesRange( Int_t runNumber, const Char_t* objName ) const
 {
   /// Check if the object contains the run
@@ -262,7 +237,7 @@ Bool_t AliMuonTrackCuts::StreamParameters( Int_t runNumber,  Int_t runMax )
 }
 
 //________________________________________________________________________
-Bool_t AliMuonTrackCuts::SetParameter(Int_t iparam, Float_t value)
+Bool_t AliMuonTrackCuts::SetParameter(Int_t iparam, Double_t value)
 {
   /// Set parameter
   if ( fUseCustomParam ) {
@@ -293,23 +268,18 @@ UInt_t AliMuonTrackCuts::GetSelectionMask( const TObject* obj )
   
   const AliVParticle* track = static_cast<const AliVParticle*> ( obj );
   
-  Bool_t isESD = IsESDTrack(track);
-
   UInt_t selectionMask = 0;
 
-  Bool_t isMuon = ( isESD ) ? ((AliESDMuonTrack*)track)->ContainTrackerData() :  ((AliAODTrack*)track)->IsMuonTrack();
-
-  if ( ! isMuon ) return selectionMask;
+  if ( ! AliAnalysisMuonUtility::IsMuonTrack(track) ) return selectionMask;
 
   Double_t eta = track->Eta();
   if ( eta > -4. && eta < -2.5 ) selectionMask |= kMuEta;
 
-  Double_t rAbsEnd =  ( isESD ) ? ((AliESDMuonTrack*)track)->GetRAtAbsorberEnd() : ((AliAODTrack*)track)->GetRAtAbsorberEnd();
-  Double_t thetaAbsEndDeg = TMath::ATan( rAbsEnd / 505. ) * TMath::RadToDeg();
+  Double_t thetaAbsEndDeg = AliAnalysisMuonUtility::GetThetaAbsDeg(track);
 
   if ( thetaAbsEndDeg > 2. && thetaAbsEndDeg < 10. ) selectionMask |= kMuThetaAbs;
 
-  Int_t matchTrig = ( isESD ) ? ((AliESDMuonTrack*)track)->GetMatchTrigger() : ((AliAODTrack*)track)->GetMatchTrigger();
+  Int_t matchTrig = AliAnalysisMuonUtility::GetMatchTrigger(track);
   Int_t cutLevel[3] = {kMuMatchApt, kMuMatchLpt, kMuMatchHpt};
   Double_t pt = track->Pt();
   for ( Int_t ilevel=0; ilevel<3; ilevel++ ) {
@@ -318,8 +288,7 @@ UInt_t AliMuonTrackCuts::GetSelectionMask( const TObject* obj )
     selectionMask |= cutLevel[ilevel];
   }
 
-  Double_t chi2norm = ( isESD ) ? ((AliESDMuonTrack*)track)->GetNormalizedChi2() : ((AliAODTrack*)track)->Chi2perNDF();
-  if ( chi2norm < GetChi2NormCut() ) selectionMask |= kMuTrackChiSquare;
+  if ( AliAnalysisMuonUtility::GetChi2perNDFtracker(track) < GetChi2NormCut() ) selectionMask |= kMuTrackChiSquare;
 
   TVector3 dcaAtVz = GetCorrectedDCA(track);
   Double_t pTotMean = GetAverageMomentum(track);
@@ -327,6 +296,7 @@ UInt_t AliMuonTrackCuts::GetSelectionMask( const TObject* obj )
   Double_t pDca = pTotMean * dcaAtVz.Mag();
     
   Double_t pTot = track->P();
+  Double_t rAbsEnd = AliAnalysisMuonUtility::GetRabs(track);
   
   // Momentum resolution only
   // The cut depends on the momentum resolution. In particular:
@@ -392,22 +362,10 @@ Int_t AliMuonTrackCuts::GetThetaAbsBin ( Double_t rAtAbsEnd ) const
 TVector3 AliMuonTrackCuts::GetCorrectedDCA ( const AliVParticle* track ) const
 {
   /// Get corrected DCA
-
-  Bool_t isESD = IsESDTrack(track);
   
-  Double_t vtxPos[3];
-  if ( isESD ) {
-    vtxPos[0] = ((AliESDMuonTrack*)track)->GetNonBendingCoor();
-    vtxPos[1] = ((AliESDMuonTrack*)track)->GetBendingCoor();
-    vtxPos[2] = ((AliESDMuonTrack*)track)->GetZ();
-  }
-  else ((AliAODTrack*)track)->GetXYZ(vtxPos);
+  TVector3 vertex(AliAnalysisMuonUtility::GetXatVertex(track), AliAnalysisMuonUtility::GetYatVertex(track), AliAnalysisMuonUtility::GetZatVertex(track));
   
-  TVector3 vertex(vtxPos);
-  
-  TVector3 dcaTrack(0.,0., vtxPos[2]);
-  dcaTrack.SetX( ( isESD ) ? ((AliESDMuonTrack*)track)->GetNonBendingCoorAtDCA() : ((AliAODTrack*)track)->XAtDCA() );
-  dcaTrack.SetY( ( isESD ) ? ((AliESDMuonTrack*)track)->GetBendingCoorAtDCA() : ((AliAODTrack*)track)->YAtDCA() );
+  TVector3 dcaTrack(AliAnalysisMuonUtility::GetXatDCA(track), AliAnalysisMuonUtility::GetYatDCA(track), AliAnalysisMuonUtility::GetZatDCA(track));
   
   TVector3 dcaAtVz = dcaTrack - vertex - GetMeanDCA();
 
@@ -418,13 +376,11 @@ TVector3 AliMuonTrackCuts::GetCorrectedDCA ( const AliVParticle* track ) const
 Double_t AliMuonTrackCuts::GetAverageMomentum ( const AliVParticle* track ) const
 {
   /// Get average momentum before and after the absorber
-
-  Bool_t isESD = IsESDTrack(track);
   
   Double_t pTotMean = 0.;
   Double_t pTot = track->P();
   //if ( isESD ) pTotMean = 0.5 * ( pTot + ((AliESDMuonTrack*)track)->PUncorrected() );
-  if ( isESD ) pTotMean = ((AliESDMuonTrack*)track)->PUncorrected(); // Increased stability if using uncorrected value
+  if ( ! AliAnalysisMuonUtility::IsAODTrack(track) ) pTotMean = ((AliESDMuonTrack*)track)->PUncorrected(); // Increased stability if using uncorrected value
   else {
     pTotMean = pTot - GetMeanPCorr(((AliAODTrack*)track)->GetRAtAbsorberEnd());
   }
@@ -564,7 +520,22 @@ void AliMuonTrackCuts::SetDefaultFilterMask ()
 {
   /// Standard cuts for single muon
   SetFilterMask ( kMuEta | kMuThetaAbs | kMuPdca | kMuMatchApt );
-}  
+}
+
+//________________________________________________________________________
+Bool_t AliMuonTrackCuts::TrackPtCutMatchTrigClass ( const AliVParticle* track, const TArrayI ptCutFromClass) const
+{
+  /// Check if track passes the trigger pt cut level used in the trigger class
+  Int_t matchTrig = AliAnalysisMuonUtility::GetMatchTrigger(track);
+  Bool_t matchTrackerPt = kTRUE;
+  if ( IsApplySharpPtCutInMatching() ) {
+    matchTrackerPt = ( track->Pt() >= GetSharpPtCut(ptCutFromClass[0]-1,kFALSE) );
+  }
+  Bool_t passCut = ( ( matchTrig >= ptCutFromClass[0] ) && matchTrackerPt );
+  AliDebug(1,Form("Class matchTrig %i %i  trackMatchTrig %i  trackPt %g (required %i)  passCut %i", ptCutFromClass[0], ptCutFromClass[1], matchTrig, track->Pt(), IsApplySharpPtCutInMatching(),passCut));
+  return passCut;
+}
+
 
 //________________________________________________________________________
 void AliMuonTrackCuts::Print(Option_t* option) const
index ed3b222..9c2e109 100644 (file)
@@ -7,6 +7,7 @@
 class AliVParticle;
 class TList;
 class TVector3;
+class TArrayI;
 
 class AliMuonTrackCuts : public AliAnalysisCuts
 {
@@ -24,7 +25,6 @@ class AliMuonTrackCuts : public AliAnalysisCuts
   
   AliMuonTrackCuts();
   AliMuonTrackCuts(const char* name, const char* title);
-  AliMuonTrackCuts(const char* name, const char* title, Bool_t isESD); // Obsolete
   AliMuonTrackCuts(const AliMuonTrackCuts& obj);
   AliMuonTrackCuts& operator=(const AliMuonTrackCuts& obj);
 
@@ -41,6 +41,8 @@ class AliMuonTrackCuts : public AliAnalysisCuts
   void SetIsMC(Bool_t isMC = kTRUE) { fIsMC = isMC; }
 
   void Print ( Option_t* option = "" ) const;
+  
+  Bool_t TrackPtCutMatchTrigClass ( const AliVParticle* track, const TArrayI ptCutFromClass) const;
 
   TVector3 GetCorrectedDCA ( const AliVParticle* track ) const;
   Double_t GetAverageMomentum ( const AliVParticle* track ) const;
@@ -104,18 +106,16 @@ class AliMuonTrackCuts : public AliAnalysisCuts
  private:
   
   Int_t GetThetaAbsBin ( Double_t rAtAbsEnd ) const;
-  Bool_t SetParameter ( Int_t iparam, Float_t value );
+  Bool_t SetParameter ( Int_t iparam, Double_t value );
   Bool_t RunMatchesRange ( Int_t runNumber, const Char_t* objName ) const;
-  Bool_t IsESDTrack ( const AliVParticle* track ) const;
 
-  Bool_t fIsESD;            ///< Event is ESD
   Bool_t fIsMC;             ///< Monte Carlo analysis
   Bool_t fUseCustomParam;   ///< Use custom parameters (do not search in OADB)
   Bool_t fSharpPtCut;       ///< Flag to apply sharp pt cut in track-trigger matching
 
   TArrayD fParameters;      ///< List of parameters
 
-  ClassDef(AliMuonTrackCuts, 2); // Class for muon track filters
+  ClassDef(AliMuonTrackCuts, 3); // Class for muon track filters
 };
  
 #endif
index cb979db..bee6e34 100644 (file)
 
 // PWG3 includes
 #include "AliMergeableCollection.h"
+#include "AliMuonEventCuts.h"
 #include "AliMuonTrackCuts.h"
 #include "AliMuonPairCuts.h"
+#include "AliAnalysisMuonUtility.h"
 
 /// \cond CLASSIMP
 ClassImp(AliVAnalysisMuon) // Class implementation in ROOT context
@@ -77,6 +79,7 @@ ClassImp(AliVAnalysisMuon) // Class implementation in ROOT context
 //________________________________________________________________________
 AliVAnalysisMuon::AliVAnalysisMuon() :
   AliAnalysisTaskSE(),
+  fMuonEventCuts(0x0),
   fMuonTrackCuts(0x0),
   fMuonPairCuts(0x0),
   fESDEvent(0x0),
@@ -85,15 +88,9 @@ AliVAnalysisMuon::AliVAnalysisMuon() :
   fChargeKeys(0x0),
   fSrcKeys(0x0),
   fPhysSelKeys(0x0),
-  fTriggerClasses(0x0),
-  fCentralityClasses(0x0),
   fEventCounters(0x0),
   fMergeableCollection(0x0),
   fOutputList(0x0),
-  fMinNvtxContirbutors(0),
-  fSelectedTrigPattern(0x0),
-  fRejectedTrigPattern(0x0),
-  fSelectedTrigLevel(0x0),
   fOutputPrototypeList(0x0)
 {
   /// Default ctor.
@@ -102,6 +99,7 @@ AliVAnalysisMuon::AliVAnalysisMuon() :
 //________________________________________________________________________
 AliVAnalysisMuon::AliVAnalysisMuon(const char *name, const AliMuonTrackCuts& trackCuts, const AliMuonPairCuts& pairCuts) :
   AliAnalysisTaskSE(name),
+  fMuonEventCuts(new AliMuonEventCuts("stdEventCuts","stdEventCuts")),
   fMuonTrackCuts(new AliMuonTrackCuts(trackCuts)),
   fMuonPairCuts(new AliMuonPairCuts(pairCuts)),
   fESDEvent(0x0),
@@ -110,25 +108,17 @@ AliVAnalysisMuon::AliVAnalysisMuon(const char *name, const AliMuonTrackCuts& tra
   fChargeKeys(0x0),
   fSrcKeys(0x0),
   fPhysSelKeys(0x0),
-  fTriggerClasses(new THashList()),
-  fCentralityClasses(0x0),
   fEventCounters(0x0),
   fMergeableCollection(0x0),
   fOutputList(0x0),
-  fMinNvtxContirbutors(1),
-  fSelectedTrigPattern(new TObjArray()),
-  fRejectedTrigPattern(new TObjArray()),
-  fSelectedTrigLevel(new TObjArray()),
   fOutputPrototypeList(0x0)
 {
   //
   /// Constructor.
   //
   
-  fTriggerClasses->SetOwner();
   InitKeys();
-  SetTrigClassPatterns();
-  SetTrigClassLevels();
+  SetTrigClassPatterns("");
   SetCentralityClasses();
 
   DefineOutput(1, TObjArray::Class());
@@ -137,6 +127,7 @@ AliVAnalysisMuon::AliVAnalysisMuon(const char *name, const AliMuonTrackCuts& tra
 //________________________________________________________________________
 AliVAnalysisMuon::AliVAnalysisMuon(const char *name, const AliMuonTrackCuts& trackCuts) :
   AliAnalysisTaskSE(name),
+  fMuonEventCuts(new AliMuonEventCuts("stdEventCuts","stdEventCuts")),
   fMuonTrackCuts(new AliMuonTrackCuts(trackCuts)),
   fMuonPairCuts(0x0),
   fESDEvent(0x0),
@@ -145,25 +136,17 @@ AliVAnalysisMuon::AliVAnalysisMuon(const char *name, const AliMuonTrackCuts& tra
   fChargeKeys(0x0),
   fSrcKeys(0x0),
   fPhysSelKeys(0x0),
-  fTriggerClasses(new THashList()),
-  fCentralityClasses(0x0),
   fEventCounters(0x0),
   fMergeableCollection(0x0),
   fOutputList(0x0),
-  fMinNvtxContirbutors(1),
-  fSelectedTrigPattern(new TObjArray()),
-  fRejectedTrigPattern(new TObjArray()),
-  fSelectedTrigLevel(new TObjArray()),
   fOutputPrototypeList(0x0)
 {
   //
   /// Constructor.
   //
   
-  fTriggerClasses->SetOwner();
   InitKeys();
-  SetTrigClassPatterns();
-  SetTrigClassLevels();
+  SetTrigClassPatterns("");
   SetCentralityClasses();
   
   DefineOutput(1, TObjArray::Class());
@@ -173,6 +156,7 @@ AliVAnalysisMuon::AliVAnalysisMuon(const char *name, const AliMuonTrackCuts& tra
 //________________________________________________________________________
 AliVAnalysisMuon::AliVAnalysisMuon(const char *name, const AliMuonPairCuts& pairCuts) :
   AliAnalysisTaskSE(name),
+  fMuonEventCuts(new AliMuonEventCuts("stdEventCuts","stdEventCuts")),
   fMuonTrackCuts(0x0),
   fMuonPairCuts(new AliMuonPairCuts(pairCuts)),
   fESDEvent(0x0),
@@ -181,24 +165,16 @@ AliVAnalysisMuon::AliVAnalysisMuon(const char *name, const AliMuonPairCuts& pair
   fChargeKeys(0x0),
   fSrcKeys(0x0),
   fPhysSelKeys(0x0),
-  fTriggerClasses(new THashList()),
-  fCentralityClasses(0x0),
   fEventCounters(0x0),
   fMergeableCollection(0x0),
   fOutputList(0x0),
-  fMinNvtxContirbutors(1),
-  fSelectedTrigPattern(new TObjArray()),
-  fRejectedTrigPattern(new TObjArray()),
-  fSelectedTrigLevel(new TObjArray()),
   fOutputPrototypeList(0x0)
 {
   //
   /// Constructor.
   //
-  fTriggerClasses->SetOwner();
   InitKeys();
-  SetTrigClassPatterns();
-  SetTrigClassLevels();
+  SetTrigClassPatterns("");
   SetCentralityClasses();
     
   DefineOutput(1, TObjArray::Class());
@@ -212,17 +188,13 @@ AliVAnalysisMuon::~AliVAnalysisMuon()
   /// Destructor
   //
 
+  delete fMuonEventCuts;
   delete fMuonTrackCuts;
   delete fMuonPairCuts;
   delete fTerminateOptions;
   delete fChargeKeys;
   delete fSrcKeys;
   delete fPhysSelKeys;
-  delete fTriggerClasses;
-  delete fCentralityClasses;
-  delete fSelectedTrigPattern;
-  delete fRejectedTrigPattern;
-  delete fSelectedTrigLevel;
   delete fOutputPrototypeList;
 
 
@@ -286,11 +258,11 @@ void AliVAnalysisMuon::UserCreateOutputObjects()
 
   fEventCounters = new AliCounterCollection("eventCounters");
 
-  if ( ! fCentralityClasses ) SetCentralityClasses();
+  if ( ! GetCentralityClasses() ) SetCentralityClasses();
   TString centralityClasses = "";
-  for ( Int_t icent=1; icent<=fCentralityClasses->GetNbins(); ++icent ) {
+  for ( Int_t icent=1; icent<=GetCentralityClasses()->GetNbins(); ++icent ) {
     if ( ! centralityClasses.IsNull() ) centralityClasses += "/";
-    centralityClasses += fCentralityClasses->GetBinLabel(icent);
+    centralityClasses += GetCentralityClasses()->GetBinLabel(icent);
   }
   fEventCounters->AddRubric("selected", "yes/no");
   fEventCounters->AddRubric("trigger", 100);
@@ -303,6 +275,8 @@ void AliVAnalysisMuon::UserCreateOutputObjects()
 
   PostData(1, fOutputList);
   
+  fMuonEventCuts->Print();
+  
   MyUserCreateOutputObjects();
 }
 
@@ -323,25 +297,19 @@ void AliVAnalysisMuon::UserExec(Option_t * /*option*/)
     AliError ("AOD or ESD event not found. Nothing done!");
     return;
   }
+  
+  if ( ! fMuonEventCuts->IsSelected(fInputHandler) ) return;
 
   Int_t physSel = ( fInputHandler->IsEventSelected() & AliVEvent::kAny ) ? kPhysSelPass : kPhysSelReject;
 
   //
   // Global event info
   //
+  TObjArray* selectTrigClasses = fMuonEventCuts->GetSelectedTrigClassesInEvent(InputEvent());
 
-  TString firedTrigClasses = ( fAODEvent ) ? fAODEvent->GetFiredTriggerClasses() : fESDEvent->GetFiredTriggerClasses();
-  firedTrigClasses.Prepend("ANY ");
-  AliDebug(2, Form("Fired classes %s", firedTrigClasses.Data()));
-  TObjArray* selectTrigClasses = BuildTriggerClasses(firedTrigClasses);
-  if ( selectTrigClasses->GetEntries() == 0 ) {
-    delete selectTrigClasses;
-    return;
-  }
-
-  Double_t centrality = InputEvent()->GetCentrality()->GetCentralityPercentile("V0M");
-  Int_t centralityBin = fCentralityClasses->FindBin(centrality);
-  TString centralityBinLabel = fCentralityClasses->GetBinLabel(centralityBin);
+  Double_t centrality = fMuonEventCuts->GetCentrality(InputEvent());
+  Int_t centralityBin = GetCentralityClasses()->FindBin(centrality);
+  TString centralityBinLabel = GetCentralityClasses()->GetBinLabel(centralityBin);
 
   TString selKey = ( physSel == kPhysSelPass ) ? "yes" : "no";
   for ( Int_t itrig=0; itrig<selectTrigClasses->GetEntries(); ++itrig ) {
@@ -349,9 +317,7 @@ void AliVAnalysisMuon::UserExec(Option_t * /*option*/)
     fEventCounters->Count(Form("trigger:%s/selected:%s/centrality:%s", trigName.Data(), selKey.Data(), centralityBinLabel.Data()));
   }
 
-  ProcessEvent(fPhysSelKeys->At(physSel)->GetName(), *selectTrigClasses, fCentralityClasses->GetBinLabel(centralityBin));
-
-  delete selectTrigClasses;
+  ProcessEvent(fPhysSelKeys->At(physSel)->GetName(), *selectTrigClasses, centralityBinLabel);
 
   // Post final data. It will be written to a file with option "RECREATE"
   PostData(1, fOutputList);
@@ -383,135 +349,6 @@ void AliVAnalysisMuon::Terminate(Option_t *)
 }
 
 
-
-//________________________________________________________________________
-Int_t AliVAnalysisMuon::GetNTracks()
-{
-  //
-  /// Return the number of tracks in event
-  //
-  return ( fAODEvent ) ? fAODEvent->GetNTracks() : fESDEvent->GetNumberOfMuonTracks();
-}
-
-
-//________________________________________________________________________
-AliVParticle* AliVAnalysisMuon::GetTrack(Int_t itrack)
-{
-  //
-  /// Get the current track
-  //
-  AliVParticle* track = 0x0;
-  if ( fAODEvent ) track = fAODEvent->GetTrack(itrack);
-  else track = fESDEvent->GetMuonTrack(itrack);
-  return track;
-}
-
-//________________________________________________________________________
-Double_t AliVAnalysisMuon::MuonMass2() const
-{
-  /// A usefull constant
-  static Double_t m2 = 1.11636129640000012e-02;
-  return m2;
-}
-
-//________________________________________________________________________
-TLorentzVector AliVAnalysisMuon::GetTrackPair(AliVParticle* track1, AliVParticle* track2) const
-{
-  //
-  /// Get track pair
-  //
-  
-  AliVParticle* tracks[2] = {track1, track2};
-  
-  TLorentzVector vec[2];
-  for ( Int_t itrack=0; itrack<2; ++itrack ) {
-    Double_t trackP = tracks[itrack]->P();
-    Double_t energy = TMath::Sqrt(trackP*trackP + MuonMass2());
-    vec[itrack].SetPxPyPzE(tracks[itrack]->Px(), tracks[itrack]->Py(), tracks[itrack]->Pz(), energy);
-  }
-  
-  TLorentzVector vecPair = vec[0] + vec[1];
-  return vecPair;
-}
-
-
-//________________________________________________________________________
-Int_t AliVAnalysisMuon::GetNMCTracks()
-{
-  //
-  /// Return the number of MC tracks in event
-  //
-  Int_t nMCtracks = 0;
-  if ( fMCEvent ) nMCtracks = fMCEvent->GetNumberOfTracks();
-  else if ( fAODEvent ) {
-    TClonesArray* mcArray = (TClonesArray*)fAODEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
-    if ( mcArray ) nMCtracks = mcArray->GetEntries();
-  }
-  return nMCtracks;
-}
-
-//________________________________________________________________________
-AliVParticle* AliVAnalysisMuon::GetMCTrack(Int_t trackLabel)
-{
-  //
-  /// MC information can be provided by the MC input handler
-  /// (mostly when analyising ESDs) or can be found inside AODs
-  /// This method returns the correct one
-  //
-  AliVParticle* mcTrack = 0x0;
-  if ( fMCEvent ) mcTrack = fMCEvent->GetTrack(trackLabel);
-  else if ( fAODEvent ) {
-    TClonesArray* mcArray = (TClonesArray*)fAODEvent->FindListObject(AliAODMCParticle::StdBranchName());
-    if ( mcArray ) mcTrack =  (AliVParticle*)mcArray->At(trackLabel);
-  }
-  if ( ! mcTrack ) AliWarning(Form("No track with label %i!", trackLabel));
-  return mcTrack;
-}
-
-//________________________________________________________________________
-Int_t AliVAnalysisMuon::GetMotherIndex(AliVParticle* mcParticle)
-{
-  //
-  /// Return the mother index
-  //
-  Int_t imother = ( fMCEvent ) ? ((AliMCParticle*)mcParticle)->GetMother() : ((AliAODMCParticle*)mcParticle)->GetMother();
-  return imother;
-}
-
-//________________________________________________________________________
-Int_t AliVAnalysisMuon::GetDaughterIndex(AliVParticle* mcParticle, Int_t idaughter)
-{
-  //
-  /// Return the daughter index
-  /// idaughter can be:
-  /// 0 -> first daughter
-  /// 1 -> last daughter
-  //
-  if ( idaughter < 0 || idaughter > 1 ) {
-    AliError(Form("Requested daughter %i Daughter index can be either 0 (first) or 1 (last)", idaughter));
-    return -1;
-  }
-  
-  if ( fMCEvent ) {
-    if ( idaughter == 0 ) return ((AliMCParticle*)mcParticle)->GetFirstDaughter();
-    else return ((AliMCParticle*)mcParticle)->GetLastDaughter();
-  }
-  
-  return ((AliAODMCParticle*)mcParticle)->GetDaughter(idaughter);
-}
-
-
-
-//________________________________________________________________________
-Bool_t AliVAnalysisMuon::IsMC()
-{
-  //
-  /// Contains MC info
-  //
-  return ( fMCEvent || ( fAODEvent && fAODEvent->FindListObject(AliAODMCParticle::StdBranchName()) ) );
-}
-
-
 //________________________________________________________________________
 Int_t AliVAnalysisMuon::GetParticleType(AliVParticle* track)
 {
@@ -522,7 +359,7 @@ Int_t AliVAnalysisMuon::GetParticleType(AliVParticle* track)
   Int_t trackSrc = kUnidentified;
   Int_t trackLabel = track->GetLabel();
   if ( trackLabel >= 0 ) {
-    AliVParticle* matchedMCTrack = GetMCTrack(trackLabel);
+    AliVParticle* matchedMCTrack = AliAnalysisMuonUtility::GetMCTrack(trackLabel,InputEvent(),MCEvent());
     if ( matchedMCTrack ) trackSrc = RecoTrackMother(matchedMCTrack);
   } // track has MC label
   return trackSrc;
@@ -541,18 +378,18 @@ Int_t AliVAnalysisMuon::RecoTrackMother(AliVParticle* mcParticle)
   // Track is not a muon
   if ( TMath::Abs(recoPdg) != 13 ) return kRecoHadron;
   
-  Int_t imother = GetMotherIndex(mcParticle);
+  Int_t imother = AliAnalysisMuonUtility::GetMotherIndex(mcParticle);
   
   Int_t den[3] = {100, 1000, 1};
   
   Int_t motherType = kDecayMu;
   while ( imother >= 0 ) {
-    AliVParticle* part = GetMCTrack(imother);
+    AliVParticle* part = AliAnalysisMuonUtility::GetMCTrack(imother,InputEvent(),MCEvent());
     //if ( ! part ) return motherType;
     
     Int_t absPdg = TMath::Abs(part->PdgCode());
     
-    Bool_t isPrimary = ( fMCEvent ) ? ( imother < fMCEvent->GetNumberOfPrimaries() ) : ((AliAODMCParticle*)part)->IsPrimary();
+    Bool_t isPrimary = AliAnalysisMuonUtility::IsPrimary(part, MCEvent());
     
     if ( isPrimary ) {
       if ( absPdg == 24 ) return kWbosonMu;
@@ -580,7 +417,7 @@ Int_t AliVAnalysisMuon::RecoTrackMother(AliVParticle* mcParticle)
       }
     } // is secondary
     
-    imother = GetMotherIndex(part);
+    imother = AliAnalysisMuonUtility::GetMotherIndex(part);
     
   } // loop on mothers
   
@@ -589,18 +426,6 @@ Int_t AliVAnalysisMuon::RecoTrackMother(AliVParticle* mcParticle)
 
 
 //________________________________________________________________________
-AliVVertex* AliVAnalysisMuon::GetVertexSPD() const
-{
-  //
-  /// Get vertex SPD
-  //
-  
-  AliVVertex* primaryVertex = ( fAODEvent ) ? (AliVVertex*)fAODEvent->GetPrimaryVertexSPD() : (AliVVertex*)fESDEvent->GetPrimaryVertexSPD();
-  return primaryVertex;
-}
-
-
-//________________________________________________________________________
 Bool_t AliVAnalysisMuon::AddObjectToCollection(TObject* object, Int_t index)
 {
   //
@@ -648,14 +473,14 @@ TObject* AliVAnalysisMuon::GetSum(TString physSel, TString trigClassNames, TStri
   /// - physSel, trigClassNames must be in the form: key1,key2
   /// - centrality must be in the form minValue_maxValue
   /// - objectPattern must be in the form match1,match2
-  ///   meaning that the object name must contain match1 and match2 (and either one of match3 and match4
+  ///   meaning that the object name must contain match1 or match2
   ///   wildcard * is allowed
   
   if ( ! fMergeableCollection ) return 0x0;
   
   // Get centrality range
   Int_t firstCentrality = 1;
-  Int_t lastCentrality = fCentralityClasses->GetNbins();
+  Int_t lastCentrality = GetCentralityClasses()->GetNbins();
   
   TObjArray* centralityRange = centrality.Tokenize("_");
   Float_t range[2] = {0., 100.};
@@ -663,15 +488,15 @@ TObject* AliVAnalysisMuon::GetSum(TString physSel, TString trigClassNames, TStri
     for ( Int_t irange=0; irange<2; ++irange ) {
       range[irange] = ((TObjString*)centralityRange->At(irange))->GetString().Atof();
     }
-    firstCentrality = fCentralityClasses->FindBin(range[0]+0.0001);
-    lastCentrality = fCentralityClasses->FindBin(range[1]-0.0001);
+    firstCentrality = GetCentralityClasses()->FindBin(range[0]+0.0001);
+    lastCentrality = GetCentralityClasses()->FindBin(range[1]-0.0001);
   }
   delete centralityRange;
   
   TString sumCentralityString = "";
   for ( Int_t icent=firstCentrality; icent<=lastCentrality; ++icent ) {
     if ( ! sumCentralityString.IsNull() ) sumCentralityString += ",";
-    sumCentralityString += fCentralityClasses->GetBinLabel(icent);
+    sumCentralityString += GetCentralityClasses()->GetBinLabel(icent);
   }
   
 //  objectPattern.ReplaceAll(" ","");
@@ -822,84 +647,30 @@ Bool_t AliVAnalysisMuon::SetSparseRange(AliCFGridSparse* gridSparse,
 }
 
 //________________________________________________________________________
-void AliVAnalysisMuon::SetTrigClassPatterns(TString pattern)
+TString AliVAnalysisMuon::GetDefaultTrigClassPatterns() const
 {
-  /// Set trigger classes
-  ///
-  /// Classes are filled dynamically according to the pattern
-  /// - if name contains ! (without spaces): reject it
-  /// - otherwise, keep it
-  /// example:
-  /// SetTrigClassPatterns("CMBAC !ALLNOTRD")
-  /// keeps classes containing CMBAC, and not containing ALLNOTRD.
-  ///
-  /// CAVEAT: if you use an fCFContainer and you want an axis to contain the trigger classes,
-  ///         please be sure that each pattern matches only 1 trigger class, or triggers will be mixed up
-  ///         when merging different chuncks.
-
-  fSelectedTrigPattern->SetOwner();
-  if ( fSelectedTrigPattern->GetEntries() > 0 ) fSelectedTrigPattern->Delete();
-  fRejectedTrigPattern->SetOwner();
-  if ( fRejectedTrigPattern->GetEntries() > 0 ) fRejectedTrigPattern->Delete();
-
-  pattern.ReplaceAll("  "," ");
-  pattern.ReplaceAll("! ","!");
-
-  TObjArray* fullList = pattern.Tokenize(" ");
-
-  for ( Int_t ipat=0; ipat<fullList->GetEntries(); ++ipat ) {
-    TString currPattern = fullList->At(ipat)->GetName();
-    if ( currPattern.Contains("!") ) {
-      currPattern.ReplaceAll("!","");
-      fRejectedTrigPattern->AddLast(new TObjString(currPattern));
-    }
-    else fSelectedTrigPattern->AddLast(new TObjString(currPattern));
-  }
-  
-  delete fullList;
+  /// Get default trigger class patterns
+  return fMuonEventCuts->GetDefaultTrigClassPatterns();
 }
 
 //________________________________________________________________________
-void AliVAnalysisMuon::SetTrigClassLevels(TString pattern)
+void AliVAnalysisMuon::SetTrigClassPatterns(const TString pattern)
 {
-  /// Set trigger cut level associated to the trigger class
-  ///
-  /// example:
-  /// SetTrigClassLevels("MSL:Lpt MSH:Hpt")
-  ///
-  /// For the trigger classes defined in SetTrigClassPatterns
-  /// it check if they contains the keywords MSL or MSH
-  /// Hence, in the analysis, the function
-  /// TrackPtCutMatchTrigClass(track, "CPBIMSL") returns true if track match Lpt
-  /// TrackPtCutMatchTrigClass(track, "CPBIMSH") returns true if track match Hpt
-  /// TrackPtCutMatchTrigClass(track, "CMBAC") always returns true
-  
-  fSelectedTrigLevel->SetOwner();
-  if ( fSelectedTrigLevel->GetEntries() > 0 ) fSelectedTrigLevel->Delete();
-  
-  pattern.ReplaceAll("  "," ");
-  pattern.ReplaceAll(" :",":");
-  
-  TObjArray* fullList = pattern.Tokenize(" ");
-  
-  for ( Int_t ipat=0; ipat<fullList->GetEntries(); ++ipat ) {
-    TString currPattern = fullList->At(ipat)->GetName();
-    TObjArray* arr = currPattern.Tokenize(":");
-    TObjString* trigClassPattern = new TObjString(arr->At(0)->GetName());
-    TString selTrigLevel = arr->At(1)->GetName();
-    selTrigLevel.ToUpper();
-    UInt_t trigLevel = 0;
-    if ( selTrigLevel.Contains("APT") ) trigLevel = 1;
-    else if ( selTrigLevel.Contains("LPT") ) trigLevel = 2;
-    else if ( selTrigLevel.Contains("HPT") ) trigLevel = 3;
-    trigClassPattern->SetUniqueID(trigLevel);
-    fSelectedTrigLevel->AddLast(trigClassPattern);
-    delete arr;
-  }
-  
-  delete fullList;
+  /// Set trigger classes
+  TString currPattern = pattern;
+  if ( currPattern.IsNull() ) { 
+    currPattern = GetDefaultTrigClassPatterns();
+    currPattern.Append(" !CMUP"); // by default do not account for UltraPeripheral events
+  } 
+  fMuonEventCuts->SetTrigClassPatterns(currPattern);
 }
 
+//________________________________________________________________________
+TList* AliVAnalysisMuon::GetAllSelectedTrigClasses() const
+{
+  /// Get trigger classes
+  return fMuonEventCuts->GetAllSelectedTrigClasses();
+}
 
 //________________________________________________________________________
 void AliVAnalysisMuon::SetCentralityClasses(Int_t nCentralityBins, Double_t* centralityBins)
@@ -907,22 +678,16 @@ void AliVAnalysisMuon::SetCentralityClasses(Int_t nCentralityBins, Double_t* cen
   //
   /// Set centrality classes
   //
-  Double_t* bins = centralityBins;
-  Int_t nbins = nCentralityBins;
-  
-  Double_t defaultCentralityBins[] = {-5., 0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 100., 105.};
-  if ( ! centralityBins ) {
-    bins = defaultCentralityBins;
-    nbins = sizeof(defaultCentralityBins)/sizeof(defaultCentralityBins[0])-1;
-  }
+  fMuonEventCuts->SetCentralityClasses(nCentralityBins, centralityBins);
+}
 
-  if ( fCentralityClasses ) delete fCentralityClasses;
-  fCentralityClasses = new TAxis(nbins, bins);
-  TString currClass = "";
-  for ( Int_t ibin=1; ibin<=fCentralityClasses->GetNbins(); ++ibin ){
-    currClass = Form("%.0f_%.0f",fCentralityClasses->GetBinLowEdge(ibin),fCentralityClasses->GetBinUpEdge(ibin));
-    fCentralityClasses->SetBinLabel(ibin, currClass.Data());
-  }
+//________________________________________________________________________
+TAxis* AliVAnalysisMuon::GetCentralityClasses() const
+{
+  //
+  /// Set centrality classes
+  //
+  return fMuonEventCuts->GetCentralityClasses();
 }
 
 //________________________________________________________________________
@@ -994,95 +759,3 @@ void AliVAnalysisMuon::InitKeys()
   TString physSelKeys = "PhysSelPass PhysSelReject";
   fPhysSelKeys = physSelKeys.Tokenize(" ");
 }
-
-//________________________________________________________________________
-TObjArray* AliVAnalysisMuon::BuildTriggerClasses(TString firedTrigClasses)
-{
-  //
-  /// Return the list of trigger classes to be considered
-  /// for current event. Update the global list if necessary
-  //
-
-  TObjArray* selectedTrigClasses = new TObjArray(0);
-  selectedTrigClasses->SetOwner();
-  
-  TObjArray* firedTrigClassesList = firedTrigClasses.Tokenize(" ");
-  
-  UInt_t trigLevel = 0;
-  for ( Int_t itrig=0; itrig<firedTrigClassesList->GetEntries(); ++itrig ) {
-    TString trigName = ((TObjString*)firedTrigClassesList->At(itrig))->GetString();
-    
-    TObject* foundTrig = fTriggerClasses->FindObject(trigName.Data());
-    if ( foundTrig ) trigLevel = foundTrig->GetUniqueID();
-    else {
-      Bool_t rejectTrig = kFALSE;
-      for ( Int_t ipat=0; ipat<fRejectedTrigPattern->GetEntries(); ++ipat ) {
-        if ( trigName.Contains(fRejectedTrigPattern->At(ipat)->GetName() ) ) {
-          rejectTrig = kTRUE;
-          break;
-        }
-      } // loop on reject pattern
-      if ( rejectTrig ) continue;
-
-      rejectTrig = kTRUE;
-      for ( Int_t ipat=0; ipat<fSelectedTrigPattern->GetEntries(); ++ipat ) {
-        if ( trigName.Contains(fSelectedTrigPattern->At(ipat)->GetName() ) ) {
-          rejectTrig = kFALSE;
-          break;
-        }
-      } // loop on keep pattern
-      if ( rejectTrig ) continue;
-    
-      trigLevel = 0;
-      for ( Int_t ipat=0; ipat<fSelectedTrigLevel->GetEntries(); ++ipat ) {
-        if ( trigName.Contains(fSelectedTrigLevel->At(ipat)->GetName() ) ) {
-          trigLevel = fSelectedTrigLevel->At(ipat)->GetUniqueID();
-          break;
-        }
-      } // loop on trig level patterns      
-    }
-    TObjString* currTrig = new TObjString(trigName);
-    currTrig->SetUniqueID(trigLevel);
-    selectedTrigClasses->AddLast(currTrig);
-    
-    if ( foundTrig ) continue;
-    AliInfo(Form("Adding %s (trig level %u) to considered trigger classes",trigName.Data(),trigLevel));
-    TObjString* addTrig = new TObjString(trigName);
-    addTrig->SetUniqueID(trigLevel);
-    fTriggerClasses->Add(addTrig);
-  } // loop on trigger classes
-
-  delete firedTrigClassesList;
-
-  return selectedTrigClasses;
-}
-
-
-//________________________________________________________________________
-Bool_t AliVAnalysisMuon::TrackPtCutMatchTrigClass(AliVParticle* track, TString trigClassName)
-{
-  /// Check if track passes the trigger pt cut level used in the trigger class
-  Int_t matchTrig = ( fAODEvent ) ? ((AliAODTrack*)track)->GetMatchTrigger() : ((AliESDMuonTrack*)track)->GetMatchTrigger();
-  Int_t classMatchLevel = GetTrigClassPtCutLevel(trigClassName);
-  Bool_t matchTrackerPt = kTRUE;
-  if ( fMuonTrackCuts && fMuonTrackCuts->IsApplySharpPtCutInMatching() ) {
-    matchTrackerPt = ( track->Pt() >= fMuonTrackCuts->GetSharpPtCut(classMatchLevel-1,kFALSE) );
-  }
-  Bool_t passCut = ( ( matchTrig >= classMatchLevel ) && matchTrackerPt );
-  AliDebug(1,Form("Class %s  matchTrig %i  trackMatchTrig %i  trackPt %g (required %i)  passCut %i", trigClassName.Data(), classMatchLevel, matchTrig, track->Pt(), fMuonTrackCuts->IsApplySharpPtCutInMatching(),passCut));
-  return passCut;
-}
-
-
-//________________________________________________________________________
-Int_t AliVAnalysisMuon::GetTrigClassPtCutLevel(TString trigClassName)
-{
-  /// Get trigger class pt cut level for tracking/trigger matching
-  TObject* obj = fTriggerClasses->FindObject(trigClassName.Data());
-  if ( ! obj ) {
-    AliWarning(Form("Class %s not in the list!", trigClassName.Data()));
-    return -1;
-  }
-  
-  return obj->GetUniqueID();
-}
index cd24b4c..ab9bd4b 100644 (file)
@@ -22,6 +22,7 @@ class AliVParticle;
 class AliAODEvent;
 class AliESDEvent;
 class AliCFGridSparse;
+class AliMuonEventCuts;
 class AliMuonTrackCuts;
 class AliMuonPairCuts;
 class AliVVertex;
@@ -42,12 +43,18 @@ class AliVAnalysisMuon : public AliAnalysisTaskSE {
   virtual void   FinishTaskOutput();
 
   void SetCentralityClasses(Int_t nCentralityBins = -1, Double_t* centralityBins = 0x0);
+  TAxis* GetCentralityClasses() const;
   Bool_t SetCentralityClassesFromOutput();
 
-  void SetTrigClassPatterns(TString pattern = "CINT CMU !CMUP CMBAC CPBI !-ACE- !-AC- !-E- !WU !EGA !EJE !PHS");
-  void SetTrigClassLevels(TString pattern = "MSL:Lpt MSH:Hpt MUL:Lpt MLL:Lpt");
+  void SetTrigClassPatterns(const TString pattern);
+  TString GetDefaultTrigClassPatterns() const;
+  /// Get trigger classes
+  TList* GetAllSelectedTrigClasses() const;
+  
   void SetTerminateOptions(TString physSel="All", TString trigClass="ANY", TString centralityRange="", TString furtherOpts="");
   
+  /// Get muon event cuts
+  AliMuonEventCuts* GetMuonEventCuts() { return fMuonEventCuts; }
   /// Get muon track cuts
   AliMuonTrackCuts* GetMuonTrackCuts() { return fMuonTrackCuts; }
   /// Get muon pair cuts
@@ -58,9 +65,6 @@ class AliVAnalysisMuon : public AliAnalysisTaskSE {
                                Int_t ivar, TString labelName,
                                Double_t varMin, Double_t varMax,
                                TString option = "");
-  
-  /// Set minimum number of vertex contributors
-  void SetMinNvtxContributors(Int_t minNvtxContributors) { fMinNvtxContirbutors = minNvtxContributors; }
 
  protected:
   
@@ -83,18 +87,8 @@ class AliVAnalysisMuon : public AliAnalysisTaskSE {
   /////////////////////
   // Utility methods //
   /////////////////////
-  
-  // Transparently handle tracks in ESD/AOD
-  Int_t GetNTracks();
-  AliVParticle* GetTrack(Int_t itrack);
-  TLorentzVector GetTrackPair(AliVParticle* track1, AliVParticle* track2) const;
-  
+    
   // Methods for MC
-  Bool_t IsMC();
-  Int_t GetNMCTracks();
-  AliVParticle* GetMCTrack(Int_t trackLabel);
-  Int_t GetMotherIndex(AliVParticle* mcParticle);
-  Int_t GetDaughterIndex(AliVParticle* mcParticle, Int_t idaughter);
   Int_t GetParticleType(AliVParticle* track);
   Int_t RecoTrackMother(AliVParticle* mcParticle);
   
@@ -103,16 +97,6 @@ class AliVAnalysisMuon : public AliAnalysisTaskSE {
   TObject* GetMergeableObject(TString physSel, TString trigClassName, TString centrality, TString objectName);
   TObject* GetSum(TString physSel, TString trigClassNames, TString centrality, TString objectPattern);
   
-  // A useful constant
-  Double_t MuonMass2() const;
-  
-  // Handle triggers
-  Bool_t TrackPtCutMatchTrigClass(AliVParticle* track, TString trigClassName);
-  Int_t GetTrigClassPtCutLevel(TString trigClassName);
-  
-  // Methods for event information
-  AliVVertex* GetVertexSPD() const;
-  
   enum {
     kPhysSelPass,    ///< Physics selected events
     kPhysSelReject,  ///< Events non-passing selection
@@ -130,7 +114,8 @@ class AliVAnalysisMuon : public AliAnalysisTaskSE {
     kUnidentified,  ///< Particle that fails matching kine
     kNtrackSources  ///< Total number of track sources
   };
-    
+  
+  AliMuonEventCuts* fMuonEventCuts; ///< Muon event cuts
   AliMuonTrackCuts* fMuonTrackCuts; ///< Muon track cuts
   AliMuonPairCuts* fMuonPairCuts;   ///< Muon pair track cuts
   AliESDEvent* fESDEvent;      //!< ESD event, not owner
@@ -139,13 +124,10 @@ class AliVAnalysisMuon : public AliAnalysisTaskSE {
   TObjArray* fChargeKeys;      ///< Muon charge keys
   TObjArray* fSrcKeys;         ///< MC sources names
   TObjArray* fPhysSelKeys;     ///< Physics selection names
-  TList* fTriggerClasses;      ///< List of trigger classes
-  TAxis* fCentralityClasses;   ///< Centrality classes
   
   AliCounterCollection* fEventCounters;  //!< event counters
   AliMergeableCollection* fMergeableCollection; //!< collection of mergeable objects
   TObjArray* fOutputList;  //!< List of outputs  
-  Int_t fMinNvtxContirbutors;  ///< Minimum number of vertex contributors
 
  private:
   AliVAnalysisMuon(const AliVAnalysisMuon&);
@@ -153,14 +135,9 @@ class AliVAnalysisMuon : public AliAnalysisTaskSE {
   
   void InitKeys();
   void CreateMergeableObjects(TString physSel, TString trigClassName, TString centrality);
-  TObjArray* BuildTriggerClasses(TString firedTrigClasses);
-  
-  TObjArray* fSelectedTrigPattern; ///< List of triggers to be kept
-  TObjArray* fRejectedTrigPattern; ///< List of triggers to be rejected
-  TObjArray* fSelectedTrigLevel;   ///< Track-trigger pt cut for selected trigger class
   TObjArray* fOutputPrototypeList; //!< List of prototype object to be used in collection
 
-  ClassDef(AliVAnalysisMuon, 3);
+  ClassDef(AliVAnalysisMuon, 4);
 };
 
 #endif
index 2b3238f..e1feda5 100644 (file)
@@ -61,7 +61,9 @@
 
 // PWG includes
 #include "AliMergeableCollection.h"
+#include "AliMuonEventCuts.h"
 #include "AliMuonTrackCuts.h"
+#include "AliAnalysisMuonUtility.h"
 
 
 /// \cond CLASSIMP
@@ -193,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);
+    track = AliAnalysisMuonUtility::GetTrack(itrack, InputEvent());
     fMuonTrackCuts->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);
@@ -214,7 +213,7 @@ 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 chi2 = pDca / fMuonTrackCuts->GetSigmaPdca(AliAnalysisMuonUtility::GetRabs(track)) ;
     chi2 *= chi2;
     Double_t chiProb = TMath::Prob(chi2, 2);
 
@@ -336,7 +335,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;
 
@@ -765,13 +764,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,12 +794,12 @@ 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[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("%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");
index b9e39ce..fd48711 100644 (file)
 
 // PWG3 includes
 #include "AliVAnalysisMuon.h"
+#include "AliMuonEventCuts.h"
 #include "AliMuonTrackCuts.h"
 #include "AliMergeableCollection.h"
+#include "AliAnalysisMuonUtility.h"
 
 
 /// \cond CLASSIMP
@@ -191,8 +193,8 @@ void AliAnalysisTaskTrigChEff::FinishTaskOutput()
 
   TString histoName = "";
   for ( Int_t isel=0; isel<kNselections; ++isel ) {
-    for ( Int_t itrig=0; itrig<fTriggerClasses->GetEntries(); ++itrig ) {
-      for ( Int_t icent=1; icent<=fCentralityClasses->GetNbins(); ++icent ) {
+    for ( Int_t itrig=0; itrig<GetAllSelectedTrigClasses()->GetEntries(); ++itrig ) {
+      for ( Int_t icent=1; icent<=GetCentralityClasses()->GetNbins(); ++icent ) {
         for ( Int_t imethod=0; imethod<kNeffMethods; ++imethod ) {
           for ( Int_t itype=0; itype<kNhistoTypes; ++itype ) {
             for ( Int_t icount=-1; icount<kNcounts; ++icount ) {
@@ -201,10 +203,10 @@ void AliAnalysisTaskTrigChEff::FinishTaskOutput()
                   TH1* histo = 0x0;
                   for ( Int_t jmatch=imatch+1; jmatch<=kMatchHpt; ++jmatch ) {
                     histoName = GetHistoName(itype, icount, ich, jmatch, imethod);
-                    TH1* histoAdd = (TH1*)fMergeableCollection->GetObject(Form("/%s/%s/%s/",fPhysSelKeys->At(isel)->GetName(), fTriggerClasses->At(itrig)->GetName(), fCentralityClasses->GetBinLabel(icent)), histoName);
+                    TH1* histoAdd = (TH1*)fMergeableCollection->GetObject(Form("/%s/%s/%s/",fPhysSelKeys->At(isel)->GetName(), GetAllSelectedTrigClasses()->At(itrig)->GetName(), GetCentralityClasses()->GetBinLabel(icent)), histoName);
                     if ( ! histoAdd ) continue;
                     histoName = GetHistoName(itype, icount, ich, imatch, imethod);
-                    if ( ! histo ) histo = (TH1*)GetMergeableObject(fPhysSelKeys->At(isel)->GetName(), fTriggerClasses->At(itrig)->GetName(), fCentralityClasses->GetBinLabel(icent), histoName);
+                    if ( ! histo ) histo = (TH1*)GetMergeableObject(fPhysSelKeys->At(isel)->GetName(), GetAllSelectedTrigClasses()->At(itrig)->GetName(), GetCentralityClasses()->GetBinLabel(icent), histoName);
                     AliDebug(2,Form("Adding %s (%g) to %s (%g)", histoAdd->GetName(), histoAdd->Integral(), histo->GetName(), histo->Integral()));
                     histo->Add(histoAdd);
                   } // loop on higher pt matching
@@ -320,7 +322,7 @@ void AliAnalysisTaskTrigChEff::MyUserCreateOutputObjects()
     } // loop on track selection
   } // loop on eff method
 
-  fMuonTrackCuts->Print("mask");
+  fMuonTrackCuts->Print();
   
   fList = new TList();
   fList->SetOwner();
@@ -358,14 +360,14 @@ void AliAnalysisTaskTrigChEff::ProcessEvent(TString physSel, const TObjArray& se
 
   Int_t itrackSel = -1;
 
-  Int_t nTracks = GetNTracks();
+  Int_t nTracks = AliAnalysisMuonUtility::GetNTracks(InputEvent());
   for ( Int_t itrack = 0; itrack < nTracks; ++itrack ) {
-    track = GetTrack(itrack);
+    track = AliAnalysisMuonUtility::GetTrack(itrack,InputEvent());
 
-    Bool_t matchTracker = ( fAODEvent && ((AliAODTrack*)track)->IsMuonTrack() ) || ((AliESDMuonTrack*)track)->ContainTrackerData();
+    Bool_t matchTracker = AliAnalysisMuonUtility::IsMuonTrack(track);
     if ( ! matchTracker && ! fUseGhosts ) continue;
 
-    Int_t matchTrig = ( fAODEvent ) ? ((AliAODTrack*)track)->GetMatchTrigger() : ((AliESDMuonTrack*)track)->GetMatchTrigger();
+    Int_t matchTrig = AliAnalysisMuonUtility::GetMatchTrigger(track);
     itrackSel = matchTrig;
     UInt_t selection = fMuonTrackCuts->GetSelectionMask(track);
     Bool_t isSelected = ( ( selection & fMuonTrackCuts->GetFilterMask() ) == fMuonTrackCuts->GetFilterMask() );
@@ -377,12 +379,12 @@ void AliAnalysisTaskTrigChEff::ProcessEvent(TString physSel, const TObjArray& se
     for ( Int_t imethod=0; imethod<kNeffMethods; ++imethod ) {
       if ( imethod == kEffFromTrack ) {
         if ( track->P() < 10. ) continue;
-        pattern = ( fAODEvent ) ? ((AliAODTrack*)track)->GetMUONTrigHitsMapTrk() :  ((AliESDMuonTrack*)track)->GetHitsPatternInTrigChTrk();
+        pattern = AliAnalysisMuonUtility::GetMUONTrigHitsMapTrk(track);
         board = AliESDMuonTrack::GetCrossedBoard(pattern);
       }
       else {
-        pattern = ( fAODEvent ) ? ((AliAODTrack*)track)->GetMUONTrigHitsMapTrg() :  ((AliESDMuonTrack*)track)->GetHitsPatternInTrigCh();
-        board = ( fAODEvent ) ? AliESDMuonTrack::GetCrossedBoard(pattern) : ((AliESDMuonTrack*)track)->LoCircuit();
+        pattern = AliAnalysisMuonUtility::GetMUONTrigHitsMapTrg(track);
+        board = ( AliAnalysisMuonUtility::IsAODEvent(InputEvent()) ) ? AliESDMuonTrack::GetCrossedBoard(pattern) : ((AliESDMuonTrack*)track)->LoCircuit();
       }
       
       Int_t effFlag = AliESDMuonTrack::GetEffFlag(pattern);
index f6473e4..7fe38ed 100644 (file)
@@ -41,7 +41,12 @@ AliAnalysisTaskTrigChEff* AddTaskMTRchamberEfficiency(Bool_t useGhosts = kFALSE,
   AliAnalysisTaskTrigChEff* taskTrigChEff = new AliAnalysisTaskTrigChEff("TriggerChamberEfficiency", *muonTrackCuts);
   taskTrigChEff->SetUseGhostTracks(useGhosts);
   if ( isMC ) taskTrigChEff->SetTrigClassPatterns("ANY");
-  else taskTrigChEff->SetTrigClassPatterns("ANY CINT CMU !CMUP CMBAC CPBI !-ACE- !-AC- !-E- !WU !EGA !EJE");
+  else {
+    TString trigClassPatterns = taskTrigChEff->GetDefaultTrigClassPatterns();
+    trigClassPatterns.Prepend("ANY !CMUP ");
+    taskTrigChEff->SetTrigClassPatterns(trigClassPatterns);
+  }
+  taskTrigChEff->GetMuonEventCuts()->SetFilterMask(AliMuonEventCuts::kSelectedTrig);
   mgr->AddTask(taskTrigChEff);
 
   // Create container