MFT/AliAnalysisTaskMFTExample.cxx

   1 #include "TH1D.h"
   2 #include "AliAODEvent.h"
   3 #include "AliAODMCParticle.h"
   4 #include "AliAODMCHeader.h"
   5 #include "AliAnalysisTaskMFTExample.h"
   6 #include "TDatabasePDG.h"
   7 #include "AliMUONTrackParam.h"
   8 #include "AliMUONTrackExtrap.h"
   9 #include "AliAODDimuon.h"
  10 #include "AliAODTrack.h"
  11 #include "AliAODHeader.h"
  12 #include "TClonesArray.h"
  13 #include "AliMFTConstants.h"
  14 #include "AliMFTSupport.h"
  15 #include "TRandom.h"
  16
  17 ClassImp(AliAnalysisTaskMFTExample)
  18
  19 //====================================================================================================================================================
  20
  21 AliAnalysisTaskMFTExample::AliAnalysisTaskMFTExample() :
  22   AliAnalysisTaskSE(),
  23   fVertexMode(0),
  24   fHistPtSingleMuons(0),
  25   fHistPtSingleMuonsFromJpsi(0),
  26   fHistPtDimuonsOS(0),
  27   fHistMassDimuonsOS(0),
  28   fHistPtDimuonsJpsi(0),
  29   fHistMassDimuonsJpsi(0),
  30   fHistResidualXVtxJpsi(0),
  31   fHistResidualYVtxJpsi(0),
  32   fHistResidualZVtxJpsi(0) {
  33
  34   // Default constructor
  35
  36   for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = 0;
  37   fVtxResolutionITS[0] = 5.e-4;
  38   fVtxResolutionITS[1] = 5.e-4;
  39   fVtxResolutionITS[2] = 4.e-4;
  40
  41 }
  42
  43 //====================================================================================================================================================
  44
  45 AliAnalysisTaskMFTExample::AliAnalysisTaskMFTExample(const Char_t *name) :
  46   AliAnalysisTaskSE(name),
  47   fVertexMode(0),
  48   fHistPtSingleMuons(0),
  49   fHistPtSingleMuonsFromJpsi(0),
  50   fHistPtDimuonsOS(0),
  51   fHistMassDimuonsOS(0),
  52   fHistPtDimuonsJpsi(0),
  53   fHistMassDimuonsJpsi(0),
  54   fHistResidualXVtxJpsi(0),
  55   fHistResidualYVtxJpsi(0),
  56   fHistResidualZVtxJpsi(0) {
  57
  58   // Constructor
  59
  60   for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = 0;
  61   fVtxResolutionITS[0] = 5.e-4;
  62   fVtxResolutionITS[1] = 5.e-4;
  63   fVtxResolutionITS[2] = 4.e-4;
  64
  65   // Define input and output slots here
  66   DefineOutput(1, TH1D::Class());
  67   DefineOutput(2, TH1D::Class());
  68   DefineOutput(3, TH1D::Class());
  69   DefineOutput(4, TH1D::Class());
  70   DefineOutput(5, TH1D::Class());
  71   DefineOutput(6, TH1D::Class());
  72   DefineOutput(7, TH1D::Class());
  73   DefineOutput(8, TH1D::Class());
  74   DefineOutput(9, TH1D::Class());
  75
  76 }
  77
  78 //====================================================================================================================================================
  79
  80 void AliAnalysisTaskMFTExample::UserCreateOutputObjects() {
  81
  82   // Called once
  83
  84   fHistPtSingleMuons = new TH1D("fHistPtSingleMuons","p_{T} of single muons (All)", 100, 0, 10);
  85   fHistPtSingleMuons -> SetXTitle("p_{T}  [GeV/c]");
  86   fHistPtSingleMuons -> Sumw2();
  87
  88   fHistPtSingleMuonsFromJpsi = new TH1D("fHistPtSingleMuonsFromJpsi", "p_{T} of single muons (from J/#psi)", 100, 0, 10);
  89   fHistPtSingleMuonsFromJpsi -> SetXTitle("p_{T}  [GeV/c]");
  90   fHistPtSingleMuonsFromJpsi -> Sumw2();
  91
  92   fHistMassDimuonsOS = new TH1D("fHistMassDimuonsOS", "Mass of OS dimuons", 500, 0, 10);
  93   fHistMassDimuonsOS -> SetXTitle("Mass  [GeV/c^{2}]");
  94   fHistMassDimuonsOS -> Sumw2();
  95
  96   fHistMassDimuonsJpsi = new TH1D("fHistMassDimuonsJpsi", "Mass of J/#psi dimuons", 500, 0, 10);
  97   fHistMassDimuonsJpsi -> SetXTitle("Mass  [GeV/c^{2}]");
  98   fHistMassDimuonsJpsi -> Sumw2();
  99
 100   fHistPtDimuonsOS = new TH1D("fHistPtDimuonsOS", "p_{T} of OS dimuons", 100, 0, 10);
 101   fHistPtDimuonsOS -> SetXTitle("p_{T}  [GeV/c]");
 102   fHistPtDimuonsOS -> Sumw2();
 103
 104   fHistPtDimuonsJpsi = new TH1D("fHistPtDimuonsJpsi", "p_{T} of J/#psi dimuons", 100, 0, 10);
 105   fHistPtDimuonsJpsi -> SetXTitle("p_{T}  [GeV/c]");
 106   fHistPtDimuonsJpsi -> Sumw2();
 107
 108   fHistResidualXVtxJpsi = new TH1D("fHistResidualXVtxJpsi", "J/#psi vertex residuals along x", 100, -100, 100);
 109   fHistResidualXVtxJpsi -> SetXTitle("Residual  [#mum]");
 110   fHistResidualXVtxJpsi -> Sumw2();
 111
 112   fHistResidualYVtxJpsi = new TH1D("fHistResidualYVtxJpsi", "J/#psi vertex residuals along y", 100, -100, 100);
 113   fHistResidualYVtxJpsi -> SetXTitle("Residual  [#mum]");
 114   fHistResidualYVtxJpsi -> Sumw2();
 115
 116   fHistResidualZVtxJpsi = new TH1D("fHistResidualZVtxJpsi", "J/#psi vertex residuals along z", 100, -1000, 1000);
 117   fHistResidualZVtxJpsi -> SetXTitle("Residual  [#mum]");
 118   fHistResidualZVtxJpsi -> Sumw2();
 119
 120   PostData(1, fHistPtSingleMuons);
 121   PostData(2, fHistPtSingleMuonsFromJpsi);
 122   PostData(3, fHistMassDimuonsOS);
 123   PostData(4, fHistMassDimuonsJpsi);
 124   PostData(5, fHistPtDimuonsOS);
 125   PostData(6, fHistPtDimuonsJpsi);
 126   PostData(7, fHistResidualXVtxJpsi);
 127   PostData(8, fHistResidualYVtxJpsi);
 128   PostData(9, fHistResidualZVtxJpsi);
 129
 130 }
 131
 132 //====================================================================================================================================================
 133
 134 void AliAnalysisTaskMFTExample::UserExec(Option_t *) {
 135
 136   // Main loop
 137   // Called for each event
 138
 139   AliAODEvent *aodEv = dynamic_cast<AliAODEvent*>(InputEvent());
 140   if (!aodEv) return;
 141
 142   aodEv->GetHeader()->InitMagneticField();
 143   AliMUONTrackExtrap::SetField();
 144
 145   TClonesArray *stackMC = (TClonesArray*) (aodEv->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
 146
 147   // Getting primary vertex, either from the generation or from the reconstruction -------------------
 148
 149   if (fVertexMode == kGenerated) {
 150     AliAODMCHeader *mcHeader = (AliAODMCHeader*) (aodEv->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 151     mcHeader->GetVertex(fPrimaryVertex);
 152     for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = gRandom->Gaus(fPrimaryVertex[i], fVtxResolutionITS[i]);
 153   }
 154   else if (fVertexMode == kReconstructed) {
 155     aodEv->GetPrimaryVertex()->GetXYZ(fPrimaryVertex);
 156   }
 157
 158   // -------------------------------------------------------------------------------------------------
 159
 160   AliAODTrack *recMuon1=0, *recMuon2=0;
 161   AliAODMCParticle *mcMuon1=0, *mcMuon2=0;
 162   AliAODMCParticle *mcMother1=0;
 163   Bool_t isMuon1FromJpsi=0;
 164
 165   // Loop over MUON+MFT muons
 166
 167   //--------------------------------------------------------------------------------------------------
 168
 169   for (Int_t iTrack=0; iTrack<aodEv->GetNumberOfTracks(); iTrack++) {
 170
 171     if (!(aodEv->GetTrack(iTrack)->IsMuonGlobalTrack())) continue;
 172
 173     recMuon1 = aodEv->GetTrack(iTrack);
 174     isMuon1FromJpsi = kFALSE;
 175
 176     // pt all muons
 177     fHistPtSingleMuons -> Fill(recMuon1->Pt());
 178
 179     // pt muons from J/psi
 180     if (recMuon1->GetLabel() >= 0) {
 181       mcMuon1 = (AliAODMCParticle*) stackMC->At(recMuon1->GetLabel());
 182       if (mcMuon1) {
 183         if (mcMuon1->GetMother()>=0) {
 184           mcMother1 = (AliAODMCParticle*) stackMC->At(mcMuon1->GetMother());
 185           if (mcMother1->PdgCode()==443) {
 186             isMuon1FromJpsi = kTRUE;
 187             fHistPtSingleMuonsFromJpsi -> Fill(recMuon1->Pt());
 188           }
 189         }
 190       }
 191     }
 192
 193     for (Int_t jTrack=0; jTrack<iTrack; jTrack++) {
 194
 195       if (!(aodEv->GetTrack(jTrack)->IsMuonGlobalTrack())) continue;
 196
 197       recMuon2 = aodEv->GetTrack(jTrack);
 198
 199       AliAODDimuon *dimuon = new AliAODDimuon(recMuon1, recMuon2);
 200       if (dimuon->Charge()) continue;
 201
 202       // pt and mass all OS dimuons
 203       fHistPtDimuonsOS   -> Fill(dimuon->Pt());
 204       fHistMassDimuonsOS -> Fill(dimuon->Mass());
 205
 206       // pt and mass J/psi dimuons
 207       if (!isMuon1FromJpsi) continue;
 208       if (recMuon2->GetLabel() >= 0) {
 209         mcMuon2 = (AliAODMCParticle*) stackMC->At(recMuon2->GetLabel());
 210         if (mcMuon2) {
 211           if (mcMuon2->GetMother() == mcMuon1->GetMother()) {
 212             Double_t vertex[3]={0};
 213             TLorentzVector kinem(0,0,0,0);
 214             if (!AliMFTSupport::RefitAODDimuonWithCommonVertex(dimuon, vertex, kinem)) continue;
 215             fHistPtDimuonsJpsi    -> Fill(kinem.Pt());
 216             fHistMassDimuonsJpsi  -> Fill(kinem.M());
 217             fHistResidualXVtxJpsi -> Fill(1.e4*(vertex[0] - fPrimaryVertex[0]));
 218             fHistResidualYVtxJpsi -> Fill(1.e4*(vertex[1] - fPrimaryVertex[1]));
 219             fHistResidualZVtxJpsi -> Fill(1.e4*(vertex[2] - fPrimaryVertex[2]));
 220           }
 221         }
 222       }
 223
 224     }   // end of loop on 2nd muon
 225
 226   }   // end of loop on 1st muon
 227
 228   //--------------------------------------------------------------------------------------------------
 229
 230   PostData(1, fHistPtSingleMuons);
 231   PostData(2, fHistPtSingleMuonsFromJpsi);
 232   PostData(3, fHistMassDimuonsOS);
 233   PostData(4, fHistMassDimuonsJpsi);
 234   PostData(5, fHistPtDimuonsOS);
 235   PostData(6, fHistPtDimuonsJpsi);
 236   PostData(7, fHistResidualXVtxJpsi);
 237   PostData(8, fHistResidualYVtxJpsi);
 238   PostData(9, fHistResidualZVtxJpsi);
 239
 240 }
 241
 242 //====================================================================================================================================================
 243
 244 void AliAnalysisTaskMFTExample::Terminate(Option_t *) {
 245
 246   // Draw result to the screen
 247   // Called once at the end of the query
 248
 249 }
 250
 251 //====================================================================================================================================================