Single Muon Analysis from generic or muon-selected AOD's (Diego)
authormartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 24 Nov 2008 17:30:49 +0000 (17:30 +0000)
committermartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 24 Nov 2008 17:30:49 +0000 (17:30 +0000)
PWG3/muon/AliAnalysisTaskSingleMu.cxx
PWG3/muon/AliAnalysisTaskSingleMu.h
PWG3/muon/RunSingleMuonAnalysisFromAOD.C

index 0262735..ca358bd 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+//-----------------------------------------------------------------------------
+/// \class AliAnalysisTaskSingleMu
+/// Analysis task for single muons in the spectrometer.
+/// The output is a tree with:
+///  - pt, y and phi of the muon
+///  - z position of primary vertex
+///  - transverse distance at vertex (DCA)
+///  - matched trigger
+///
+/// \author Diego Stocco
+//-----------------------------------------------------------------------------
+
 //----------------------------------------------------------------------------
 //    Implementation of the single muon analysis class
-// An example of usage can be found in the macro runSingleMuAnalysis.C.
+// An example of usage can be found in the macro RunSingleMuAnalysisFromAOD.C.
 //----------------------------------------------------------------------------
 
 #define AliAnalysisTaskSingleMu_cxx
 
 // ROOT includes
-#include "Riostream.h"
 #include "TChain.h"
-#include "TH1.h"
-#include "TCanvas.h"
-#include "TSystem.h"
 #include "TROOT.h"
-#include "TParticle.h"
 #include "TLorentzVector.h"
+#include "TCanvas.h"
 
 // STEER includes
 #include "AliLog.h"
 #include "AliAnalysisManager.h"
 #include "AliAnalysisTaskSingleMu.h"
 
-ClassImp(AliAnalysisTaskSingleMu)
+/// \cond CLASSIMP
+ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
+/// \endcond
 
 //________________________________________________________________________
 AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name) :
   AliAnalysisTask(name,""), 
   fAOD(0),
-  fOutputContainer(0)
+  fResults(0),
+  fVarFloat(0),
+  fVarInt(0),
+  fFloatVarName(0),
+  fIntVarName(0)
 {
   //
   /// Constructor.
   //
-  ResetHistos();
+  InitVariables();
   // Input slot #0 works with an Ntuple
   DefineInput(0, TChain::Class());
-  // Output slot #0 writes into a TObjArray container
-  DefineOutput(0,  TObjArray::Class());
+  // Output slot #0 writes into a TTree container
+  DefineOutput(0,  TTree::Class());
 }
 
 //___________________________________________________________________________
@@ -86,7 +100,7 @@ void AliAnalysisTaskSingleMu::ConnectInputData(Option_t *)
       Printf("ERROR: Could not get AODInputHandler");
     } else
       printf("   ConnectInputData of task %s\n", GetName());
-      fAOD = aodH->GetEvent();
+    fAOD = aodH->GetEvent();
   }
 }
 
@@ -98,34 +112,23 @@ void AliAnalysisTaskSingleMu::CreateOutputObjects()
   //
   printf("   CreateOutputObjects of task %s\n", GetName());
 
-  Int_t ptBins = 60;
-  Float_t ptLow = 0., ptHigh = 30.;
-  const Char_t *ptName = "P_{t} (GeV/c)";
-
-  Int_t vzBins = 40;
-  Float_t vzLow = -20., vzHigh = 20.;
-  const Char_t *vzName = "Vz (cm)";
-
-  TString baseName, histoName;
-  fOutputContainer = new TObjArray(fgkNhistos*fgkNTrigCuts);
-  fOutputContainer->SetName("SingleMuAnalysisContainer");
-  Int_t iHisto = 0;
-
-  for(Int_t iTrig=0; iTrig<fgkNTrigCuts; iTrig++){
-    
-    // 2D histos
-    if(!fVzVsPt[iTrig]){
-      baseName = "fVzVsPt";
-      histoName = baseName + trigName[iTrig];
-      fVzVsPt[iTrig] = new TH2F(histoName, histoName,
-                               ptBins, ptLow, ptHigh,
-                               vzBins, vzLow, vzHigh);
-      fVzVsPt[iTrig]->GetXaxis()->SetTitle(ptName);
-      fVzVsPt[iTrig]->GetYaxis()->SetTitle(vzName);
-
-      fOutputContainer->AddAt(fVzVsPt[iTrig], iHisto);
-      iHisto++;
-    }
+  // initialize tree
+  if(!fResults) fResults = new TTree("Results", "Single mu selection results");
+
+  TString baseName, suffixName;
+
+  suffixName="/F";
+  for(Int_t iVar=0; iVar<kNfloatVars; iVar++){
+    baseName = fFloatVarName[iVar];
+    if(iVar==0) baseName += suffixName;
+    fResults->Branch(fFloatVarName[iVar].Data(), &fVarFloat[iVar], baseName.Data());
+  }
+
+  suffixName="/I";
+  for(Int_t iVar=0; iVar<kNintVars; iVar++){
+    baseName = fIntVarName[iVar];
+    if(iVar==0) baseName += suffixName;
+    fResults->Branch(fIntVarName[iVar].Data(), &fVarInt[iVar], baseName.Data());
   }
 }
 
