PAR: includes from previously enabled PARfiles
[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()) {
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 //====================================================================================================================================================