Update of the Xiaoming code for pp900 first muon analysis. Fixing wanirngs and violti...
authormartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 10 Mar 2010 20:44:35 +0000 (20:44 +0000)
committermartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 10 Mar 2010 20:44:35 +0000 (20:44 +0000)
23 files changed:
PWG3/PWG3muonLinkDef.h
PWG3/libPWG3muon.pkg
PWG3/muon/AddTaskMuonsHF.C
PWG3/muon/AliAODMuonPair.cxx [deleted file]
PWG3/muon/AliAODMuonPair.h [deleted file]
PWG3/muon/AliAODMuonTrack.cxx [deleted file]
PWG3/muon/AliAODMuonTrack.h [deleted file]
PWG3/muon/AliAnalysisTaskSEMuonsHF.cxx
PWG3/muon/AliAnalysisTaskSEMuonsHF.h
PWG3/muon/AliDimuInfoStoreMC.cxx [new file with mode: 0644]
PWG3/muon/AliDimuInfoStoreMC.h [new file with mode: 0644]
PWG3/muon/AliDimuInfoStoreRD.cxx [new file with mode: 0644]
PWG3/muon/AliDimuInfoStoreRD.h [new file with mode: 0644]
PWG3/muon/AliMCMuonPair.cxx [deleted file]
PWG3/muon/AliMCMuonPair.h [deleted file]
PWG3/muon/AliMCMuonTrack.cxx [deleted file]
PWG3/muon/AliMCMuonTrack.h [deleted file]
PWG3/muon/AliMuonInfoStoreMC.cxx [new file with mode: 0644]
PWG3/muon/AliMuonInfoStoreMC.h [new file with mode: 0644]
PWG3/muon/AliMuonInfoStoreRD.cxx [new file with mode: 0644]
PWG3/muon/AliMuonInfoStoreRD.h [new file with mode: 0644]
PWG3/muon/AliMuonsHFHeader.cxx
PWG3/muon/AliMuonsHFHeader.h

index 2dbac69..0394fd7 100644 (file)
 #pragma link C++ class AliAnalysisTaskCreateMixedDimuons+;
 #pragma link C++ class AliAnalysisTaskMuonAODCreation+;
 #pragma link C++ class AliAnalysisTaskMuonDistributions+;
+#pragma link C++ class AliMuonInfoStoreRD+;
+#pragma link C++ class AliDimuInfoStoreRD+;
+#pragma link C++ class AliMuonInfoStoreMC+;
+#pragma link C++ class AliDimuInfoStoreMC+;
 #pragma link C++ class AliMuonsHFHeader+;
-#pragma link C++ class AliMCMuonTrack+;
-#pragma link C++ class AliMCMuonPair+;
-#pragma link C++ class AliAODMuonTrack+;
-#pragma link C++ class AliAODMuonPair+;
 #pragma link C++ class AliAnalysisTaskSEMuonsHF+;
 #pragma link C++ class AliAnalysisTaskDimuonCFContainerBuilder+;
 #endif
index 639a22e..867e28d 100644 (file)
@@ -16,12 +16,12 @@ SRCS:=       muon/AliAnalysisTaskESDMuonFilter.cxx \
                  muon/AliAnalysisTaskCreateMixedDimuons.cxx \
                 muon/AliAnalysisTaskMuonAODCreation.cxx \
                 muon/AliAnalysisTaskMuonDistributions.cxx \
-                muon/AliMuonsHFHeader.cxx \
-                muon/AliMCMuonTrack.cxx \
-                muon/AliMCMuonPair.cxx \
-                muon/AliAODMuonTrack.cxx \
-                muon/AliAODMuonPair.cxx \
-                muon/AliAnalysisTaskSEMuonsHF.cxx \
+                 muon/AliMuonInfoStoreRD.cxx \
+                 muon/AliDimuInfoStoreRD.cxx \
+                 muon/AliMuonInfoStoreMC.cxx \
+                 muon/AliDimuInfoStoreMC.cxx \
+                 muon/AliMuonsHFHeader.cxx \
+                 muon/AliAnalysisTaskSEMuonsHF.cxx \
                 muon/AliAnalysisTaskDimuonCFContainerBuilder.cxx
 
 
index 9f2acf8..9c56da8 100644 (file)
@@ -28,38 +28,36 @@ AliAnalysisTaskSEMuonsHF* AddTaskMuonsHF(Int_t mode=0, Bool_t isMC=kFALSE, Bool_
     }
   }
 
-  // set cuts for single muon track selection
-  Double_t cuts[10]={-1.,      // 0, min of 3-momentum
-                     999999.,  // 1, max of 3-momnentum
-                     1.,       // 2, PtMin
-                     999999.,  // 3, PtMax
-                     -4.,      // 4, EtaMin
-                     -2.5,     // 5, EtaMax
-                     -1.,      // 6, DCAmin
-                     10.,      // 7, DCAmax
-                     0.5,      // 8, for trigger
-                     3.5.      // 9, for trigger
-                    };
-
+  // set cuts for events or muons selection
+  Double_t cutsEvsH[3] ={-999999.0,   // low limit of Ncontrs
+                          999999.0,   // up limit of |vz|
+                          999999.0};  // up limit of vt
+  Double_t cutsMuon[10]={-999999.0,   // 0, min of 3-momentum
+                          999999.0,   // 1, max of 3-momnentum
+                         -999999.0,   // 2, PtMin
+                          999999.0,   // 3, PtMax
+                         -999999.0,   // 4, EtaMin
+                          999999.0,   // 5, EtaMax
+                         -999999.0,   // 6, DCAmin
+                          999999.0,   // 7, DCAmax
+                         -999999.0,   // 8, for trigger
+                          999999.0};  // 9, for trigger
+  Double_t cutsDimu[10]={-999999.0, 999999.0,  // single muon cuts used for dimuon selection
+                         -999999.0, 999999.0,
+                         -999999.0, 999999.0,
+                         -999999.0, 999999.0,
+                         -999999.0, 999999.0};
   AliAnalysisTaskSEMuonsHF *taskMuonsHF = new AliAnalysisTaskSEMuonsHF("MuonsHF Analysis Task");
   taskMuonsHF->SetAnaMode(mode);
   taskMuonsHF->SetIsUseMC(isMC);
   taskMuonsHF->SetIsOutputTree(isTree);
-  taskMuonsHF->SetSingleMuonCuts(cuts);
+  taskMuonsHF->SetEvsHCuts(cutsEvsH);
+  taskMuonsHF->SetMuonCuts(cutsMuon);
+  taskMuonsHF->SetDimuCuts(cutsDimu);
   mgr->AddTask(taskMuonsHF);
-  
-  TString outputfile = AliAnalysisManager::GetCommonFileName();
-  outputfile += ":PWG3Muon_MuonHF";
-
-  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("listHisEventH",TList::Class(),AliAnalysisManager::kOutputContainer,outputfile);
-  AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("listHisSingleMuon",TList::Class(),AliAnalysisManager::kOutputContainer,outputfile);
-  AliAnalysisDataContainer *coutput4 = mgr->CreateContainer("listHisDimuon",TList::Class(),AliAnalysisManager::kOutputContainer,outputfile);
 
   mgr->ConnectInput(taskMuonsHF,0,mgr->GetCommonInputContainer());
   if (isTree) mgr->ConnectOutput(taskMuonsHF,0,mgr->GetCommonOutputContainer());
-  mgr->ConnectOutput(taskMuonsHF,1,coutput2);
-  mgr->ConnectOutput(taskMuonsHF,2,coutput3);
-  mgr->ConnectOutput(taskMuonsHF,3,coutput4);
 
   return taskMuonsHF;
 }
diff --git a/PWG3/muon/AliAODMuonPair.cxx b/PWG3/muon/AliAODMuonPair.cxx
deleted file mode 100644 (file)
index 7b294dc..0000000
+++ /dev/null
@@ -1,68 +0,0 @@
-#include "AliAODMuonTrack.h"
-#include "AliAODMuonPair.h"
-
-ClassImp(AliAODMuonPair)
-
-//-----------------------------------------------------------------------------
-AliAODMuonPair::AliAODMuonPair() :
-TObject(),
-fP(),
-fCharge(0)
-{
-  //
-  // default constructor
-  //
-  for (Int_t i=0; i<2; i++) {
-    fTrk[i] = 0;
-    fTrigger[i] = 0;
-    fDca[i] = 0.;
-    fChi2[i] = 0.;
-  }
-}
-
-//-----------------------------------------------------------------------------
-AliAODMuonPair::AliAODMuonPair(AliAODMuonTrack *trk0, AliAODMuonTrack *trk1) :
-TObject(),
-fP(),
-fCharge(0)
-{
-  //
-  // default constructor
-  //
-  for (Int_t i=0; i<2; i++) {
-    fTrigger[i]=0; fDca[i]=0.; fChi2[i]=0.;
-  }
-
-  fTrk[0] = trk0;
-  fTrk[1] = trk1;
-  FillPairInfo();
-}
-
-//-----------------------------------------------------------------------------
-AliAODMuonPair::~AliAODMuonPair()
-{
-  //
-  // destructor
-  //
-}
-
-//-----------------------------------------------------------------------------
-void AliAODMuonPair::FillPairInfo()
-{
-  AliAODMuonTrack *trk0 = (AliAODMuonTrack*)fTrk[0].GetObject();
-  AliAODMuonTrack *trk1 = (AliAODMuonTrack*)fTrk[1].GetObject();
-
-  fP = trk0->GetP() + trk1->GetP();
-  fCharge = trk0->GetCharge() + trk1->GetCharge();
-
-  fTrigger[0] = trk0->GetTrigger();
-  fTrigger[1] = trk1->GetTrigger();
-
-  fDca[0] = trk0->GetDCA();
-  fDca[1] = trk1->GetDCA();
-
-  fChi2[0] = trk0->GetChi2();
-  fChi2[1] = trk1->GetChi2();
-
-  return;
-}
diff --git a/PWG3/muon/AliAODMuonPair.h b/PWG3/muon/AliAODMuonPair.h
deleted file mode 100644 (file)
index 85c1fe7..0000000
+++ /dev/null
@@ -1,42 +0,0 @@
-#ifndef ALIAODMUONPAIR_H
-#define ALIAODMUONPAIR_H
-
-#include <TObject.h>
-#include <TRef.h>
-#include <TLorentzVector.h>
-
-#include "AliAODMuonTrack.h"
-
-class AliAODMuonPair : public TObject {
- public:
-
-  AliAODMuonPair();
-  AliAODMuonPair(AliAODMuonTrack *trk0, AliAODMuonTrack *trk1);
-  virtual ~AliAODMuonPair();
-
-  AliAODMuonTrack* GetTrack(Int_t i) const { return (i<2 ? (AliAODMuonTrack*)(fTrk[i].GetObject()) : 0x0); }
-
-  TLorentzVector GetP()                     const { return fP;                           }
-  Int_t          GetCharge()                const { return fCharge;                      }
-  Int_t          GetMuMatchTrigger(Int_t i) const { return (i<2 ? fTrigger[i] : -9999 ); }
-  Double_t       GetMuDCA(Int_t i)          const { return (i<2 ? fDca[i]     : -9999.); }
-  Double_t       GetMuChi2(Int_t i)         const { return (i<2 ? fChi2[i]    : -9999.); }
-
- protected:
-
-  TRef fTrk[2];
-  void FillPairInfo();
-
- private:
-
-
-  TLorentzVector fP;
-  Int_t fCharge;
-  Int_t fTrigger[2];
-  Double_t fDca[2];
-  Double_t fChi2[2];
-
-  ClassDef(AliAODMuonPair, 3);
-};
-
-#endif
diff --git a/PWG3/muon/AliAODMuonTrack.cxx b/PWG3/muon/AliAODMuonTrack.cxx
deleted file mode 100644 (file)
index dd3be40..0000000
+++ /dev/null
@@ -1,115 +0,0 @@
-#include <TObject.h>
-#include <TDatabasePDG.h>
-#include <TLorentzVector.h>
-#include <TClonesArray.h>
-
-#include "AliAODTrack.h"
-#include "AliESDMuonTrack.h"
-#include "AliMCMuonTrack.h"
-#include "AliAODMuonTrack.h"
-
-ClassImp(AliAODMuonTrack)
-
-//-----------------------------------------------------------------------------
-AliAODMuonTrack::AliAODMuonTrack() :
-TObject(),
-fP(),
-fCharge(0),
-fTrigger(0),
-fDca(0.),
-fChi2(0.),
-fCentr(0.)
-{
-  //
-  // default constructor
-  //
-}
-
-//-----------------------------------------------------------------------------
-AliAODMuonTrack::AliAODMuonTrack(AliAODTrack *trk) :
-TObject(),
-fP(),
-fCharge(0),
-fTrigger(0),
-fDca(0.),
-fChi2(0.),
-fCentr(0.)
-{
-  //
-  // default constructor
-  //
- this->FillTrackInfo(trk);
-}
-
-//-----------------------------------------------------------------------------
-AliAODMuonTrack::AliAODMuonTrack(AliESDMuonTrack *trk) :
-TObject(),
-fP(),
-fCharge(0),
-fTrigger(0),
-fDca(0.),
-fChi2(0.),
-fCentr(0.)
-{
-  //
-  // default constructor
-  //
- this->FillTrackInfo(trk);
-}
-
-//-----------------------------------------------------------------------------
-AliAODMuonTrack::~AliAODMuonTrack()
-{
-  //
-  // destructor
-  //
-}
-
-//-----------------------------------------------------------------------------
-void AliAODMuonTrack::FillTrackInfo(AliAODTrack *trk)
-{
-  Double_t mMu    = TDatabasePDG::Instance()->GetParticle(13)->Mass();
-  Double_t px     = trk->Px();
-  Double_t py     = trk->Py();
-  Double_t pz     = trk->Pz();
-  Double_t energy = TMath::Sqrt(mMu*mMu + px*px + py*py + pz*pz);
-  fP.SetPxPyPzE(px,py,pz,energy);
-  fCharge = trk->Charge();
-  fTrigger = trk->GetMatchTrigger();
-  fDca = trk->DCA();
-  fChi2 = trk->Chi2perNDF();
-  return;
-}
-
-//-----------------------------------------------------------------------------
-void AliAODMuonTrack::FillTrackInfo(AliESDMuonTrack *trk)
-{
-  trk->LorentzP(fP);
-  fCharge = trk->Charge();
-  fTrigger = trk->GetMatchTrigger();
-  fDca = trk->GetDCA();
-  fChi2 = trk->GetChi2()/(2.*trk->GetNHit()-5.);
-  return;
-}
-
-//-----------------------------------------------------------------------------
-Bool_t AliAODMuonTrack::SelectSingleMuon(Double_t cuts[10])
-{
-  TLorentzVector lorentzP = this->GetP();
-  Double_t p = lorentzP.P();
-  if (p<cuts[0] || (p>cuts[1])) return kFALSE;
-
-  Double_t pt = lorentzP.Pt();
-  if (pt<cuts[2] || pt>cuts[3]) return kFALSE;
-
-  Double_t eta = lorentzP.Eta();
-  if (eta<cuts[4] || eta>cuts[5]) return kFALSE;
-
-  Double_t dca = this->GetDCA();
-  if (dca<cuts[6] || dca>cuts[7]) return kFALSE;
-
-  Int_t trigger = this->GetTrigger();
-  if (trigger<cuts[8] || trigger>cuts[9]) return kFALSE;
-
-  return kTRUE;
-}
diff --git a/PWG3/muon/AliAODMuonTrack.h b/PWG3/muon/AliAODMuonTrack.h
deleted file mode 100644 (file)
index 371b84a..0000000
+++ /dev/null
@@ -1,43 +0,0 @@
-#ifndef ALIAODMUONTRACK_H
-#define ALIAODMUONTRACK_H
-
-#include <TObject.h>
-#include <TLorentzVector.h>
-#include <TClonesArray.h>
-
-#include "AliAODTrack.h"
-#include "AliESDMuonTrack.h"
-
-class AliAODMuonTrack : public TObject {
- public:
-
-  AliAODMuonTrack();
-  AliAODMuonTrack(AliAODTrack *trk);
-  AliAODMuonTrack(AliESDMuonTrack *trk);
-  virtual ~AliAODMuonTrack();
-
-  Bool_t SelectSingleMuon(Double_t cuts[10]);
-
-  TLorentzVector GetP()       const { return fP;       }
-  Int_t          GetCharge()  const { return fCharge;  }
-  Int_t          GetTrigger() const { return fTrigger; }
-  Double_t       GetDCA()     const { return fDca;     }
-  Double_t       GetChi2()    const { return fChi2;    }
-  Double_t       GetCentr()   const { return fCentr;   }
-
- private:
-
-  void FillTrackInfo(AliAODTrack *trk);
-  void FillTrackInfo(AliESDMuonTrack *trk);
-
-  TLorentzVector fP;
-  Short_t fCharge;
-  Int_t fTrigger;
-  Double_t fDca;
-  Double_t fChi2;
-  Double_t fCentr;  // used for PbPb conllisions
-
-  ClassDef(AliAODMuonTrack, 3);
-};
-
-#endif
index 09a0a3e..9063a6b 100644 (file)
@@ -1,26 +1,52 @@
-#include <TH1F.h>
-#include <TH2F.h>
-#include <TString.h>
+/**************************************************************************
+ * Copyright(c) 1998-2008, 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.                  *
+ **************************************************************************/
+
+/////////////////////////////////////////////////////////////
+//
+// AliAnalysisTaskSE for the single muon and dimuon from HF analysis,
+// using the classes AliMuonsHFHeader,
+//                   AliMuonInfoStoreRD,
+//                   AliDimuInfoStoreRD,
+//                   AliMuonInfoStoreMC,
+//                   AliDimuInfoStoreMC.
+//
+// Author: X-M. Zhang, zhang@clermont.in2p3.fr
+//                     zhangxm@iopp.ccnu.edu.cn
+/////////////////////////////////////////////////////////////
+
+#include <TList.h>
 #include <TClonesArray.h>
+#include <TFile.h>
 
 #include "AliAnalysisManager.h"
 #include "AliMCEventHandler.h"
 #include "AliAODMCParticle.h"
-#include "AliVVertex.h"
-#include "AliVHeader.h"
 #include "AliVEvent.h"
 #include "AliESDEvent.h"
 #include "AliAODEvent.h"
 #include "AliESDMuonTrack.h"
 #include "AliAODTrack.h"
 #include "AliMuonsHFHeader.h"
-#include "AliAODMuonTrack.h"
-#include "AliAODMuonPair.h"
-#include "AliMCMuonTrack.h"
-#include "AliMCMuonPair.h"
-#include "AliAnalysisTaskSE.h"
+#include "AliMuonInfoStoreRD.h"
+#include "AliMuonInfoStoreMC.h"
+#include "AliDimuInfoStoreRD.h"
+#include "AliDimuInfoStoreMC.h"
 #include "AliAnalysisTaskSEMuonsHF.h"
 
+class AliAnalysisTaskSE;
+
 ClassImp(AliAnalysisTaskSEMuonsHF)
 
 //_____________________________________________________________________________
@@ -30,25 +56,15 @@ fAnaMode(0),
 fIsOutputTree(kFALSE),
 fIsUseMC(kFALSE),
 fHeader(0),
-fMuTrkClArr(0),
-fMuPairClArr(0),
-fListHisAtEvLevel(0),
-fListHisSingleMuon(0),
-fListHisDimuon(0)
+fMuonClArr(0),
+fDimuClArr(0),
+fListHisHeader(0),
+fListHisMuon(),
+fListHisDimu(0)
 {
   //
   // Default constructor
   //
-  fSingleMuonCuts[0] = -1.;
-  fSingleMuonCuts[1] = 999999.;
-  fSingleMuonCuts[2] = -1.;
-  fSingleMuonCuts[3] = 999999.;
-  fSingleMuonCuts[4] = -1.;
-  fSingleMuonCuts[5] = 999999.;
-  fSingleMuonCuts[6] = -1.;
-  fSingleMuonCuts[7] = 999999.;
-  fSingleMuonCuts[8] = -1.;
-  fSingleMuonCuts[9] = 4.;
 }
 
 //_____________________________________________________________________________
@@ -58,29 +74,15 @@ fAnaMode(0),
 fIsOutputTree(kFALSE),
 fIsUseMC(kFALSE),
 fHeader(0),
