Dimuon branch in AODEvent
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 25 Feb 2010 12:05:26 +0000 (12:05 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 25 Feb 2010 12:05:26 +0000 (12:05 +0000)
(Roberta Arnaldi)

15 files changed:
ANALYSIS/AliAnalysisTaskSE.cxx
ANALYSIS/AliAnalysisTaskSE.h
PWG3/PWG3muonLinkDef.h
PWG3/libPWG3muon.pkg
PWG3/muon/AliAnalysisTaskESDMuonFilter.cxx
PWG3/muon/AliAnalysisTaskESDMuonFilter.h
PWG3/muon/AliAnalysisTaskMuonAODfromGeneral.cxx
STEER/AODLinkDef.h
STEER/AliAODDimuon.cxx [new file with mode: 0644]
STEER/AliAODDimuon.h [new file with mode: 0644]
STEER/AliAODEvent.cxx
STEER/AliAODEvent.h
STEER/AliAODHandler.cxx
STEER/AliAODHandler.h
STEER/libAOD.pkg

index e8c906b..da11178 100644 (file)
@@ -44,6 +44,7 @@
 #include "AliMCEvent.h"
 #include "AliStack.h"
 #include "AliLog.h"
+#include "AliAODDimuon.h"
 
 
 ClassImp(AliAnalysisTaskSE)
@@ -61,6 +62,7 @@ TClonesArray*    AliAnalysisTaskSE::fgAODMCParticles    = NULL;
 AliAODTracklets* AliAnalysisTaskSE::fgAODTracklets      = NULL;
 AliAODCaloCells* AliAnalysisTaskSE::fgAODEmcalCells     = NULL;
 AliAODCaloCells* AliAnalysisTaskSE::fgAODPhosCells      = NULL;
+TClonesArray*    AliAnalysisTaskSE::fgAODDimuons        = NULL;
 
 AliAnalysisTaskSE::AliAnalysisTaskSE():
     AliAnalysisTask(),
@@ -276,6 +278,14 @@ void AliAnalysisTaskSE::CreateOutputObjects()
                fgAODMCParticles->SetName("mcparticles");
                handler->AddBranch("TClonesArray", &fgAODMCParticles);
            }
+           if ((handler->NeedsDimuonsBranchReplication() || merging) && !(fgAODDimuons))      
+           {   
+               if (fDebug > 1) AliInfo("Replicating dimuon branch\n");
+               fgAODDimuons = new TClonesArray("AliAODDimuon",0);
+               fgAODDimuons->SetName("dimuons");
+               handler->AddBranch("TClonesArray", &fgAODDimuons);
+           }    
+
            // cache the pointerd in the AODEvent
            fOutputAOD->GetStdContent();
        }
@@ -375,7 +385,13 @@ void AliAnalysisTaskSE::Exec(Option_t* option)
                TClonesArray* mcParticles = (TClonesArray*) ((dynamic_cast<AliAODEvent*>(InputEvent()))->FindListObject("mcparticles"));
                new (fgAODMCParticles) TClonesArray(*mcParticles);
            }
-       
+           
+           if ((handler->NeedsDimuonsBranchReplication() || merging) && (fgAODDimuons))
+           {
+               TClonesArray* dimuons = (dynamic_cast<AliAODEvent*>(InputEvent()))->GetDimuons();
+               new (fgAODDimuons) TClonesArray(*dimuons);
+           }
+
            // Additional merging if needed
            if (merging) {
                // mcParticles
index 355a979..62d706c 100644 (file)
@@ -83,10 +83,11 @@ class AliAnalysisTaskSE : public AliAnalysisTask
     static AliAODTracklets* fgAODTracklets;     //! Tracklets for replication
     static AliAODCaloCells* fgAODEmcalCells;    //! Emcal Cell replication
     static AliAODCaloCells* fgAODPhosCells;     //! Phos  Cell replication
+    static TClonesArray*    fgAODDimuons;       //! Dimuons replication
     // Event Selection
     Bool_t                  fSelectCollisions;   //  Task processes collision candidates only
      
-    ClassDef(AliAnalysisTaskSE, 2); // Analysis task for standard jet analysis
+    ClassDef(AliAnalysisTaskSE, 3); // Analysis task for standard jet analysis
 };
  
 #endif
index c0caf94..793deed 100644 (file)
@@ -4,7 +4,6 @@
 #pragma link off all classes;
 #pragma link off all functions;
 
-#pragma link C++ class AliAODDimuon+;
 #pragma link C++ class AliAODEventInfo+;
 #pragma link C++ class AliAnalysisTaskMuonAODfromGeneral+;
 #pragma link C++ class AliAnalysisTaskFromStandardToMuonAOD+;
@@ -29,4 +28,3 @@
 #pragma link C++ class AliAnalysisTaskSEMuonsHF+;
 #endif
 
-
index 2db0c79..6754e4b 100644 (file)
@@ -7,7 +7,6 @@ SRCS:=   muon/AliAnalysisTaskESDMuonFilter.cxx \
                  muon/AliAnalysisTaskLUT.cxx  \
                  muon/AliAnalysisTaskTrigChEff.cxx \
                  muon/AliAnalysisTaskLinkToMC.cxx \
-                 muon/AliAODDimuon.cxx \
                  muon/AliAODEventInfo.cxx \
                  muon/AliESDMuonTrackCuts.cxx \
                  muon/AliAnalysisTaskSingleMuESD.cxx \
@@ -35,3 +34,4 @@ ifeq (win32gcc,$(ALICE_TARGET))
 PACKSOFLAGS:= $(SOFLAGS) -L$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET) -lSTEERBase \
                          -lESD -lAOD -lANALYSIS -lANALYSISalice
 endif
+
index 339f08f..ea0b9da 100644 (file)
@@ -36,6 +36,7 @@
 #include "AliMCEvent.h"
 #include "AliMCEventHandler.h"
 #include "AliAODMCParticle.h"