@@ -148,23 +151,18 @@ void AliAnalysisTaskSingleMu::Exec(Option_t *)
 
   // Object declaration
   AliAODTrack *muonTrack = 0x0;
-  TLorentzVector lorVec;
-  Int_t trigMatch = -1;
 
   Int_t nTracks = fAOD->GetNumberOfTracks();
   for (Int_t itrack = 0; itrack < nTracks; itrack++) {
     muonTrack = fAOD->GetTrack(itrack);
 
     // Apply cuts
-    if(!MuonPassesCuts(*muonTrack, lorVec, trigMatch)) continue;
-
-    for(Int_t iTrig=0; iTrig<=trigMatch; iTrig++){
-      fVzVsPt[iTrig]->Fill(lorVec.Pt(), fAOD->GetPrimaryVertex()->GetZ());
-    }
+    if(!FillTrackVariables(*muonTrack)) continue;
+    fResults->Fill();
   }
 
   // Post final data. It will be written to a file with option "RECREATE"
-  PostData(0, fOutputContainer);
+  PostData(0, fResults);
 }
 
 //________________________________________________________________________
@@ -176,51 +174,46 @@ void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
     TCanvas *c1 = new TCanvas("c1","Vz vs Pt",10,10,310,310);
     c1->SetFillColor(10); c1->SetHighLightColor(10);
     c1->SetLeftMargin(0.15); c1->SetBottomMargin(0.15);  
-    c1->Divide(2,2);
-    for(Int_t iTrig=0; iTrig<fgkNTrigCuts; iTrig++){
-      c1->cd(iTrig+1);
-      fVzVsPt[iTrig]->DrawCopy("COLZ");
-    }
+    fResults->Draw("pt:vz","","COLZ");
   }
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskSingleMu::ResetHistos() 
+void AliAnalysisTaskSingleMu::InitVariables() 
 {
   //
   /// Reset histograms
   //
-  for(Int_t iTrig=0; iTrig<fgkNTrigCuts; iTrig++){
-    fVzVsPt[iTrig] = 0x0;
-  }
-  trigName[kNoMatchTrig] = "NoMatchTrig";
-  trigName[kAllPtTrig]   = "AllPtTrig";
-  trigName[kLowPtTrig]   = "LowPtTrig";
-  trigName[kHighPtTrig]  = "HighPtTrig";
+
+  fVarFloat = new Float_t[kNfloatVars];
+  fVarInt = new Int_t[kNintVars];
+
+  fFloatVarName = new TString[kNfloatVars];
+  fFloatVarName[kVarPt]     = "pt";
+  fFloatVarName[kVarY]      = "y";
+  fFloatVarName[kVarPhi]    = "phi";
+  fFloatVarName[kVarVz]     = "vz";
+  fFloatVarName[kVarDCA]    = "dca";
+
+  fIntVarName = new TString[kNintVars];
+  fIntVarName[kVarTrig]     = "matchTrig";
 }
 
 
 //________________________________________________________________________