-fMuTrkClArr(0),
-fMuPairClArr(0),
-fListHisAtEvLevel(0),
-fListHisSingleMuon(0),
-fListHisDimuon(0)
+fMuonClArr(0),
+fDimuClArr(0),
+fListHisHeader(0),
+fListHisMuon(0),
+fListHisDimu(0)
 {
   //
   // Constructor
   //
-  fSingleMuonCuts[0] = -1.;
-  fSingleMuonCuts[1] = 999999.;
-  fSingleMuonCuts[2] = -1.;
-  fSingleMuonCuts[3] = 999999.;
-  fSingleMuonCuts[4] = -1.;
-  fSingleMuonCuts[5] = 999999.;
-  fSingleMuonCuts[6] = -1.;
-  fSingleMuonCuts[7] = 999999.;
-  fSingleMuonCuts[8] = -1.;
-  fSingleMuonCuts[9] = 4.;
-
-  DefineOutput(1, TList::Class());
-  DefineOutput(2, TList::Class());
-  DefineOutput(3, TList::Class());
 }
 
 //_____________________________________________________________________________
@@ -89,41 +91,71 @@ AliAnalysisTaskSEMuonsHF::~AliAnalysisTaskSEMuonsHF()
   //
   // Default destructor
   //
-  if (fHeader) { delete fHeader; fHeader=NULL; }
-  if (fMuTrkClArr)  { delete fMuTrkClArr;  fMuTrkClArr =NULL; }
-  if (fMuPairClArr) { delete fMuPairClArr; fMuPairClArr=NULL; }
+  if (fHeader)    { delete fHeader;     fHeader    =NULL; }
+  if (fMuonClArr) { delete fMuonClArr;  fMuonClArr =NULL; }
+  if (fDimuClArr) { delete fDimuClArr;  fDimuClArr =NULL; }
 
-  if (fListHisSingleMuon) { delete fListHisSingleMuon; fListHisSingleMuon=NULL; }
-  if (fListHisDimuon)     { delete fListHisDimuon;     fListHisDimuon=NULL;     }
+  if (fListHisHeader) { delete    fListHisHeader; fListHisHeader=NULL; }
+  if (fListHisMuon)   { delete [] fListHisMuon;   fListHisMuon  =NULL; }
+  if (fListHisDimu)   { delete [] fListHisDimu;   fListHisDimu  =NULL; }
 }
 
 //_____________________________________________________________________________
