]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG3/muon/AliAnalysisTaskTrigChEff.cxx
Updates in macro for 2d v2 analysis
[u/mrichter/AliRoot.git] / PWG3 / muon / AliAnalysisTaskTrigChEff.cxx
index 4862d3652234e1fd225e0edbad2d28a74fc7898e..37cf83ae51bb2254ef115985f24ff2d4f0c64518 100644 (file)
@@ -1,96 +1,98 @@
+/**************************************************************************
+ * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//-----------------------------------------------------------------------------
+/// \class AliAnalysisTaskSingleMu
+/// Analysis task for single muons in the spectrometer.
+/// The output is a list of histograms.
+/// The macro class can run on AOD or ESDs.
+/// If Monte Carlo information is present, some basics checks are performed.
+///
+/// \author Diego Stocco
+//-----------------------------------------------------------------------------
+
+//----------------------------------------------------------------------------
+//    Implementation of the class for trigger chamber efficiency determinaltion
+//----------------------------------------------------------------------------
+
+
 #define AliAnalysisTaskTrigChEff_cxx
 
 // ROOT includes
-#include "TChain.h"
 #include "TH1.h"
+#include "TH2.h"
 #include "TCanvas.h"
 #include "TROOT.h"
 #include "TString.h"
 #include "TList.h"
+#include "TGraphAsymmErrors.h"
 
 // STEER includes
 #include "AliLog.h"
-
 #include "AliESDEvent.h"
 #include "AliESDMuonTrack.h"
-#include "AliESDInputHandler.h"
-
-#include "AliAODEvent.h"
-#include "AliAODTrack.h"
-#include "AliAODInputHandler.h"
+#include "AliAnalysisManager.h"
 
 // ANALYSIS includes
-#include "AliAnalysisTask.h"
-#include "AliAnalysisDataSlot.h"
-#include "AliAnalysisManager.h"
+#include "AliAnalysisTaskSE.h"
+
 #include "AliAnalysisTaskTrigChEff.h"
 
 ClassImp(AliAnalysisTaskTrigChEff)
 
 //________________________________________________________________________
-AliAnalysisTaskTrigChEff::AliAnalysisTaskTrigChEff(const char *name) :
-  AliAnalysisTask(name,""), 
-  fESD(0),
-  fAOD(0),
-  fAnalysisType("ESD"),
-  fList(0)
+AliAnalysisTaskTrigChEff::AliAnalysisTaskTrigChEff() :
+AliAnalysisTaskSE(), 
+fUseGhosts(kFALSE),
+fList(0)
 {
   //
-  /// Constructor.
-  //
-  // Input slot #0 works with an Ntuple
-  DefineInput(0, TChain::Class());
-  // Output slot #0 writes into a TObjArray container
-  DefineOutput(0,  TList::Class());
+  /// Default constructor
 }
 
-//___________________________________________________________________________
-void AliAnalysisTaskTrigChEff::ConnectInputData(Option_t *) {
+
+//________________________________________________________________________
+AliAnalysisTaskTrigChEff::AliAnalysisTaskTrigChEff(const char *name) :
+  AliAnalysisTaskSE(name), 
+  fUseGhosts(kFALSE),
+  fList(0)
+{
   //
-  /// Connect ESD or AOD here
-  /// Called once
+  /// Constructor
   //
 
-  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
-  if (!tree) {
-    Printf("ERROR: Could not read chain from input slot 0");
-  } else {
-    // Disable all branches and enable only the needed ones
-    // The next two lines are different when data produced as AliESDEvent is read
-    if(fAnalysisType == "ESD") {
-      tree->SetBranchStatus("*", kFALSE);
-      tree->SetBranchStatus("MuonTracks.*", kTRUE);
-
-      AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
-      
-      if (!esdH) {
-       Printf("ERROR: Could not get ESDInputHandler");
-      } else
-       fESD = esdH->GetEvent();
-    }
-    else if(fAnalysisType == "AOD") {
-      AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
-      
-      if (!aodH) {
-       Printf("ERROR: Could not get AODInputHandler");
-      } else
-       fAOD = aodH->GetEvent();
-    }
-    else 
-      Printf("Wrong analysis type: Only ESD and AOD types are allowed!");
-  }
+  DefineOutput(1,  TList::Class());
+}
+
+//________________________________________________________________________
+AliAnalysisTaskTrigChEff::~AliAnalysisTaskTrigChEff()
+{
+  if ( ! AliAnalysisManager::GetAnalysisManager()->IsProofMode() )
+    delete fList;
 }
 
-//___________________________________________________________________________
-void AliAnalysisTaskTrigChEff::CreateOutputObjects() {
+//_________________________________________________________________________
+void AliAnalysisTaskTrigChEff::UserCreateOutputObjects() {
   //
   /// Create histograms
   /// Called once
   //
-  Printf("   CreateOutputObjects of task %s\n", GetName());
 
-  TString cathCode[2] = {"bendPlane", "nonBendPlane"};
-  TString countTypeName[2] = {"CountInCh", "NonCountInCh"};
+  TString countTypeName[kNcounts] = {"bendPlane", "nonBendPlane","bothPlanes", "allTracks"};
 
   const Char_t* yAxisTitle = "counts";
 
@@ -109,152 +111,164 @@ void AliAnalysisTaskTrigChEff::CreateOutputObjects() {
   Float_t boardLow = 1-0.5, boardHigh = kNboards+1.-0.5;
   const Char_t* boardName = "board";
 
-  TString baseName, histoName;
+  TString baseName, histoName, histoTitle;
   fList = new TList();
+  fList->SetOwner();
 
   TH1F* histo;
+  TH2F* histo2D;
+
+  Int_t histoIndex = -1;
+
+  for(Int_t icount=0; icount<kNcounts; icount++){
+    histoName = Form("%sCountChamber", countTypeName[icount].Data());
+    histo = new TH1F(histoName, histoName,
+                    chamberBins, chamberLow, chamberHigh);
+    histo->GetXaxis()->SetTitle(chamberName);
+    histo->GetYaxis()->SetTitle(yAxisTitle);
+    histoIndex = GetHistoIndex(kHchamberEff, icount);
+    fList->AddAt(histo, histoIndex);
+  } // loop on counts
 
-  histo = new TH1F("nTracksInSlat", "Num. of tracks used for efficiency calculation", 
-                  slatBins, slatLow, slatHigh);
-  histo->GetXaxis()->SetTitle(slatName);
-  histo->GetYaxis()->SetTitle("num of used tracks");
-
-  fList->AddAt(histo, kHtracksInSlat);
-
-  histo = new TH1F("nTracksInBoard", "Num. of tracks used for efficiency calculation", 
-                  boardBins, boardLow, boardHigh);
-  histo->GetXaxis()->SetTitle(boardName);
-  histo->GetYaxis()->SetTitle("num of used tracks");
-
-  fList->AddAt(histo, kHtracksInBoard);
-
-  for(Int_t hType=0; hType<kNcounts; hType++){
-    Int_t hindex = (hType==0) ? kHchamberAllEff : kHchamberNonEff;
-    for(Int_t cath=0; cath<kNcathodes; cath++){
-      histoName = Form("%sChamber%s", cathCode[cath].Data(), countTypeName[hType].Data());
+  for(Int_t icount=0; icount<kNcounts; icount++){
+    for(Int_t ch=0; ch<kNchambers; ch++){
+      histoName = Form("%sCountSlatCh%i", countTypeName[icount].Data(), kFirstTrigCh+ch);
       histo = new TH1F(histoName, histoName,
-                      chamberBins, chamberLow, chamberHigh);
-      histo->GetXaxis()->SetTitle(chamberName);
+                      slatBins, slatLow, slatHigh);
+      histo->GetXaxis()->SetTitle(slatName);
       histo->GetYaxis()->SetTitle(yAxisTitle);
-       
-      fList->AddAt(histo, hindex + cath);
-    } // loop on cath
+      histoIndex = GetHistoIndex(kHslatEff, icount, ch);
+      fList->AddAt(histo, histoIndex);
+    } // loop on chamber
   } // loop on counts
 
-  for(Int_t hType=0; hType<kNcounts; hType++){
-    Int_t hindex = (hType==0) ? kHslatAllEff : kHslatNonEff;
-    for(Int_t cath=0; cath<kNcathodes; cath++){
-      for(Int_t ch=0; ch<kNchambers; ch++){
-       Int_t chCath = GetPlane(cath, ch);
-       histoName = Form("%sSlat%s%i", cathCode[cath].Data(), countTypeName[hType].Data(), kFirstTrigCh+ch);
-       histo = new TH1F(histoName, histoName,
-                        slatBins, slatLow, slatHigh);
-       histo->GetXaxis()->SetTitle(slatName);
-       histo->GetYaxis()->SetTitle(yAxisTitle);
-
-       fList->AddAt(histo, hindex + chCath);
-      } // loop on chamber
-    } // loop on cath
+  for(Int_t icount=0; icount<kNcounts; icount++){
+    for(Int_t ch=0; ch<kNchambers; ch++){
+      histoName = Form("%sCountBoardCh%i", countTypeName[icount].Data(), kFirstTrigCh+ch);
+      histo = new TH1F(histoName, histoName,
+                      boardBins, boardLow, boardHigh);
+      histo->GetXaxis()->SetTitle(boardName);
+      histo->GetYaxis()->SetTitle(yAxisTitle);
+      histoIndex = GetHistoIndex(kHboardEff, icount, ch);
+      fList->AddAt(histo, histoIndex);
+    } // loop on chamber
   } // loop on counts
 
-  for(Int_t hType=0; hType<kNcounts; hType++){
-    Int_t hindex = (hType==0) ? kHboardAllEff : kHboardNonEff;
-    for(Int_t cath=0; cath<kNcathodes; cath++){
-      for(Int_t ch=0; ch<kNchambers; ch++){
-       Int_t chCath = GetPlane(cath, ch);
-       histoName = Form("%sBoard%s%i", cathCode[cath].Data(), countTypeName[hType].Data(), kFirstTrigCh+ch);
-       histo = new TH1F(histoName, histoName,
-                        boardBins, boardLow, boardHigh);
-       histo->GetXaxis()->SetTitle(boardName);
-       histo->GetYaxis()->SetTitle(yAxisTitle);
-
-       fList->AddAt(histo, hindex + chCath);
-      } // loop on chamber
-    } // loop on cath
-  } // loop on counts
+  histo2D = new TH2F("checkRejectedBoard", "Rejected tracks motivation", 
+                    4, 20.5, 24.5, boardBins, boardLow, boardHigh);
+  histo2D->GetXaxis()->SetBinLabel(1,"Many pads");
+  histo2D->GetXaxis()->SetBinLabel(2,"Few pads");
+  histo2D->GetXaxis()->SetBinLabel(3,"Outside geom");
+  histo2D->GetXaxis()->SetBinLabel(4,"Tracker track");
+  histo2D->GetYaxis()->SetTitle(boardName);
+  histoIndex = GetHistoIndex(kHcheckBoard);
+  fList->AddAt(histo2D, histoIndex);
+  
+  PostData(1, fList);
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskTrigChEff::Exec(Option_t *) {
+void AliAnalysisTaskTrigChEff::UserExec(Option_t *) {
   //
   /// Main loop
   /// Called for each event
   //
-  Int_t nTracks = 0, board = 0;
+  AliESDEvent* esdEvent = dynamic_cast<AliESDEvent*> (InputEvent());
+
+  if (!esdEvent) {
+    Printf("ERROR: esdEvent not available\n");
+    return;
+  }
+
+  Int_t slat = 0, board = 0;
   UShort_t pattern = 0;
   AliESDMuonTrack *esdTrack = 0x0;
-  AliAODTrack* aodTrack = 0x0;
 
-  if(fAnalysisType == "ESD") {
-    if (!fESD) {
-      Printf("ERROR: fESD not available");
-      return;
-    }
-    nTracks = fESD->GetNumberOfMuonTracks(); 
-  }
-  else if(fAnalysisType == "AOD") {
-    if (!fAOD) {
-      Printf("ERROR: fAOD not available");
-      return;
-    }
-    nTracks = fAOD->GetNumberOfTracks();
-  }
+  Int_t nTracks = esdEvent->GetNumberOfMuonTracks();
 
-  // Object declaration
   const Int_t kFirstTrigCh = 11; //AliMpConstants::NofTrackingChambers()+1;
 
+  TArrayI othersEfficient(kNchambers);
+  Int_t histoIndex = -1;
+
   for (Int_t itrack = 0; itrack < nTracks; itrack++) {
-    if(fAnalysisType == "ESD") {
-      esdTrack = fESD->GetMuonTrack(itrack);
-      pattern =  esdTrack->GetHitsPatternInTrigCh();
-      board = esdTrack->LoCircuit();
-    }
-    else if(fAnalysisType == "AOD") {
-      aodTrack = fAOD->GetTrack(itrack);
-      if(!aodTrack->IsMuonTrack()) continue;
-      pattern =  aodTrack->GetHitsPatternInTrigCh();
-      board = 0; // aodTrack->LoCircuit(); Lo Circuit not implemented in AOD
-    }
+    esdTrack = esdEvent->GetMuonTrack(itrack);
 
-    Int_t effFlag = GetEffFlag(pattern);
+    if ( ! esdTrack->ContainTrackerData() && ! fUseGhosts ) continue;
 
-    if(effFlag < kChEff) continue; // Track not good for efficiency calculation
+    pattern =  esdTrack->GetHitsPatternInTrigCh();
+    Int_t effFlag = AliESDMuonTrack::GetEffFlag(pattern);
 
-    Int_t slat = GetSlat(pattern);
+    board = esdTrack->LoCircuit();
 
-    if(effFlag >= kSlatEff) ((TH1F*)fList->At(kHtracksInSlat))->Fill(slat);
-    if(effFlag >= kBoardEff) ((TH1F*)fList->At(kHtracksInBoard))->Fill(board);
+    if(effFlag < AliESDMuonTrack::kChEff) {
+      histoIndex = GetHistoIndex(kHcheckBoard);
+      ((TH2F*)fList->At(histoIndex))->Fill(AliESDMuonTrack::GetSlatOrInfo(pattern), board);
+      continue; // Track not good for efficiency calculation
+    }
 
+    othersEfficient.Reset(1);
     for(Int_t cath=0; cath<kNcathodes; cath++){
-      Int_t ineffCh = IsChInefficient(pattern, cath);
-      Int_t nChambers = kNchambers;
-      for(Int_t ch=0; ch<nChambers; ch++){
-       Int_t whichType = kAllChEff;
-       Int_t currCh = ch;
-       if(ineffCh>=0){
-         whichType = kChNonEff;
-         currCh = ineffCh;
-         nChambers = -1;
-       }
-
-       Int_t iChamber = kFirstTrigCh + currCh;
-       Int_t hindex = (whichType==kAllChEff) ? kHchamberAllEff : kHchamberNonEff;
-       ((TH1F*)fList->At(hindex + cath))->Fill(iChamber);
-
-       if(effFlag < kSlatEff) continue; // Track crossed different slats
-       Int_t chCath = GetPlane(cath, currCh);
-       hindex = (whichType==kAllChEff) ? kHslatAllEff : kHslatNonEff;
-       ((TH1F*)fList->At(hindex + chCath))->Fill(slat);
-
-       if(effFlag < kBoardEff) continue; // Track crossed different boards
-       hindex = (whichType==kAllChEff) ? kHboardAllEff : kHboardNonEff;
-       ((TH1F*)fList->At(hindex + chCath))->Fill(board);
+      for(Int_t ich=0; ich<kNchambers; ich++){
+       if( ! AliESDMuonTrack::IsChamberHit(pattern, cath, ich)){
+         for(Int_t jch=0; jch<kNchambers; jch++){
+           if ( jch != ich) {
+             othersEfficient[jch] = 0;
+           }
+         } // loop on other chambers
+         break;
+       } // if chamber not efficient
       } // loop on chambers
     } // loop on cathodes
-  }
+
+    Bool_t rejectTrack = kTRUE;
+    for (Int_t ich=0; ich<kNchambers; ich++){
+      if ( othersEfficient[ich] > 0 ){
+       rejectTrack = kFALSE;
+       break;
+      }
+    }
+
+    if ( rejectTrack ) continue;
+
+    slat = AliESDMuonTrack::GetSlatOrInfo(pattern);
+
+    for(Int_t ch=0; ch<kNchambers; ch++){
+      if ( ! othersEfficient[ch] )
+       continue; // Reject track if the info of the chamber under study 
+                  // is necessary to create the track itself
+
+      Int_t iChamber = kFirstTrigCh + ch;
+
+      Bool_t hitsBend = AliESDMuonTrack::IsChamberHit(pattern, 0, ch);
+      Bool_t hitsNonBend = AliESDMuonTrack::IsChamberHit(pattern, 1, ch);
+
+      Bool_t fillHisto[kNcounts] = {
+       hitsBend,
+       hitsNonBend,
+       ( hitsBend && hitsNonBend ),
+       kTRUE
+      };
+
+      for (Int_t icount=0; icount<kNcounts; icount++){
+       if ( ! fillHisto[icount] ) continue;
+
+       histoIndex = GetHistoIndex(kHchamberEff, icount);
+       ((TH1F*)fList->At(histoIndex))->Fill(iChamber);
+
+       if(effFlag < AliESDMuonTrack::kSlatEff) continue; // Track crossed different slats
+       histoIndex = GetHistoIndex(kHslatEff, icount, ch);
+       ((TH1F*)fList->At(histoIndex))->Fill(slat);
+
+       if(effFlag < AliESDMuonTrack::kBoardEff) continue; // Track crossed different boards
+       histoIndex = GetHistoIndex(kHboardEff, icount, ch);
+       ((TH1F*)fList->At(histoIndex))->Fill(board);
+      } // loop on chambers
+    } // loop on count types
+  } // loop on tracks
 
   // Post final data. It will be written to a file with option "RECREATE"
-  PostData(0, fList);
+  PostData(1, fList);
 }
 
 //________________________________________________________________________
@@ -263,42 +277,86 @@ void AliAnalysisTaskTrigChEff::Terminate(Option_t *) {
   /// Draw result to the screen
   /// Called once at the end of the query.
   //
-  if (!gROOT->IsBatch()) {
-    TCanvas *can[kNcathodes];
-    TH1F *num = 0x0;
-    TH1F *den = 0x0;
-    for(Int_t cath=0; cath<kNcathodes; cath++){
-      TString canName = Form("can%i",cath);
-      can[cath] = new TCanvas(canName.Data(),canName.Data(),10*(1+cath),10*(1+cath),310,310);
-      can[cath]->SetFillColor(10); can[cath]->SetHighLightColor(10);
-      can[cath]->SetLeftMargin(0.15); can[cath]->SetBottomMargin(0.15);  
-      can[cath]->Divide(2,2);
+  if ( gROOT->IsBatch() ) return;
+  fList = dynamic_cast<TList*> (GetOutputData(1));
+
+  if (!fList) return;
+
+  TCanvas *can;
+  TH1F *num = 0x0;
+  TH1F *den = 0x0;
+  TGraphAsymmErrors* effGraph = 0x0;
+  TString baseName[3] = {"Chamber", "RPC", "Board"};
+  Int_t baseIndex[3] = {kHchamberEff, kHslatEff, kHboardEff};
+  TString effName[kNcounts-1] = {"BendPlane", "NonBendPlane", "BothPlanes"};
+  Int_t histoIndexNum = -1, histoIndexDen = -1;
+  for (Int_t itype=0; itype<3; itype++) {
+    for(Int_t icount=0; icount<kNcounts-1; icount++){
+      TString canName = Form("efficiencyPer%s_%s",baseName[itype].Data(),effName[icount].Data());
+      can = new TCanvas(canName.Data(),canName.Data(),10*(1+kNcounts*itype+icount),10*(1+kNcounts*itype+icount),310,310);
+      can->SetFillColor(10); can->SetHighLightColor(10);
+      can->SetLeftMargin(0.15); can->SetBottomMargin(0.15);  
+      if ( itype > 0 )
+       can->Divide(2,2);
+
       for(Int_t ch=0; ch<kNchambers; ch++){
-       Int_t chCath = GetPlane(cath, ch);
-       num = (TH1F*)(fList->At(kHboardAllEff + chCath)->Clone());
-       den = (TH1F*)(fList->At(kHboardNonEff + chCath)->Clone());
-       den->Add(num);
-       num->Divide(den);
-       can[cath]->cd(ch+1);
-       num->DrawCopy("E");
-      }
-    }
-  }
+       histoIndexNum = GetHistoIndex(baseIndex[itype], icount, ch);
+       histoIndexDen = GetHistoIndex(baseIndex[itype], kAllTracks, ch);
+       num = (TH1F*)(fList->At(histoIndexNum));
+       den = (TH1F*)(fList->At(histoIndexDen));
+       effGraph = new TGraphAsymmErrors(num, den);
+       effGraph->GetYaxis()->SetRangeUser(0., 1.1);
+       effGraph->GetYaxis()->SetTitle("Efficiency");
+       effGraph->GetXaxis()->SetTitle(baseName[itype].Data());
+       can->cd(ch+1);
+       effGraph->Draw("AP");
+       if ( itype == 0 ) break;
+      } // loop on chamber
+    } // loop on count types
+  } // loop on histo
 }
 
 //________________________________________________________________________
-Int_t AliAnalysisTaskTrigChEff::IsChInefficient(UShort_t pattern, 
-                                               Int_t cathode)
+Int_t
+AliAnalysisTaskTrigChEff::GetHistoIndex(Int_t histoType, Int_t countType, 
+                                       Int_t chamber)
 {
   //
-  /// Check which chamber was inefficient.
+  /// Return the index of the histogram in the list
   //
-  Int_t ineffCh = -999;
-  for(Int_t ch=0; ch<kNchambers; ch++){
-    Int_t chCath = GetPlane(cathode, ch);
-    Int_t invert = kNplanes - chCath - 1;
-    Int_t response = (pattern >> invert) & 0x01;
-    if(!response) ineffCh = ch;
+  switch ( histoType ) {
+  case kHchamberEff:
+    return 0 + countType;
+  case kHslatEff:
+    return 4 + kNchambers*countType + chamber;
+  case kHboardEff:
+    return 20 + kNchambers*countType + chamber;
+  case kHcheckBoard:
+    return 36;
+  }
+
+  /*
+  const Int_t kNhistosPlaneCorr = 38;
+
+  switch ( histoType ){
+  case kHtracksInSlat:
+    return 0 + planeCorrelation*kNhistosPlaneCorr;
+  case kHtracksInBoard:
+    return 1 + planeCorrelation*kNhistosPlaneCorr;
+  case kHchamberEff:
+    return 2 + kNcathodes*countType + cathode
+      + planeCorrelation*kNhistosPlaneCorr;
+  case kHslatEff:
+    return 6 + kNchambers*kNcathodes*countType 
+      + kNchambers*cathode + chamber
+      + planeCorrelation*kNhistosPlaneCorr;
+  case kHboardEff:
+    return 22 + kNchambers*kNcathodes*countType 
+      + kNchambers*cathode + chamber
+      + planeCorrelation*kNhistosPlaneCorr;
+  case kHcheckBoard:
+    return 0 + 2*kNhistosPlaneCorr;
   }
-  return ineffCh;
+  */
+  return -1;
 }