+#include "AliAODDimuon.h"
 
 ClassImp(AliAnalysisTaskESDMuonFilter)
 
@@ -71,6 +72,7 @@ void AliAnalysisTaskESDMuonFilter::Init()
   AliAODHandler *aodH = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
   if (!aodH) Fatal("UserCreateOutputObjects", "No AOD handler. Aborting.");
   if(fEnableMuonAOD)aodH->AddFilteredAOD("AliAOD.Muons.root", "MuonEvents");
+  if(fEnableDimuonAOD)aodH->AddFilteredAOD("AliAOD.Dimuons.root", "DimuonEvents");
 }
 
 
@@ -119,12 +121,21 @@ void AliAnalysisTaskESDMuonFilter::ConvertESDtoAOD()
   
   // Loop on muon tracks to fill the AOD track branch
   Int_t nMuTracks = esd->GetNumberOfMuonTracks();
+
   
   // Update number of positive and negative tracks from AOD event (M.G.)
   Int_t nPosTracks = header->GetRefMultiplicityPos();
   Int_t nNegTracks = header->GetRefMultiplicityNeg();
   
+  // Access to the AOD container of dimuons
+  TClonesArray &dimuons = *(AODEvent()->GetDimuons());
+  AliAODDimuon *aodDimuon = 0x0;
+  
   Bool_t MuonsExist = kFALSE;
+  Bool_t DimuonsExist = kFALSE;
+  Int_t firstMuonTrack=0;
+  Int_t nMuons=0;
+  Int_t jDimuons=0;
 
   for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack) {
     esdMuTrack = esd->GetMuonTrack(nMuTrack);
@@ -140,8 +151,6 @@ void AliAnalysisTaskESDMuonFilter::ConvertESDtoAOD()
        }  
      }
 
-    if(!MuonsExist) MuonsExist=kTRUE;
-
     p[0] = esdMuTrack->Px(); 
     p[1] = esdMuTrack->Py(); 
     p[2] = esdMuTrack->Pz();
@@ -177,17 +186,37 @@ void AliAnalysisTaskESDMuonFilter::ConvertESDtoAOD()
     aodTrack->SetMuonClusterMap(esdMuTrack->GetMuonClusterMap());
     aodTrack->SetMatchTrigger(esdMuTrack->GetMatchTrigger());
     aodTrack->Connected(esdMuTrack->IsConnected());
-    
     primary->AddDaughter(aodTrack);
     
     if (esdMuTrack->Charge() > 0) nPosTracks++;
     else nNegTracks++;
+    
+    // fill dimuon branch
+    if(!MuonsExist) {
+      MuonsExist=kTRUE;
+      firstMuonTrack=jTracks-1.;
+    }  
+    nMuons++;
+    if(nMuons==2) DimuonsExist = kTRUE;   
+    if(DimuonsExist) { 
+      AliAODTrack *track0 = (AliAODTrack*)tracks.At(firstMuonTrack);
+      AliAODTrack *track1 = (AliAODTrack*)tracks.At(jTracks-1);
+      aodDimuon = new(dimuons[jDimuons++]) AliAODDimuon(tracks.At(jTracks-1),tracks.At(firstMuonTrack));
+      //AliAODDimuon *dimuon0 = (AliAODDimuon*)dimuons.At(0);
+      //printf("Dimuon: mass = %f, px=%f, py=%f, pz=%f\n",dimuon0->M(),dimuon0->Px(),dimuon0->Py(),dimuon0->Pz());  
+      //AliAODTrack  *mu0 = (AliAODTrack*) dimuon0->GetMu(0);
+      //AliAODTrack  *mu1 = (AliAODTrack*) dimuon0->GetMu(1);
+      //printf("Muon0 px=%f py=%f pz=%f\n",mu0->Px(),mu0->Py(),mu0->Pz());
+      //printf("Muon1 px=%f py=%f pz=%f\n",mu1->Px(),mu1->Py(),mu1->Pz());
+      break;
+    }
+
   }
   
   header->SetRefMultiplicity(jTracks); 
   header->SetRefMultiplicityPos(nPosTracks);
   header->SetRefMultiplicityNeg(nNegTracks);
-
+  
   // From Andrei
   if(fEnableMuonAOD && MuonsExist){
     AliAODExtension *extMuons = dynamic_cast<AliAODHandler*>
@@ -195,6 +224,12 @@ void AliAnalysisTaskESDMuonFilter::ConvertESDtoAOD()
     extMuons->SelectEvent();
   }
 
+  if(fEnableDimuonAOD && DimuonsExist){
+    AliAODExtension *extDimuons = dynamic_cast<AliAODHandler*>
+   ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler())->GetFilteredAOD("AliAOD.Dimuons.root");
+    extDimuons->SelectEvent();
+  }
+
 }
 
 void AliAnalysisTaskESDMuonFilter::Terminate(Option_t */*option*/)
index 76289ed..b4723c5 100644 (file)
@@ -28,6 +28,7 @@ class AliAnalysisTaskESDMuonFilter : public AliAnalysisTaskSE
     // Setters\r
     virtual void SetTrackFilter(AliAnalysisFilter* trackF) {fTrackFilter = trackF;}\r
     void SetWriteMuonAOD(Bool_t enableMuonAOD){fEnableMuonAOD = enableMuonAOD;}\r
+    void SetWriteDimuonAOD(Bool_t enableDimuonAOD){fEnableDimuonAOD = enableDimuonAOD;}\r
 \r
  private:\r
     AliAnalysisTaskESDMuonFilter(const AliAnalysisTaskESDMuonFilter&);\r
@@ -35,6 +36,7 @@ class AliAnalysisTaskESDMuonFilter : public AliAnalysisTaskSE
     void PrintMCInfo(AliStack *pStack,Int_t label); // for debugging\r
     AliAnalysisFilter* fTrackFilter; //  Track Filter\r
     Bool_t fEnableMuonAOD; // flag for enabling Muon AOD production\r