-void AliAnalysisTaskSEMuonsHF::UserCreateOutputObjects()
+void AliAnalysisTaskSEMuonsHF::Init()
 {
-  CreateOutoutHistosAtEvnetLevel();
-  if (fAnaMode!=2) CreateOutputHistosSingleMuon();
-  if (fAnaMode!=1) CreateOutputHistosDimuon();
+  // Initialization
+  // Setting and initializing the running mode and status
 
+  AliMuonsHFHeader::SetAnaMode(fAnaMode);
+  AliMuonsHFHeader::SetIsMC(fIsUseMC);
   if (!fHeader) {
     fHeader = new AliMuonsHFHeader();
-    fHeader->SetName("MuonsHFHeader");
+    fHeader->SetName(AliMuonsHFHeader::StdBranchName());
   }
-  if (!fMuTrkClArr) {
-    if (fIsUseMC) fMuTrkClArr = new TClonesArray("AliMCMuonTrack", 0);
-    else fMuTrkClArr = new TClonesArray("AliAODMuonTrack", 0);
-    fMuTrkClArr->SetName("MuonTrack");
+
+  if (!fMuonClArr) {
+    if (fIsUseMC) { 
+      fMuonClArr = new TClonesArray("AliMuonInfoStoreMC", 0);
+      fMuonClArr->SetName(AliMuonInfoStoreMC::StdBranchName());
+    } else {
+      fMuonClArr = new TClonesArray("AliMuonInfoStoreRD", 0);
+      fMuonClArr->SetName(AliMuonInfoStoreRD::StdBranchName());
+    }
   }
 
-  if (fAnaMode!=1 && !fMuPairClArr) {
-    if (fIsUseMC) fMuPairClArr = new TClonesArray("AliMCMuonPair", 0);
-    else fMuPairClArr = new TClonesArray("AliAODMuonPair", 0);
-    fMuPairClArr->SetName("MuonPair");
+  if (fAnaMode!=1 && !fDimuClArr) {
+    if (fIsUseMC) {
+      fDimuClArr = new TClonesArray("AliDimuInfoStoreMC", 0);
+      fDimuClArr->SetName(AliDimuInfoStoreMC::StdBranchName());
+    } else {
+      fDimuClArr = new TClonesArray("AliDimuInfoStoreRD", 0);
+      fDimuClArr->SetName(AliDimuInfoStoreRD::StdBranchName());
+    }
+  }
+
+  return;
+}
+
+//_____________________________________________________________________________
+void AliAnalysisTaskSEMuonsHF::UserCreateOutputObjects()
+{
+  // Create the output container
+
+  if (!fListHisHeader) fListHisHeader = new TList();
+  if (fAnaMode!=2 && !fListHisMuon) {
+    if (fIsUseMC) fListHisMuon = new TList[AliMuonInfoStoreMC::NSources()];
+    else fListHisMuon = new TList();
   }
+  if (fAnaMode!=1 && !fListHisDimu) {
+    if (fIsUseMC) fListHisDimu = new TList[AliDimuInfoStoreMC::NSources()];
+    else fListHisDimu = new TList();
+  }
+  fHeader->CreateHistograms(fListHisHeader, fListHisMuon, fListHisDimu);
 
   if (fIsOutputTree) {
     AddAODBranch("AliMuonsHFHeader", &fHeader);
-    AddAODBranch("TClonesArray", &fMuTrkClArr);
-    if (fAnaMode!=1) AddAODBranch("TClonesArray", &fMuPairClArr);
+    AddAODBranch("TClonesArray", &fMuonClArr);
+    if (fAnaMode!=1) AddAODBranch("TClonesArray", &fDimuClArr);
   }
 
   return;
@@ -132,12 +164,13 @@ void AliAnalysisTaskSEMuonsHF::UserCreateOutputObjects()
 //_____________________________________________________________________________
 void AliAnalysisTaskSEMuonsHF::UserExec(Option_t *)
 {
+  // Execute analysis for current event:
+  // muon event header & (di)muon info store
+
   AliVEvent *event = dynamic_cast<AliVEvent*>(InputEvent());
   //if (AODEvent() && IsStandardAOD()) event = dynamic_cast<AliVEvent*>(AODEvent());
   TString evName = event->IsA()->GetName();
   Bool_t isAOD = ((evName=="AliAODEvent") ? kTRUE : kFALSE);
-  fHeader->SetHeader(((AliVHeader*)event->GetHeader()));
-  fHeader->SetVertex(((AliVVertex*)event->GetPrimaryVertex()));
   event = 0x0;
 
   Int_t ntrks = 0;
@@ -148,373 +181,136 @@ void AliAnalysisTaskSEMuonsHF::UserExec(Option_t *)
   if (isAOD) {
     aod = dynamic_cast<AliAODEvent*>(InputEvent());
     //if (AODEvent() && IsStandardAOD()) aod = dynamic_cast<AliAODEvent*>(AODEvent());
-    if (!aod) {
-      AliError("AOD event not found. Nothing done!");
-      return;
-    }
+    if (!aod) { AliError("AOD event not found. Nothing done!"); return; }
     if (fIsUseMC) {
       mcClArr = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
-      if (!mcClArr) {
-        AliError("MC Array not found. Nothing done!");
-        return;
-      }
+      if (!mcClArr) { AliError("MC Array not found. Nothing done!"); return; }
     }
+    fHeader->SetEvent(aod);
     ntrks = aod->GetNTracks();
   } else {
     esd = dynamic_cast<AliESDEvent*>(InputEvent());
-    if (!esd) {
-      AliError("ESD event not found. Nothing done!");
-      return;
-    }
+    if (!esd) { AliError("ESD event not found. Nothing done!"); return; }
     if (fIsUseMC) {
-      mcH = (AliMCEventHandler*)
-              ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
-      if (!mcH) {
-        AliError("MC Handler not found. Nothing done!");
-        return;
-      }
+      mcH = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
+      if (!mcH) { AliError("MC Handler not found. Nothing done!"); return; }
     }
+    fHeader->SetEvent(esd);
     ntrks = esd->GetNumberOfMuonTracks();
   }
 
-  fMuTrkClArr->Delete();
-  TClonesArray &muTrkRef = *fMuTrkClArr;
-  Int_t countN = fMuTrkClArr->GetEntriesFast();
+  fHeader->FillHistosEventH(fListHisHeader);
+
+  fMuonClArr->Delete();
+  TClonesArray &muonRef = *fMuonClArr;
+  Int_t countN = fMuonClArr->GetEntriesFast();
 
-  Int_t theMult = 0;
-  AliAODTrack *trkAOD = 0;
-  AliESDMuonTrack *trkESD = 0;
-  AliAODMuonTrack *track = 0;
-  AliMCMuonTrack *trkMC = 0;
+  AliAODTrack        *trkAOD = 0;
+  AliESDMuonTrack    *trkESD = 0;
+  AliMuonInfoStoreRD *trkRD  = 0;
+  AliMuonInfoStoreMC *trkMC  = 0;
   for (Int_t itrk=0; itrk<ntrks; itrk++) {  // loop over all tracks
-    track=0; trkMC=0;
     if (isAOD) {
       trkAOD = (AliAODTrack*)aod->GetTrack(itrk);
-      if (!trkAOD->IsMuonTrack()) { trkAOD=0; continue; }
-      if (fIsUseMC) trkMC = new AliMCMuonTrack(trkAOD, mcClArr);
-      else track = new AliAODMuonTrack(trkAOD);
+      if (!trkAOD->IsMuonTrack()) { trkAOD=0; trkRD=0; trkMC=0; continue; }
+      if (fIsUseMC) trkMC = new AliMuonInfoStoreMC(trkAOD, mcClArr);
+      else trkRD = new AliMuonInfoStoreRD(trkAOD);
+      trkAOD = 0;
     } else {
       trkESD = (AliESDMuonTrack*)esd->GetMuonTrack(itrk);
-      if (!trkESD->ContainTrackerData()) { trkESD=0; continue; }
-      if (fIsUseMC) trkMC = new AliMCMuonTrack(trkESD, esd, mcH);
-      else track = new AliAODMuonTrack(trkESD);
+      if (!trkESD->ContainTrackerData()) { trkESD=0; trkRD=0; trkMC=0; continue; }
+      if (fIsUseMC) trkMC = new AliMuonInfoStoreMC(trkESD, esd, mcH);
+      else trkRD = new AliMuonInfoStoreRD(trkESD);
+      trkESD = 0;
     }
 
-    if (fAnaMode!=2) {
-      if (track) FillDistributionsSingleMuon(track);
-      else FillDistributionsSingleMuon((AliAODMuonTrack*)trkMC, trkMC->GetSource());
+    if (trkRD) {
+      if (fAnaMode!=2) fHeader->FillHistosMuonRD(fListHisMuon, trkRD);
+      new(muonRef[countN++]) AliMuonInfoStoreRD(*trkRD);
     }
-
-    if (fIsOutputTree || fAnaMode!=1) {
-      if (track) new(muTrkRef[countN++]) AliAODMuonTrack(*track);
-      else new(muTrkRef[countN++]) AliMCMuonTrack(*trkMC);
+    if (trkMC) {
+      if (fAnaMode!=2) fHeader->FillHistosMuonMC(fListHisMuon, trkMC);
+      new(muonRef[countN++]) AliMuonInfoStoreMC(*trkMC);
     }
 
-    theMult++;
-    trkESD = 0;
-    trkAOD = 0;
-    if (track) { delete track; track=0; }
+    if (trkRD) { delete trkRD; trkRD=0; }
     if (trkMC) { delete trkMC; trkMC=0; }
   }  // end loop of all tracks
-  fHeader->SetMultSingleMuon(theMult);
-  FillDistributionsAtEventLeavel();
-  PostData(1, fListHisAtEvLevel);
-  if (fAnaMode!=2) PostData(2, fListHisSingleMuon);
+
   if (fAnaMode==1) return;
+  fDimuClArr->Delete();
+  countN = fDimuClArr->GetEntriesFast();
+  TClonesArray &dimuRef = *fDimuClArr;
 
-  theMult = 0;
-  fMuPairClArr->Delete();
-  TClonesArray &muPairRef = *fMuPairClArr;
-  countN = fMuPairClArr->GetEntriesFast();
-  AliAODMuonPair *pair = 0;
-  AliMCMuonPair *pairMC = 0;
-  ntrks = fMuTrkClArr->GetEntriesFast();
+  AliDimuInfoStoreRD *dimuRD = 0;
+  AliDimuInfoStoreMC *dimuMC = 0;
+  ntrks = fMuonClArr->GetEntriesFast();
   for (Int_t itrk=0; itrk<ntrks-1; itrk++) {  // 1st loop over muon tracks
-    for (Int_t jtrk=itrk+1; jtrk<ntrks; jtrk++) {  // 2nd loop over muon tracks
-      pair=0; pairMC=0;
+    for (Int_t jtrk=itrk+1; jtrk<ntrks; jtrk++) {  // 2nd loop ofver muon tracks
       if (fIsUseMC)
-        pairMC = new AliMCMuonPair((AliMCMuonTrack*)fMuTrkClArr->At(itrk), (AliMCMuonTrack*)fMuTrkClArr->At(jtrk));
+        dimuMC = new AliDimuInfoStoreMC((AliMuonInfoStoreMC*)fMuonClArr->At(itrk), (AliMuonInfoStoreMC*)fMuonClArr->At(jtrk));
       else
-        pair = new AliAODMuonPair((AliAODMuonTrack*)fMuTrkClArr->At(itrk), (AliAODMuonTrack*)fMuTrkClArr->At(jtrk));
-
-      if (pair) FillDistributionsDimuon(pair);
-      else FillDistributionsDimuon(pairMC, pairMC->GetSource());
+        dimuRD = new AliDimuInfoStoreRD((AliMuonInfoStoreRD*)fMuonClArr->At(itrk), (AliMuonInfoStoreRD*)fMuonClArr->At(jtrk));
 
-      if (fIsOutputTree) {
-        if (pair) new(muPairRef[countN++]) AliAODMuonPair(*pair);
-        else new(muPairRef[countN++]) AliMCMuonPair(*pairMC);
+      if (dimuRD) {
+        fHeader->FillHistosDimuRD(fListHisDimu, dimuRD);
+        if (fIsOutputTree) new(dimuRef[countN++]) AliDimuInfoStoreRD(*dimuRD);
+      }
+      if (dimuMC) {
+        fHeader->FillHistosDimuMC(fListHisDimu, dimuMC);
+        if (fIsOutputTree) new(dimuRef[countN++]) AliDimuInfoStoreMC(*dimuMC);
       }
 
-      theMult++;
-      if (pair)   { delete pair;   pair  =0; }
-      if (pairMC) { delete pairMC; pairMC=0; }
-    }  // end 2nd loop of muon traks
+      if (dimuRD) { delete dimuRD; dimuRD=0; }
+      if (dimuMC) { delete dimuMC; dimuMC=0; }
+    }  // end 2nd loop of muon tracks
   }  // end 1st loop of muon tracks
-  fHeader->SetMultDimuon(theMult);
 
-  PostData(3, fListHisDimuon);
   return;
 }
 
 //_____________________________________________________________________________
 void AliAnalysisTaskSEMuonsHF::Terminate(Option_t *)
 {
-  printf("This is Terminate!\n");
-  // adding the fitting
-  return;
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskSEMuonsHF::CreateOutoutHistosAtEvnetLevel()
-{
-  if (!fListHisAtEvLevel) fListHisAtEvLevel = new TList();
-  Int_t nbins[kNHistEv];
-  Double_t xlow[kNHistEv], xup[kNHistEv];
-  TString name[kNHistEv], axis[kNHistEv], unit[kNHistEv];
-  nbins[kHistVx]=200; xlow[kHistVx]  =-50.; xup[kHistVx]  =50.; name[kHistVx]  =  "vx"; axis[kHistVx] ="v_{x}"; unit[kHistVx]="cm";
-  nbins[kHistVy]=200; xlow[kHistVy]  =-50.; xup[kHistVy]  =50.; name[kHistVy]  =  "vy"; axis[kHistVy] ="v_{y}"; unit[kHistVy]="cm";
-  nbins[kHistVz]=200; xlow[kHistVz]  =-50.; xup[kHistVz]  =50.; name[kHistVz]  =  "vz"; axis[kHistVz] ="v_{z}"; unit[kHistVz]="cm";
-  nbins[kHistMult]=9; xlow[kHistMult]=-0.5; xup[kHistMult]=8.5; name[kHistMult]="mult"; axis[kHistMult]="mult"; unit[kHistMult]="";
-
-  TH1F *histo = 0;
-  TString histName, histTitle;
-  for (Int_t ihist=0; ihist<kNHistEv; ihist++) {  // loop over histos
-    histName  = Form("hsMu%s", name[ihist].Data());
-    histTitle = Form("%s of event", axis[ihist].Data());
-    histo = new TH1F(histName.Data(), histTitle.Data(), nbins[ihist], xlow[ihist], xup[ihist]);
-    histo->GetXaxis()->SetTitle(Form("%s %s", axis[ihist].Data(), unit[ihist].Data()));
-    histo->Sumw2();
-    fListHisAtEvLevel->AddAt(histo, ihist);
-  }  // end loop over histos
-
-  return;
-}
-
-
-//_____________________________________________________________________________
-void AliAnalysisTaskSEMuonsHF::CreateOutputHistosSingleMuon()
-{
-  if (!fListHisSingleMuon) fListHisSingleMuon = new TList();
-
-  Int_t nbins[kNHistMu];
-  Double_t xlow[kNHistMu], xup[kNHistMu];
-  TString name[kNHistMu], axis[kNHistMu], unit[kNHistMu];
-  nbins[kHistP]  =300 ; xlow[kHistP]=   0. ; xup[kHistP]  =150.; name[kHistP]  ="P"  ; axis[kHistP]  ="P"      ; unit[kHistP]  ="[GeV/c]";
-  nbins[kHistPt] =120 ; xlow[kHistPt]=  0. ; xup[kHistPt] = 30.; name[kHistPt] ="Pt" ; axis[kHistPt] ="p_{t}"  ; unit[kHistPt] ="[GeV/c]";
-  nbins[kHistEta]=50  ; xlow[kHistEta]=-4. ; xup[kHistEta]=-2.5; name[kHistEta]="Eta"; axis[kHistEta]="#eta"   ; unit[kHistEta]="";
-  nbins[kHistDca]=1000; xlow[kHistDca]= 0. ; xup[kHistDca]=500.; name[kHistDca]="Dca"; axis[kHistDca]="DCA"    ; unit[kHistDca]="cm";
-  nbins[kHistTrg]=   4; xlow[kHistTrg]=-0.5; xup[kHistTrg]= 3.5; name[kHistTrg]="trg"; axis[kHistTrg]="Trigger"; unit[kHistTrg]="";
-
-  TH1F *histo = 0;
-  TString histName, histTitle;
-  Int_t index=-1, lowBound=0;
-  for (Int_t ihist=0; ihist<kNHistMu; ihist++) {  // loop over histos
-    index = lowBound + ihist;
-    histName  = Form("hsMu%s", name[ihist].Data());
-    histTitle = Form("%s of muon candidates", axis[ihist].Data());
-    histo = new TH1F(histName.Data(), histTitle.Data(), nbins[ihist], xlow[ihist], xup[ihist]);
-    histo->GetXaxis()->SetTitle(Form("%s %s", axis[ihist].Data(), unit[ihist].Data()));
-    histo->Sumw2();
-    fListHisSingleMuon->AddAt(histo, index);
-  }  // end loop over histos
-
-  if (fIsUseMC) CreateOutputHistosSingleMuonMC(nbins, xlow, xup, name, axis, unit);
-  TH2F *hVzPt = new TH2F("hVzPt", "v_{z} vs. p_{t}", 300, -30., 30., 120, 0., 30);
-  hVzPt->GetXaxis()->SetTitle("v_{z} [cm]");
-  hVzPt->GetYaxis()->SetTitle("p_{t} [GeV/c]");
-  fListHisSingleMuon->AddLast(hVzPt);
-  return;
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskSEMuonsHF::CreateOutputHistosDimuon()
-{
-  if (!fListHisDimuon) fListHisDimuon = new TList();
-
-  TString dimuName[kNMuMus], dimuTitle[kNMuMus];
-  dimuName[kMuNMuP]="MuNMuP"; dimuTitle[kMuNMuP]="#mu^{-}#mu^{+}";
-  dimuName[kMuNMuN]="MuNMuN"; dimuTitle[kMuNMuN]="#mu^{-}#mu^{-}";
-  dimuName[kMuPMuP]="MuPMuP"; dimuTitle[kMuPMuP]="#mu^{+}#mu^{+}";
-
-  Int_t nbins[kNHistDimu];
-  Double_t xlow[kNHistDimu], xup[kNHistDimu];
-  TString name[kNHistDimu], axis[kNHistDimu], unit[kNHistDimu];
-  nbins[kInvM]  =300; xlow[kInvM]  =0.; xup[kInvM]  =30.; name[kInvM]  ="InvM"; axis[kInvM]  ="M";     unit[kInvM]  ="[GeV/c^{2}]";
-  nbins[kPtPair]=120; xlow[kPtPair]=0.; xup[kPtPair]=60.; name[kPtPair]="Pt";   axis[kPtPair]="p_{t}"; unit[kPtPair]="[GeV/c]";
-
-  TH1F *histo = 0;
-  TString histName, histTitle;
-  Int_t index=-1, lowBound=0;
-  for (Int_t idmu=0; idmu<kNMuMus; idmu++) {  // loop over different kinds of dimuon
-    for (Int_t ihist=0; ihist<kNHistDimu; ihist++) {  // loop over histos
-      index = lowBound + ihist + idmu*kNHistDimu;
-      histName  = Form("h%s%s", dimuName[idmu].Data(), name[ihist].Data());
-      histTitle = Form("%s of %s candidates", axis[ihist].Data(), dimuTitle[idmu].Data());
-      histo = new TH1F(histName.Data(), histTitle.Data(), nbins[ihist], xlow[ihist], xup[ihist]);
-      histo->GetXaxis()->SetTitle(Form("%s(%s) %s", axis[ihist].Data(), dimuTitle[idmu].Data(), unit[ihist].Data()));
-      histo->Sumw2();
-      fListHisDimuon->AddAt(histo, index);
-    }  // end loop over histos
-  }  // end loop of different kinds of dimuon
-
-  if (fIsUseMC) CreateOutputHistosDimuonMC(nbins, xlow, xup, name, axis, unit, dimuName, dimuTitle);
-  return;
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskSEMuonsHF::CreateOutputHistosSingleMuonMC(Int_t *nbins, Double_t *xlow, Double_t *xup,
-                                                              TString *name, TString *axis, TString *unit)
-{
-  // define the histos of single muons
-  TString srcName[kNSingleMuSrcs];
-  srcName[kBeautyMu]    = "Bottom";
-  srcName[kCharmMu]     = "Charm";
-  srcName[kPrimaryMu]   = "Primary";
-  srcName[kSecondaryMu] = "Secondary";
-  srcName[kNotMu]       = "NotMu";
-
-  TH1F *histo = 0;
-  TString histName, histTitle;
-  Int_t index=-1, lowBound=kNHistMu;
-  for (Int_t src=0; src<kNSingleMuSrcs; src++) {  // loop over single muon src
-    for (Int_t ihist=0; ihist<kNHistMu; ihist++) {  // loop over histos
-      index = lowBound + ihist + src*kNHistMu;
-      histName  = Form("hsMu%s_%s", srcName[src].Data(), name[ihist].Data());
-      histTitle = Form("%s of muon candidates<-%s", axis[ihist].Data(), srcName[src].Data());
-      histo = new TH1F(histName.Data(), histTitle.Data(), nbins[ihist], xlow[ihist], xup[ihist]);
-      histo->GetXaxis()->SetTitle(Form("%s %s", axis[ihist].Data(), unit[ihist].Data()));
-      histo->Sumw2();
-      fListHisSingleMuon->AddAt(histo, index);
-    }  // end loop over histos
-  }  // // end loop of single muon src
-
-  return;
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskSEMuonsHF::CreateOutputHistosDimuonMC(Int_t *nbins, Double_t *xlow, Double_t *xup,
-                                                          TString *name, TString *axis, TString *unit,
-                                                          TString *dimuName, TString *dimuTitle)
-{
-  TString srcName[kNDimuSrcs];
-  srcName[kBBdiff]    = "BBdiff";
-  srcName[kBchain]    = "Bchain";
-  srcName[kDDdiff]    = "DDdiff";
-  srcName[kDchain]    = "Dchain";
-  srcName[kResonance] = "Resonance";
-  srcName[kUncorr]    = "Uncorr";
-
-  TH1F *histo = 0;
-  TString histName, histTitle;
-  Int_t index=-1, lowBound=kNMuMus*kNHistDimu;
-  for (Int_t idmu=0; idmu<kNMuMus; idmu++) {  // loop over different kinds of dimuon
-    for (Int_t ihist=0; ihist<kNHistDimu; ihist++) {  // loop over histos
-      for (Int_t src=0; src<kNDimuSrcs; src++) {  // loop over dimuon src
-        index = lowBound + src + ihist*kNDimuSrcs + idmu*kNHistDimu*kNDimuSrcs;
-        histName  = Form("hs%s%s_%s", dimuName[idmu].Data(), srcName[src].Data(), name[ihist].Data());
-        histTitle = Form("%s of %s<-%s", name[ihist].Data(), dimuName[idmu].Data(), srcName[src].Data());
-        histo = new TH1F(histName.Data(), histTitle.Data(), nbins[ihist], xlow[ihist], xup[ihist]);
-        histo->GetXaxis()->SetTitle(Form("%s(%s) %s", axis[ihist].Data(), dimuTitle[idmu].Data(), unit[ihist].Data()));
-        fListHisDimuon->AddAt(histo, index);
-      }  // end loop of dimuon src
-    }  // end loop over histos
-  }  // end loop of different kinds of dimuon
-
-  return;
-}
+  // Terminate analysis
 
-//_____________________________________________________________________________
-void AliAnalysisTaskSEMuonsHF::FillDistributionsAtEventLeavel()
-{
-  Double_t disEv[kNHistEv];
-  disEv[kHistVx]   = fHeader->GetXv();
-  disEv[kHistVy]   = fHeader->GetYv();
-  disEv[kHistVz]   = fHeader->GetZv();
-  disEv[kHistMult] = (Double_t)fHeader->GetMultSingleMuon();
-
-  Int_t index=-1, lowBound=0;
-  for (Int_t ihist=0; ihist<kNHistEv; ihist++) {
-    index = lowBound + ihist;
-    ((TH1F*)fListHisAtEvLevel->At(index))->Fill(disEv[ihist]);
-  }
-
-  return;
-}
+  TFile *fout = new TFile("muonsHF.root", "RECREATE");
+  gDirectory->mkdir("header");
+  gDirectory->cd("header");
+  fListHisHeader->Write();
+  gDirectory->cd("..");
 
-//_____________________________________________________________________________
-void AliAnalysisTaskSEMuonsHF::FillDistributionsSingleMuon(AliAODMuonTrack *track, Int_t src)
-{
-  if (fHeader->IsUnrecoVertex()) return;
-  if (!(track->SelectSingleMuon(fSingleMuonCuts))) return;
-
-  TLorentzVector lorentzP = track->GetP();
-  Double_t disMu[kNHistMu];
-  disMu[kHistP] = lorentzP.P();
-  disMu[kHistPt] = lorentzP.Pt();
-  disMu[kHistEta] = lorentzP.Eta();
-  disMu[kHistDca] = track->GetDCA();
-  disMu[kHistTrg] = (Double_t)track->GetTrigger();
-
-  Int_t index=-1, lowBound=0;
-  for (Int_t ihist=0; ihist<kNHistMu; ihist++) {
-    index = lowBound + ihist;
-    ((TH1F*)fListHisSingleMuon->At(index))->Fill(disMu[ihist]);
-  } 
-
-  if (fIsUseMC && src>=0) FillDistributionsSingleMuonMC(disMu, src);
-
-  Double_t vtx[3];
-  fHeader->GetXYZ(vtx);
-  ((TH2F*)fListHisSingleMuon->FindObject("hVzPt"))->Fill(vtx[2], disMu[kHistPt]);
-  return;
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskSEMuonsHF::FillDistributionsSingleMuonMC(Double_t *disMu, Int_t src)
-{
-  if (src<0) return;
-  Int_t index=-1, lowBound=kNHistMu;
-  for (Int_t ihist=0; ihist<kNHistMu; ihist++) {
-    index = lowBound + ihist + src*kNHistMu;
-    ((TH1F*)fListHisSingleMuon->At(index))->Fill(disMu[ihist]);
+  if (fAnaMode!=2) {
+    gDirectory->mkdir("muon");
+    gDirectory->cd("muon");
+    if (fIsUseMC) {
+      char *muonName[] = {"BottomMu", "CharmMu", "PrimaryMu", "SecondaryMu", "NotMu", "Undentified", "All"};
+      for (Int_t i=AliMuonInfoStoreMC::NSources(); i--;) {
+        gDirectory->mkdir(muonName[i]);
+        gDirectory->cd(muonName[i]);
+        fListHisMuon[i].Write();
+        gDirectory->cd("..");
+      }
+    } else fListHisMuon->Write();
+    gDirectory->cd("..");
   }
-  return;
-}
 
-//_____________________________________________________________________________
-void AliAnalysisTaskSEMuonsHF::FillDistributionsDimuon(AliAODMuonPair *pair, Int_t src)
-{
-  if (!(pair->GetTrack(0)->SelectSingleMuon(fSingleMuonCuts)) ||
-      !(pair->GetTrack(1)->SelectSingleMuon(fSingleMuonCuts))) return;
-
-  TLorentzVector lorentzP = pair->GetP();
-  Double_t disDimu[kNHistDimu];
-  disDimu[kInvM]   = lorentzP.Mag();
-  disDimu[kPtPair] = lorentzP.Pt();
-
-  Int_t dimuK=-1, charge=pair->GetCharge();
-  if (charge==0) dimuK = kMuNMuP;
-  if (charge<0)  dimuK = kMuNMuN;
-  if (charge>0)  dimuK = kMuPMuP;
-
-  Int_t index=-1, lowBound=0;
-  for (Int_t ihist=0; ihist<kNHistDimu; ihist++) {
-    index = lowBound + ihist + dimuK*kNHistDimu;
-    ((TH1F*)fListHisDimuon->At(index))->Fill(disDimu[ihist]);
+  if (fAnaMode!=1) {
+    gDirectory->mkdir("dimu");
+    gDirectory->cd("dimu");
+    if (fIsUseMC) {
+      char *dimuName[] = {"BBdiff", "Bchain", "DDdiff", "Dchain", "Resonance", "Background", "All"};
+      for (Int_t i=AliDimuInfoStoreMC::NSources(); i--;) {
+        gDirectory->mkdir(dimuName[i]);
+        gDirectory->cd(dimuName[i]);
+        fListHisDimu[i].Write();
+        gDirectory->cd("..");
+      }
+    } else fListHisDimu->Write();
+    gDirectory->cd("..");
   }
 
-  if (fIsUseMC && src>=0) FillDistributionsDimuonMC(disDimu, dimuK, src);
-  return;
-}
-
-//_____________________________________________________________________________
-void AliAnalysisTaskSEMuonsHF::FillDistributionsDimuonMC(Double_t *disDimu, Int_t dimuK, Int_t src)
-{
-  if (src<0) return;
-  Int_t index=-1, lowBound=kNMuMus*kNHistDimu;
-  for (Int_t ihist=0; ihist<kNHistDimu; ihist++) {
-    index = lowBound + src + ihist*kNDimuSrcs + dimuK*kNHistDimu*kNDimuSrcs;
-    ((TH1F*)fListHisDimuon->At(index))->Fill(disDimu[ihist]);
-  }
+  fout->Write();
+  fout->Close();
   return;
 }
index b8095b4..c48c8aa 100644 (file)
@@ -1,15 +1,23 @@
 #ifndef ALIANALYSISTASKSEMUONSHF_H
 #define ALIANALYSISTASKSEMUONSHF_H
 
-#include <TString.h>
-#include <TList.h>
-#include <TClonesArray.h>
+/* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//*************************************************************************
+// Class AliAnalysisTaskSEMuonsHF
+// AliAnalysisTaskSE for the single muon and dimuon from HF analysis
+// Author: X-M. Zhang, zhang@clermont.in2p3.fr
+//                     zhangxm@iopp.ccnu.edu.cn
+//*************************************************************************
 
-#include "AliMuonsHFHeader.h"
-#include "AliAODMuonTrack.h"
-#include "AliAODMuonPair.h"
 #include "AliAnalysisTaskSE.h"
 
+class TString;
+class TList;
+class TClonesArray;
+class AliMuonsHFHeader;
+
 class AliAnalysisTaskSEMuonsHF : public AliAnalysisTaskSE {
  public:
 
@@ -17,106 +25,40 @@ class AliAnalysisTaskSEMuonsHF : public AliAnalysisTaskSE {
   AliAnalysisTaskSEMuonsHF(const char *name);
   virtual ~AliAnalysisTaskSEMuonsHF();
 
+  virtual void Init();
+  virtual void LocalInit() { Init(); }
   virtual void UserCreateOutputObjects();
   virtual void UserExec(Option_t *opt);
   virtual void Terminate(Option_t *opt);
 
-  void SetAnaMode(Int_t mode) { fAnaMode = (mode<3 ? mode: 0); }
+  void SetAnaMode(Int_t mode) { fAnaMode = (mode<3 ? mode : 0); }
   void SetIsOutputTree(Bool_t ist) { fIsOutputTree = ist; }
   void SetIsUseMC(Bool_t isMC) { fIsUseMC = isMC; }
 
-  void SetSingleMuonCuts(Double_t cuts[10]) {
-    for (Int_t i=0; i<10; i++) fSingleMuonCuts[i]=cuts[i];
-  }
+  void SetEvsHCuts(Double_t cuts[3])  const { AliMuonsHFHeader::SetSelectionCuts(cuts);   }
+  void SetMuonCuts(Double_t cuts[10]) const { AliMuonInfoStoreRD::SetSelectionCuts(cuts); }
+  void SetDimuCuts(Double_t cuts[10]) const { AliDimuInfoStoreRD::SetSelectionCuts(cuts); }
 
  private:
 
-  void CreateOutoutHistosAtEvnetLevel();
-  void CreateOutputHistosSingleMuon();
-  void CreateOutputHistosDimuon();
-  void CreateOutputHistosSingleMuonMC(Int_t *nbins, Double_t *xlow, Double_t *xup,
-                                      TString *name, TString *axis, TString *unit);
-  void CreateOutputHistosDimuonMC(Int_t *nbins, Double_t *xlow, Double_t *xup,
-                                  TString *name, TString *axis, TString *unit,
-                                  TString *dimuName, TString *dimuTitle);
-
-  void FillDistributionsAtEventLeavel();
-  void FillDistributionsSingleMuon(AliAODMuonTrack *track, Int_t src=-1);
-  void FillDistributionsDimuon(AliAODMuonPair *pair, Int_t src=-1);
-  void FillDistributionsSingleMuonMC(Double_t *disMu, Int_t src);
-  void FillDistributionsDimuonMC(Double_t *disDimu, Int_t dimuK, Int_t src);
-
-  enum {  // histos at event level
-    kHistVx,    // x position of vtx
-    kHistVy,    // y position of vtx
-    kHistVz,    // z position of vtx
-    kHistMult,  // multiplicity of single muon track
-    kNHistEv    // number of histos at event level
-  };
-  enum {  // distributions of single mu
-    kHistP,    // histo of single mu p
-    kHistPt,   // histo of single mu pt
-    kHistEta,  // histo of single mu eta
-    kHistDca,  // histo of single mu DCA
-    kHistTrg,  // histo of single mu pass trigger
-    kNHistMu   // number of histos of single mu
-  };
-  enum {  // distribution of dimuon
-    kInvM,      // invariance mass if dimuon
-    kPtPair,    // pt of dimu pair
-    kNHistDimu  // number of histos of dimuon
-  };
-
-  enum {  // MC sources of single muon
-    kBeautyMu,      // mu<-b
-    kCharmMu,       // mu<-c
-    kPrimaryMu,     // primary muon
-    kSecondaryMu,   // mu<-secondary paritcle
-    kNotMu,         // not muon
-    kNSingleMuSrcs  // number of single muon sources
-  };
-  enum {  // MC soureces of dimuon
-    kBBdiff,     // dimuon<-BB_diff
-    kBchain,     // dimuon<-B-chain
-    kDDdiff,     // dimuon<-DD_diff
-    kDchain,     // dimuon<-D-chain
-    kResonance,  // dimuon<-resonances
-    kUncorr,     // uncorr dimuon
-    kNDimuSrcs   // number of dimuon sources
-  };
-  enum {  // different kinds of correlated dimu
-    kMuNMuP,  // mu-mu+
-    kMuNMuN,  // mu-mu-
-    kMuPMuP,  // mu+mu+
-    kNMuMus   // number of corr dimu kinds
-  };
-
-  Int_t fAnaMode;  // = 0, ana both single muon and dimuon
-                   // = 1, ana single muon
-                   // = 2, ana dimuon
+  AliAnalysisTaskSEMuonsHF(const AliAnalysisTaskSEMuonsHF&);
+  AliAnalysisTaskSEMuonsHF& operator=(const AliAnalysisTaskSEMuonsHF&);
+
+  Int_t fAnaMode;        // = 0, ana both single muon and dimuon
+                         // = 1, ana single muon
+                         // = 2, ana dimuon
   Bool_t fIsOutputTree;  // flag used to switch on/off tree output
   Bool_t fIsUseMC;       // flag used to switch on/off MC ana
 
-  Double_t fSingleMuonCuts[10];  // 0, max of 3-momentum
-                                 // 1, min of 3-momentum
-                                 // 2, pt_Min
-                                 // 3, pt_Max
-                                 // 4, eta_Min
-                                 // 5, eta_Max
-                                 // 6, dca_Min
-                                 // 7, dca_Max
-                                 // 8, about trigger
-                                 // 9, about trigger
-
-  AliMuonsHFHeader *fHeader;   // output clones array for info at ev level
-  TClonesArray *fMuTrkClArr;   // output clones array for single mu
-  TClonesArray *fMuPairClArr;  // output clones array for dimu
-
-  TList *fListHisAtEvLevel;   // output list of histos at event level
-  TList *fListHisSingleMuon;  // output list of histos for single mu
-  TList *fListHisDimuon;      // output list of histos for dimuon
-
-  ClassDef(AliAnalysisTaskSEMuonsHF, 5);
+  AliMuonsHFHeader *fHeader;  // output for info at ev level
+  TClonesArray *fMuonClArr;   // output clones array for single mu
+  TClonesArray *fDimuClArr;   // output clones array for dimu
+
+  TList *fListHisHeader;  // output list of histos at event level
+  TList *fListHisMuon;    // output list of histos for single mu
+  TList *fListHisDimu;    // output list of histos for dimuon
+
+  ClassDef(AliAnalysisTaskSEMuonsHF, 6);
 };
 
 #endif
diff --git a/PWG3/muon/AliDimuInfoStoreMC.cxx b/PWG3/muon/AliDimuInfoStoreMC.cxx
new file mode 100644 (file)
index 0000000..5ba0e2c
--- /dev/null
@@ -0,0 +1,150 @@
+/**************************************************************************
+ * Copyright(c) 1998-2006, 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.                  *
+ **************************************************************************/
+
+/////////////////////////////////////////////////////////////
+//
+// class used to extract and store info of MC particles
+//
+// Author: X-M. Zhang, zhang@clermont.in2p3.fr
+//                     zhangxm@iopp.ccnu.edu.cn
+/////////////////////////////////////////////////////////////
+
+#include "AliMuonInfoStoreMC.h"
+#include "AliDimuInfoStoreRD.h"
+#include "AliDimuInfoStoreMC.h"
+
+ClassImp(AliDimuInfoStoreMC)
+
+const TString AliDimuInfoStoreMC::fgkStdBranchName("DimuMC");
+const Int_t   AliDimuInfoStoreMC::fgkNSources = 7;
+
+//-----------------------------------------------------------------------------
+AliDimuInfoStoreMC::AliDimuInfoStoreMC() :
+AliDimuInfoStoreRD(),
+fIsFull(kFALSE),
+fLorentzP(),
+fSource(-1)
+{
+  //
+  // default constructor
+  //
+}
+
+//-----------------------------------------------------------------------------
+AliDimuInfoStoreMC::AliDimuInfoStoreMC(AliMuonInfoStoreMC* const trk0, AliMuonInfoStoreMC* const trk1, Bool_t full) :
+AliDimuInfoStoreRD(),
+fIsFull(full),
+fLorentzP(),
+fSource(-1)
+{
+  //
+  // default constructor
+  //
+  fMuonRef[0] = trk0;
+  fMuonRef[1] = trk1;
+  fLorentzP = trk0->LorentzP() + trk1->LorentzP();
+  AliDimuInfoStoreRD::FillDimuInfo();
+  if (fIsFull) this->FindDimuonSourceFull();
+  else this->FindDimuonSourceFast();
+}
+
+//-----------------------------------------------------------------------------
+AliDimuInfoStoreMC::AliDimuInfoStoreMC(const AliDimuInfoStoreMC &src) :
+AliDimuInfoStoreRD(src),
+fIsFull(src.fIsFull),
+fLorentzP(src.fLorentzP),
+fSource(src.fSource)
+{
+  //
+  // copy constructor
+  //
+}
+
+//-----------------------------------------------------------------------------
+AliDimuInfoStoreMC& AliDimuInfoStoreMC::operator=(const AliDimuInfoStoreMC &src)
+{
+  //
+  // assignment constructor
+  //
+  if(&src==this) return *this;
+
+  fIsFull   = src.fIsFull;
+  fLorentzP = src.fLorentzP;
+  fSource   = src.fSource;
+
+  return *this;
+}
+
+
+//-----------------------------------------------------------------------------
+AliDimuInfoStoreMC::~AliDimuInfoStoreMC()
+{
+  //
+  // destructor
+  //
+}
+
+//-----------------------------------------------------------------------------
+void AliDimuInfoStoreMC::FindDimuonSourceFast()
+{
+  // find corr relation of two particles (fast for p-p)
+
+  AliMuonInfoStoreMC *trk0 = (AliMuonInfoStoreMC*)fMuonRef[0].GetObject();
+  Int_t src0 = trk0->MuonSource();
+  if (src0<0 || src0==4 || src0==3) {
+    fSource=5; return;
+  }
+
+  AliMuonInfoStoreMC *trk1 = (AliMuonInfoStoreMC*)fMuonRef[1].GetObject();
+  Int_t src1 = trk1->MuonSource();
+  if (src1<0 || src1==4 || src1==3) {
+    fSource=5; return;
+  }
+
+  // Drell-Yan is expected very small at LHC, we ingore it
+  Int_t np0 = trk0->NParents() - 1;
+  if (np0<0) {
+    fSource=5; return;
+  }
+  Int_t np1 = trk1->NParents() - 1;
+  if (np1<0) {
+    fSource=5; return;
+  }
+
+  if (trk0->IsMotherAResonance(np0) &&
+      trk1->IsMotherAResonance(np1) &&
+      trk0->ParentIndex(np0)==trk1->ParentIndex(np1)) {
+    fSource=4; return;
+  }
+
+  if (src0==0 && src1==0) {
+    if ((trk0->ParentIndex(0))==(trk1->ParentIndex(0)))
+      fSource = 1;
+    else
+      fSource = 0;
+    return;
+  }
+
+  if (src0==1 && src1==1) {
+    if ((trk0->ParentIndex(0))==(trk1->ParentIndex(0)))
+      fSource = 3;
+    else
+      fSource = 2;
+    return;
+  }
+
+  fSource = 5;
+  return;
+}
diff --git a/PWG3/muon/AliDimuInfoStoreMC.h b/PWG3/muon/AliDimuInfoStoreMC.h
new file mode 100644 (file)
index 0000000..2a3cd2f
--- /dev/null
@@ -0,0 +1,56 @@
+#ifndef ALIDIMUINFOSTOREMC_H
+#define ALIDIMUINFOSTOREMC_H
+
+/* Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//***********************************************************
+// Class AliDimuInfoStoreMC
+// class used to extract and store info of MC particles
+// Author: X-M. Zhang, zhang@clermont.in2p3.fr
+//                     zhangxm@iopp.ccnu.edu.cn
+//***********************************************************
+
+#include <TLorentzVector.h>
+#include "AliDimuInfoStoreRD.h"
+
+class AliMuonInfoStoreMC;
+class AliDimuInfoStoreMC : public AliDimuInfoStoreRD {
+ public:
+
+  AliDimuInfoStoreMC();
+  AliDimuInfoStoreMC(AliMuonInfoStoreMC* const trk0, AliMuonInfoStoreMC* const trk1, Bool_t full=kFALSE);
+  AliDimuInfoStoreMC(const AliDimuInfoStoreMC &src);
+  AliDimuInfoStoreMC& operator=(const AliDimuInfoStoreMC &src);
+  virtual ~AliDimuInfoStoreMC();
+
+  AliMuonInfoStoreMC* Muon(Int_t i) const { return (i<2 ? (AliMuonInfoStoreMC*)(fMuonRef[i].GetObject()) : 0x0); }
+
+  TLorentzVector LorentzP() const { return fLorentzP; }
+  Int_t DimuSource() const { return fSource; }
+
+  static const char* StdBranchName() { return fgkStdBranchName.Data(); }
+  static const Int_t NSources()      { return fgkNSources;             }
+
+
+ private:
+
+  void FindDimuonSourceFast();
+  void FindDimuonSourceFull();
+
+  static const TString fgkStdBranchName;  // Standard branch name
+  static const Int_t   fgkNSources;       // num. of dimuon sources
+
+  Bool_t fIsFull;  // flag of using full analysis (for Pb-Pb)
+  TLorentzVector fLorentzP;  // lorentz momentum of MC particle
+  Int_t fSource;  // = 0, BBdiff
+                  // = 1, Bchain
+                  // = 2, DDdiff
+                  // = 3, Dchain
+                  // = 4, Resonance
+                  // = 5, UnCorr bkg
+
+  ClassDef(AliDimuInfoStoreMC, 1);
+};
+
+#endif
diff --git a/PWG3/muon/AliDimuInfoStoreRD.cxx b/PWG3/muon/AliDimuInfoStoreRD.cxx
new file mode 100644 (file)
index 0000000..8dbdaa9
--- /dev/null
@@ -0,0 +1,142 @@
+/**************************************************************************
+ * Copyright(c) 1998-2006, 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.                  *
+ **************************************************************************/
+
+/////////////////////////////////////////////////////////////
+//
+// class used to extract and store reco info of dimu candidate
+//
+// Author: X-M. Zhang, zhang@clermont.in2p3.fr
+//                     zhangxm@iopp.ccnu.edu.cn
+/////////////////////////////////////////////////////////////
+
+#include <TDatabasePDG.h>
+#include <TLorentzVector.h>
+
+#include "AliMuonInfoStoreRD.h"
+#include "AliDimuInfoStoreRD.h"
+
+ClassImp(AliDimuInfoStoreRD)
+
+const TString AliDimuInfoStoreRD::fgkStdBranchName("DimuRD");
+Double_t      AliDimuInfoStoreRD::fgCutd[10] = {-999999., 999999.,
+                                                -999999., 999999.,
+                                                -999999., 999999.,
+                                                -999999., 999999.,
+                                                -999999., 999999.};
+
+//-----------------------------------------------------------------------------
+AliDimuInfoStoreRD::AliDimuInfoStoreRD() :
+TObject(),
+fMomentum(),
+fCharge(0),
+fInvM(0.)
+{
+  //
+  // default constructor
+  //
+  for (Int_t i=2; i--;) fMuonRef[i] = 0;
+}
+
+//-----------------------------------------------------------------------------
+AliDimuInfoStoreRD::AliDimuInfoStoreRD(AliMuonInfoStoreRD* const trk0, AliMuonInfoStoreRD* const trk1) :
+TObject(),
+fMomentum(),
+fCharge(0),
+fInvM(0.)
+{
+  //
+  // default constructor
+  //
+  fMuonRef[0] = trk0;
+  fMuonRef[1] = trk1;
+  FillDimuInfo();
+}
+
+//-----------------------------------------------------------------------------
+AliDimuInfoStoreRD::AliDimuInfoStoreRD(const AliDimuInfoStoreRD &src) :
+TObject(src),
+fMomentum(src.fMomentum),
+fCharge(src.fCharge),
+fInvM(src.fInvM)
+{
+  //
+  // copy constructor
+  //
+  for (Int_t i=2; i--;) fMuonRef[i] = src.fMuonRef[i];
+}
+
+//-----------------------------------------------------------------------------
+AliDimuInfoStoreRD& AliDimuInfoStoreRD::operator=(const AliDimuInfoStoreRD &src)
+{
+  //
+  // assignment constructor
+  //
+  if(&src==this) return *this;
+
+  fMomentum = src.fMomentum;
+  fCharge   = src.fCharge;
+  fInvM     = src.fInvM;
+  for (Int_t i=2; i--;) fMuonRef[i] = src.fMuonRef[i];
+
+  return *this;
+}
+
+//-----------------------------------------------------------------------------
+AliDimuInfoStoreRD::~AliDimuInfoStoreRD()
+{
+  //
+  // destructor
+  //
+}
+
+//-----------------------------------------------------------------------------
+void AliDimuInfoStoreRD::FillDimuInfo()
+{
+  // fill dimuon candidate info from the corresponding two muon tracks
+
+  AliMuonInfoStoreRD *trk0 = (AliMuonInfoStoreRD*)fMuonRef[0].GetObject();
+  AliMuonInfoStoreRD *trk1 = (AliMuonInfoStoreRD*)fMuonRef[1].GetObject();
+
+  fMomentum = trk0->Momentum() + trk1->Momentum();
+  fCharge   = trk0->Charge()   + trk1->Charge();
+
+  Double_t mMu = TDatabasePDG::Instance()->GetParticle(13)->Mass();
+  TLorentzVector lorentzP0, lorentzP1, lorentzP;
+  lorentzP0.SetVectM(trk0->Momentum(), mMu);
+  lorentzP1.SetVectM(trk1->Momentum(), mMu);
+  lorentzP = lorentzP0 + lorentzP1;
+  fInvM = lorentzP.Mag();
+
+  trk0 = 0;
+  trk1 = 0;
+  return;
+}
+
+//-----------------------------------------------------------------------------
+Bool_t AliDimuInfoStoreRD::DimuSelection()
+{
+  // select dimuon candidates according to the corresponding two muon tracks cuts
+
+  Double_t cutsOld[10];
+  AliMuonInfoStoreRD::SelectionCust(cutsOld);
+  AliMuonInfoStoreRD::SetSelectionCuts(AliDimuInfoStoreRD::fgCutd);
+  AliMuonInfoStoreRD *trk0 = (AliMuonInfoStoreRD*)fMuonRef[0].GetObject();
+  AliMuonInfoStoreRD *trk1 = (AliMuonInfoStoreRD*)fMuonRef[1].GetObject(); 
+  if (!trk0->MuonSelection()) { AliMuonInfoStoreRD::SetSelectionCuts(cutsOld); return kFALSE; }
+  if (!trk1->MuonSelection()) { AliMuonInfoStoreRD::SetSelectionCuts(cutsOld); return kFALSE; }
+
+  AliMuonInfoStoreRD::SetSelectionCuts(cutsOld);
+  return kTRUE;
+}
diff --git a/PWG3/muon/AliDimuInfoStoreRD.h b/PWG3/muon/AliDimuInfoStoreRD.h
new file mode 100644 (file)
index 0000000..fec5df8
--- /dev/null
@@ -0,0 +1,58 @@
+#ifndef ALIDIMUINFOSTORERD_H
+#define ALIDIMUINFOSTORERD_H
+
+/* Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//***********************************************************
+// Class AliDimuInfoStoreRD
+// class used to extract and store reco info of dimu candidate
+// Author: X-M. Zhang, zhang@clermont.in2p3.fr
+//                     zhangxm@iopp.ccnu.edu.cn
+//***********************************************************
+
+#include <TObject.h>
+#include <TRef.h>
+#include <TVector3.h>
+#include <TString.h>
+
+#include "AliMuonInfoStoreRD.h"
+
+class AliDimuInfoStoreRD : public TObject {
+ public:
+
+  AliDimuInfoStoreRD();
+  AliDimuInfoStoreRD(AliMuonInfoStoreRD* const trk0, AliMuonInfoStoreRD* const trk1);
+  AliDimuInfoStoreRD(const AliDimuInfoStoreRD &src);
+  AliDimuInfoStoreRD& operator=(const AliDimuInfoStoreRD &src);
+  virtual ~AliDimuInfoStoreRD();
+
+  AliMuonInfoStoreRD* Muon(Int_t i) const { return (i<2 ? (AliMuonInfoStoreRD*)(fMuonRef[i].GetObject()) : 0x0); }
+
+  TVector3 Momentum() const { return fMomentum; }
+  Short_t  Charge()   const { return fCharge;   }
+  Double_t InvM()     const { return fInvM;     }
+
+  Bool_t DimuSelection();
+
+  static const char* StdBranchName() { return fgkStdBranchName.Data(); }
+  static void SetSelectionCuts(Double_t cuts[10]) { for (Int_t i=10; i--;) fgCutd[i]=cuts[i]; }
+
+ protected:
+
+  void FillDimuInfo();
+  TRef fMuonRef[2];  // ref to the two corresponding muon tracks
+
+ private:
+
+  static const TString fgkStdBranchName;  // Standard branch name
+  static Double_t fgCutd[10];             // single muon cuts for dimuon selection
+
+  TVector3 fMomentum;  // 3-momentum of dimuon
+  Short_t  fCharge;    // charge of dimuon
+  Double_t fInvM;      // invariance mass of dimuon
+
+  ClassDef(AliDimuInfoStoreRD, 3);
+};
+
+#endif
diff --git a/PWG3/muon/AliMCMuonPair.cxx b/PWG3/muon/AliMCMuonPair.cxx
deleted file mode 100644 (file)
index 8fb7ba8..0000000
+++ /dev/null
@@ -1,96 +0,0 @@
-#include "AliAODMuonPair.h"
-#include "AliMCMuonTrack.h"
-#include "AliMCMuonPair.h"
-
-ClassImp(AliAODMuonPair)
-
-//-----------------------------------------------------------------------------
-AliMCMuonPair::AliMCMuonPair() :
-AliAODMuonPair(),
-fIsFull(kFALSE),
-fPGen(),
-fSource(-1)
-{
-  //
-  // default constructor
-  //
-  for (Int_t i=0; i<2; i++) fTrk[i] = 0;
-
-}
-
-//-----------------------------------------------------------------------------
-AliMCMuonPair::AliMCMuonPair(AliMCMuonTrack *trk0, AliMCMuonTrack *trk1, Bool_t full) :
-AliAODMuonPair(),
-fIsFull(full),
-fPGen(),
-fSource(-1)
-{
-  //
-  // default constructor
-  //
-  fTrk[0] = trk0;
-  fTrk[1] = trk1;
-  fPGen = trk0->GetPGen() + trk1->GetPGen();
-  AliAODMuonPair::FillPairInfo();
-  if (fIsFull) this->FindDimuonSourceFull();
-  else this->FindDimuonSourceFast();
-}
-
-//-----------------------------------------------------------------------------
-AliMCMuonPair::~AliMCMuonPair()
-{
-  //
-  // destructor
-  //
-}
-
-//-----------------------------------------------------------------------------
-void AliMCMuonPair::FindDimuonSourceFast()
-{
-  AliMCMuonTrack *trk0 = (AliMCMuonTrack*)fTrk[0].GetObject();
-  Int_t src0 = trk0->GetSource();
-  if (src0<0 || src0==4 || src0==3) {
-    fSource=5; return;
-  }
-
-  AliMCMuonTrack *trk1 = (AliMCMuonTrack*)fTrk[1].GetObject();
-  Int_t src1 = trk1->GetSource();
-  if (src1<0 || src1==4 || src1==3) {
-    fSource=5; return;
-  }
-
-  // Drell-Yan is expected very small at LHC, we ingore it
-  Int_t np0 = trk0->GetNParents() - 1;
-  if (np0<0) {
-    fSource=5; return;
-  }
-
-  Int_t np1 = trk1->GetNParents() - 1;
-  if (np1<0) {
-    fSource=5; return;
-  }
-
-  if (trk0->IsMotherAResonance(np0) && trk1->IsMotherAResonance(np1) &&
-      (trk0->GetParentIndex(np0))==(trk1->GetParentIndex(np1))) {
-    fSource=4; return;
-  }
-
-  if (src0==0 && src1==0) {
-    if ((trk0->GetParentIndex(0))==(trk1->GetParentIndex(0)))
-      fSource = 1;
-    else
-      fSource = 0;
-    return;
-  }
-
-  if (src0==1 && src1==1) {
-    if ((trk0->GetParentIndex(0))==(trk1->GetParentIndex(0)))
-      fSource = 3;
-    else
-      fSource = 2;
-    return;
-  }
-
-  fSource = 5;
-  return;
-}
diff --git a/PWG3/muon/AliMCMuonPair.h b/PWG3/muon/AliMCMuonPair.h
deleted file mode 100644 (file)
index 89003e5..0000000
+++ /dev/null
@@ -1,38 +0,0 @@
-#ifndef ALIMCMUONPAIR_H
-#define ALIMCMUONPAIR_H
-
-#include <TLorentzVector.h>
-
-#include "AliMCMuonTrack.h"
-#include "AliAODMuonPair.h"
-
-class AliMCMuonPair : public AliAODMuonPair {
- public:
-
-  AliMCMuonPair();
-  AliMCMuonPair(AliMCMuonTrack *trk0, AliMCMuonTrack *trk1, Bool_t full=kFALSE);
-  virtual ~AliMCMuonPair();
-
-  AliMCMuonTrack* GetTrack(Int_t i) const { return (i<2 ? (AliMCMuonTrack*)(fTrk[i].GetObject()) : 0x0); }
-
-  TLorentzVector GetPGen() const { return fPGen; }
-  Int_t GetSource() const { return fSource; }
-
- private:
-
-  void FindDimuonSourceFast();
-  void FindDimuonSourceFull();
-
-  Bool_t fIsFull;
-  TLorentzVector fPGen;
-  Int_t fSource;  // = 0, BBdiff
-                  // = 1, Bchain
-                  // = 2, DDdiff
-                  // = 3, Dchain
-                  // = 4, Resonance
-                  // = 5, UnCorr bkg
-
-  ClassDef(AliMCMuonPair, 1);
-};
-
-#endif
diff --git a/PWG3/muon/AliMCMuonTrack.cxx b/PWG3/muon/AliMCMuonTrack.cxx
deleted file mode 100644 (file)
index 6f646b5..0000000
+++ /dev/null
@@ -1,407 +0,0 @@
-#include <TParticle.h>
-#include <TLorentzVector.h>
-#include <TClonesArray.h>
-
-#include "AliMCEventHandler.h"
-#include "AliMCEvent.h"
-#include "AliStack.h"
-#include "AliAODMCParticle.h"
-#include "AliESDEvent.h"
-#include "AliESDMuonTrack.h"
-#include "AliESDMuonCluster.h"
-#include "AliAODTrack.h"
-#include "AliAODMuonTrack.h"
-#include "AliMCMuonTrack.h"
-
-ClassImp(AliMCMuonTrack)
-
-const Double_t AliMCMuonTrack::fgkSigmaCut = 10.;
-
-//-----------------------------------------------------------------------------
-AliMCMuonTrack::AliMCMuonTrack() :
-AliAODMuonTrack(),
-fIsFull(kFALSE),
-fPGen(),
-fTrackIndex(-1),
-fTrackPDGCode(0),
-fSource(-1),
-fNParents(0),
-fOscillation(kFALSE),
-fWeight(0.)
-{
-  //
-  // default constructor
-  //
-  for (Int_t i=0; i<fgkNParentsMax; i++) {
-    fParentIndex[i] = -1;
-    fParentPDGCode[i] = 0;
-  }
-  for (Int_t i=0; i<4; i++) {
-    fQuarkIndex[i] = -1;
-    fQuarkPDGCode[i] = 0;
-  }
-}
-
-//-----------------------------------------------------------------------------
-AliMCMuonTrack::AliMCMuonTrack(AliAODTrack *trkAOD, TClonesArray *mcClArr, Bool_t full) :
-AliAODMuonTrack(trkAOD),
-fIsFull(full),
-fPGen(),
-fTrackIndex(-1),
-fTrackPDGCode(0),
-fSource(-1),
-fNParents(0),
-fOscillation(kFALSE),
-fWeight(0.)
-{
-  //
-  // default constructor
-  //
-  AliAODMCParticle *pMC = this->FindTrackRef(trkAOD, mcClArr);
-  if (pMC) this->SetMCInfo(pMC, mcClArr);
-}
-
-//-----------------------------------------------------------------------------
-AliMCMuonTrack::AliMCMuonTrack(AliESDMuonTrack *trkESD, AliESDEvent *esd, AliMCEventHandler *mcH, Bool_t full) :
-AliAODMuonTrack(trkESD),
-fIsFull(full),
-fPGen(),
-fTrackIndex(-1),
-fTrackPDGCode(0),
-fSource(-1),
-fNParents(0),
-fOscillation(kFALSE),
-fWeight(0.)
-{
-  //
-  // default constructor
-  //
-  TParticle *pMC = this->FindTrackRef(trkESD, esd, mcH);  // do nothing when running on official train
-  if (pMC) this->SetMCInfo(pMC, mcH);
-  // MC infomation is not implement when running on the official train for they depend on the MUON module.
-  // In the case of private running, one can get the information from ESD by uncommenting the definitions
-  // of FindTrackRef() and CovESDtoMuonTrack()
-  // and including the following head files in both AliMCMuonTrack.h and this file,
-  // #include "AliMUONRecoCheck.h"
-  // #include "AliMUONVTrackStore.h"
-  // #include "AliMUONVClusterStore.h"
-  // #include "AliMUONTrackParam.h"
-  // #include "AliMUONTrack.h"
-  // #include "AliMUONVCluster.h"
-  // #include "AliMUONESDInterface.h"
-}
-
-
-
-//-----------------------------------------------------------------------------
-AliMCMuonTrack::~AliMCMuonTrack()
-{
-  //
-  // destructor
-  //
-}
-
-//-----------------------------------------------------------------------------
-AliAODMCParticle* AliMCMuonTrack::FindTrackRef(AliAODTrack *trkAOD, TClonesArray *mcClArr)
-{
-  AliAODMCParticle *pMC = 0;
-  fTrackIndex = trkAOD->GetLabel();
-  if (fTrackIndex>=0)
-    pMC = (AliAODMCParticle*)mcClArr->At(fTrackIndex);
-  return pMC;
-}
-
-//-----------------------------------------------------------------------------
-TParticle* AliMCMuonTrack::FindTrackRef(AliESDMuonTrack *trkESD, AliESDEvent *esd, AliMCEventHandler *mcH)
-{
-  TParticle *pMCRef = 0;
-
-  /*AliMUONTrack *trkMuon = CovESDtoMuonTrack(*trkESD);
-  AliMUONRecoCheck rc(esd,mcH);
-  AliMUONVTrackStore* trkRefArr = rc.TrackRefs(-1);
-  TIter next(trkRefArr->CreateIterator());
-  AliMUONTrack* trkRef = 0;
-  while ((trkRef=static_cast<AliMUONTrack*>(next()))) {
-    Bool_t trkCompArr[10];
-    Int_t nMatchClusters = trkMuon->CompatibleTrack(trkRef, fgkSigmaCut, trkCompArr);
-    Double_t matchClusterFrac = ((Double_t)nMatchClusters) / ((Double_t)trkMuon->GetNClusters());
-    if ((trkCompArr[0] || trkCompArr[1] || trkCompArr[2] || trkCompArr[3]) &&
-        (trkCompArr[6] || trkCompArr[7] || trkCompArr[8] || trkCompArr[9]) && matchClusterFrac>0.5) {
-      fTrackIndex = trkRef->GetUniqueID();
-      trkRefArr->Remove(*trkRef);
-      break;
-    }
-  }
-  if (fTrackIndex>=0) pMCRef = mcH->MCEvent()->Stack()->Particle(fTrackIndex);*/
-
-  return pMCRef;
-}
-
-//-----------------------------------------------------------------------------
-void AliMCMuonTrack::SetMCInfo(AliAODMCParticle *pMC, TClonesArray *mcClArr)
-{
-  fPGen.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->E());
-  fTrackPDGCode = pMC->GetPdgCode();
-  if (TMath::Abs(fTrackPDGCode)!=13) {
-    fSource = 4;
-    return;
-  } 
-
-  Int_t lineM = pMC->GetMother();
-  if (lineM<0) {
-    fSource = 2;
-    return;
-  }
-
-  Bool_t isPrimary =
-    ((AliAODMCParticle*)mcClArr->At(lineM))->IsPrimary();
-  if (!isPrimary) {
-    fSource = 3;
-    return;
-  }
-
-  this->FillHistoryParents(lineM, mcClArr);
-  fSource = this->SelectHFMuon();
-  return;
-}
-
-//-----------------------------------------------------------------------------
-void AliMCMuonTrack::SetMCInfo(TParticle *pMC, AliMCEventHandler *mcH)
-{
-  fPGen.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->Energy());
-  fTrackPDGCode = pMC->GetPdgCode();
-  if (TMath::Abs(fTrackPDGCode)!=13) {
-    fSource = 4;
-    return;
-  }
-
-  Int_t lineM = pMC->GetFirstMother();
-  if (lineM<0) {
-    fSource = 2;
-    return;
-  }
-
-  AliStack *stack = mcH->MCEvent()->Stack();
-  if (lineM>=stack->GetNprimary()) {
-    fSource = 3;
-    return;
-  }
-
-  this->FillHistoryParents(lineM, stack);
-  fSource = this->SelectHFMuon();
-  return;
-}
-
-//-----------------------------------------------------------------------------
-void AliMCMuonTrack::FillHistoryParents(Int_t lineM, TClonesArray *mcClArr)
-{
-  Int_t countP=-1, pdg=0;
-  Int_t parents[10], parLine[10];
-  AliAODMCParticle *mother = 0;
-  while(lineM>=0){
-    mother = (AliAODMCParticle*)mcClArr->At(lineM);
-    pdg = mother->GetPdgCode();
-    if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
-    parents[++countP] = pdg;
-    parLine[countP] = lineM;
-    lineM = mother->GetMother();
-  }
-  for(Int_t i=0; i<=countP; i++){
-    fParentIndex[i] = parLine[countP-i];
-    fParentPDGCode[i] = parents[countP-i];
-  }
-  fNParents = countP + 1;
-
-  if (fIsFull && lineM>=0) this->FillHistoryQuarks(lineM, mcClArr);
-  return;
-}
-
-//-----------------------------------------------------------------------------
-void AliMCMuonTrack::FillHistoryParents(Int_t lineM, AliStack *stack)
-{
-  Int_t countP=-1, pdg=0;
-  Int_t parents[10], parLine[10];
-  TParticle *mother = 0;
-  while(lineM>=0){
-    mother = stack->Particle(lineM);
-    pdg = mother->GetPdgCode();
-    if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
-    parents[++countP] = pdg;
-    parLine[countP] = lineM;
-    lineM = mother->GetFirstMother();
-  }
-  for(Int_t i=0; i<=countP; i++){
-    fParentIndex[i] = parLine[countP-i];
-    fParentPDGCode[i] = parents[countP-i];
-  }
-  fNParents = countP + 1;
-
-  if (fIsFull && lineM>=0) this->FillHistoryQuarks(lineM, stack);
-  return;
-}
-
-//-----------------------------------------------------------------------------
-void AliMCMuonTrack::FillHistoryQuarks(Int_t lineM, TClonesArray *mcClArr)
-{
-  // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx 
-  if (lineM<0) return;
-  Int_t countP=-1, pdg=0;
-  AliAODMCParticle *mother = 0;
-  while(lineM>=0){
-    mother = (AliAODMCParticle*)mcClArr->At(lineM);
-    pdg = mother->GetPdgCode();
-    fQuarkIndex[++countP] = lineM;
-    fQuarkPDGCode[countP] = pdg;
-    lineM = mother->GetMother();
-  }
-
-  // for PYTHIA checking
-  countP = 1;
-  for(Int_t par=0; par<4; par++) {
-    if(TMath::Abs(this->GetQuarkPDGCode(par))<6) { countP=par; break; }
-  }
-  if(this->GetQuarkIndex(countP)>-1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)) {
-    if(this->GetParentFlavour(0)!=TMath::Abs(this->GetQuarkPDGCode(countP))) {
-      AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
-                 this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP))));
-
-      pdg = this->GetQuarkPDGCode(countP);
-      Int_t line = this->GetQuarkIndex(countP);
-      this->ResetQuarkInfo();
-      while(TMath::Abs(pdg)!=this->GetParentFlavour(0)) {
-        pdg = ((AliAODMCParticle*)mcClArr->At(++line))->GetPdgCode();
-      }
-      while(line>=0){
-        mother = (AliAODMCParticle*)mcClArr->At(line);
-        pdg = mother->GetPdgCode();
-        fQuarkIndex[countP] = line;
-        fQuarkPDGCode[countP++] = pdg;
-        line = mother->GetMother();
-      }
-    }
-  }
-  return;
-}
-
-//-----------------------------------------------------------------------------
-void AliMCMuonTrack::FillHistoryQuarks(Int_t lineM, AliStack *stack)
-{
-  // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx 
-  if (lineM<0) return;
-  Int_t countP=-1, pdg=0;
-  TParticle *mother = 0;
-  while(lineM>=0){
-    mother = stack->Particle(lineM);
-    pdg = mother->GetPdgCode();
-    fQuarkIndex[++countP] = lineM;
-    fQuarkPDGCode[countP] = pdg;
-    lineM = mother->GetFirstMother();
-  }
-
-  // for PYTHIA checking
-  countP = 1;
-  for(Int_t par=0; par<4; par++) {
-    if(TMath::Abs(this->GetQuarkPDGCode(par))<6) { countP=par; break; }
-  }
-  if(this->GetQuarkIndex(countP)>-1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)) {
-    if(this->GetParentFlavour(0)!=TMath::Abs(this->GetQuarkPDGCode(countP))) {
-      AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
-                 this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP))));
-
-      pdg = this->GetQuarkPDGCode(countP);
-      Int_t line = this->GetQuarkIndex(countP);
-      this->ResetQuarkInfo();
-      while(TMath::Abs(pdg)!=this->GetParentFlavour(0)) {
-        pdg = stack->Particle(++line)->GetPdgCode();
-      }
-      while(line>=0){
-        mother = stack->Particle(line);
-        pdg = mother->GetPdgCode();
-        fQuarkIndex[countP] = line;
-        fQuarkPDGCode[countP++] = pdg;
-        line = mother->GetFirstMother();
-      }
-    }
-  }
-  return;
-}
-
-//-----------------------------------------------------------------------------
-/*AliMUONTrack* AliMCMuonTrack::CovESDtoMuonTrack(AliESDMuonTrack &trkESD)
-{
-  // method in $ALICE_ROOT/PWG3/muondep/AliAnalysisTaskESDMCLabelAddition.cxx
-  AliMUONTrack *trkMuon = new AliMUONTrack();
-  if(!trkESD.ClustersStored()) return trkMuon;
-
-  AliMUONTrackParam param;
-  AliMUONESDInterface::GetParamAtFirstCluster(trkESD, param);
-  AliMUONESDInterface::GetParamCov(trkESD, param);
-  AliMUONVClusterStore *cStore = AliMUONESDInterface::NewClusterStore();
-  AliMUONVCluster *cluster = cStore->CreateCluster(0,0,0);
-
-  AliESDMuonCluster *esdCluster = (AliESDMuonCluster*)trkESD.GetClusters().First();
-  while (esdCluster) {
-    AliMUONESDInterface::ESDToMUON(*esdCluster, *cluster);
-    param.SetZ(cluster->GetZ());
-    trkMuon->AddTrackParamAtCluster(param, *cluster, kTRUE);
-    esdCluster = (AliESDMuonCluster*)trkESD.GetClusters().After(esdCluster);
-  }
-
-  delete cluster; cluster = 0;
-  delete cStore;  cStore = 0;
-  return trkMuon;
-}*/
-
-//-----------------------------------------------------------------------------
-Int_t AliMCMuonTrack::SelectHFMuon()
-{
-  Int_t flv = GetParentFlavour(0);
-  if (flv!=4 && flv!=5) return 2;
-
-  Bool_t isRes = kFALSE;
-  Int_t i=0, nparents=this->GetNParents();
-  while (i<nparents && !isRes) {
-    isRes = IsMotherAResonance(i++);
-  }
-
-  if (isRes) return 2;
-  if (flv==5) return 0;
-  else return 1;
-}
-
-//-----------------------------------------------------------------------------
-Bool_t AliMCMuonTrack::IsDiquark(Int_t pdg)
-{
-  // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
-  pdg = TMath::Abs(pdg);
-  if(pdg>1000 && (pdg%100)<10) return kTRUE;
-  else return kFALSE;
-}
-
-//-----------------------------------------------------------------------------
-void AliMCMuonTrack::ResetQuarkInfo()
-{
-  for(Int_t pos=1; pos<4; pos++) {
-    fQuarkIndex[pos] = -1;
-    fQuarkPDGCode[pos] = 0;
-  }
-  return;
-}
-
-//-----------------------------------------------------------------------------
-Int_t AliMCMuonTrack::GetParentFlavour(Int_t i) const
-{
-  Int_t pdg = GetParentPDGCode(i);
-  pdg = TMath::Abs(pdg/100);
-  if(pdg>9) pdg /= 10;
-  return pdg;
-}
-
-//-----------------------------------------------------------------------------
-Bool_t AliMCMuonTrack::IsMotherAResonance(Int_t i) const
-{
-  // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
-  Int_t pdg = GetParentPDGCode(i);
-  Int_t id=pdg%100000;
-  return (!((id-id%10)%110));
-}
diff --git a/PWG3/muon/AliMCMuonTrack.h b/PWG3/muon/AliMCMuonTrack.h
deleted file mode 100644 (file)
index 6e2a97f..0000000
+++ /dev/null
@@ -1,83 +0,0 @@
-#ifndef ALIMCMUONTRACK_H
-#define ALIMCMUONTRACK_H
-
-#include <TParticle.h>
-#include <TLorentzVector.h>
-#include <TClonesArray.h>
-
-#include "AliMCEventHandler.h"
-#include "AliStack.h"
-#include "AliAODMCParticle.h"
-#include "AliESDEvent.h"
-#include "AliESDMuonTrack.h"
-#include "AliAODTrack.h"
-#include "AliAODMuonTrack.h"
-
-class AliMCMuonTrack : public AliAODMuonTrack {
- public:
-
-  AliMCMuonTrack();
-  AliMCMuonTrack(AliAODTrack *trkAOD, TClonesArray *mcClArr, Bool_t full=kFALSE);
-  AliMCMuonTrack(AliESDMuonTrack *trkESD, AliESDEvent *esd, AliMCEventHandler *mcH, Bool_t full=kFALSE);
-  AliMCMuonTrack(AliESDMuonTrack *trkESD, AliStack *stack);
-  virtual ~AliMCMuonTrack();
-
-  Int_t GetParentFlavour(Int_t i=0) const;
-  Bool_t IsMotherAResonance(Int_t i) const;
-  Bool_t IsAMuon() const { return (fSource>=0 && fSource!=4); }
-
-  TLorentzVector GetPGen() const { return fPGen; }
-  Int_t GetSource() const { return fSource; }
-  Int_t GetTrackIndex() const { return fTrackIndex; }
-  Int_t GetTrackPDGCode() const { return fTrackPDGCode; }
-  Int_t GetNParents() const { return fNParents; }
-  Int_t GetParentIndex(Int_t i=0) const { return (i<fNParents ? fParentIndex[i] : -1); }
-  Int_t GetParentPDGCode(Int_t i=0) const { return (i<fNParents ? fParentPDGCode[i] : 0); }
-  Int_t GetQuarkIndex(Int_t i=0) const { return (i<4 ? fQuarkIndex[i] : -1); }
-  Int_t GetQuarkPDGCode(Int_t i=0) const { return (i<4 ? fQuarkPDGCode[i] : 0); }
-  Bool_t IsOscillation() const { return fOscillation; }
-  Double_t GetWeight() const { return fWeight; }
-
- private:
-
-  AliAODMCParticle* FindTrackRef(AliAODTrack *trkAOD, TClonesArray *mcClArr);
-  TParticle* FindTrackRef(AliESDMuonTrack *trkESD, AliESDEvent *esd, AliMCEventHandler *mcH);  // do not implement on official train
-  void SetMCInfo(AliAODMCParticle *pMC, TClonesArray *mcClArr);
-  void SetMCInfo(TParticle *pMC, AliMCEventHandler *mcH);
-  void FillHistoryParents(Int_t lineM, TClonesArray *mcClArr);
-  void FillHistoryParents(Int_t lineM, AliStack *stack);
-  void FillHistoryQuarks(Int_t lineM, TClonesArray *mcClArr);
-  void FillHistoryQuarks(Int_t lineM, AliStack *stack);
-  //AliMUONTrack* CovESDtoMuonTrack(AliESDMuonTrack &trkESD);  // do not implement on official train
-  Int_t SelectHFMuon();
-
-  Bool_t IsDiquark(Int_t pdg);
-  void ResetQuarkInfo();
-
-  static const Double_t fgkSigmaCut;
-
-  Bool_t fIsFull;
-  static const Int_t fgkNParentsMax = 5;
-  TLorentzVector fPGen;
-  Int_t fTrackIndex;
-  Int_t fTrackPDGCode;
-  Int_t fSource;  // = 0, mu<-b 
-                  // = 1, mu<-c 
-                  // = 2, primary mu
-                  // = 3, secondary mu
-                  // = 4, not mu
-
-  Int_t fParentIndex[fgkNParentsMax];
-  Int_t fParentPDGCode[fgkNParentsMax];
-  Int_t fNParents;
-  Bool_t fOscillation;
-
-  Int_t fQuarkIndex[4];
-  Int_t fQuarkPDGCode[4];
-
-  Double_t fWeight;  // for PbPb collisoions
-
-  ClassDef(AliMCMuonTrack, 2);
-};
-
-#endif
diff --git a/PWG3/muon/AliMuonInfoStoreMC.cxx b/PWG3/muon/AliMuonInfoStoreMC.cxx
new file mode 100644 (file)
index 0000000..85b398c
--- /dev/null
@@ -0,0 +1,490 @@
+/**************************************************************************
+ * Copyright(c) 1998-2006, 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.                  *
+ **************************************************************************/
+
+/////////////////////////////////////////////////////////////
+//
+// class used to extract and store info of MC particle
+//
+// Author: X-M. Zhang, zhang@clermont.in2p3.fr
+//                     zhangxm@iopp.ccnu.edu.cn
+/////////////////////////////////////////////////////////////
+
+#include <TParticle.h>
+#include <TClonesArray.h>
+
+#include "AliMCEvent.h"
+#include "AliMCEventHandler.h"
+#include "AliStack.h"
+#include "AliAODMCParticle.h"
+#include "AliESDMuonTrack.h"
+#include "AliAODTrack.h"
+#include "AliMuonInfoStoreRD.h"
+#include "AliMuonInfoStoreMC.h"
+
+class AliESDEvent;
+
+ClassImp(AliMuonInfoStoreMC)
+
+const TString AliMuonInfoStoreMC::fgkStdBranchName("MuonMC");
+const Int_t   AliMuonInfoStoreMC::fgkNSources = 7;
+
+//-----------------------------------------------------------------------------
+AliMuonInfoStoreMC::AliMuonInfoStoreMC() :
+AliMuonInfoStoreRD(),
+fIsFull(kFALSE),
+fLorentzP(),
+fTrackIndex(-1),
+fTrackPDGCode(0),
+fSource(-1),
+fNParents(0),
+fOscillation(kFALSE),
+fWeight(0.)
+{
+  //
+  // default constructor
+  //
+  for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
+  for (Int_t i=4; i--;) { fQuarkIndex[i]  = -1; fQuarkPDGCode[i]  = 0; }
+}
+
+//-----------------------------------------------------------------------------
+AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliAODTrack *trkAOD, TClonesArray *mcClArr, Bool_t full) :
+AliMuonInfoStoreRD(trkAOD),
+fIsFull(full),
+fLorentzP(),
+fTrackIndex(-1),
+fTrackPDGCode(0),
+fSource(-1),
+fNParents(0),
+fOscillation(kFALSE),
+fWeight(0.)
+{
+  //
+  // default constructor
+  //
+  for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
+  for (Int_t i=4; i--;) { fQuarkIndex[i]  = -1; fQuarkPDGCode[i]  = 0; }
+
+  AliAODMCParticle *pMC = this->FindTrackRef(trkAOD, mcClArr);
+  if (pMC) this->SetMCInfo(pMC, mcClArr);
+}
+
+//-----------------------------------------------------------------------------
+AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliESDMuonTrack *trkESD, AliMCEventHandler *mcH, Bool_t full) :
+AliMuonInfoStoreRD(trkESD),
+fIsFull(full),
+fLorentzP(),
+fTrackIndex(-1),
+fTrackPDGCode(0),
+fSource(-1),
+fNParents(0),
+fOscillation(kFALSE),
+fWeight(0.)
+{
+  //
+  // default constructor
+  //
+  for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
+  for (Int_t i=4; i--;) { fQuarkIndex[i]  = -1; fQuarkPDGCode[i]  = 0; }
+
+  TParticle *pMC = this->FindTrackRef(trkESD, mcH);
+  if (pMC) this->SetMCInfo(pMC, mcH);
+}
+
+//-----------------------------------------------------------------------------
+/*AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliESDMuonTrack *trkESD, AliESDEvent *esd, AliMCEventHandler *mcH, Bool_t full) :
+AliMuonInfoStoreRD(trkESD),
+fIsFull(full),
+fLorentzP(),
+fTrackIndex(-1),
+fTrackPDGCode(0),
+fSource(-1),
+fNParents(0),
+fOscillation(kFALSE),
+fWeight(0.)
+{
+#include "AliMUONRecoCheck.h"
+#include "AliMUONVTrackStore.h"
+#include "AliMUONTrack.h"
+#include "AliMUONESDInterface.h"
+  //
+  // default constructor
+  //
+  for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
+  for (Int_t i=4; i--;) { fQuarkIndex[i]  = -1; fQuarkPDGCode[i]  = 0; }
+
+  TParticle *pMC = this->FindTrackRef(trkESD, esd, mcH);
+  if (pMC) this->SetMCInfo(pMC, mcH);
+}*/
+
+//-----------------------------------------------------------------------------
+AliMuonInfoStoreMC::AliMuonInfoStoreMC(const AliMuonInfoStoreMC &src) :
+AliMuonInfoStoreRD(src),
+fIsFull(src.fIsFull),
+fLorentzP(src.fLorentzP),
+fTrackIndex(src.fTrackIndex),
+fTrackPDGCode(src.fTrackPDGCode),
+fSource(src.fSource),
+fNParents(src.fNParents),
+fOscillation(src.fOscillation),
+fWeight(src.fWeight)
+{
+  //
+  // copy constructor
+  //
+  for (Int_t i=5; i--;) {
+    fParentIndex[i]   = src.fParentIndex[i];
+    fParentPDGCode[i] = src.fParentPDGCode[i];
+  }
+  for (Int_t i=4; i--;) {
+    fQuarkIndex[i]    = src.fQuarkIndex[i];
+    fQuarkPDGCode[i]  = src.fQuarkPDGCode[i];
+  }
+}
+
+//-----------------------------------------------------------------------------
+AliMuonInfoStoreMC& AliMuonInfoStoreMC::operator=(const AliMuonInfoStoreMC &src)
+{
+  //
+  // assignment constructor
+  //
+  if(&src==this) return *this;
+  AliMuonInfoStoreRD::operator=(src);
+
+  fIsFull       = src.fIsFull;
+  fLorentzP     = src.fLorentzP;
+  fTrackIndex   = src.fTrackIndex;
+  fTrackPDGCode = src.fTrackPDGCode;
+  fSource       = src.fSource;
+  fNParents     = src.fNParents;
+  fOscillation  = src.fOscillation;
+  fWeight       = src.fWeight;
+  for (Int_t i=5; i--;) {
+    fParentIndex[i]   = src.fParentIndex[i];
+    fParentPDGCode[i] = src.fParentPDGCode[i];
+  }
+  for (Int_t i=4; i--;) {
+    fQuarkIndex[i]    = src.fQuarkIndex[i];
+    fQuarkPDGCode[i]  = src.fQuarkPDGCode[i];
+  }
+
+  return *this;
+}
+
+//-----------------------------------------------------------------------------
+AliMuonInfoStoreMC::~AliMuonInfoStoreMC()
+{
+  //
+  // destructor
+  //
+}
+
+//-----------------------------------------------------------------------------
+AliAODMCParticle* AliMuonInfoStoreMC::FindTrackRef(AliAODTrack* const trkAOD, TClonesArray* const mcClArr)
+{
+  // find MC track ref with AOD base
+
+  AliAODMCParticle *pMC = 0;
+  fTrackIndex = trkAOD->GetLabel();
+  if (fTrackIndex>=0)
+    pMC = (AliAODMCParticle*)mcClArr->At(fTrackIndex);
+  return pMC;
+}
+
+//-----------------------------------------------------------------------------
+TParticle* AliMuonInfoStoreMC::FindTrackRef(AliESDMuonTrack* const trkESD, AliMCEventHandler* const mcH)
+{
+  // find MC track ref with ESD base
+
+  TParticle *pMCRef = 0;
+  fTrackIndex = trkESD->GetLabel();
+  if (fTrackIndex>=0) pMCRef = mcH->MCEvent()->Stack()->Particle(fTrackIndex);
+  return pMCRef;
+}
+
+//-----------------------------------------------------------------------------
+/*TParticle* AliMuonInfoStoreMC::FindTrackRef(AliESDMuonTrack* const trkESD, AliESDEvent* const esd, AliMCEventHandler* const mcH)
+{
+  // find MC track ref with ESD trackRef base
+
+  TParticle *pMCRef = 0;
+  AliMUONRecoCheck rc(esd,mcH);
+  AliMUONVTrackStore *trkRefArr = rc.TrackRefs(-1);
+
+  AliMUONTrack trkMuon;
+  AliMUONESDInterface::ESDToMUON(*trkESD, trkMuon, kFALSE);
+
+  Int_t nMatchClusters = 0;
+  AliMUONTrack *trkRef = rc.FindCompatibleTrack(trkMuon, *trkRefArr, nMatchClusters, kFALSE, 10.);
+  if (trkRef) fTrackIndex = trkRef->GetUniqueID();
+  if (fTrackIndex>=0) pMCRef = mcH->MCEvent()->Stack()->Particle(fTrackIndex);
+  return pMCRef;
+}*/
+
+//-----------------------------------------------------------------------------
+void AliMuonInfoStoreMC::SetMCInfo(AliAODMCParticle *pMC, TClonesArray *mcClArr)
+{
+  // fill track MC info with AOD base
+
+  fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->E());
+  fTrackPDGCode = pMC->GetPdgCode();
+  if (TMath::Abs(fTrackPDGCode)!=13) {
+    fSource = 4;
+    return;
+  } 
+
+  Int_t lineM = pMC->GetMother();
+  if (lineM<0) {
+    fSource = 2;
+    return;
+  }
+
+  Bool_t isPrimary = ((AliAODMCParticle*)mcClArr->At(lineM))->IsPrimary();
+  if (!isPrimary) {
+    fSource = 3;
+    return;
+  }
+
+  this->FillHistoryParents(lineM, mcClArr);
+  fSource = this->SelectHFMuon();
+  return;
+}
+
+//-----------------------------------------------------------------------------
+void AliMuonInfoStoreMC::SetMCInfo(TParticle *pMC, AliMCEventHandler* const mcH)
+{
+  // fill track MC info with ESD base
+
+  fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->Energy());
+  fTrackPDGCode = pMC->GetPdgCode();
+  if (TMath::Abs(fTrackPDGCode)!=13) {
+    fSource = 4;
+    return;
+  }
+
+  Int_t lineM = pMC->GetFirstMother();
+  if (lineM<0) {
+    fSource = 2;
+    return;
+  }
+
+  AliStack *stack = mcH->MCEvent()->Stack();
+  if (lineM>=stack->GetNprimary()) {
+    fSource = 3;
+    return;
+  }
+
+  this->FillHistoryParents(lineM, stack);
+  fSource = this->SelectHFMuon();
+  return;
+}
+
+//-----------------------------------------------------------------------------
+void AliMuonInfoStoreMC::FillHistoryParents(Int_t lineM, TClonesArray *mcClArr)
+{
+  // find track hadron parents with AOD base
+
+  Int_t countP=-1, pdg=0;
+  Int_t parents[10], parLine[10];
+  AliAODMCParticle *mother = 0;
+  while(lineM>=0){
+    mother = (AliAODMCParticle*)mcClArr->At(lineM);
+    pdg = mother->GetPdgCode();
+    if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
+    parents[++countP] = pdg;
+    parLine[countP] = lineM;
+    lineM = mother->GetMother();
+  }
+  for(Int_t i=0; i<=countP; i++){
+    fParentIndex[i] = parLine[countP-i];
+    fParentPDGCode[i] = parents[countP-i];
+  }
+  fNParents = countP + 1;
+
+  if (fIsFull && lineM>=0) this->FillHistoryQuarks(lineM, mcClArr);
+  return;
+}
+
+//-----------------------------------------------------------------------------
+void AliMuonInfoStoreMC::FillHistoryParents(Int_t lineM, AliStack *stack)
+{
+  // find track hadron parents with ESD base
+
+  Int_t countP=-1, pdg=0;
+  Int_t parents[10], parLine[10];
+  TParticle *mother = 0;
+  while(lineM>=0){
+    mother = stack->Particle(lineM);
+    pdg = mother->GetPdgCode();
+    if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
+    parents[++countP] = pdg;
+    parLine[countP] = lineM;
+    lineM = mother->GetFirstMother();
+  }
+  for(Int_t i=0; i<=countP; i++){
+    fParentIndex[i] = parLine[countP-i];
+    fParentPDGCode[i] = parents[countP-i];
+  }
+  fNParents = countP + 1;
+
+  if (fIsFull && lineM>=0) this->FillHistoryQuarks(lineM, stack);
+  return;
+}
+
+//-----------------------------------------------------------------------------
+void AliMuonInfoStoreMC::FillHistoryQuarks(Int_t lineM, TClonesArray* const mcClArr)
+{
+  // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx 
+
+  if (lineM<0) return;
+  Int_t countP=-1, pdg=0;
+  AliAODMCParticle *mother = 0;
+  while(lineM>=0){
+    mother = (AliAODMCParticle*)mcClArr->At(lineM);
+    pdg = mother->GetPdgCode();
+    fQuarkIndex[++countP] = lineM;
+    fQuarkPDGCode[countP] = pdg;
+    lineM = mother->GetMother();
+  }
+
+  // for PYTHIA checking
+  countP = 1;
+  for(Int_t par=0; par<4; par++) {
+    if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
+  }
+  if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
+    if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
+      AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
+                 this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));
+
+      pdg = this->QuarkPDGCode(countP);
+      Int_t line = this->QuarkIndex(countP);
+      this->ResetQuarkInfo();
+      while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
+        pdg = ((AliAODMCParticle*)mcClArr->At(++line))->GetPdgCode();
+      }
+      while(line>=0){
+        mother = (AliAODMCParticle*)mcClArr->At(line);
+        pdg = mother->GetPdgCode();
+        fQuarkIndex[countP] = line;
+        fQuarkPDGCode[countP++] = pdg;
+        line = mother->GetMother();
+      }
+    }
+  }
+  return;
+}
+
+//-----------------------------------------------------------------------------
+void AliMuonInfoStoreMC::FillHistoryQuarks(Int_t lineM, AliStack* const stack)
+{
+  // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx 
+
+  if (lineM<0) return;
+  Int_t countP=-1, pdg=0;
+  TParticle *mother = 0;
+  while(lineM>=0){
+    mother = stack->Particle(lineM);
+    pdg = mother->GetPdgCode();
+    fQuarkIndex[++countP] = lineM;
+    fQuarkPDGCode[countP] = pdg;
+    lineM = mother->GetFirstMother();
+  }
+
+  // for PYTHIA checking
+  countP = 1;
+  for(Int_t par=0; par<4; par++) {
+    if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
+  }
+  if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
+    if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
+      AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
+                 this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));
+
+      pdg = this->QuarkPDGCode(countP);
+      Int_t line = this->QuarkIndex(countP);
+      this->ResetQuarkInfo();
+      while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
+        pdg = stack->Particle(++line)->GetPdgCode();
+      }
+      while(line>=0){
+        mother = stack->Particle(line);
+        pdg = mother->GetPdgCode();
+        fQuarkIndex[countP] = line;
+        fQuarkPDGCode[countP++] = pdg;
+        line = mother->GetFirstMother();
+      }
+    }
+  }
+  return;
+}
+
+//-----------------------------------------------------------------------------
+Int_t AliMuonInfoStoreMC::SelectHFMuon()
+{
+  // set info of muon from HF
+
+  Int_t flv = ParentFlavour(0);
+  if (flv!=4 && flv!=5) return 2;
+
+  Bool_t isRes = kFALSE;
+  Int_t i=0, nparents=this->NParents();
+  while (i<nparents && !isRes) {
+    isRes = IsMotherAResonance(i++);
+  }
+
+  if (isRes) return 2;
+  if (flv==5) return 0;
+  else return 1;
+}
+
+//-----------------------------------------------------------------------------
+Bool_t AliMuonInfoStoreMC::IsDiquark(Int_t pdg)
+{
+  // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
+  pdg = TMath::Abs(pdg);
+  if(pdg>1000 && (pdg%100)<10) return kTRUE;
+  else return kFALSE;
+}
+
+//-----------------------------------------------------------------------------
+void AliMuonInfoStoreMC::ResetQuarkInfo()
+{
+  // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
+  for(Int_t pos=1; pos<4; pos++) {
+    fQuarkIndex[pos] = -1;
+    fQuarkPDGCode[pos] = 0;
+  }
+  return;
+}
+
+//-----------------------------------------------------------------------------
+Int_t AliMuonInfoStoreMC::ParentFlavour(Int_t i) const
+{
+  // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
+  Int_t pdg = ParentPDGCode(i);
+  pdg = TMath::Abs(pdg/100);
+  if(pdg>9) pdg /= 10;
+  return pdg;
+}
+
+//-----------------------------------------------------------------------------
+Bool_t AliMuonInfoStoreMC::IsMotherAResonance(Int_t i) const
+{
+  // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
+  Int_t pdg = ParentPDGCode(i);
+  Int_t id=pdg%100000;
+  return (!((id-id%10)%110));
+}
diff --git a/PWG3/muon/AliMuonInfoStoreMC.h b/PWG3/muon/AliMuonInfoStoreMC.h
new file mode 100644 (file)
index 0000000..b8a2ddd
--- /dev/null
@@ -0,0 +1,101 @@
+#ifndef ALIMUONINFOSTOREMC_H
+#define ALIMUONINFOSTOREMC_H
+
+/* Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//***********************************************************
+// Class AliMuonInfoStoreRD
+// class used to extract and store info of MC particle
+// Author: X-M. Zhang, zhang@clermont.in2p3.fr
+//                     zhangxm@iopp.ccnu.edu.cn
+//***********************************************************
+
+#include <TString.h>
+#include <TParticle.h>
+#include <TLorentzVector.h>
+#include <TClonesArray.h>
+
+#include "AliMCEventHandler.h"
+#include "AliStack.h"
+#include "AliAODMCParticle.h"
+#include "AliESDEvent.h"
+#include "AliESDMuonTrack.h"
+#include "AliAODTrack.h"
+#include "AliMuonInfoStoreRD.h"
+
+class AliMuonInfoStoreMC : public AliMuonInfoStoreRD {
+ public:
+
+  AliMuonInfoStoreMC();
+  AliMuonInfoStoreMC(AliAODTrack *trkAOD, TClonesArray *mcClArr, Bool_t full=kFALSE);
+  AliMuonInfoStoreMC(AliESDMuonTrack *trkESD, AliESDEvent *esd, AliMCEventHandler *mcH, Bool_t full=kFALSE);
+  AliMuonInfoStoreMC(AliESDMuonTrack *trkESD, AliMCEventHandler *mcH, Bool_t full=kFALSE);
+  AliMuonInfoStoreMC(AliESDMuonTrack *trkESD, AliStack *stack);
+  AliMuonInfoStoreMC(const AliMuonInfoStoreMC &src);
+  AliMuonInfoStoreMC& operator=(const AliMuonInfoStoreMC &src);
+  virtual ~AliMuonInfoStoreMC();
+
+  Int_t ParentFlavour(Int_t i=0) const;
+  Bool_t IsMotherAResonance(Int_t i) const;
+  Bool_t IsAMuon() const { return (fSource>=0 && fSource!=4); }
+
+  TLorentzVector LorentzP()      const { return fLorentzP; }
+  Int_t MuonSource()             const { return fSource; }
+  Int_t TrackIndex()             const { return fTrackIndex; }
+  Int_t TrackPDGCode()           const { return fTrackPDGCode; }
+  Int_t NParents()               const { return fNParents; }
+  Int_t ParentIndex(Int_t i=0)   const { return (i<fNParents ? fParentIndex[i] : -1); }
+  Int_t ParentPDGCode(Int_t i=0) const { return (i<fNParents ? fParentPDGCode[i] : 0); }
+  Int_t QuarkIndex(Int_t i=0)    const { return (i<4 ? fQuarkIndex[i] : -1); }
+  Int_t QuarkPDGCode(Int_t i=0)  const { return (i<4 ? fQuarkPDGCode[i] : 0); }
+  Bool_t IsOscillation()         const { return fOscillation; }
+  Double_t Weight()              const { return fWeight; }
+
+  static const char* StdBranchName() { return fgkStdBranchName.Data(); }
+  static const Int_t NSources()      { return fgkNSources;             }
+
+ private:
+
+  AliAODMCParticle* FindTrackRef(AliAODTrack*     const trkAOD, TClonesArray* const mcClArr);
+  TParticle*        FindTrackRef(AliESDMuonTrack* const trkESD, AliMCEventHandler* const mcH);
+  TParticle*        FindTrackRef(AliESDMuonTrack* const trkESD, AliESDEvent* const esd, AliMCEventHandler* const mcH);
+  void SetMCInfo(AliAODMCParticle *pMC, TClonesArray* const mcClArr);
+  void SetMCInfo(TParticle *pMC, AliMCEventHandler *mcH);
+  void FillHistoryParents(Int_t lineM, TClonesArray* const mcClArr);
+  void FillHistoryParents(Int_t lineM, AliStack* const stack);
+  void FillHistoryQuarks(Int_t lineM, TClonesArray *mcClArr);
+  void FillHistoryQuarks(Int_t lineM, AliStack *stack);
+  Int_t SelectHFMuon();
+
+  Bool_t IsDiquark(Int_t pdg);
+  void ResetQuarkInfo();
+
+  static const TString fgkStdBranchName;  // Standard branch name
+  static const Int_t   fgkNSources;       // num. of muon sources
+
+  Bool_t fIsFull;            // whether to use full mode (Pb-Pb)
+  TLorentzVector fLorentzP;  // lorentz momentum of particle
+  Int_t fTrackIndex;         // index of the MC particle
+  Int_t fTrackPDGCode;       // PDG code of the MC particle
+  Int_t fSource;  // = 0, mu<-b 
+                  // = 1, mu<-c 
+                  // = 2, primary mu
+                  // = 3, secondary mu
+                  // = 4, not mu
+                  // = 5, unidentified track
+
+  Int_t fParentIndex[5];    // index of parents
+  Int_t fParentPDGCode[5];  // PDG code of parents
+  Int_t fNParents;          // num. of parents
+  Bool_t fOscillation;      // flag of oscillation
+
+  Int_t fQuarkIndex[4];    // index of quarks
+  Int_t fQuarkPDGCode[4];  // PDG code of quarks
+
+  Double_t fWeight;  // for PbPb collisoions
+
+  ClassDef(AliMuonInfoStoreMC, 3);
+};
+
+#endif
diff --git a/PWG3/muon/AliMuonInfoStoreRD.cxx b/PWG3/muon/AliMuonInfoStoreRD.cxx
new file mode 100644 (file)
index 0000000..04f7e9d
--- /dev/null
@@ -0,0 +1,218 @@
+/**************************************************************************
+ * Copyright(c) 1998-2006, 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.                  *
+ **************************************************************************/
+
+/////////////////////////////////////////////////////////////
+//
+// class used to extract and store reco info of muon track
+//
+// Author: X-M. Zhang, zhang@clermont.in2p3.fr
+//                     zhangxm@iopp.ccnu.edu.cn
+/////////////////////////////////////////////////////////////
+
+#include "AliAODTrack.h"
+#include "AliESDMuonTrack.h"
+#include "AliMuonInfoStoreRD.h"
+
+class TObject;
+
+ClassImp(AliMuonInfoStoreRD)
+
+const TString AliMuonInfoStoreRD::fgkStdBranchName("MuonRD");
+Double_t      AliMuonInfoStoreRD::fgCuts[10] = {-999999., 999999.,
+                                                -999999., 999999.,
+                                                -999999., 999999.,
+                                                -999999., 999999.,
+                                                -999999., 999999.};
+
+//-----------------------------------------------------------------------------
+AliMuonInfoStoreRD::AliMuonInfoStoreRD() :
+TObject(),
+fMomentum(),
+fMomentumAtDCA(),
+fCharge(0),
+fMatchTrigger(-1),
+fNClusters(0),
+fMUONClusterMap(0),
+fChi2FitMomentum(0.),
+fChi2MatchTrigger(0.)
+{
+  //
+  // default constructor
+  //
+  for (Int_t i=3; i--;) fDCA[i] = 0.;
+}
+
+//-----------------------------------------------------------------------------
+AliMuonInfoStoreRD::AliMuonInfoStoreRD(AliAODTrack *trk) :
+TObject(),
+fMomentum(),
+fMomentumAtDCA(),
+fCharge(0),
+fMatchTrigger(-1),
+fNClusters(0),
+fMUONClusterMap(0),
+fChi2FitMomentum(0.),
+fChi2MatchTrigger(0.)
+{
+  //
+  // AOD-base constructor
+  //
+  for (Int_t i=3; i--;) fDCA[i] = 0.;
+  this->FillMuonInfo(trk);
+}
+
+//-----------------------------------------------------------------------------
+AliMuonInfoStoreRD::AliMuonInfoStoreRD(AliESDMuonTrack *trk) :
+TObject(),
+fMomentum(),
+fMomentumAtDCA(),
+fCharge(0),
+fMatchTrigger(-1),
+fNClusters(0),
+fMUONClusterMap(0),
+fChi2FitMomentum(0.),
+fChi2MatchTrigger(0.)
+{
+  //
+  // ESD-base constructor
+  //
+  for (Int_t i=3; i--;) fDCA[i] = 0.;
+  this->FillMuonInfo(trk);
+}
+
+//-----------------------------------------------------------------------------
+AliMuonInfoStoreRD::AliMuonInfoStoreRD(const AliMuonInfoStoreRD &src) :
+TObject(src),
+fMomentum(src.fMomentum),
+fMomentumAtDCA(src.fMomentumAtDCA),
+fCharge(src.fCharge),
+fMatchTrigger(src.fMatchTrigger),
+fNClusters(src.fNClusters),
+fMUONClusterMap(src.fMUONClusterMap),
+fChi2FitMomentum(src.fChi2FitMomentum),
+fChi2MatchTrigger(src.fChi2MatchTrigger)
+{
+  //
+  // copy constructor
+  //
+  for (Int_t i=3; i--;) fDCA[i] = src.fDCA[i];
+}
+
+//-----------------------------------------------------------------------------
+AliMuonInfoStoreRD& AliMuonInfoStoreRD::operator=(const AliMuonInfoStoreRD &src)
+{
+  //
+  // assignment constructor
+  //
+  if(&src==this) return *this;
+
+  fMomentum         = src.fMomentum;
+  fMomentumAtDCA    = src.fMomentumAtDCA;
+
+  fCharge           = src.fCharge;
+  fMatchTrigger     = src.fMatchTrigger;
+  fNClusters        = src.fNClusters;
+  fMUONClusterMap   = src.fMUONClusterMap;
+  fChi2FitMomentum  = src.fChi2FitMomentum;
+  fChi2MatchTrigger = src.fChi2MatchTrigger;
+
+  for (Int_t i=3; i--;) fDCA[i] = src.fDCA[i];
+
+  return *this;
+}
+
+//-----------------------------------------------------------------------------
+AliMuonInfoStoreRD::~AliMuonInfoStoreRD()
+{
+  //
+  // destructor
+  //
+}
+
+//-----------------------------------------------------------------------------
+void AliMuonInfoStoreRD::FillMuonInfo(AliAODTrack *trk)
+{
+  // extract reco info of muon track from AOD
+
+  Double_t arr[3];
+  trk->PxPyPz(arr);
+  this->SetMomentum(arr);
+
+  trk->PxPyPzAtDCA(arr);
+  this->SetMomentumAtDCA(arr);
+
+  trk->XYZAtDCA(arr);
+  this->SetDCA(arr);
+
+  this->SetCharge(trk->Charge());
+  this->SetMatchTrigger(trk->GetMatchTrigger());
+  this->SetMUONClusterMap(trk->GetMUONClusterMap());
+  this->SetChi2FitMomentum(trk->Chi2perNDF());
+  this->SetChi2MatchTrigger(trk->GetChi2MatchTrigger());
+
+  return;
+}
+
+
+//-----------------------------------------------------------------------------
+void AliMuonInfoStoreRD::FillMuonInfo(AliESDMuonTrack *trk)
+{
+  // extract reco info of muon track from ESD
+  // tract params before related to vertex are extracted
+
+  Double_t arr[3];
+  trk->PxPyPz(arr);
+  this->SetMomentum(arr);
+
+  trk->PxPyPzAtDCA(arr);
+  this->SetMomentumAtDCA(arr);
+
+  arr[0] = trk->GetNonBendingCoorAtDCA();
+  arr[1] = trk->GetBendingCoorAtDCA();
+  arr[2] = trk->GetZ();
+  this->SetDCA(arr);
+
+  this->SetCharge(trk->Charge());
+  this->SetMatchTrigger(trk->GetMatchTrigger());
+  this->SetMUONClusterMap(trk->GetMuonClusterMap());
+  this->SetChi2FitMomentum(trk->GetChi2()/(2.*trk->GetNHit()-5.));
+  this->SetChi2MatchTrigger(trk->GetChi2MatchTrigger());
+
+  this->SetNClusters(trk->GetNClusters());
+  return;
+}
+
+//-----------------------------------------------------------------------------
+Bool_t AliMuonInfoStoreRD::MuonSelection()
+{
+  // select muon tracks according to the selection cuts
+
+  Double_t p = Momentum().Mag();
+  if (p<AliMuonInfoStoreRD::fgCuts[0] || p>AliMuonInfoStoreRD::fgCuts[1])             return kFALSE;
+
+  Double_t pt = Momentum().Pt();
+  if (pt<AliMuonInfoStoreRD::fgCuts[2] || pt>AliMuonInfoStoreRD::fgCuts[3])           return kFALSE;
+
+  Double_t eta = Momentum().Eta();
+  if (eta<AliMuonInfoStoreRD::fgCuts[4] || eta>AliMuonInfoStoreRD::fgCuts[5])         return kFALSE;
+
+  Double_t dca = this->DCA();
+  if (dca<AliMuonInfoStoreRD::fgCuts[6] || dca>AliMuonInfoStoreRD::fgCuts[7])         return kFALSE;
+
+  Int_t trigger = this->MatchTrigger();
+  if (trigger<AliMuonInfoStoreRD::fgCuts[8] || trigger>AliMuonInfoStoreRD::fgCuts[9]) return kFALSE;
+
+  return kTRUE;
+}
diff --git a/PWG3/muon/AliMuonInfoStoreRD.h b/PWG3/muon/AliMuonInfoStoreRD.h
new file mode 100644 (file)
index 0000000..365b15c
--- /dev/null
@@ -0,0 +1,92 @@
+#ifndef ALIMUONINFOSTORERD_H
+#define ALIMUONINFOSTORERD_H
+
+/* Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//***********************************************************
+// Class AliMuonInfoStoreRD
+// class used to extract and store reco info of muon track
+// Author: X-M. Zhang, zhang@clermont.in2p3.fr
+//                     zhangxm@iopp.ccnu.edu.cn
+//***********************************************************
+
+#include <TObject.h>
+#include <TVector3.h>
+#include <TString.h>
+
+class AliAODTrack;
+class AliESDMuonTrack;
+
+class AliMuonInfoStoreRD : public TObject {
+ public:
+
+  AliMuonInfoStoreRD();
+  AliMuonInfoStoreRD(AliAODTrack *trk);
+  AliMuonInfoStoreRD(AliESDMuonTrack *trk);
+  AliMuonInfoStoreRD(const AliMuonInfoStoreRD &src);
+  AliMuonInfoStoreRD& operator=(const AliMuonInfoStoreRD &src);
+  virtual ~AliMuonInfoStoreRD();
+
+  TVector3 Momentum()      const { return fMomentum; }
+  TVector3 MomentumAtDCA() const { return fMomentumAtDCA; }
+
+  void XYZAtDCA(Double_t dca[3]) const { for (Int_t i=3; i--;) dca[i]=fDCA[i]; }
+  Double_t DCA() const  { return TMath::Sqrt(fDCA[0]*fDCA[0]+fDCA[1]*fDCA[1]); }
+
+  Short_t  Charge()           const { return fCharge;           }
+  Int_t    MatchTrigger()     const { return fMatchTrigger;     }
+  Int_t    NClusters()        const { return fNClusters;        }
+  UInt_t   MUONClusterMap()   const { return fMUONClusterMap;   }
+  Double_t Chi2FitMomentum()  const { return fChi2FitMomentum;  }
+  Double_t Chi2MatchTrigger() const { return fChi2MatchTrigger; }
+
+  Bool_t MuonSelection();
+
+  static const char* StdBranchName()                  { return fgkStdBranchName.Data(); }
+  static const void SelectionCust(Double_t cuts[10])  { for (Int_t i=10; i--;) cuts[i]=fgCuts[i]; }
+  static void SetSelectionCuts(Double_t cuts[10]) { for (Int_t i=10; i--;) fgCuts[i]=cuts[i]; }
+
+ private:
+
+  void FillMuonInfo(AliAODTrack *trk);
+  void FillMuonInfo(AliESDMuonTrack *trk);
+
+  void SetMomentum(Double_t p[3])      { fMomentum.SetXYZ(p[0],p[1],p[2]);      }
+  void SetMomentumAtDCA(Double_t p[3]) { fMomentumAtDCA.SetXYZ(p[0],p[1],p[2]); }
+
+  void SetDCA(Double_t dca[3]) { for (Int_t i=3; i--;) fDCA[i]=dca[i]; }
+  void SetCharge(Short_t charge)           { fCharge        = charge;  }
+  void SetNClusters(Int_t ncls)            { fNClusters     = ncls;    }
+  void SetMatchTrigger(Int_t trigger)      { fMatchTrigger  = trigger; }
+  void SetMUONClusterMap(UInt_t clMap)     { fMUONClusterMap  = clMap; }
+  void SetChi2FitMomentum(Double_t chi2)   { fChi2FitMomentum  = chi2; }
+  void SetChi2MatchTrigger(Double_t chi2)  { fChi2MatchTrigger = chi2; }
+
+  static const TString fgkStdBranchName;  // Standard branch name
+  static Double_t fgCuts[10];  // 0, min of 3-momentum
+                               // 1, max of 3-momentum
+                               // 2, pt_Min
+                               // 3, pt_Max
+                               // 4, eta_Min
+                               // 5, eta_Max
+                               // 6, dca_Min
+                               // 7, dca_Max
+                               // 8, about trigger matching
+                               // 9, about trigger matching
+
+  TVector3 fMomentum;       // momentum corrected w vtx
+  TVector3 fMomentumAtDCA;  // momentum at DCA in vtx plane
+
+  Double_t fDCA[3];            // distance of closet approach
+  Short_t  fCharge;            // track charge
+  Int_t    fMatchTrigger;      // type of match trigger
+  Int_t    fNClusters;         // number of clusters in the track
+  UInt_t   fMUONClusterMap;    // map of MUON clusters
+  Double_t fChi2FitMomentum;   // chi2/NDF of momentum fit
+  Double_t fChi2MatchTrigger;  // chi2 of trigger matching
+
+  ClassDef(AliMuonInfoStoreRD, 4);
+};
+
+#endif
index df87c52..515f81b 100644 (file)
-#include <TNamed.h>
+/**************************************************************************
+ * Copyright(c) 1998-2006, 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.                  *
+ **************************************************************************/
+
+/////////////////////////////////////////////////////////////
+//
+// class used to extract and store info at event level
+//
+// Author: X-M. Zhang, zhang@clermont.in2p3.fr
+//                     zhangxm@iopp.ccnu.edu.cn
+/////////////////////////////////////////////////////////////
+
 #include <TMath.h>
