Analysis code updated
[u/mrichter/AliRoot.git] / 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()) 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 pca[3]={0};
213             Double_t pcaQuality=0;
214             TClonesArray *muons = new TClonesArray("AliAODTrack",2);
215             muons -> Add(recMuon1);
216             muons -> Add(recMuon2);
217             TLorentzVector kinem(0,0,0,0);
218             if (!AliMFTAnalysisTools::RefitAODDimuonWithCommonVertex(muons, pca, pcaQuality, kinem)) continue;
219             fHistPtDimuonsJpsi    -> Fill(kinem.Pt());
220             fHistMassDimuonsJpsi  -> Fill(kinem.M());
221             fHistResidualXVtxJpsi -> Fill(1.e4*(pca[0] - fPrimaryVertex[0]));
222             fHistResidualYVtxJpsi -> Fill(1.e4*(pca[1] - fPrimaryVertex[1]));
223             fHistResidualZVtxJpsi -> Fill(1.e4*(pca[2] - fPrimaryVertex[2]));
224           }
225         }
226       }
227    
228     }   // end of loop on 2nd muon
229
230   }   // end of loop on 1st muon
231  
232   //--------------------------------------------------------------------------------------------------
233    
234   PostData(1, fHistogramList);
235
236 }
237
238 //====================================================================================================================================================
239
240 void AliAnalysisTaskMFTExample::Terminate(Option_t *) {
241
242   // Draw result to the screen
243   // Called once at the end of the query
244   
245 }
246
247 //====================================================================================================================================================