-Bool_t AliAnalysisTaskSingleMu::MuonPassesCuts(AliAODTrack &muonTrack,
-                                              TLorentzVector &lorVec,
-                                              Int_t &trigMatch)
+Bool_t AliAnalysisTaskSingleMu::FillTrackVariables(AliAODTrack &muonTrack)
 {
   //
   /// Fill lorentz vector and check cuts
   //
 
+  TLorentzVector lorVec;
+
   // Check if track is a muon
   if(muonTrack.GetMostProbablePID()!=AliAODTrack::kMuon) return kFALSE;
 
   // Check if track is triggered
-  trigMatch = kNoMatchTrig;
-  if (muonTrack.MatchTriggerHighPt()) {
-    trigMatch = kHighPtTrig;
-  } else if (muonTrack.MatchTriggerLowPt()) {
-    trigMatch = kLowPtTrig;
-  } else if (muonTrack.MatchTriggerAnyPt()){
-    trigMatch = kAllPtTrig;
-  }
+  fVarInt[kVarTrig] = (muonTrack.GetMatchTrigger() && 0x3);
   
   // Fill track parameters
   Double_t px = muonTrack.Px();
@@ -233,5 +226,16 @@ Bool_t AliAnalysisTaskSingleMu::MuonPassesCuts(AliAODTrack &muonTrack,
   Double_t energy = TMath::Sqrt(p*p + kMuonMass*kMuonMass);
   lorVec.SetPxPyPzE(px,py,pz,energy);
 
+  fVarFloat[kVarPt]  = lorVec.Pt();
+  fVarFloat[kVarY]   = lorVec.Rapidity();
+  fVarFloat[kVarPhi] = lorVec.Phi();
+
+  fVarFloat[kVarVz] = fAOD->GetPrimaryVertex()->GetZ();
+
+  Double_t xDca = muonTrack.XAtDCA();
+  Double_t yDca = muonTrack.YAtDCA();
+
+  fVarFloat[kVarDCA] = TMath::Sqrt(xDca*xDca + yDca*yDca);
+
   return kTRUE;
 }
index d1c66df..fd37d27 100644 (file)
@@ -1,8 +1,10 @@
-#include "TH2.h"
-#include "TObjArray.h"
+/// \ingroup "PWG3muon"
+/// \class AliAnalysisTaskSingleMu
+/// \brief Analysis task for single muons in the spectrometer
+///
+//  Author Diego Stocco
 
 #include "AliAODEvent.h"
-#include "AliAODVertex.h"
 #include "AliAODTrack.h"
 
 class AliAnalysisTaskSingleMu : public AliAnalysisTask {
@@ -15,37 +17,40 @@ class AliAnalysisTaskSingleMu : public AliAnalysisTask {
   virtual void   Exec(Option_t *option);
   virtual void   Terminate(Option_t *);
 
-protected:
-  Bool_t MuonPassesCuts(AliAODTrack &muonTrack,
-                       TLorentzVector &lorVec,
-                       Int_t &trigMatch);
-
-  const AliAODVertex* GetVertex();
-  void ResetHistos();
+ protected:
+  Bool_t FillTrackVariables(AliAODTrack &muonTrack);
+  
+  void InitVariables();
 
-private:
+ private:
   AliAnalysisTaskSingleMu(const AliAnalysisTaskSingleMu&);
   AliAnalysisTaskSingleMu& operator=(const AliAnalysisTaskSingleMu&);
 
   AliAODEvent *fAOD; //!< ESDevent object
 
-  static const Int_t fgkNhistos = 1;
-  static const Int_t fgkNTrigCuts = 4;
+  TTree *fResults; //!< Tree with results
 
-  enum
-  {
-    kNoMatchTrig,
-    kAllPtTrig,
-    kLowPtTrig,
-    kHighPtTrig
+  enum {
+    kVarPt,     //!< Muon pt
+    kVarY,      //!< Muon rapidity
+    kVarPhi,    //!< Muon phi
+    kVarVz,     //!< Primary vertex longitudinal position
+    kVarDCA,    //!< Transverse distance at vertex
+    kNfloatVars
   };
 
-  TString trigName[fgkNTrigCuts]; //!< trigger cut names 
+  enum {
+    kVarTrig,   //!< Matched trigger
+    kNintVars
+  };
+  
 
-  TObjArray * fOutputContainer; //!< output data container
+  Float_t* fVarFloat; //!< Array of float variables
+  Int_t* fVarInt;     //!< Array of int variables
+  
+  TString* fFloatVarName; //!< Float variable names for branches
+  TString* fIntVarName;   //!< Intt variable names for branches
 
-  TH2F *fVzVsPt[fgkNTrigCuts]; //!< Single muon spectrum
-     
   ClassDef(AliAnalysisTaskSingleMu, 0); // Single muon analysis
 };
 
index 64bb450..8a83ef0 100644 (file)
@@ -3,50 +3,77 @@
 // 
 // In case it is not run with full aliroot, it needs to have in the working directory:
 //  - STEERBase.par
+//  - ESD.par
 //  - AOD.par
 //  - ANALYSIS.par
+//  - ANALYSISalice.par
+//  - PWG3muon.par
+//
+// The inputPath is either:
+//  - The directory containing the AOD file in local mode
+//  - The xml file with the list AODs in the alien catalogue in grid mode
+//  - The proof dataset in proof mode
 // 
 // The macro reads AODs and outputs file:
-// - outputDir/singleMuAnalysis.root
+//  - outputDir/singleMuAnalysis.root
 //--------------------------------------------------------------------------
 
-void runSingleMuAnalysis(Char_t *inputDir=".", Char_t *outputDir=".") {
+enum analysisMode {kMlocal, kMgridInteractive, kMgridBatch, kMproof};
+TString modeName[4] = {"local", "local", "grid", "proof"};
+
+void RunSingleMuonAnalysisFromAOD(Int_t mode=kMlocal, Char_t *inputPath=".", Char_t *outputDir=".", Char_t *aodFilename = "AliAODs.root", Long64_t nRuns = -1, Long64_t offset = 0) {
   TStopwatch timer;
   timer.Start();
+
+  // Check if user is running root or aliroot
+  TString checkString = gSystem->Getenv("ALICE_ROOT");
+  checkString.Append("/lib/tgt_linux/libSTEERBase.so");
+  TString foundLib = gSystem->GetLibraries(checkString.Data());
+  Bool_t isFullAliroot = (foundLib.Length()==0) ? kFALSE : kTRUE;
+
+  // Load libraries
   gSystem->Load("libTree");
   gSystem->Load("libGeom");
   gSystem->Load("libVMC");
-  gSystem->Load("libANALYSIS");
-  gSystem->Load("libSTEERBase");
-  gSystem->Load("libAOD");
-  gSystem->Load("libESD");  
-  gSystem->Load("libPWG3base.so");
 
-  TString outFileName("singleMuAnalysis.root");
-  outFileName.Prepend(Form("%s/",outputDir));
+  if(mode==kMproof)
+    TProof::Open("alicecaf.cern.ch");
 
-  //____________________________________________//
-  AliTagAnalysis *TagAna = new AliTagAnalysis("AOD"); 
+  if(isFullAliroot){
+    gSystem->Load("libANALYSIS");
+    gSystem->Load("libANALYSISalice");
+    gSystem->Load("libSTEERBase");
+    gSystem->Load("libAOD");
+    gSystem->Load("libESD");  
+    gSystem->Load("${ALICE_ROOT}/lib/tgt_${ALICE_TARGET}/libPWG3muon.so");
+  }
+  else {
+    const Int_t nNeededPar = 6;
+    TString parList[nNeededPar] = {"STEERBase", "ESD", "AOD", "ANALYSIS", "ANALYSISalice", "PWG3muon"};
+    if(mode==kMproof){
+      gProof->UploadPackage("AF-v4-15");
+      gProof->EnablePackage("AF-v4-15");
+      if(!SetupPar("PWG3muon")) return;
+    }
+    else {
+      for(Int_t ipar=0; ipar<nNeededPar; ipar++){
+       if(!SetupPar(parList[ipar].Data())) return;
+      }
+    }
+  }
 
-  AliRunTagCuts *runCuts = new AliRunTagCuts();
-  AliLHCTagCuts *lhcCuts = new AliLHCTagCuts();
-  AliDetectorTagCuts *detCuts = new AliDetectorTagCuts();
-  AliEventTagCuts *evCuts = new AliEventTagCuts();
+  // Connect to alien
+  //
+  if(mode==kMgridInteractive || mode==kMgridBatch)
+    TGrid::Connect("alien://"); 
 
-  TagAna->ChainLocalTags(inputDir);
-  
-
-  // Temporary workaround to avoid problems with AOD tags.
-  TChain* chain = new TChain("aodTree");
-  TString inFileName("AliAOD.root");
-  inFileName.Prepend(Form("%s/",inputDir));
-  chain->Add(inFileName);
-  // When problems will be solved and/or you manage in getting a
-  // Run*.Merged.AOD.tag.root, substitute with the following lines:
 
-  //TChain* chain = 0x0;
-  //chain = TagAna->QueryTags(runCuts,lhcCuts,detCuts,evCuts);
+  TString outFileName("singleMuAnalysis.root");
+  outFileName.Prepend(Form("%s/",outputDir));  
 
+  // Get the chain.
+  TChain* chain = 0x0;
+  if(mode!=kMproof) chain = CreateChain(mode, inputPath, aodFilename);
 
   //____________________________________________//
   // Make the analysis manager
@@ -54,20 +81,90 @@ void runSingleMuAnalysis(Char_t *inputDir=".", Char_t *outputDir=".") {
   AliVEventHandler* aodH = new AliAODInputHandler;
   mgr->SetInputEventHandler(aodH);
   //____________________________________________//
-  // 1st Pt task
+  // Single muon task
   AliAnalysisTaskSingleMu *task1 = new AliAnalysisTaskSingleMu("SingleMu");
   mgr->AddTask(task1);
   // Create containers for input/output
   AliAnalysisDataContainer *cinput1 = mgr->CreateContainer("cchain1",TChain::Class(),AliAnalysisManager::kInputContainer);
-  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("cobjArray1", TObjArray::Class(),AliAnalysisManager::kOutputContainer,outFileName.Data());
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("tree1", TTree::Class(),AliAnalysisManager::kOutputContainer,outFileName.Data());
   
   //____________________________________________//
   mgr->ConnectInput(task1,0,cinput1);
   mgr->ConnectOutput(task1,0,coutput1);
   if (!mgr->InitAnalysis()) return;
   mgr->PrintStatus();
-  mgr->StartAnalysis("local",chain);
+
+  if(mode==kMproof)
+    mgr->StartAnalysis(modeName[mode].Data(), inputPath, nRuns, offset);
+  else 
+    mgr->StartAnalysis(modeName[mode].Data(),chain);
 
   timer.Stop();
   timer.Print();
 }
+
+
+//______________________________________________________________________________
+Bool_t SetupPar(char* pararchivename)
+{
+  if (pararchivename) {
+    FileStat_t fs;
+    char pararchivenameFull[1024];
+    sprintf(pararchivenameFull, "%s.par", pararchivename);
+    if(gSystem->GetPathInfo(pararchivenameFull, fs)){
+      Error("SetupPar", "PAR Archive %s not found!\nPlease either copy it in the current directory\nor run full aliroot", pararchivenameFull);
+      return kFALSE;
+    }
+    char processline[1024];
+    sprintf(processline,".! tar xvzf %s.par",pararchivename);
+    gROOT->ProcessLine(processline);
+    const char* ocwd = gSystem->WorkingDirectory();
+    gSystem->ChangeDirectory(pararchivename);
+    printf("Current directory = %s\n",gSystem->pwd());
+
+    // check for BUILD.sh and execute
+    if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
+      printf("*******************************\n");
+      printf("*** Building PAR archive    ***\n");
+      printf("*******************************\n");
+
+      if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
+       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
+       return kFALSE;
+      }
+    }
+    // check for SETUP.C and execute
+    if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
+      printf("*******************************\n");
+      printf("*** Setup PAR archive       ***\n");
+      printf("*******************************\n");
+      gROOT->Macro("PROOF-INF/SETUP.C");
+    }
+       
+    gSystem->ChangeDirectory("../");
+    return kTRUE;
+  }
+  return kFALSE;
+}
+
+
+//______________________________________________________________________________
+TChain* CreateChain(Int_t mode, Char_t* inputPath, Char_t* aodFilename = "AliAOD.root")
+{
+  printf("*******************************\n");
+  printf("*** Getting the Chain       ***\n");
+  printf("*******************************\n");
+  TChain *chain = 0x0;
+  if(mode==kMgridInteractive || mode==kMgridBatch){
+    AliTagAnalysis *analysis = new AliTagAnalysis();
+    chain = analysis->GetChainFromCollection(inputPath,"aodTree");
+  }
+  else{
+    chain = new TChain("aodTree");
+    TString inFileName(aodFilename);
+    inFileName.Prepend(Form("%s/",inputPath));
+    chain->Add(inFileName);
+  }
+  //if (chain) chain->ls();
+  return chain;
+}