+#include <TH1.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TList.h>
+#include <TObjArray.h>
+#include <TObjString.h>
 
-#include "AliVHeader.h"
-#include "AliVVertex.h"
+#include "AliTriggerAnalysis.h"
+#include "AliBackgroundSelection.h"
+#include "AliAODVertex.h"
+#include "AliAODEvent.h"
+#include "AliESDEvent.h"
+#include "AliMultiplicity.h"
+#include "AliMuonInfoStoreRD.h"
+#include "AliMuonInfoStoreMC.h"
+#include "AliDimuInfoStoreRD.h"
+#include "AliDimuInfoStoreMC.h"
 #include "AliMuonsHFHeader.h"
 
+class TNamed;
+class AliESDVertex;
+
 ClassImp(AliMuonsHFHeader)
 
+const TString AliMuonsHFHeader::fgkStdBranchName("MuEvsH");
+Int_t         AliMuonsHFHeader::fgAnaMode         = 0;
+Bool_t        AliMuonsHFHeader::fgIsMC            = kFALSE;
+Bool_t        AliMuonsHFHeader::fgIsEventSelected = kFALSE;
+Double_t      AliMuonsHFHeader::fgCuts[3] = { -999999., 999999., 999999.};
+
 //_____________________________________________________________________________
 AliMuonsHFHeader::AliMuonsHFHeader() :
 TNamed(),
 fTriggerMask(0),
