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 "AliMFTAnalysisTools.h"
  15 #include "TRandom.h"
  16 #include "TList.h"
  17
  18 ClassImp(AliAnalysisTaskMFTExample)
  19
  20 //====================================================================================================================================================
  21
  22 AliAnalysisTaskMFTExample::AliAnalysisTaskMFTExample() :
  23   AliAnalysisTaskSE(),
  24   fVertexMode(0),
  25   fHistogramList(0),
  26   fHistPtSingleMuons(0),
  27   fHistPtSingleMuonsFromJpsi(0),
  28   fHistPtDimuonsOS(0),
  29   fHistMassDimuonsOS(0),
  30   fHistPtDimuonsJpsi(0),
  31   fHistMassDimuonsJpsi(0),
  32   fHistResidualXVtxJpsi(0),
  33   fHistResidualYVtxJpsi(0),
  34   fHistResidualZVtxJpsi(0) {
  35
  36   // Default constructor
  37
  38   for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = 0;
  39   fVtxResolutionITS[0] = 5.e-4;
  40   fVtxResolutionITS[1] = 5.e-4;
  41   fVtxResolutionITS[2] = 4.e-4;
  42
  43 }
  44
  45 //====================================================================================================================================================
  46
  47 AliAnalysisTaskMFTExample::AliAnalysisTaskMFTExample(const Char_t *name) :
  48   AliAnalysisTaskSE(name),
  49   fVertexMode(0),
  50   fHistogramList(0),
  51   fHistPtSingleMuons(0),
  52   fHistPtSingleMuonsFromJpsi(0),
  53   fHistPtDimuonsOS(0),
  54   fHistMassDimuonsOS(0),
  55   fHistPtDimuonsJpsi(0),
  56   fHistMassDimuonsJpsi(0),
  57   fHistResidualXVtxJpsi(0),
  58   fHistResidualYVtxJpsi(0),
  59   fHistResidualZVtxJpsi(0) {
  60
  61   // Constructor
  62
  63   for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = 0;
  64   fVtxResolutionITS[0] = 5.e-4;
  65   fVtxResolutionITS[1] = 5.e-4;
  66   fVtxResolutionITS[2] = 4.e-4;
  67
  68   // Define input and output slots here
  69   DefineOutput(1, TList::Class());
  70
  71 }
  72
  73 //====================================================================================================================================================
  74
  75 void AliAnalysisTaskMFTExample::UserCreateOutputObjects() {
  76
  77   // Called once
  78
  79   fHistogramList = new TList();
  80   fHistogramList->SetOwner(kTRUE);
  81
  82   fHistPtSingleMuons = new TH1D("fHistPtSingleMuons","p_{T} of single muons (All)", 100, 0, 10);
  83   fHistPtSingleMuons -> SetXTitle("p_{T}  [GeV/c]");
  84   fHistPtSingleMuons -> Sumw2();
  85
  86   fHistPtSingleMuonsFromJpsi = new TH1D("fHistPtSingleMuonsFromJpsi", "p_{T} of single muons (from J/#psi)", 100, 0, 10);
  87   fHistPtSingleMuonsFromJpsi -> SetXTitle("p_{T}  [GeV/c]");
  88   fHistPtSingleMuonsFromJpsi -> Sumw2();
  89
  90   fHistMassDimuonsOS = new TH1D("fHistMassDimuonsOS", "Mass of OS dimuons", 500, 0, 10);
  91   fHistMassDimuonsOS -> SetXTitle("Mass  [GeV/c^{2}]");
  92   fHistMassDimuonsOS -> Sumw2();
  93
  94   fHistMassDimuonsJpsi = new TH1D("fHistMassDimuonsJpsi", "Mass of J/#psi dimuons", 500, 0, 10);
  95   fHistMassDimuonsJpsi -> SetXTitle("Mass  [GeV/c^{2}]");
  96   fHistMassDimuonsJpsi -> Sumw2();
  97
  98   fHistPtDimuonsOS = new TH1D("fHistPtDimuonsOS", "p_{T} of OS dimuons", 100, 0, 10);
  99   fHistPtDimuonsOS -> SetXTitle("p_{T}  [GeV/c]");
 100   fHistPtDimuonsOS -> Sumw2();
 101
 102   fHistPtDimuonsJpsi = new TH1D("fHistPtDimuonsJpsi", "p_{T} of J/#psi dimuons", 100, 0, 10);
 103   fHistPtDimuonsJpsi -> SetXTitle("p_{T}  [GeV/c]");
 104   fHistPtDimuonsJpsi -> Sumw2();
 105
 106   fHistResidualXVtxJpsi = new TH1D("fHistResidualXVtxJpsi", "J/#psi vertex residuals along x", 100, -100, 100);
 107   fHistResidualXVtxJpsi -> SetXTitle("Residual  [#mum]");
 108   fHistResidualXVtxJpsi -> Sumw2();
 109
 110   fHistResidualYVtxJpsi = new TH1D("fHistResidualYVtxJpsi", "J/#psi vertex residuals along y", 100, -100, 100);
 111   fHistResidualYVtxJpsi -> SetXTitle("Residual  [#mum]");
 112   fHistResidualYVtxJpsi -> Sumw2();
 113
 114   fHistResidualZVtxJpsi = new TH1D("fHistResidualZVtxJpsi", "J/#psi vertex residuals along z", 100, -1000, 1000);
 115   fHistResidualZVtxJpsi -> SetXTitle("Residual  [#mum]");
 116   fHistResidualZVtxJpsi -> Sumw2();
 117
 118   fHistogramList->Add(fHistPtSingleMuons);
 119   fHistogramList->Add(fHistPtSingleMuonsFromJpsi);
 120   fHistogramList->Add(fHistMassDimuonsOS);
 121   fHistogramList->Add(fHistMassDimuonsJpsi);
 122   fHistogramList->Add(fHistPtDimuonsOS);
 123   fHistogramList->Add(fHistPtDimuonsJpsi);
 124   fHistogramList->Add(fHistResidualXVtxJpsi);
 125   fHistogramList->Add(fHistResidualYVtxJpsi);
 126   fHistogramList->Add(fHistResidualZVtxJpsi);
 127
 128   PostData(1, fHistogramList);
 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()) {
 201         delete dimuon;
 202         continue;
 203       }
 204
 205       // pt and mass all OS dimuons
 206       fHistPtDimuonsOS   -> Fill(dimuon->Pt());
 207       fHistMassDimuonsOS -> Fill(dimuon->Mass());
 208
 209       delete dimuon;
 210
 211       // pt and mass J/psi dimuons
 212       if (!isMuon1FromJpsi) continue;
 213       if (recMuon2->GetLabel() >= 0) {
 214         mcMuon2 = (AliAODMCParticle*) stackMC->At(recMuon2->GetLabel());
 215         if (mcMuon2) {
 216           if (mcMuon2->GetMother() == mcMuon1->GetMother()) {
 217             AliAODDimuon *dimuonJpsi = new AliAODDimuon;
 218             dimuonJpsi->SetMuons(recMuon1,recMuon2);
 219             Double_t pca[3]={0};
 220             Double_t pcaQuality=0;
 221             TLorentzVector kinem(0,0,0,0);
 222             if (!AliMFTAnalysisTools::CalculatePCA(dimuonJpsi, pca, pcaQuality, kinem)) {
 223               delete dimuonJpsi;
 224               continue;
 225             }
 226             fHistPtDimuonsJpsi    -> Fill(kinem.Pt());
 227             fHistMassDimuonsJpsi  -> Fill(kinem.M());
 228             fHistResidualXVtxJpsi -> Fill(1.e4*(pca[0] - fPrimaryVertex[0]));
 229             fHistResidualYVtxJpsi -> Fill(1.e4*(pca[1] - fPrimaryVertex[1]));
 230             fHistResidualZVtxJpsi -> Fill(1.e4*(pca[2] - fPrimaryVertex[2]));
 231             delete dimuonJpsi;
 232           }
 233         }
 234       }
 235
 236     }   // end of loop on 2nd muon
 237
 238   }   // end of loop on 1st muon
 239
 240   //--------------------------------------------------------------------------------------------------
 241
 242   PostData(1, fHistogramList);
 243
 244 }
 245
 246 //====================================================================================================================================================
 247
 248 void AliAnalysisTaskMFTExample::Terminate(Option_t *) {
 249
 250   // Draw result to the screen
 251   // Called once at the end of the query
 252
 253 }
 254
 255 //====================================================================================================================================================