AliCFMuonResTask1.C :
[u/mrichter/AliRoot.git] / CORRFW / test / muons / AliCFMuonResTask1.cxx
CommitLineData
7441207c 1/**************************************************************************\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3 * *\r
4 * Author: The ALICE Off-line Project. *\r
5 * Contributors are mentioned in the code where appropriate. *\r
6 * *\r
7 * Permission to use, copy, modify and distribute this software and its *\r
8 * documentation strictly for non-commercial purposes is hereby granted *\r
9 * without fee, provided that the above copyright notice appears in all *\r
10 * copies and that both the copyright notice and this permission notice *\r
11 * appear in the supporting documentation. The authors make no claims *\r
12 * about the suitability of this software for any purpose. It is *\r
13 * provided "as is" without express or implied warranty. *\r
14 **************************************************************************/\r
15\r
16//-----------------------------------------------------------------------\r
17// Example of task (running locally, on AliEn and CAF),\r
18// which provides standard way of calculating acceptance and efficiency\r
19// between different steps of the procedure.\r
20// The ouptut of the task is a AliCFContainer from which the efficiencies\r
21// can be calculated\r
22//-----------------------------------------------------------------------\r
23// Author : R. Vernet, Consorzio Cometa - Catania (it)\r
24//-----------------------------------------------------------------------\r
25// Modification done by X. Lopez - LPC Clermont (fr)\r
26//-----------------------------------------------------------------------\r
27\r
28\r
29#ifndef ALICFMUONRESTASK1_CXX\r
30#define ALICFMUONRESTASK1_CXX\r
31\r
32#include "AliCFMuonResTask1.h"\r
33#include "AliHeader.h"\r
34#include "AliESDHeader.h"\r
35#include "AliStack.h"\r
36#include "TParticle.h"\r
37#include "TH1I.h"\r
38#include "AliMCEvent.h"\r
39#include "AliAnalysisManager.h"\r
40#include "AliESDEvent.h"\r
41#include "AliAODEvent.h"\r
42#include "AliCFManager.h"\r
43#include "AliCFCutBase.h"\r
44#include "AliCFContainer.h"\r
45#include "TChain.h"\r
46#include "AliESDtrack.h"\r
47#include "AliLog.h"\r
48#include "AliESDMuonTrack.h"\r
49#include "AliESDtrack.h"\r
50#include "AliESDInputHandler.h"\r
51#include "TCanvas.h"\r
52\r
53ClassImp(AliCFMuonResTask1)\r
54\r
55//__________________________________________________________________________\r
56AliCFMuonResTask1::AliCFMuonResTask1() :\r
57 fReadAODData(0),\r
58 fCFManager(0x0),\r
59 fQAHistList(0x0),\r
60 fHistEventsProcessed(0x0),\r
61 fNevt(0)\r
62{\r
63 //\r
64 //Default ctor\r
65 //\r
66}\r
67//___________________________________________________________________________\r
68AliCFMuonResTask1::AliCFMuonResTask1(const Char_t* name) :\r
69 AliAnalysisTaskSE(name),\r
70 fReadAODData(0),\r
71 fCFManager(0x0),\r
72 fQAHistList(0x0),\r
73 fHistEventsProcessed(0x0),\r
74 fNevt(0)\r
75{\r
76 //\r
77 // Constructor. Initialization of Inputs and Outputs\r
78 //\r
79 Info("AliCFMuonResTask1","Calling Constructor");\r
80\r
81 fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;\r
82\r
83 DefineOutput(1,TH1I::Class());\r
84 DefineOutput(2,AliCFContainer::Class());\r
85\r
86}\r
87\r
88//___________________________________________________________________________\r
89AliCFMuonResTask1& AliCFMuonResTask1::operator=(const AliCFMuonResTask1& c) \r
90{\r
91 //\r
92 // Assignment operator\r
93 //\r
94 if (this!=&c) {\r
95 AliAnalysisTaskSE::operator=(c) ;\r
96 fReadAODData = c.fReadAODData ;\r
97 fCFManager = c.fCFManager;\r
98 fQAHistList = c.fQAHistList ;\r
99 fHistEventsProcessed = c.fHistEventsProcessed;\r
100 fNevt = c.fNevt ;\r
101 }\r
102 return *this;\r
103}\r
104\r
105//___________________________________________________________________________\r
106AliCFMuonResTask1::AliCFMuonResTask1(const AliCFMuonResTask1& c) :\r
107 AliAnalysisTaskSE(c),\r
108 fReadAODData(c.fReadAODData),\r
109 fCFManager(c.fCFManager),\r
110 fQAHistList(c.fQAHistList),\r
111 fHistEventsProcessed(c.fHistEventsProcessed),\r
112 fNevt(c.fNevt) \r
113{\r
114 //\r
115 // Copy Constructor\r
116 //\r
117}\r
118\r
119//___________________________________________________________________________\r
120AliCFMuonResTask1::~AliCFMuonResTask1() {\r
121 //\r
122 //destructor\r
123 //\r
124 Info("~AliCFMuonResTask1","Calling Destructor");\r
125 if (fCFManager) delete fCFManager ;\r
126 if (fHistEventsProcessed) delete fHistEventsProcessed ;\r
127 if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}\r
128}\r
129\r
130//_________________________________________________\r
131void AliCFMuonResTask1::UserExec(Option_t *)\r
132{\r
133 //\r
134 // Main loop function\r
135 // \r
136 \r
137 Info("UserExec","") ;\r
138 if (!fMCEvent) {\r
139 Error("UserExec","NO MC EVENT FOUND!");\r
140 return;\r
141 }\r
142\r
143 fNevt++; \r
144 Double_t containerInput[9] ;\r
145 \r
146////////\r
147//// MC\r
148//////// \r
149\r
150 fCFManager->SetEventInfo(fMCEvent);\r
151 AliStack* stack = fMCEvent->Stack();\r
152\r
153 // loop on the MC event\r
154 for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) { \r
155 AliMCParticle *mcPart = fMCEvent->GetTrack(ipart);\r
156 \r
157 TParticle *part = mcPart->Particle(); \r
158 TParticle *part0 = mcPart->Particle();\r
159 TParticle *part1 = mcPart->Particle();\r
160 \r
161 // Selection of the resonance\r
162 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;\r
163\r
164 // Mother kinematics\r
165 Float_t e = part->Energy();\r
166 Float_t pz = part->Pz(); \r
167 Float_t py = part->Py();\r
168 Float_t px = part->Px();\r
169 Float_t rapmc = Rap(e,pz);\r
170 Float_t mass = part->GetCalcMass();\r
171\r
172 // Decays kinematics\r
173\r
174 Int_t p0 = part->GetDaughter(0);\r
175 part0 = stack->Particle(p0); \r
176 // selection of the rapidity for first muon\r
177 AliMCParticle *mcpart0 = new AliMCParticle(part0);\r
0dbc8694 178 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart0)) continue;\r
7441207c 179 Int_t pdg0 = part0->GetPdgCode();\r
180\r
181 Int_t p1 = part->GetDaughter(1);\r
182 part1 = stack->Particle(p1);\r
183 // selection of the rapidity for second muon\r
184 AliMCParticle *mcpart1 = new AliMCParticle(part1);\r
0dbc8694 185 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart1)) continue;\r
7441207c 186 Int_t pdg1 = part1->GetPdgCode();\r
187 \r
188 // 0 mu- 1 mu+\r
189 Float_t e0, pz0, py0, px0, phi0, rapmc0;\r
190 Float_t e1, pz1, py1, px1, phi1, rapmc1;\r
191\r
192 // ordering sign: first = mu-\r
193 if(pdg0==13){\r
194 e0 = part0->Energy();\r
195 pz0 = part0->Pz();\r
196 py0 = part0->Py();\r
197 px0 = part0->Px();\r
0dbc8694 198 phi0 = part0->Phi(); \r
7441207c 199 phi0 = Phideg(phi0); \r
200 rapmc0=Rap(e0,pz0);\r
201\r
202 e1 = part1->Energy();\r
203 pz1 = part1->Pz();\r
204 py1 = part1->Py();\r
205 px1 = part1->Px();\r
206 phi1 = part1->Phi();\r
207 phi1 = Phideg(phi1);\r
208 rapmc1=Rap(e1,pz1);\r
209 }\r
210 else{\r
790b26af 211 if(pdg0==-13){\r
212 e1 = part0->Energy();\r
213 pz1 = part0->Pz();\r
214 py1 = part0->Py();\r
215 px1 = part0->Px();\r
216 phi1 = part0->Phi();\r
217 phi1 = Phideg(phi1); \r
218 rapmc1=Rap(e1,pz1);\r
219 \r
220 e0 = part1->Energy();\r
221 pz0 = part1->Pz();\r
222 py0 = part1->Py();\r
223 px0 = part1->Px();\r
224 phi0 = part1->Phi();\r
225 phi0 = Phideg(phi0);\r
226 rapmc0=Rap(e0,pz0); \r
227 }\r
7441207c 228 }\r
790b26af 229 \r
230 if(pdg0==13 || pdg1==13) { \r
7441207c 231\r
232 Float_t pmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1)+\r
233 (pz0+pz1)*(pz0+pz1));\r
234 Float_t ptmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1));\r
235 Float_t imassmc = Imass(e0,px0,py0,pz0,e1,px1,py1,pz1);\r
236 Float_t rapmc_check=Rap((e0+e1),(pz0+pz1));\r
237\r
238 containerInput[0] = fNevt ; \r
239 containerInput[1] = rapmc0 ; \r
240 containerInput[2] = phi0 ; \r
241 containerInput[3] = rapmc1 ; \r
242 containerInput[4] = phi1 ; \r
243 containerInput[5] = imassmc ; \r
244 containerInput[6] = rapmc ; \r
245 containerInput[7] = ptmc;\r
246 containerInput[8] = pmc; \r
247\r
248 // fill the container at the first step\r
249 fCFManager->GetParticleContainer()->Fill(containerInput,0);\r
790b26af 250 }\r
7441207c 251 } \r
252\r
253////////\r
254//// ESD\r
255////////\r
256 \r
257 AliESDEvent *fESD; \r
258 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>\r
259 (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
260 fESD = esdH->GetEvent();\r
7441207c 261 Int_t mult1 = fESD->GetNumberOfMuonTracks() ;\r
262\r
0dbc8694 263 \r
7441207c 264 for (Int_t j = 0; j<mult1; j++) { \r
265 AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));\r
266 Float_t zr1 = mu1->Charge();\r
267// Select mu-\r
268 if(zr1<0){\r
269 Float_t pxr1 = mu1->Px();\r
270 Float_t pyr1 = mu1->Py();\r
271 Float_t pzr1 = mu1->Pz();\r
0dbc8694 272 Float_t phir1 = mu1->Phi(); \r
273 phir1=Phideg(phir1);\r
7441207c 274 Float_t er1 = mu1->E();\r
275 Float_t rap1=Rap(er1,pzr1);\r
276 // selection of the rapidity for first muon\r
0dbc8694 277 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mu1)) continue;\r
7441207c 278\r
279 for (Int_t jj = 0; jj<mult1; jj++) {\r
280 AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));\r
281 Float_t zr2 = mu2->Charge();\r
282// Select mu+\r
283 if(zr2>0 && jj !=j){\r
284 Float_t pxr2 = mu2->Px();\r
285 Float_t pyr2 = mu2->Py();\r
286 Float_t pzr2 = mu2->Pz();\r
287 Float_t phir2 = mu2->Phi();\r
0dbc8694 288 phir2 = Phideg(phir2);\r
7441207c 289 Float_t er2 = mu2->E();\r
290 Float_t rap2=Rap(er2,pzr2);\r
291 // selection of the rapidity for second muon\r
0dbc8694 292 if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mu2)) continue;\r
7441207c 293\r
294 Float_t prec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2)+\r
295 (pzr1+pzr2)*(pzr1+pzr2));\r
296 Float_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));\r
297 Float_t raprec=Rap((er1+er2),(pzr1+pzr2));\r
298 Float_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);\r
299\r
300 containerInput[0] = fNevt ; \r
301 containerInput[1] = rap1 ; \r
302 containerInput[2] = phir1 ; \r
303 containerInput[3] = rap2 ;\r
304 containerInput[4] = phir2 ; \r
305 containerInput[5] = imassrec ; \r
306 containerInput[6] = raprec ; \r
307 containerInput[7] = ptrec;\r
308 containerInput[8] = prec;\r
309 \r
310 // fill the container at the second step\r
311 fCFManager->GetParticleContainer()->Fill(containerInput,1);\r
312\r
313 } // mu+ Selection\r
314 } // second mu Loop\r
315 } // mu- Selection\r
316 } \r
0dbc8694 317 \r
7441207c 318// ----------\r
319 fHistEventsProcessed->Fill(0);\r
320 PostData(1,fHistEventsProcessed) ;\r
321 PostData(2,fCFManager->GetParticleContainer()) ;\r
322}\r
323//________________________________________________________________________\r
324const Float_t AliCFMuonResTask1::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,\r
325 Float_t e2, Float_t px2, Float_t py2, Float_t pz2) \r
326{\r
327// invariant mass calculation\r
328 Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+\r
329 (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));\r
330 return imassrec;\r
331}\r
332//________________________________________________________________________\r
333const Float_t AliCFMuonResTask1::Rap(Float_t e, Float_t pz) \r
334{\r
335// calculate rapidity\r
336 Float_t rap;\r
337 if(e!=pz){\r
338 rap = 0.5*TMath::Log((e+pz)/(e-pz));\r
339 return rap;\r
340 }\r
341 else{\r
342 rap = -200;\r
343 return rap;\r
344 }\r
345}\r
346//________________________________________________________________________\r
347const Float_t AliCFMuonResTask1::Phideg(Float_t phi) \r
348{\r
0dbc8694 349// calculate Phi in range [-180,180] \r
7441207c 350 Float_t phideg;\r
351 \r
352 phideg = phi-TMath::Pi();\r
353 phideg = phideg*57.296;\r
354 return phideg;\r
355}\r
356//________________________________________________________________________\r
357void AliCFMuonResTask1::Terminate(Option_t *) \r
358{\r
359 // draw result of the Invariant mass MC and ESD\r
360\r
361 AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2)); \r
362\r
363 TH1D *kmass = cont->ShowProjection(5,0);\r
364 TH1D *rmass = cont->ShowProjection(5,1);\r
365\r
366 TCanvas *c1 = new TCanvas("AliCFMuonResTask1","JPSI MC & ESD",10,10,510,510);\r
367 c1->Divide(1,2);\r
368 c1->cd(1);\r
369 kmass->Draw("HIST");\r
370 c1->cd(2);\r
371 rmass->Draw("HIST");\r
372}\r
373//________________________________________________________________________\r
374\r
375#endif\r