+fFiredTrigger(0),
+fNFiredTrigger(0),
+fIsPhysicsTriggered(kFALSE),
+fIsPhysicsAccepted(kFALSE),
+fEventType(0),
+fUnrecoVertex(kFALSE),
 fNContributors(0),
-fMultMuon(0),
-fMultDimuon(0),
-fCentrality(0.),
-fUnrecoVertex(kFALSE)
+fUnrecoVtxSPD(kFALSE),
+fNContributorsSPD(0),
+fNTrackletsSPD(0),
+fCentrality(0.)
 {
   //
   // default constructor
   //
-  for (Int_t i=0; i<3; i++) fPosition[i]=0.;
+  for (Int_t i=3; i--;) fVtx[i] = 0.;
+  for (Int_t i=3; i--;) fVtxSPD[i] = 0.;
+}
+
+//_____________________________________________________________________________
+AliMuonsHFHeader::AliMuonsHFHeader(const AliMuonsHFHeader &src) :
+TNamed(),
+fTriggerMask(src.fTriggerMask),
+fFiredTrigger(src.fFiredTrigger),
+fNFiredTrigger(src.fNFiredTrigger),
+fIsPhysicsTriggered(src.fIsPhysicsTriggered),
+fIsPhysicsAccepted(src.fIsPhysicsAccepted),
+fEventType(src.fEventType),
+fUnrecoVertex(src.fUnrecoVertex),
+fNContributors(src.fNContributors),
+fUnrecoVtxSPD(src.fUnrecoVtxSPD),
+fNContributorsSPD(src.fNContributorsSPD),
+fNTrackletsSPD(src.fNTrackletsSPD),
+fCentrality(src.fCentrality)
+{
+  //
+  // copy constructor
+  //
+  for (Int_t i=3; i--;) fVtx[i]    = src.fVtx[i];
+  for (Int_t i=3; i--;) fVtxSPD[i] = src.fVtxSPD[i];
+}
+
+//_____________________________________________________________________________
+AliMuonsHFHeader& AliMuonsHFHeader::operator=(const AliMuonsHFHeader &src)
+{
+  //
+  // assignment constructor
+  //
+  fTriggerMask        = src.fTriggerMask;
+  fFiredTrigger       = src.fFiredTrigger;
+  fNFiredTrigger      = src.fNFiredTrigger; 
+  fIsPhysicsTriggered = src.fIsPhysicsTriggered;
+  fIsPhysicsAccepted  = src.fIsPhysicsAccepted;
+  fEventType          = src.fEventType;
+  fUnrecoVertex       = src.fUnrecoVertex;
+  fNContributors      = src.fNContributors;
+  fUnrecoVtxSPD       = src.fUnrecoVtxSPD;
+  fNContributorsSPD   = src.fNContributorsSPD;
+  fNTrackletsSPD      = src.fNTrackletsSPD;
+  fCentrality         = src.fCentrality;
+
+  for (Int_t i=3; i--;) fVtx[i]      = src.fVtx[i];
+  for (Int_t i=3; i--;) fVtxSPD[i]   = src.fVtxSPD[i];
+
+  return *this;
 }
 
 //_____________________________________________________________________________