+    Bool_t fEnableDimuonAOD; // flag for enabling Dimuon AOD production\r
     ClassDef(AliAnalysisTaskESDMuonFilter, 1); // Analysis task for standard ESD filtering\r
 \r
 };\r
index 637d8f7..9686194 100644 (file)
@@ -171,13 +171,12 @@ void AliAnalysisTaskMuonAODfromGeneral::Exec(Option_t *) {
   Int_t jDimuons=0;
   for(Int_t i=0; i<nMuTracks; i++){
      for(Int_t j=i+1; j<nMuTracks; j++){
-       new(rDimuons[jDimuons++]) AliAODDimuon(tracks[i],tracks[j],fInfos);
+       new(rDimuons[jDimuons++]) AliAODDimuon(tracks[i],tracks[j]);
      } 
   }
   
   fInfos->SetBeamEnergy(fBeamEnergy);
   fInfos->SetEv(fNewAOD);
-  fInfos->SetEi(fInfos);
   fInfos->SetHe(header);
   fInfos->SetTr(fNewAOD->GetTracks());
   fInfos->SetDi(fDimuons);
index fe765bf..3f74305 100644 (file)
@@ -43,6 +43,7 @@
 #pragma link C++ class AliAODMCHeader+;
 #pragma link C++ class AliAODPWG4Particle+;
 #pragma link C++ class AliAODPWG4ParticleCorrelation+;
+#pragma link C++ class AliAODDimuon+;
 
 
 #endif
diff --git a/STEER/AliAODDimuon.cxx b/STEER/AliAODDimuon.cxx
new file mode 100644 (file)
index 0000000..c9b9180
--- /dev/null
@@ -0,0 +1,456 @@
+// AliAODDimuon: a class for AODs for the MUON Arm of the ALICE Experiment
+// Author: P. Cortese, Universita' del Piemonte Orientale in Alessandria and
+// INFN of Torino - Italy
+//
+// The class defines a dimuon pair object from two AliAODTrack objects.
+// AliAODDimuon objects are supposed to be added to the AliAODEvent structure
+// during analysis. They would then allow to calculate the dimuon-related
+// kinematic variables with a minimal disk occupancy.
+// The payload of the class has been reduced to two pointers to the two
+// tracks. An instance of this class has also to be added to the AliAODEvent 
+// structure to provide additional information that is specific to MUON and 
+// therefore has not been included into the AOD header.
+// Two transient data members are not stored on file as they can be recomputed
+// at runtime.
+//
+
+#include "AliAODDimuon.h"
+#include "TLorentzVector.h"
+#define AliAODDimuon_CXX
+
+ClassImp(AliAODDimuon)
+
+//______________________________________________________________________________
+AliAODDimuon::AliAODDimuon():AliVParticle(),fP(0),fMProton(0.93827231)
+{
+  // default constructor
+  fMu[0]=0;
+  fMu[1]=0;
+}
+
+//______________________________________________________________________________
+AliAODDimuon::AliAODDimuon(const AliAODDimuon& dimu):
+  AliVParticle(dimu),
+  fP(0),fMProton(0.93827231)
+{
+  // copy constructor
+  fMu[0]=dimu.Mu(0);
+  fMu[1]=dimu.Mu(1);
+}
+
+//______________________________________________________________________________
+AliAODDimuon &AliAODDimuon::operator=(const AliAODDimuon& dimu)
+{
+  // assignment operator
+  if(&dimu != this){
+    fP=0;
+    fMProton=0.93827231;
+    fMu[0]=dimu.Mu(0);
+    fMu[1]=dimu.Mu(1);
+  }
+  return *this;
+}
+
+//______________________________________________________________________________
+AliAODDimuon::AliAODDimuon(TObject *mu0, TObject *mu1):
+  fP(0),fMProton(0.93827231)
+{
+  // Creates a dimuon pair from two tracks
+  
+  //printf("Creating dimuon from %p %p\n",mu0,mu1);
+  fMu[0]=mu0;
+  fMu[1]=mu1;
+}
+
+//______________________________________________________________________________
+AliAODDimuon::~AliAODDimuon()
+{
+  // destructor
+  if(fP)delete fP;
+  fP=0;
+}
+
+//______________________________________________________________________________
+void AliAODDimuon::BookP(){
+  // Fills the dimuon momentum if not filled yet
+  static UInt_t unID[2]={0,0};
+  if(!fP){
+    fP=new TLorentzVector(Px(),Py(),Pz(),E());
+    unID[0]=fMu[0].GetUniqueID();
+    unID[1]=fMu[1].GetUniqueID();
+  }
+  // For efficiency reasons
+  if((unID[0]!=fMu[0].GetUniqueID())||(unID[1]!=fMu[1].GetUniqueID())){
+    fP->SetPxPyPzE(Px(),Py(),Pz(),E());
+    unID[0]=fMu[0].GetUniqueID();
+    unID[1]=fMu[1].GetUniqueID();
+  }
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::Px() const {
+  // Px of the dimuon
+  if(this->CheckPointers())return -999999999;
+  return ((AliAODTrack*)fMu[0].GetObject())->Px()+
+         ((AliAODTrack*)fMu[1].GetObject())->Px();
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::Py() const {
+  // Py of the dimuon
+  if(this->CheckPointers())return -999999999;
+  return ((AliAODTrack*)fMu[0].GetObject())->Py()+
+         ((AliAODTrack*)fMu[1].GetObject())->Py();
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::Pz() const {
+  // Pz of the dimuon
+  if(this->CheckPointers())return -999999999;
+  return ((AliAODTrack*)fMu[0].GetObject())->Pz()+
+         ((AliAODTrack*)fMu[1].GetObject())->Pz();
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::Pt() const {
+  // Pt of the dimuon
+  if(this->CheckPointers())return -999999999;
+  Double_t px=Px();
+  Double_t py=Py();
+  return TMath::Sqrt(px*px+py*py);
+  return -999999999;
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::E() const {
+  // Dimuon energy
+  if(this->CheckPointers())return -999999999;
+  return ((AliAODTrack*)fMu[0].GetObject())->E()+
+         ((AliAODTrack*)fMu[1].GetObject())->E();
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::P() const {
+  // This is just to override the virtual function
+  printf("You should never call: Double_t AliAODDimuon::P() const\n");
+  return -999999999;
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::P() {
+  // Dimuon momentum
+  if(this->CheckPointers())return -999999999;
+  BookP();
+  return fP->P();
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::M() const {
+  // This is just to override the virtual function
+  printf("You should never call: Double_t AliAODDimuon::M() const\n");
+  return -999999999;
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::M() {
+  // Dimuon invariant mass
+  if(this->CheckPointers())return -999999999;
+  BookP();
+  return fP->M();
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::Mass() {
+  // Dimuon invariant mass
+  if(this->CheckPointers())return -999999999;
+  BookP();
+  return fP->M();
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::Eta() const {
+  // This is just to override the virtual function
+  printf("You should never call: Double_t AliAODDimuon::Eta() const\n");
+  return -999999999;
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::Eta() {
+  // Dimuon pseudorapidity
+  if(this->CheckPointers())return -999999999;
+  BookP();
+  return fP->Eta();
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::Phi() const {
+  // This is just to override the virtual function
+  printf("You should never call: Double_t AliAODDimuon::Phi() const\n");
+  return -999999999;
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::Phi() {
+  // Dimuon asimuthal angle
+  if(this->CheckPointers())return -999999999;
+  BookP();
+  return fP->Phi();
+}
+//______________________________________________________________________________
+Double_t AliAODDimuon::Theta() const {
+  // This is just to override the virtual function
+  printf("You should never call: Double_t AliAODDimuon::Theta() const\n");
+  return -999999999;
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::Theta() {
+  // Dimuon polar angle
+  if(this->CheckPointers())return -999999999;
+  BookP();
+  return fP->Theta();
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::Y() const {
+  // This is just to override the virtual function
+  printf("You should never call: Double_t AliAODDimuon::Y() const\n");
+  return -999999999;
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::Y() {
+  // Dimuon rapidity
+  if(this->CheckPointers())return -999999999;
+  BookP();
+  return fP->Rapidity();
+}
+
+//______________________________________________________________________________
+Short_t AliAODDimuon::Charge() const {
+  // Dimuon charge
+  if(this->CheckPointers())return -999;
+  return ((AliAODTrack*)fMu[0].GetObject())->Charge()+
+         ((AliAODTrack*)fMu[1].GetObject())->Charge();
+}
+
+//______________________________________________________________________________
+Int_t AliAODDimuon::CheckPointers() const{
+  // Checks if the track pointers have been initialized
+  if(fMu[0]==0||fMu[1]==0){
+    printf("Dimuon not initialized\n");
+    return -999;
+  }
+  if((fMu[0].GetObject())==0||(fMu[1].GetObject())==0){
+    printf("Can not get objects. Got: %p %p\n",fMu[0].GetObject(),fMu[1].GetObject());
+    return -999;
+  }
+  return 0;
+}
+
+//______________________________________________________________________________
+void AliAODDimuon::SetMu(Int_t imu, AliAODTrack *mu){
+  // Assign a track pointer
+  if (imu==0||imu==1){
+    fMu[imu]=mu;
+  }
+}
+
+//______________________________________________________________________________
+void AliAODDimuon::SetMuons(AliAODTrack *mu0, AliAODTrack *mu1){
+  // Assign the track pointers
+  fMu[0]=mu0;
+  fMu[1]=mu1;
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::XF() {
+  // Dimuon Feynman x
+  //Double_t ebeam=((AliAODEventInfo*)fEi.GetObject())->EBeam();
+  Double_t ebeam = 3500.; // temporary
+  if(ebeam<=0){
+    printf("AliAODDimuon::xf: can not compute xf with EBeam=%f\n",ebeam);
+    return -999999999;
+  }
+  if(this->CheckPointers())return -999999999;
+  BookP();
+  Double_t mDimu=M();
+  Double_t pMax=TMath::Sqrt(ebeam*ebeam-mDimu*mDimu);
+  return Pz()/pMax;
+}
+
+//______________________________________________________________________________
+// Calculation the Collins-Soper angle (adapted from code by R. Arnaldi)
+Double_t AliAODDimuon::CostCS(){
+  // Cosinus of the Collins-Soper polar decay angle
+  if(CheckPointers())return -999999999;
+  Double_t ebeam=3500.;  //temporary
+  if(ebeam<=0){
+    printf("Can not compute costCS with EBeam=%f\n",ebeam);
+    return -999999999;
+  }
+  Double_t mp=fMProton;
+  Double_t pbeam=TMath::Sqrt(ebeam*ebeam-mp*mp);
+  Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
+  Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
+  Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
+  Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
+  Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
+  Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
+  Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
+  Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
+  Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
+  Double_t mu2Charge=((AliAODTrack*)fMu[1].GetObject())->Charge();
+
+  // Fill the Lorentz vector for projectile and target
+  // For the moment we do not consider the crossing angle
+  // Projectile runs towards the MUON arm
+  TLorentzVector pProjLab(0.,0.,-pbeam,ebeam); // projectile
+  TLorentzVector pTargLab(0.,0., pbeam,ebeam); // target
+  //
+  // --- Get the muons parameters in the LAB frame
+  //
+  TLorentzVector pMu1Lab(pla10,pla11,pla12,e1);
+  TLorentzVector pMu2Lab(pla20,pla21,pla22,e2);
+  //
+  // --- Obtain the dimuon parameters in the LAB frame
+  //
+  TLorentzVector pDimuLab=pMu1Lab+pMu2Lab;
+  //
+  // --- Translate the dimuon parameters in the dimuon rest frame
+  //
+  TVector3 beta=(-1./pDimuLab.E())*pDimuLab.Vect();
+  TLorentzVector pMu1Dimu=pMu1Lab;
+  TLorentzVector pMu2Dimu=pMu2Lab;
+  TLorentzVector pProjDimu=pProjLab;
+  TLorentzVector pTargDimu=pTargLab;
+  pMu1Dimu.Boost(beta);
+  pMu2Dimu.Boost(beta);
+  pProjDimu.Boost(beta);
+  pTargDimu.Boost(beta);
+  //
+  // --- Determine the z axis for the CS angle 
+  //
+  TVector3 zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
+  //
+  // --- Determine the CS angle (angle between mu+ and the z axis defined above)
+  //
+  Double_t cost;
+  if(mu1Charge > 0) {
+    cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());
+    // Theta CS is not properly defined for Like-Sign muons
+    if(mu2Charge > 0 && cost<0) cost=-cost;
+  } else { 
+    // Theta CS is not properly defined for Like-Sign muons
+    cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());
+    if(mu2Charge < 0 && cost<0) cost=-cost;
+  }
+  return cost;
+}
+
+//______________________________________________________________________________
+// Calculation the Helicity polarization angle (adapted from code by R. Arnaldi)
+Double_t AliAODDimuon::CostHe(){
+  // Cosinus of the polar decay angle in the Helicity reference frame
+  if(CheckPointers())return -999999999;
+  Double_t ebeam=3500; //temporary
+  if(ebeam<=0){
+    printf("Can not compute costCS with EBeam=%f\n",ebeam);
+    return -999999999;
+  }
+  Double_t pbeam=TMath::Sqrt(ebeam*ebeam-fMProton*fMProton);
+  Double_t pla10=((AliAODTrack*)fMu[0].GetObject())->Px();
+  Double_t pla11=((AliAODTrack*)fMu[0].GetObject())->Py();
+  Double_t pla12=((AliAODTrack*)fMu[0].GetObject())->Pz();
+  Double_t e1=((AliAODTrack*)fMu[0].GetObject())->E();
+  Double_t mu1Charge=((AliAODTrack*)fMu[0].GetObject())->Charge();
+  Double_t pla20=((AliAODTrack*)fMu[1].GetObject())->Px();
+  Double_t pla21=((AliAODTrack*)fMu[1].GetObject())->Py();
+  Double_t pla22=((AliAODTrack*)fMu[1].GetObject())->Pz();
+  Double_t e2=((AliAODTrack*)fMu[1].GetObject())->E();
+  Double_t mu2Charge=((AliAODTrack*)fMu[1].GetObject())->Charge();
+
+  // Fill the Lorentz vector for projectile and target
+  // For the moment we consider no crossing angle
+  // Projectile runs towards the MUON arm
+  TLorentzVector pProjLab(0.,0.,-pbeam,ebeam); // projectile
+  TLorentzVector pTargLab(0.,0., pbeam,ebeam); // target
+  //
+  // --- Get the muons parameters in the LAB frame
+  //
+  TLorentzVector pMu1Lab(pla10,pla11,pla12,e1);
+  TLorentzVector pMu2Lab(pla20,pla21,pla22,e2);
+  //
+  // --- Obtain the dimuon parameters in the LAB frame
+  //
+  TLorentzVector pDimuLab=pMu1Lab+pMu2Lab;
+  //
+  // --- Translate the dimuon parameters in the dimuon rest frame
+  //
+  TVector3 beta=(-1./pDimuLab.E())*pDimuLab.Vect();
+  TLorentzVector pMu1Dimu=pMu1Lab;
+  TLorentzVector pMu2Dimu=pMu2Lab;
+  pMu1Dimu.Boost(beta);
+  pMu2Dimu.Boost(beta);
+  //
+  // --- Translate the dimuon parameters in the CM frame
+  //
+  TLorentzVector pDimuCM; //CM frame
+  TVector3 beta2;
+  beta2=(-1./(fMProton+pProjLab.E()))*pProjLab.Vect();
+  pDimuCM=pDimuLab;
+  pDimuCM.Boost(beta2);
+  //
+  // --- Determine the z axis for the calculation of the polarization angle
+  // (i.e. the direction of the dimuon in the CM system)
+  //
+  TVector3 zaxis;
+  zaxis=(pDimuCM.Vect()).Unit();
+  //
+  // --- Calculation of the polarization angle (Helicity)
+  // (angle between mu+ and the z axis defined above)
+  //
+  Double_t cost;
+  if(mu1Charge > 0) {
+    cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());
+    // Theta Helicity is not properly defined for Like-Sign muons
+    if(mu2Charge > 0 && cost<0) cost=-cost;
+  } else { 
+    cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());
+    // Theta Helicity is not properly defined for Like-Sign muons
+    if(mu2Charge < 0 && cost<0) cost=-cost;
+  }  
+  return cost;
+}
+
+//______________________________________________________________________________
+Int_t AliAODDimuon::AnyPt(){
+  // Test if the two muons match two trigger tracks
+  if(this->CheckPointers())return 0;
+  return (((AliAODTrack*)fMu[0].GetObject())->MatchTrigger())&&
+         (((AliAODTrack*)fMu[0].GetObject())->MatchTrigger());
+}
+
+//______________________________________________________________________________
+Int_t AliAODDimuon::LowPt(){
+  // Test if the two muons match two trigger tracks with a "Low Pt" cut
+  if(this->CheckPointers())return 0;
+  return (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerLowPt())&&
+         (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerLowPt());
+}
+
+//______________________________________________________________________________
+Int_t AliAODDimuon::HighPt(){
+  // Test if the two muons match two trigger tracks with a "High Pt" cut
+  if(this->CheckPointers())return 0;
+  return (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerHighPt())&&
+         (((AliAODTrack*)fMu[0].GetObject())->MatchTriggerHighPt());
+}
+
+//______________________________________________________________________________
+Double_t AliAODDimuon::MaxChi2Match(){
+  // Maximum matching Chi2 between track and trigger track
+  if(this->CheckPointers())return -999999999;
+  return TMath::Max((((AliAODTrack*)fMu[0].GetObject())->GetChi2MatchTrigger()),
+                    (((AliAODTrack*)fMu[0].GetObject())->GetChi2MatchTrigger()));
+}
diff --git a/STEER/AliAODDimuon.h b/STEER/AliAODDimuon.h
new file mode 100644 (file)
index 0000000..21d05ef
--- /dev/null
@@ -0,0 +1,109 @@
+#ifndef AliAODDimuon_H
+#define AliAODDimuon_H
+
+// AliAODDimuon: a class for AODs for the MUON Arm of the ALICE Experiment
+// Author: P. Cortese, Universita' del Piemonte Orientale in Alessandria and
+// INFN of Torino - Italy
+//
+// The class defines a dimuon pair object from two AliAODTrack objects.
+// AliAODDimuon objects are supposed to be added to the AliAODEvent structure
+// during analysis. They would then allow to calculate the dimuon-related
+// kinematic variables with a minimal disk occupancy.
+// The payload of the class has been reduced to two pointers to the two
+// tracks. An instance of this class has also to be added to the AliAODEvent 
+// structure to provide additional information that is specific to MUON and 
+// therefore has not been included into the AOD header.
+// Two transient data members are not stored on file as they can be recomputed
+// at runtime.
+//
+
+// 2007/07/07 v1.00 Initial version
+// 2007/12/06 v1.01 Introduction of AliAODEventInfo
+// 2007/12/18 v1.02 Corrected CostCS for Like-Sign, added CostKh, CostHe and xf
+// 2008/02/01 v1.03 Apply coding conventions
+
+#include "TRef.h"
+#include "AliVParticle.h"
+#include "AliAODTrack.h"
+
+class TLorentzVector;
+
+class AliAODDimuon: public AliVParticle {
+public:
+  AliAODDimuon();
+  AliAODDimuon(const AliAODDimuon& dimu);
+  AliAODDimuon &operator=(const AliAODDimuon& dimu);
+  AliAODDimuon(TObject *mu0, TObject *mu1);
+  virtual ~AliAODDimuon();
+
+  // Methods to access kinematics
+  virtual Double_t Px() const;
+  virtual Double_t Py() const;
+  virtual Double_t Pz() const;
+  virtual Bool_t PxPyPz(Double_t* p) const { p[0]=Px(); p[1]=Py(); p[2]=Pz(); return 1;}
+  virtual Double_t Pt() const;
+  virtual Double_t P() const;
+
+  virtual Double_t OneOverPt() const {return Pt()>0 ? 1./Pt() : -999999999;}
+  virtual Double_t Phi() const;
+  virtual Double_t Theta() const;
+
+  virtual Double_t E() const;
+  virtual Double_t M() const;
+  
+  virtual Double_t Eta() const;
+  virtual Double_t Y() const;
+  
+  virtual Short_t Charge() const;
+
+  // Dimuon vertex will be implemented when the muon track covariance matrix 
+  // at vertex will be included in the ESD (and AOD)
+  // It would require also the information about magnetic field when filling AOD
+  virtual Double_t Xv() const {return -999999999;}
+  virtual Double_t Yv() const {return -999999999;}
+  virtual Double_t Zv() const {return -999999999;}
+  virtual Bool_t XvYvZv(Double_t* v) const { v[0]=-999999999; v[1]=-999999999; v[2]=-999999999; return 0;}
+
+  Double_t P();
+  Double_t Phi();
+  Double_t Theta();
+  Double_t M();
+  Double_t Mass();
+  Double_t Eta();
+  Double_t Y();
+
+  // Added functions
+  Double_t XF();     // Feynman x
+  Double_t CostCS(); // Cosinus of the Collins-Soper polar decay angle
+  Double_t CostHe(); // Cosinus of the Helicity polar decay angle
+  Int_t AnyPt();
+  Int_t LowPt();
+  Int_t HighPt();
+  Double_t MaxChi2Match();
+  // PID
+  virtual const Double_t *PID() const {return 0;} // return PID object (to be defined, still)
+  //
+  Int_t GetLabel() const {return -1;}
+
+  // Additional getters and setters
+  AliAODTrack* GetMu(Int_t imu=0) const {return (imu==0||imu==1)&&(fMu[imu]!=0) ? (AliAODTrack*)fMu[imu].GetObject() : 0; } // Get a pointer to a muon
+  AliAODTrack* Mu(Int_t imu=0) const {return (imu==0||imu==1)&&(fMu[imu]!=0) ? (AliAODTrack*)fMu[imu].GetObject() : 0; } // Get a pointer to a muon
+
+  void SetMu(Int_t imu=0, AliAODTrack *mu=0);
+  void SetMuons(AliAODTrack *mu0=0, AliAODTrack *mu1=0);
+
+private:
+  Int_t CheckPointers() const;
+  void BookP();
+
+  // Data members
+  TRef fMu[2]; // Pointers to the reconstructed muons
+  TLorentzVector *fP; //! TLorentzVector of dimuon momentum (not stored into file)
+
+  // Useful constants
+  Double_t fMProton; //! Proton mass (not stored into file)
+
+  ClassDef(AliAODDimuon,1)  // AliAODDimuon track
+};
+
+#endif
index d0d1497..d159ace 100644 (file)
@@ -30,6 +30,7 @@
 #include "AliAODEvent.h"
 #include "AliAODHeader.h"
 #include "AliAODTrack.h"
+#include "AliAODDimuon.h"
 
 ClassImp(AliAODEvent)
 
@@ -45,7 +46,8 @@ ClassImp(AliAODEvent)
                                                      "phosCells",
                                                      "caloClusters",
                                                      "fmdClusters",
-                                                     "pmdClusters"
+                                                     "pmdClusters",
+                                                     "dimuons"
                                                      
 };
 //______________________________________________________________________________
@@ -65,7 +67,8 @@ AliAODEvent::AliAODEvent() :
   fPhosCells(0),
   fCaloClusters(0),
   fFmdClusters(0),
-  fPmdClusters(0)
+  fPmdClusters(0),
+  fDimuons(0)
 {
   // default constructor
 }
@@ -87,7 +90,8 @@ AliAODEvent::AliAODEvent(const AliAODEvent& aod):
   fPhosCells(new AliAODCaloCells(*aod.fPhosCells)),
   fCaloClusters(new TClonesArray(*aod.fCaloClusters)),
   fFmdClusters(new TClonesArray(*aod.fFmdClusters)),
-  fPmdClusters(new TClonesArray(*aod.fPmdClusters))
+  fPmdClusters(new TClonesArray(*aod.fPmdClusters)),
+  fDimuons(new TClonesArray(*aod.fDimuons))
 {
   // Copy constructor
   AddObject(fHeader);
@@ -102,6 +106,7 @@ AliAODEvent::AliAODEvent(const AliAODEvent& aod):
   AddObject(fCaloClusters);
   AddObject(fFmdClusters);
   AddObject(fPmdClusters);
+  AddObject(fDimuons);
   fConnected = aod.fConnected;
   GetStdContent();
   CreateStdFolders();
@@ -260,6 +265,7 @@ void AliAODEvent::CreateStdContent()
   AddObject(new TClonesArray("AliAODCaloCluster", 0));
   AddObject(new TClonesArray("AliAODFmdCluster", 0));
   AddObject(new TClonesArray("AliAODPmdCluster", 0));
+  AddObject(new TClonesArray("AliAODDimuon", 0));
   // set names
   SetStdNames();
 
@@ -341,6 +347,7 @@ void AliAODEvent::GetStdContent()
   fCaloClusters  = (TClonesArray*)fAODObjects->FindObject("caloClusters");
   fFmdClusters   = (TClonesArray*)fAODObjects->FindObject("fmdClusters");
   fPmdClusters   = (TClonesArray*)fAODObjects->FindObject("pmdClusters");
+  fDimuons       = (TClonesArray*)fAODObjects->FindObject("dimuons");
 }
 
 //______________________________________________________________________________
@@ -351,7 +358,8 @@ void AliAODEvent::ResetStd(Int_t trkArrSize,
                           Int_t jetSize, 
                           Int_t caloClusSize, 
                           Int_t fmdClusSize, 
-                          Int_t pmdClusSize
+                          Int_t pmdClusSize,
+                          Int_t dimuonArrSize
                           )
 {
   // deletes content of standard arrays and resets size 
@@ -387,6 +395,10 @@ void AliAODEvent::ResetStd(Int_t trkArrSize,
   fPmdClusters->Delete();
   if (pmdClusSize > fPmdClusters->GetSize()) 
     fPmdClusters->Expand(pmdClusSize);
+    
+  fDimuons->Delete();
+  if (dimuonArrSize > fDimuons->GetSize()) 
+    fDimuons->Expand(dimuonArrSize);
 
   // Reset the tracklets
   fTracklets->DeleteContainer();
@@ -410,6 +422,7 @@ void AliAODEvent::ClearStd()
   fCaloClusters  ->Delete();
   fFmdClusters   ->Clear();
   fPmdClusters   ->Clear();
+  fDimuons       ->Clear();
 }
 
 //_________________________________________________________________
index 24b386e..3fe0f68 100644 (file)
@@ -29,6 +29,7 @@
 #include "AliAODCaloCluster.h"
 #include "AliAODPmdCluster.h"
 #include "AliAODFmdCluster.h"
+#include "AliAODDimuon.h"
 
 class TTree;
 class TFolder;
@@ -48,6 +49,7 @@ class AliAODEvent : public AliVEvent {
                       kAODCaloClusters,
                       kAODFmdClusters,
                       kAODPmdClusters,
+                      kAODDimuons,
                       kAODListN
   };
 
@@ -175,6 +177,14 @@ class AliAODEvent : public AliVEvent {
   AliAODCaloCells *GetEMCALCells() const { return fEmcalCells; }
   AliAODCaloCells *GetPHOSCells() const { return fPhosCells; }
 
+  // -- Dimuons
+  TClonesArray *GetDimuons()              const { return fDimuons; }
+  Int_t         GetNDimuons()             const { return fDimuons->GetEntriesFast(); }
+  Int_t         GetNumberOfDimuons()      const { return GetNDimuons(); }
+  AliAODDimuon *GetDimuon(Int_t nDimu)    const { return (AliAODDimuon*)fDimuons->UncheckedAt(nDimu); }
+  Int_t         AddDimuon(const AliAODDimuon* dimu)
+    {new((*fDimuons)[fDimuons->GetEntriesFast()]) AliAODDimuon(*dimu); return fDimuons->GetEntriesFast()-1;}
+  
   // -- Services
   void    CreateStdContent();
   void    SetStdNames();
@@ -187,7 +197,8 @@ class AliAODEvent : public AliVEvent {
                   Int_t jetSize = 0, 
                   Int_t caloClusSize = 0, 
                   Int_t fmdClusSize = 0, 
-                  Int_t pmdClusSize = 0
+                  Int_t pmdClusSize = 0,
+                  Int_t dimuonArrsize =0
                   );
   void    ClearStd();
   void    Reset() {ClearStd();} 
@@ -215,10 +226,11 @@ class AliAODEvent : public AliVEvent {
   TClonesArray    *fCaloClusters; //! calorimeter clusters
   TClonesArray    *fFmdClusters;  //! FMDclusters
   TClonesArray    *fPmdClusters;  //! PMDclusters
+  TClonesArray    *fDimuons;      //! dimuons
 
   static const char* fAODListName[kAODListN]; //!
 
-  ClassDef(AliAODEvent, 5);
+  ClassDef(AliAODEvent, 6);
 };
 
 #endif
index b277c96..def0805 100644 (file)
@@ -62,6 +62,7 @@ AliAODHandler::AliAODHandler() :
     fNeedsFMDClustersBranchReplication(kFALSE),
     fNeedsCaloClustersBranchReplication(kFALSE),
     fNeedsMCParticlesBranchReplication(kFALSE),
+    fNeedsDimuonsBranchReplication(kFALSE),
     fAODIsReplicated(kFALSE),
     fAODEvent(NULL),
     fMCEventH(NULL),
@@ -89,6 +90,7 @@ AliAODHandler::AliAODHandler(const char* name, const char* title):
     fNeedsFMDClustersBranchReplication(kFALSE),
     fNeedsCaloClustersBranchReplication(kFALSE),
     fNeedsMCParticlesBranchReplication(kFALSE),
+    fNeedsDimuonsBranchReplication(kFALSE),
     fAODIsReplicated(kFALSE),
     fAODEvent(NULL),
     fMCEventH(NULL),
@@ -456,12 +458,14 @@ void AliAODHandler::CreateTree(Int_t flag)
     // Creates the AOD Tree
     fTreeA = new TTree("aodTree", "AliAOD tree");
     fTreeA->Branch(fAODEvent->GetList());
+    fTreeA->BranchRef();    
     if (flag == 0) fTreeA->SetDirectory(0);
 }
 
 //______________________________________________________________________________
 void AliAODHandler::FillTree()
 {
     // Fill the AOD Tree
   fTreeA->Fill();
 }
@@ -720,6 +724,7 @@ Bool_t AliAODExtension::Init(Option_t *option)
   }  
   fTreeE = new TTree("aodTree", "AliAOD tree");
   fTreeE->Branch(fAODEvent->GetList());
+  fTreeE->BranchRef();
   owd->cd();
   return kTRUE;
 }
index e9f155e..e209db2 100644 (file)
@@ -52,6 +52,7 @@ class AliAODHandler : public AliVEventHandler {
     virtual void         SetNeedsFMDClustersBranchReplication()  {fNeedsFMDClustersBranchReplication  = kTRUE;}
     virtual void         SetNeedsCaloClustersBranchReplication() {fNeedsCaloClustersBranchReplication = kTRUE;}
     virtual void         SetNeedsMCParticlesBranchReplication()  {fNeedsCaloClustersBranchReplication = kTRUE;}
+    virtual void         SetNeedsDimuonsBranchReplication()      {fNeedsDimuonsBranchReplication      = kTRUE;}
     virtual void         SetAODIsReplicated() {fAODIsReplicated = kTRUE;}
     //
     AliAODEvent*         GetAOD()  {return fAODEvent;}
@@ -78,6 +79,7 @@ class AliAODHandler : public AliVEventHandler {
     Bool_t               NeedsFMDClustersBranchReplication()  const {return  fNeedsFMDClustersBranchReplication;}
     Bool_t               NeedsCaloClustersBranchReplication() const {return  fNeedsCaloClustersBranchReplication;}
     Bool_t               NeedsMCParticlesBranchReplication()  const {return  fNeedsMCParticlesBranchReplication;}
+    Bool_t               NeedsDimuonsBranchReplication()      const {return  fNeedsDimuonsBranchReplication;}
     Bool_t               AODIsReplicated()                    const {return  fAODIsReplicated;}
     //
     void                 SetInputTree(TTree* /*tree*/) {;}
@@ -101,6 +103,7 @@ class AliAODHandler : public AliVEventHandler {
     Bool_t                   fNeedsFMDClustersBranchReplication;  // Flag for FMDClusters replication
     Bool_t                   fNeedsCaloClustersBranchReplication; // Flag for CaloClusters replication
     Bool_t                   fNeedsMCParticlesBranchReplication;  // Flag for MCParticles replication
+    Bool_t                   fNeedsDimuonsBranchReplication;      // Flag for Dimuons replication
     Bool_t                   fAODIsReplicated;                    // Flag true if replication as been executed
     AliAODEvent             *fAODEvent;               //! Pointer to the AOD event
     AliMCEventHandler       *fMCEventH;               //! Pointer to mc event handler needed not to depend on the manager
@@ -109,7 +112,7 @@ class AliAODHandler : public AliVEventHandler {
     TString                  fFileName;               //  Output file name
     TObjArray               *fExtensions;             //  List of extensions
     TObjArray               *fFilters;                // List of filtered AOD's
-    ClassDef(AliAODHandler, 5)
+    ClassDef(AliAODHandler, 6)
 };
 
 //-------------------------------------------------------------------------
index d153635..fb0f0dc 100644 (file)
@@ -7,7 +7,8 @@ SRCS = AliAODEvent.cxx AliAODHeader.cxx \
        AliAODHandler.cxx AliAODTracklets.cxx AliAODTagCreator.cxx \
        AliAODv0.cxx AliAODcascade.cxx AliAODCaloCells.cxx AliAODInputHandler.cxx \
        AliAODDiJet.cxx AliAODMCParticle.cxx AliAODMCHeader.cxx \
-        AliAODPWG4Particle.cxx AliAODPWG4ParticleCorrelation.cxx 
+        AliAODPWG4Particle.cxx AliAODPWG4ParticleCorrelation.cxx \
+       AliAODDimuon.cxx
 
 HDRS:= $(SRCS:.cxx=.h)