@@ -32,18 +132,341 @@ AliMuonsHFHeader::~AliMuonsHFHeader()
 }
 
 //_____________________________________________________________________________
-void AliMuonsHFHeader::SetHeader(AliVHeader *header)
+void AliMuonsHFHeader::SetEvent(AliAODEvent *event)
 {
-  fTriggerMask = header->GetTriggerMask();
+  // extract event info from AOD event
+
+  fTriggerMask = event->GetTriggerMask();
+
+  AliAODVertex *vertex = event->GetPrimaryVertex(); 
+  vertex->GetXYZ(fVtx);
+  fNContributors = vertex->GetNContributors();
+  fUnrecoVertex = (TMath::Abs(fVtx[0])<1e-6 && TMath::Abs(fVtx[1])<1e-6 &&
+                   TMath::Abs(fVtx[2])<1e-6);
+  this->SetTitle(vertex->GetTitle());
+
+  this->SetFiredTrigger(event->GetFiredTriggerClasses());
+  this->EventSelection();
   return;
 }
 
 //_____________________________________________________________________________
-void AliMuonsHFHeader::SetVertex(AliVVertex *vertex)
+void AliMuonsHFHeader::SetEvent(AliESDEvent *event)
 {
-  vertex->GetXYZ(fPosition);
+  // extract event info from ESD event
+  // SPD vertex and Physics event selection are implimented
+
+  fTriggerMask = event->GetTriggerMask();
+
+  const AliESDVertex *vertex = event->GetPrimaryVertex(); 
+  vertex->GetXYZ(fVtx);
   fNContributors = vertex->GetNContributors();
-  fUnrecoVertex = (TMath::Abs(fPosition[0])<1e-6 && TMath::Abs(fPosition[1])<1e-6 &&
-                   TMath::Abs(fPosition[2])<1e-6);
+  fUnrecoVertex = (TMath::Abs(fVtx[0])<1e-6 && TMath::Abs(fVtx[1])<1e-6 &&
+                   TMath::Abs(fVtx[2])<1e-6);
+  this->SetTitle(vertex->GetTitle());
+
+  const AliESDVertex *vtxSPD = event->GetPrimaryVertexSPD();
+  vtxSPD->GetXYZ(fVtxSPD);
+  fNContributorsSPD = vtxSPD->GetNContributors();
+  fUnrecoVtxSPD = (TMath::Abs(fVtxSPD[0])<1e-6 && TMath::Abs(fVtxSPD[1])<1e-6 &&
+                   TMath::Abs(fVtxSPD[2])<1e-6);
+  fNTrackletsSPD = event->GetMultiplicity()->GetNumberOfTracklets();
+
+  this->SetFiredTrigger(event->GetFiredTriggerClasses());
+  this->PhysicsTriggerAna(event);
+  this->EventSelection();
+  return;
+}
+
+//_____________________________________________________________________________
+void AliMuonsHFHeader::PhysicsTriggerAna(const AliESDEvent *esd)
+{
+  // ESD event trigger condition analysis
+  // according to the method in $ALICE_ROOT/ANALYSIS/AliPhysicsSelection.cxx
+
+  fEventType = esd->GetHeader()->GetEventType();
+  fIsPhysicsTriggered = kFALSE;
+  fIsPhysicsAccepted  = kFALSE;
+
+  AliTriggerAnalysis *triggerAna = new AliTriggerAnalysis();
+  triggerAna->SetAnalyzeMC(fgIsMC);
+  triggerAna->SetSPDGFOThreshhold(1);
+
+  Int_t  triggerHW  = triggerAna->SPDFiredChips(esd, 1);
+  Bool_t isFiredv0A = triggerAna->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0A);
+  Bool_t isFiredv0C = triggerAna->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0C);
+  if (triggerHW==0 && !isFiredv0A && !isFiredv0C) {
+    delete triggerAna;
+    triggerAna = 0;
+    return;
+  }
+
+  Bool_t triggerBG = (triggerAna->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0ABG) ||
+                      triggerAna->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0CBG));
+  Int_t  isFiredSPD = triggerAna->SPDFiredChips(esd, 0);
+  Bool_t triggerFD = ((isFiredSPD>1) || (isFiredSPD>0 && (isFiredv0A || isFiredv0C)) || (isFiredv0A && isFiredv0C));
+  if ((!triggerBG) && triggerFD) fIsPhysicsTriggered = kTRUE; 
+  delete triggerAna;
+  triggerAna = 0;
+
+  if (fIsPhysicsTriggered) {
+    AliBackgroundSelection *bkgId = new AliBackgroundSelection();
+    fIsPhysicsAccepted = bkgId->IsSelected(const_cast<AliESDEvent*>(esd));
+    delete bkgId; bkgId=0;
+  }
+  return;
+}
+
+//_____________________________________________________________________________
+void AliMuonsHFHeader::SetFiredTrigger(TString trigger)
+{
+  // get info of fired trigger classes of event
+
+  fFiredTrigger = trigger;
+  TObjArray *tokens = trigger.Tokenize(" ");
+  fNFiredTrigger = tokens->GetEntries();
+  return;
+}
+
+//_____________________________________________________________________________
+Bool_t AliMuonsHFHeader::IsTriggerFired(TString trigger)
+{
+  // check whether the trigger class "trigger" is fired in this event
+
+  if (fNFiredTrigger<=0) return kFALSE;
+
+  TObjArray *tokens = fFiredTrigger.Tokenize(" "); 
+  for (Int_t i=fNFiredTrigger; i--;) {
+    TString fired = ((TObjString*)tokens->At(i))->String(); 
+    if (fired.CompareTo(trigger)==0) return kTRUE;
+  }
+  return kFALSE;
+}
+
+//_____________________________________________________________________________
+void AliMuonsHFHeader::EventSelection(TString triggerName)
+{
+  // select event according to the "triggerName" & event selection cuts
+
+  fgIsEventSelected = kFALSE;
+  if (!this->IsPhysicsAccepted()) return;
+  if (!fgIsMC && !this->IsTriggerFired(triggerName)) return;
+  this->EventSelection();
+  return; 
+}
+
+//_____________________________________________________________________________
+void AliMuonsHFHeader::EventSelection()
+{
+  // select event according to the event selection cuts
+
+  fgIsEventSelected = kFALSE;
+  if (this->NVtxContributorsSPD()<fgCuts[0]) return;
+  if (TMath::Abs(this->VzSPD())>fgCuts[1]) return;
+  if (this->VtSPD()>fgCuts[2]) return;
+  fgIsEventSelected = kTRUE;
+  return;
+}
+
+//_____________________________________________________________________________
+void AliMuonsHFHeader::CreateHistograms(TList *listEvent, TList *listMuon, TList *listDimu)
+{
+  // create output histos of muon analysis according to the analysis mode & MC flag
+
+  this->CreateHistosEventH(listEvent);
+  if (fgIsMC) {
+    if (fgAnaMode!=2) this->CreateHistosMuonMC(listMuon);
+    if (fgAnaMode!=1) this->CreateHistosDimuMC(listDimu);
+  } else {
+    if (fgAnaMode!=2) this->CreateHistosMuonRD(listMuon);
+    if (fgAnaMode!=1) this->CreateHistosDimuRD(listDimu);
+  }
+  return;
+}
+
+//_____________________________________________________________________________
+void AliMuonsHFHeader::CreateHistosEventH(TList *list)
+{
+  // create histograms at event level
+
+  if (!list) list = new TList();
+  list->SetOwner();
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
+  const Int_t nHistos = 5;
+  char    *name[nHistos] = { "hVztSPD", "hVzxSPD", "hVzySPD", "hVzNcontr", "hVtNcontr"  };
+  Int_t  nbinsX[nHistos] = {      800 ,      800 ,      800 ,       800  ,        400   };
+  Double_t xlow[nHistos] = {      -40.,      -40.,      -40.,       -40. ,          0.  };
+  Double_t  xup[nHistos] = {       40.,       40.,       40.,        40. ,          4.  };
+  Int_t  nbinsY[nHistos] = {      400 ,      600 ,      600 ,       202  ,        202   };
+  Double_t ylow[nHistos] = {        0.,       -3.,       -3.,        -2.5,         -2.5 };
+  Double_t  yup[nHistos] = {        4.,        3.,        3.,       199.5,        199.5 };
+
+  TH2F *histo = 0;
+  for (Int_t i=0; i<nHistos; i++) {
+    histo = new TH2F(name[i], name[i], nbinsX[i], xlow[i], xup[i], nbinsY[i], ylow[i], yup[i]);
+    histo->Sumw2(); list->AddAt(histo, i); histo = 0;
+  }
+
+  TH1::AddDirectory(oldStatus);
+  return;
+}
+
+//_____________________________________________________________________________
+void AliMuonsHFHeader::CreateHistosMuonRD(TList *list)
+{
+  // create histograms for single muon
+
+  if (!list) list = new TList();
+  list->SetOwner();
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+  
+  const Int_t nHistos = 8; 
+  char    *name[nHistos] = {  "hP", "hPt", "hEta", "hDCA", "hTrg", "hCharge", "hUnfVtx" , "hEtaTimesDCA"};
+  Int_t   nbins[nHistos] = { 1500 ,  300 ,   100 ,   500 ,    4  ,       3  ,       80  ,         10000 };
+  Double_t xlow[nHistos] = {    0.,    0.,   -10.,     0.,   -0.5,      -1.5,      -40. ,             0.};
+  Double_t  xup[nHistos] = {  150.,   30.,     0.,   500.,    3.5,       1.5,       40. ,          1000.};
+
+  TH1F *histo = 0;
+  for (Int_t i=0; i<nHistos; i++) {
+    histo = new TH1F(name[i], name[i], nbins[i], xlow[i], xup[i]);
+    histo->Sumw2(); list->AddAt(histo, i); histo = 0;
+  }
+
+  TH1::AddDirectory(oldStatus);
+  return;
+}
+
+//_____________________________________________________________________________
+void AliMuonsHFHeader::CreateHistosDimuRD(TList *list)
+{
+  // create histograms for dimuon
+
+  if (!list) list = new TList();
+  list->SetOwner();
+  Bool_t oldStatus = TH1::AddDirectoryStatus();
+  TH1::AddDirectory(kFALSE);
+
+  TH1F *histo = 0;
+  const Int_t nHistos = 3;
+  char    *name[nHistos] = {  "hP", "hPt", "hInvM" };
+  Int_t   nbins[nHistos] = { 1500 ,  300 ,    300  };
+  Double_t xlow[nHistos] = {    0.,    0.,      0. };
+  Double_t  xup[nHistos] = {  150.,   30.,     30. };
+  char *dimuName[3] = {"DimuNN", "DimuNP", "DimuPP"};
+  for (Int_t i=0; i<3; i++) {
+    for (Int_t j=0; j<nHistos; j++) {
+      histo = new TH1F(Form("%s_%s",name[j],dimuName[i]), name[j], nbins[j], xlow[j], xup[j]);
+      histo->Sumw2(); list->AddAt(histo,i*nHistos+j); histo = 0;
+    }
+  }
+
+  TH1::AddDirectory(oldStatus);
+  return;
+}
+
+//_____________________________________________________________________________
+void AliMuonsHFHeader::CreateHistosMuonMC(TList *list)
+{
+  // create histograms for single muon with MC info
+
+  if (!list) list = new TList[AliMuonInfoStoreMC::NSources()];
+  for (Int_t i=AliMuonInfoStoreMC::NSources(); i--;) this->CreateHistosMuonRD(&list[i]);
+  return;
+}
+
+//_____________________________________________________________________________
+void AliMuonsHFHeader::CreateHistosDimuMC(TList *list)
+{
+  // create histograms for dimuon with MC info
+
+  if (!list) list = new TList[AliDimuInfoStoreMC::NSources()];
+  for (Int_t i=AliDimuInfoStoreMC::NSources(); i--;) this->CreateHistosDimuRD(&list[i]);
+  return;
+}
+
+//_____________________________________________________________________________
+void AliMuonsHFHeader::FillHistosEventH(TList *list)
+{
+  // fill histograms at event level according to event selection cuts
+
+  if (!list) return;
+  if (!AliMuonsHFHeader::IsEventSelected()) return;
+
+  const Int_t nHistos = 5;
+  Double_t vz = this->VzSPD();
+  Double_t vt = this->VtSPD();
+  Int_t    nc = this->NVtxContributorsSPD();
+  Double_t distX[nHistos] = { vz,            vz,            vz, vz, vt };
+  Double_t distY[nHistos] = { vt, this->VxSPD(), this->VySPD(), nc, nc };
+  for (Int_t i=nHistos; i--;) ((TH2F*)list->At(i))->Fill(distX[i], distY[i]);
+  return;
+}
+
+//_____________________________________________________________________________
+void AliMuonsHFHeader::FillHistosMuonRD(TList *list, AliMuonInfoStoreRD* const muonStoreRD)
+{
+  // fill histograms for single muon according to event & muon track selection cuts
+
+  if (!list) return;
+  if (!AliMuonsHFHeader::IsEventSelected()) return;
+  if (!muonStoreRD->MuonSelection()) return;
+
+  const Int_t nHistos = 8;
+  Double_t dist[nHistos] = {muonStoreRD->Momentum().Mag(),
+                            muonStoreRD->Momentum().Pt(),
+                            muonStoreRD->Momentum().Eta(),
+                            muonStoreRD->DCA(),
+                            muonStoreRD->MatchTrigger(),
+                            muonStoreRD->Charge(),
+                            this->VzSPD(),
+                            muonStoreRD->DCA()*TMath::Exp(-2.*muonStoreRD->Momentum().Eta())/1000.};
+  for (Int_t i=nHistos; i--;) ((TH1F*)list->At(i))->Fill(dist[i]);
+  return; 
+}
+
+//_____________________________________________________________________________
+void AliMuonsHFHeader::FillHistosDimuRD(TList *list, AliDimuInfoStoreRD* const dimuStoreRD)
+{
+  // fill histograms for dimuon according to evnet & dimuon candidates selection cuts
+
+  if (!list) return;
+  if (!AliMuonsHFHeader::IsEventSelected()) return;
+  if (!dimuStoreRD->DimuSelection()) return;
+
+  Int_t theDimu = 0;
+  if (dimuStoreRD->Charge()==0)     theDimu = 1;
+  else if (dimuStoreRD->Charge()>0) theDimu = 2;
+
+  const Int_t nHistos = 3;
+  Double_t dist[nHistos] = {dimuStoreRD->Momentum().Mag(),
+                            dimuStoreRD->Momentum().Pt(),
+                            dimuStoreRD->InvM()};
+  for (Int_t i=nHistos; i--;) ((TH1F*)list->At(theDimu*nHistos+i))->Fill(dist[i]);
+  return;
+}
+
+//_____________________________________________________________________________
+void AliMuonsHFHeader::FillHistosMuonMC(TList *list, AliMuonInfoStoreMC* const muonStoreMC)
+{
+  // fill histograms for single muon with MC info
+
+  Int_t srcMuon = muonStoreMC->MuonSource(); 
+  this->FillHistosMuonRD(&list[6], (AliMuonInfoStoreRD*)muonStoreMC);
+  if (srcMuon<0) this->FillHistosMuonRD(&list[5], (AliMuonInfoStoreRD*)muonStoreMC);
+  else this->FillHistosMuonRD(&list[srcMuon], (AliMuonInfoStoreRD*)muonStoreMC);
+  return;
+}
+
+//_____________________________________________________________________________
+void AliMuonsHFHeader::FillHistosDimuMC(TList *list, AliDimuInfoStoreMC* const dimuStoreMC)
+{
+  // fill histograms for dimuon with MC info
+
+  Int_t srcDimu = dimuStoreMC->DimuSource();
+  this->FillHistosDimuRD(&list[6], (AliDimuInfoStoreRD*)dimuStoreMC);
+  this->FillHistosDimuRD(&list[srcDimu], (AliDimuInfoStoreRD*)dimuStoreMC);
   return;
 }
index 56120e5..10c66a7 100644 (file)
 #ifndef ALIMUONSHFHEADER_H
 #define ALIMUONSHFHEADER_H
 
+/* Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//***********************************************************
+// Class AliMuonsHFHeader
+// class used to extract and store info at event level
+// Author: X-M. Zhang, zhang@clermont.in2p3.fr
+//                     zhangxm@iopp.ccnu.edu.cn
+//***********************************************************
+
 #include <TNamed.h>
+#include <TString.h>
+
+class TList;
+class AliAODEvent;
+class AliESDEvent;
 
-#include "AliVHeader.h"
-#include "AliVVertex.h"
+class AliMuonInfoStoreRD;
+class AliDimuInfoStoreRD;
+class AliMuonInfoStoreMC;
+class AliDimuInfoStoreMC;
 
 class AliMuonsHFHeader : public TNamed {
  public :
 
   AliMuonsHFHeader();
+  AliMuonsHFHeader(const AliMuonsHFHeader &src);
+  AliMuonsHFHeader& operator=(const AliMuonsHFHeader &src);
   ~AliMuonsHFHeader();
 
-  void SetHeader(AliVHeader *header);
-  void SetVertex(AliVVertex *vertex);
-
-  void SetMultSingleMuon(Int_t mul) {  fMultMuon = mul; }
-  void SetMultDimuon(Int_t mul)  {   fMultDimuon = mul; }
-  void SetCentrality(Double_t cen) { fCentrality = cen; }
-
-  void GetXYZ(Double_t *pos) const
-    { for (Int_t i=0; i<3; i++) pos[i]=fPosition[i]; }
-  Double_t GetXv() const { return fPosition[0]; }
-  Double_t GetYv() const { return fPosition[1]; }
-  Double_t GetZv() const { return fPosition[2]; }
-  Int_t GetNContributors() const { return fNContributors; }
-
-  ULong64_t GetTriggerMask() const { return fTriggerMask; }
-
-  Int_t    GetMultSingleMuon() const { return fMultMuon;   }
-  Int_t    GetMultDimuon()     const { return fMultDimuon; }
-  Double_t GetCentrality()     const { return fCentrality; }
-
-  Bool_t IsUnrecoVertex() const { return fUnrecoVertex; }
+  ULong64_t TriggerMask()     const { return fTriggerMask;        }
+  TString FiredTrigger()      const { return fFiredTrigger;       }
+  Int_t  NFiredTrigger()      const { return fNFiredTrigger;      }
+  Int_t  EventType()          const { return fEventType;          }
+  Bool_t IsPhysicsTriggered() const { return fIsPhysicsTriggered; }
+  Bool_t IsPhysicsAccepted()  const { return fIsPhysicsAccepted;  }
+  Bool_t IsTriggerFired(TString trigger);
+
+  void GetXYZ(Double_t *vtx) const { for (Int_t i=3; i--;) vtx[i]=fVtx[i]; }
+  Double_t Vx() const { return fVtx[0]; }
+  Double_t Vy() const { return fVtx[1]; }
+  Double_t Vz() const { return fVtx[2]; }
+  Double_t Vt() const { return TMath::Sqrt(fVtx[0]*fVtx[0] + fVtx[1]*fVtx[1]); }
+  Bool_t  IsUnrecoVertex() const { return fUnrecoVertex;  }
+  Int_t NVtxContributors() const { return fNContributors; }
+
+  void GetXYZSPD(Double_t *vtx) const { for (Int_t i=3; i--;) vtx[i]=fVtxSPD[i]; }
+  Double_t VxSPD() const { return fVtxSPD[0]; }
+  Double_t VySPD() const { return fVtxSPD[1]; }
+  Double_t VzSPD() const { return fVtxSPD[2]; }
+  Double_t VtSPD() const { return TMath::Sqrt(fVtxSPD[0]*fVtxSPD[0] + fVtxSPD[1]*fVtxSPD[1]); }
+  Bool_t  IsUnrecoVtxSPD()    const { return fUnrecoVtxSPD;     }
+  Int_t NVtxContributorsSPD() const { return fNContributorsSPD; }
+  Int_t NTrackletsSPD()       const { return fNTrackletsSPD;    }
+
+  Double_t Centrality() const { return fCentrality; }
+
+  void SetEvent(AliAODEvent *event);
+  void SetEvent(AliESDEvent *event);
+
+  void EventSelection(TString triggerName);
+  void CreateHistograms(TList *listEvent=0, TList *listMuon=0, TList *listDimu=0);
+  void CreateHistosEventH(TList *list);
+  void CreateHistosMuonRD(TList *list);
+  void CreateHistosDimuRD(TList *list);
+  void CreateHistosMuonMC(TList *list);
+  void CreateHistosDimuMC(TList *list);
+
+  void FillHistosEventH(TList *list);
+  void FillHistosMuonRD(TList *list, AliMuonInfoStoreRD* const muonStoreRD);
+  void FillHistosDimuRD(TList *list, AliDimuInfoStoreRD* const dimuStoreRD);
+  void FillHistosMuonMC(TList *list, AliMuonInfoStoreMC* const muonStoreMC);
+  void FillHistosDimuMC(TList *list, AliDimuInfoStoreMC* const dimuStoreMC);
+
+  static const char* StdBranchName()    { return fgkStdBranchName.Data(); }
+  static const Bool_t IsEventSelected() { return fgIsEventSelected; }
+  static void SetAnaMode(Int_t anaMode=0)        { fgAnaMode=anaMode; }
+  static void SetIsMC(Int_t isMC=kFALSE)         { fgIsMC   =isMC;    }
+  static void SetSelectionCuts(Double_t cuts[3]) { for (Int_t i=3; i--;) fgCuts[i]=cuts[i]; }
 
  private :
 
-  ULong64_t fTriggerMask;   // trigger mask
-
-  Double32_t fPosition[3];  // position of vtx
-  Int_t fNContributors;     // number of contributor for vtx
-
-  Int_t fMultMuon;      // event multiplicity
-  Int_t fMultDimuon;
-  Double_t fCentrality;     // event centrality class
-
-  Bool_t fUnrecoVertex;
-
-  ClassDef(AliMuonsHFHeader, 1)
+  void SetFiredTrigger(TString str);
+  void PhysicsTriggerAna(const AliESDEvent *esd);
+  void EventSelection();
+
+  static const TString fgkStdBranchName;  // Standard branch name
+  static Bool_t fgIsEventSelected;        // flag for event selection
+  static Int_t  fgAnaMode;                // analysis mode
+  static Bool_t fgIsMC;                   // flag to use MC
+  static Double_t fgCuts[3];  // 0, low limit of num. of vtx contributors
+                              // 1, up limit of vz
+                              // 2, up limit of vt
+
+  ULong64_t fTriggerMask;      // trigger mask
+  TString fFiredTrigger;       // fired of trigger class of event
+  Int_t fNFiredTrigger;        // num. of fired trigger class
+  Bool_t fIsPhysicsTriggered;  // flag of final physics trigger from AliPhysicsSelection
+  Bool_t fIsPhysicsAccepted;   // flag of physiscs selection w/ BKG Id
+  Int_t fEventType;            // event type
+
+  Double32_t fVtx[3];      // position of vtx
+  Bool_t fUnrecoVertex;    // flag for unreco vtx
+  Int_t fNContributors;    // num. of contributors of vtx rec
+
+  Double32_t fVtxSPD[3];   // position of vtx
+  Bool_t fUnrecoVtxSPD;    // flag for unreco vtx
+  Int_t fNContributorsSPD; // num. of contributors of vtx rec
+  Int_t fNTrackletsSPD;  // num. of SPD tracklets
+
+  Double_t fCentrality;    // event centrality class
+
+  ClassDef(AliMuonsHFHeader, 2)
 };
 
 #endif