]>
Commit | Line | Data |
---|---|---|
7cd8a4ce | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | //----------------------------------------------------------------------- | |
17 | // Example of task (running locally, on AliEn and CAF), | |
18 | // which provides standard way of calculating acceptance and efficiency | |
19 | // between different steps of the procedure. | |
20 | // The ouptut of the task is a AliCFContainer from which the efficiencies | |
21 | // can be calculated | |
22 | //----------------------------------------------------------------------- | |
23 | // Author : R. Vernet, Consorzio Cometa - Catania (it) | |
24 | //----------------------------------------------------------------------- | |
25 | // Modification done by X. Lopez - LPC Clermont (fr) | |
26 | //----------------------------------------------------------------------- | |
27 | ||
28 | ||
29 | #ifndef ALICFMUONRESTASK1_CXX | |
30 | #define ALICFMUONRESTASK1_CXX | |
31 | ||
32 | #include "AliCFMuonResTask1.h" | |
33 | #include "AliHeader.h" | |
34 | #include "AliESDHeader.h" | |
35 | #include "AliStack.h" | |
36 | #include "TParticle.h" | |
37 | #include "TLorentzVector.h" | |
38 | #include "TH1I.h" | |
39 | #include "AliMCEvent.h" | |
40 | #include "AliAnalysisManager.h" | |
41 | #include "AliESDEvent.h" | |
42 | #include "AliAODEvent.h" | |
43 | #include "AliCFManager.h" | |
44 | #include "AliCFCutBase.h" | |
45 | #include "AliCFContainer.h" | |
46 | #include "TChain.h" | |
47 | #include "AliESDtrack.h" | |
48 | #include "AliLog.h" | |
49 | #include "AliESDMuonTrack.h" | |
50 | #include "AliESDtrack.h" | |
51 | #include "AliESDInputHandler.h" | |
52 | #include "TCanvas.h" | |
53 | ||
54 | ClassImp(AliCFMuonResTask1) | |
55 | ||
56 | //__________________________________________________________________________ | |
57 | AliCFMuonResTask1::AliCFMuonResTask1() : | |
58 | fReadAODData(0), | |
59 | fCFManager(0x0), | |
60 | fQAHistList(0x0), | |
61 | fHistEventsProcessed(0x0), | |
62 | fNevt(0) | |
63 | { | |
64 | // | |
65 | //Default ctor | |
66 | // | |
67 | } | |
68 | //___________________________________________________________________________ | |
69 | AliCFMuonResTask1::AliCFMuonResTask1(const Char_t* name) : | |
70 | AliAnalysisTaskSE(name), | |
71 | fReadAODData(0), | |
72 | fCFManager(0x0), | |
73 | fQAHistList(0x0), | |
74 | fHistEventsProcessed(0x0), | |
75 | fNevt(0) | |
76 | { | |
77 | // | |
78 | // Constructor. Initialization of Inputs and Outputs | |
79 | // | |
80 | Info("AliCFMuonResTask1","Calling Constructor"); | |
81 | ||
82 | fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ; | |
83 | ||
84 | DefineOutput(1,TH1I::Class()); | |
85 | DefineOutput(2,AliCFContainer::Class()); | |
86 | ||
87 | } | |
88 | ||
89 | //___________________________________________________________________________ | |
90 | AliCFMuonResTask1& AliCFMuonResTask1::operator=(const AliCFMuonResTask1& c) | |
91 | { | |
92 | // | |
93 | // Assignment operator | |
94 | // | |
95 | if (this!=&c) { | |
96 | AliAnalysisTaskSE::operator=(c) ; | |
97 | fReadAODData = c.fReadAODData ; | |
98 | fCFManager = c.fCFManager; | |
99 | fQAHistList = c.fQAHistList ; | |
100 | fHistEventsProcessed = c.fHistEventsProcessed; | |
101 | fNevt = c.fNevt ; | |
102 | } | |
103 | return *this; | |
104 | } | |
105 | ||
106 | //___________________________________________________________________________ | |
107 | AliCFMuonResTask1::AliCFMuonResTask1(const AliCFMuonResTask1& c) : | |
108 | AliAnalysisTaskSE(c), | |
109 | fReadAODData(c.fReadAODData), | |
110 | fCFManager(c.fCFManager), | |
111 | fQAHistList(c.fQAHistList), | |
112 | fHistEventsProcessed(c.fHistEventsProcessed), | |
113 | fNevt(c.fNevt) | |
114 | { | |
115 | // | |
116 | // Copy Constructor | |
117 | // | |
118 | } | |
119 | ||
120 | //___________________________________________________________________________ | |
121 | AliCFMuonResTask1::~AliCFMuonResTask1() { | |
122 | // | |
123 | //destructor | |
124 | // | |
125 | Info("~AliCFMuonResTask1","Calling Destructor"); | |
126 | if (fCFManager) delete fCFManager ; | |
127 | if (fHistEventsProcessed) delete fHistEventsProcessed ; | |
128 | if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;} | |
129 | } | |
130 | ||
131 | //_________________________________________________ | |
132 | void AliCFMuonResTask1::UserExec(Option_t *) | |
133 | { | |
134 | // | |
135 | // Main loop function | |
136 | // | |
137 | ||
dfbceb5b | 138 | Info("UserExec"," ") ; |
7cd8a4ce | 139 | if (!fMCEvent) { |
140 | Error("UserExec","NO MC EVENT FOUND!"); | |
141 | return; | |
142 | } | |
143 | ||
144 | fNevt++; | |
145 | Double_t containerInput[13] ; | |
77743d11 | 146 | Double_t beamEnergy=7000; |
7cd8a4ce | 147 | |
148 | //////// | |
149 | //// MC | |
150 | //////// | |
151 | ||
297f59ea | 152 | fCFManager->SetMCEventInfo(fMCEvent); |
7cd8a4ce | 153 | AliStack* stack = fMCEvent->Stack(); |
154 | ||
155 | // loop on the MC event | |
156 | for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) { | |
7aad0c47 | 157 | AliMCParticle *mcPart = (AliMCParticle*) fMCEvent->GetTrack(ipart); |
7cd8a4ce | 158 | |
159 | TParticle *part = mcPart->Particle(); | |
160 | TParticle *part0 = mcPart->Particle(); | |
161 | TParticle *part1 = mcPart->Particle(); | |
162 | ||
163 | // Selection of the resonance | |
164 | if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue; | |
165 | ||
166 | // Mother kinematics | |
167 | Float_t e = part->Energy(); | |
168 | Float_t pz = part->Pz(); | |
169 | //Float_t py = part->Py(); | |
170 | //Float_t px = part->Px(); | |
171 | Float_t rapmc = Rap(e,pz); | |
172 | //Float_t mass = part->GetCalcMass(); | |
173 | ||
174 | // Decays kinematics | |
175 | ||
176 | Int_t p0 = part->GetDaughter(0); | |
177 | part0 = stack->Particle(p0); | |
178 | // selection of the rapidity for first muon | |
179 | AliMCParticle *mcpart0 = new AliMCParticle(part0); | |
180 | if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart0)) continue; | |
181 | Int_t pdg0 = part0->GetPdgCode(); | |
182 | ||
183 | Int_t p1 = part->GetDaughter(1); | |
184 | part1 = stack->Particle(p1); | |
185 | // selection of the rapidity for second muon | |
186 | AliMCParticle *mcpart1 = new AliMCParticle(part1); | |
187 | if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart1)) continue; | |
188 | Int_t pdg1 = part1->GetPdgCode(); | |
189 | ||
190 | // 0 mu- 1 mu+ | |
191 | Float_t e0=-999, pz0=-999, py0=-999, px0=-999, phi0=-999, rapmc0=-999; | |
192 | Float_t e1=-999, pz1=-999, py1=-999, px1=-999, phi1=-999, rapmc1=-999; | |
193 | Double_t charge0=-999, charge1=-999; | |
194 | ||
195 | // ordering sign: first = mu- | |
196 | if(pdg0==13){ | |
197 | e0 = part0->Energy(); | |
198 | pz0 = part0->Pz(); | |
199 | py0 = part0->Py(); | |
200 | px0 = part0->Px(); | |
201 | phi0 = part0->Phi(); | |
202 | phi0 = Phideg(phi0); | |
203 | rapmc0=Rap(e0,pz0); | |
204 | charge0 = part0->GetPDG()->Charge()/3; | |
205 | ||
206 | e1 = part1->Energy(); | |
207 | pz1 = part1->Pz(); | |
208 | py1 = part1->Py(); | |
209 | px1 = part1->Px(); | |
210 | phi1 = part1->Phi(); | |
211 | phi1 = Phideg(phi1); | |
212 | rapmc1=Rap(e1,pz1); | |
213 | charge1 = part1->GetPDG()->Charge()/3; | |
214 | } | |
215 | else{ | |
216 | if(pdg0==-13){ | |
217 | e1 = part0->Energy(); | |
218 | pz1 = part0->Pz(); | |
219 | py1 = part0->Py(); | |
220 | px1 = part0->Px(); | |
221 | phi1 = part0->Phi(); | |
222 | phi1 = Phideg(phi1); | |
223 | rapmc1=Rap(e1,pz1); | |
224 | charge1 = part0->GetPDG()->Charge()/3; | |
225 | ||
226 | e0 = part1->Energy(); | |
227 | pz0 = part1->Pz(); | |
228 | py0 = part1->Py(); | |
229 | px0 = part1->Px(); | |
230 | phi0 = part1->Phi(); | |
231 | phi0 = Phideg(phi0); | |
232 | rapmc0=Rap(e0,pz0); | |
233 | charge0 = part1->GetPDG()->Charge()/3; | |
234 | } | |
235 | } | |
236 | ||
237 | if(pdg0==13 || pdg1==13) { | |
238 | ||
239 | Float_t pmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1)+ | |
240 | (pz0+pz1)*(pz0+pz1)); | |
241 | Float_t ptmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1)); | |
242 | Float_t imassmc = Imass(e0,px0,py0,pz0,e1,px1,py1,pz1); | |
243 | //Float_t rapmc_check=Rap((e0+e1),(pz0+pz1)); | |
244 | ||
77743d11 | 245 | Double_t costCS = CostCS((Double_t) px0,(Double_t) py0,(Double_t) pz0,(Double_t) e0, (Double_t)charge0,(Double_t) px1,(Double_t) py1,(Double_t) pz1,(Double_t) e1, beamEnergy); |
7cd8a4ce | 246 | Double_t costHE = CostHE((Double_t) px0,(Double_t) py0,(Double_t) pz0,(Double_t) e0, (Double_t)charge0,(Double_t) px1,(Double_t) py1,(Double_t) pz1,(Double_t) e1); |
77743d11 | 247 | Double_t phiCS = PhiCS((Double_t) px0,(Double_t) py0,(Double_t) pz0,(Double_t) e0, (Double_t)charge0,(Double_t) px1,(Double_t) py1,(Double_t) pz1,(Double_t) e1, beamEnergy); |
248 | Double_t phiHE = PhiHE((Double_t) px0,(Double_t) py0,(Double_t) pz0,(Double_t) e0, (Double_t)charge0,(Double_t) px1,(Double_t) py1,(Double_t) pz1,(Double_t) e1, beamEnergy); | |
7cd8a4ce | 249 | |
250 | containerInput[0] = fNevt+0.5 ; | |
251 | containerInput[1] = rapmc0 ; | |
252 | containerInput[2] = phi0 ; | |
253 | containerInput[3] = rapmc1 ; | |
254 | containerInput[4] = phi1 ; | |
255 | containerInput[5] = imassmc ; | |
256 | containerInput[6] = rapmc ; | |
257 | containerInput[7] = ptmc; | |
258 | containerInput[8] = pmc; | |
259 | containerInput[9] = costCS; | |
260 | containerInput[10] = phiCS; | |
261 | containerInput[11] = costHE; | |
262 | containerInput[12] = phiHE; | |
263 | ||
264 | // fill the container at the first step | |
265 | fCFManager->GetParticleContainer()->Fill(containerInput,0); | |
266 | } | |
267 | } | |
268 | ||
269 | //////// | |
270 | //// ESD | |
271 | //////// | |
272 | ||
273 | AliESDEvent *fESD; | |
274 | AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> | |
275 | (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); | |
276 | fESD = esdH->GetEvent(); | |
277 | Int_t mult1 = fESD->GetNumberOfMuonTracks() ; | |
278 | ||
279 | ||
280 | for (Int_t j = 0; j<mult1; j++) { | |
281 | AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j))); | |
282 | Float_t zr1 = mu1->Charge(); | |
283 | // Select mu- | |
284 | if(zr1<0){ | |
285 | Float_t pxr1 = mu1->Px(); | |
286 | Float_t pyr1 = mu1->Py(); | |
287 | Float_t pzr1 = mu1->Pz(); | |
288 | Float_t phir1 = mu1->Phi(); | |
289 | phir1=Phideg(phir1); | |
290 | Float_t er1 = mu1->E(); | |
291 | Float_t rap1=Rap(er1,pzr1); | |
292 | // selection of the rapidity for first muon | |
293 | if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mu1)) continue; | |
294 | ||
295 | for (Int_t jj = 0; jj<mult1; jj++) { | |
296 | AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj))); | |
297 | Float_t zr2 = mu2->Charge(); | |
298 | // Select mu+ | |
299 | if(zr2>0 && jj !=j){ | |
300 | Float_t pxr2 = mu2->Px(); | |
301 | Float_t pyr2 = mu2->Py(); | |
302 | Float_t pzr2 = mu2->Pz(); | |
303 | Float_t phir2 = mu2->Phi(); | |
304 | phir2 = Phideg(phir2); | |
305 | Float_t er2 = mu2->E(); | |
306 | Float_t rap2=Rap(er2,pzr2); | |
307 | // selection of the rapidity for second muon | |
308 | if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mu2)) continue; | |
309 | ||
310 | Float_t prec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2)+ | |
311 | (pzr1+pzr2)*(pzr1+pzr2)); | |
312 | Float_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2)); | |
313 | Float_t raprec=Rap((er1+er2),(pzr1+pzr2)); | |
314 | Float_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2); | |
315 | ||
77743d11 | 316 | Double_t costCSrec = CostCS((Double_t) pxr1,(Double_t) pyr1,(Double_t)pzr1,(Double_t) er1, (Double_t)zr1,(Double_t) pxr2,(Double_t) pyr2,(Double_t)pzr2,(Double_t) er2, beamEnergy); |
7cd8a4ce | 317 | Double_t costHErec = CostHE((Double_t) pxr1,(Double_t) pyr1,(Double_t)pzr1,(Double_t) er1, (Double_t)zr1,(Double_t) pxr2,(Double_t) pyr2,(Double_t)pzr2,(Double_t) er2); |
77743d11 | 318 | Double_t phiCSrec = PhiCS((Double_t) pxr1,(Double_t) pyr1,(Double_t)pzr1,(Double_t) er1, (Double_t)zr1,(Double_t) pxr2,(Double_t) pyr2,(Double_t)pzr2,(Double_t) er2, beamEnergy); |
319 | Double_t phiHErec = PhiHE((Double_t) pxr1,(Double_t) pyr1,(Double_t)pzr1,(Double_t) er1, (Double_t)zr1,(Double_t) pxr2,(Double_t) pyr2,(Double_t)pzr2,(Double_t) er2, beamEnergy); | |
7cd8a4ce | 320 | |
321 | containerInput[0] = fNevt+0.5 ; | |
322 | containerInput[1] = rap1 ; | |
323 | containerInput[2] = phir1 ; | |
324 | containerInput[3] = rap2 ; | |
325 | containerInput[4] = phir2 ; | |
326 | containerInput[5] = imassrec ; | |
327 | containerInput[6] = raprec ; | |
328 | containerInput[7] = ptrec; | |
329 | containerInput[8] = prec; | |
330 | containerInput[9] = costCSrec; | |
331 | containerInput[10] = phiCSrec; | |
332 | containerInput[11] = costHErec; | |
333 | containerInput[12] = phiHErec; | |
334 | ||
335 | // fill the container at the second step | |
336 | fCFManager->GetParticleContainer()->Fill(containerInput,1); | |
337 | ||
338 | } // mu+ Selection | |
339 | } // second mu Loop | |
340 | } // mu- Selection | |
341 | } | |
342 | ||
343 | // ---------- | |
344 | fHistEventsProcessed->Fill(0); | |
345 | PostData(1,fHistEventsProcessed) ; | |
346 | PostData(2,fCFManager->GetParticleContainer()) ; | |
347 | } | |
348 | //________________________________________________________________________ | |
eedcd663 | 349 | Float_t AliCFMuonResTask1::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1, |
77743d11 | 350 | Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const |
7cd8a4ce | 351 | { |
352 | // invariant mass calculation | |
353 | Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+ | |
354 | (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2))); | |
355 | return imassrec; | |
356 | } | |
357 | //________________________________________________________________________ | |
77743d11 | 358 | Float_t AliCFMuonResTask1::Rap(Float_t e, Float_t pz) const |
7cd8a4ce | 359 | { |
360 | // calculate rapidity | |
361 | Float_t rap; | |
77743d11 | 362 | if(TMath::Abs(e-pz)>1e-7){ |
7cd8a4ce | 363 | rap = 0.5*TMath::Log((e+pz)/(e-pz)); |
364 | return rap; | |
365 | } | |
366 | else{ | |
367 | rap = -200; | |
368 | return rap; | |
369 | } | |
370 | } | |
371 | //________________________________________________________________________ | |
77743d11 | 372 | Float_t AliCFMuonResTask1::Phideg(Float_t phi) const |
7cd8a4ce | 373 | { |
374 | // calculate Phi in range [-180,180] | |
375 | Float_t phideg; | |
376 | ||
377 | phideg = phi-TMath::Pi(); | |
378 | phideg = phideg*57.296; | |
379 | return phideg; | |
380 | } | |
381 | //________________________________________________________________________ | |
eedcd663 | 382 | Double_t AliCFMuonResTask1::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, |
7cd8a4ce | 383 | Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, |
384 | Double_t Energy) | |
385 | { | |
77743d11 | 386 | // calculate costCS |
387 | TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame | |
388 | TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame | |
7cd8a4ce | 389 | TVector3 beta,zaxisCS; |
390 | Double_t mp=0.93827231; | |
391 | // | |
392 | // --- Fill the Lorentz vector for projectile and target in the CM frame | |
393 | // | |
77743d11 | 394 | pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp)); |
395 | pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp)); | |
7cd8a4ce | 396 | // |
397 | // --- Get the muons parameters in the CM frame | |
398 | // | |
77743d11 | 399 | pMu1CM.SetPxPyPzE(px1,py1,pz1,e1); |
400 | pMu2CM.SetPxPyPzE(px2,py2,pz2,e2); | |
7cd8a4ce | 401 | // |
402 | // --- Obtain the dimuon parameters in the CM frame | |
403 | // | |
77743d11 | 404 | pDimuCM=pMu1CM+pMu2CM; |
7cd8a4ce | 405 | // |
406 | // --- Translate the dimuon parameters in the dimuon rest frame | |
407 | // | |
77743d11 | 408 | beta=(-1./pDimuCM.E())*pDimuCM.Vect(); |
409 | pMu1Dimu=pMu1CM; | |
410 | pMu2Dimu=pMu2CM; | |
411 | pProjDimu=pProjCM; | |
412 | pTargDimu=pTargCM; | |
413 | pMu1Dimu.Boost(beta); | |
414 | pMu2Dimu.Boost(beta); | |
415 | pProjDimu.Boost(beta); | |
416 | pTargDimu.Boost(beta); | |
7cd8a4ce | 417 | // |
418 | // --- Determine the z axis for the CS angle | |
419 | // | |
77743d11 | 420 | zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit(); |
7cd8a4ce | 421 | // |
422 | // --- Determine the CS angle (angle between mu+ and the z axis defined above) | |
423 | Double_t cost; | |
424 | ||
77743d11 | 425 | if(charge1>0) cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit()); |
426 | else cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit()); | |
7cd8a4ce | 427 | return cost; |
428 | } | |
429 | //________________________________________________________________________ | |
430 | ||
431 | //________________________________________________________________________ | |
eedcd663 | 432 | Double_t AliCFMuonResTask1::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, |
7cd8a4ce | 433 | Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2) |
434 | { | |
77743d11 | 435 | // calculate costHE |
436 | TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame | |
437 | TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame | |
7cd8a4ce | 438 | TVector3 beta,zaxisCS; |
439 | // | |
440 | // --- Get the muons parameters in the CM frame | |
441 | // | |
77743d11 | 442 | pMu1CM.SetPxPyPzE(px1,py1,pz1,e1); |
443 | pMu2CM.SetPxPyPzE(px2,py2,pz2,e2); | |
7cd8a4ce | 444 | // |
445 | // --- Obtain the dimuon parameters in the CM frame | |
446 | // | |
77743d11 | 447 | pDimuCM=pMu1CM+pMu2CM; |
7cd8a4ce | 448 | // |
449 | // --- Translate the muon parameters in the dimuon rest frame | |
450 | // | |
77743d11 | 451 | beta=(-1./pDimuCM.E())*pDimuCM.Vect(); |
452 | pMu1Dimu=pMu1CM; | |
453 | pMu2Dimu=pMu2CM; | |
454 | pMu1Dimu.Boost(beta); | |
455 | pMu2Dimu.Boost(beta); | |
7cd8a4ce | 456 | // |
457 | // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system) | |
458 | // | |
459 | TVector3 zaxis; | |
77743d11 | 460 | zaxis=(pDimuCM.Vect()).Unit(); |
7cd8a4ce | 461 | |
462 | // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above) | |
463 | Double_t cost; | |
77743d11 | 464 | if(charge1>0) cost = zaxis.Dot((pMu1Dimu.Vect()).Unit()); |
465 | else cost = zaxis.Dot((pMu2Dimu.Vect()).Unit()); | |
7cd8a4ce | 466 | |
467 | return cost; | |
468 | } | |
469 | //________________________________________________________________________ | |
470 | ||
471 | //________________________________________________________________________ | |
eedcd663 | 472 | Double_t AliCFMuonResTask1::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, |
7cd8a4ce | 473 | Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, |
474 | Double_t Energy) | |
475 | { | |
77743d11 | 476 | // calculate phiCS |
477 | TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame | |
478 | TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame | |
7cd8a4ce | 479 | TVector3 beta,yaxisCS, xaxisCS, zaxisCS; |
480 | Double_t mp=0.93827231; | |
481 | // | |
482 | // --- Fill the Lorentz vector for projectile and target in the CM frame | |
483 | // | |
77743d11 | 484 | pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp)); |
485 | pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp)); | |
7cd8a4ce | 486 | // |
487 | // --- Get the muons parameters in the CM frame | |
488 | // | |
77743d11 | 489 | pMu1CM.SetPxPyPzE(px1,py1,pz1,e1); |
490 | pMu2CM.SetPxPyPzE(px2,py2,pz2,e2); | |
7cd8a4ce | 491 | // |
492 | // --- Obtain the dimuon parameters in the CM frame | |
493 | // | |
77743d11 | 494 | pDimuCM=pMu1CM+pMu2CM; |
7cd8a4ce | 495 | // |
496 | // --- Translate the dimuon parameters in the dimuon rest frame | |
497 | // | |
77743d11 | 498 | beta=(-1./pDimuCM.E())*pDimuCM.Vect(); |
499 | pMu1Dimu=pMu1CM; | |
500 | pMu2Dimu=pMu2CM; | |
501 | pProjDimu=pProjCM; | |
502 | pTargDimu=pTargCM; | |
503 | pMu1Dimu.Boost(beta); | |
504 | pMu2Dimu.Boost(beta); | |
505 | pProjDimu.Boost(beta); | |
506 | pTargDimu.Boost(beta); | |
7cd8a4ce | 507 | // |
508 | // --- Determine the z axis for the CS angle | |
509 | // | |
77743d11 | 510 | zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit(); |
511 | yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit(); | |
7cd8a4ce | 512 | xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit(); |
513 | ||
514 | Double_t phi; | |
77743d11 | 515 | if(charge1>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS))); |
516 | else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS))); | |
7cd8a4ce | 517 | |
518 | return phi; | |
519 | } | |
520 | //________________________________________________________________________ | |
521 | ||
522 | //________________________________________________________________________ | |
eedcd663 | 523 | Double_t AliCFMuonResTask1::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, |
7cd8a4ce | 524 | Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t Energy) |
525 | { | |
77743d11 | 526 | // calculate phiHE |
527 | TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame | |
528 | TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame | |
7cd8a4ce | 529 | TVector3 beta,xaxis,yaxis,zaxis; |
530 | Double_t mp=0.93827231; | |
531 | ||
532 | // | |
533 | // --- Get the muons parameters in the CM frame | |
534 | // | |
77743d11 | 535 | pMu1CM.SetPxPyPzE(px1,py1,pz1,e1); |
536 | pMu2CM.SetPxPyPzE(px2,py2,pz2,e2); | |
7cd8a4ce | 537 | // |
538 | // --- Obtain the dimuon parameters in the CM frame | |
539 | // | |
77743d11 | 540 | pDimuCM=pMu1CM+pMu2CM; |
7cd8a4ce | 541 | // |
542 | // --- Translate the muon parameters in the dimuon rest frame | |
543 | // | |
77743d11 | 544 | zaxis=(pDimuCM.Vect()).Unit(); |
7cd8a4ce | 545 | |
77743d11 | 546 | beta=(-1./pDimuCM.E())*pDimuCM.Vect(); |
7cd8a4ce | 547 | |
77743d11 | 548 | pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp)); |
549 | pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp)); | |
7cd8a4ce | 550 | |
77743d11 | 551 | pProjDimu=pProjCM; |
552 | pTargDimu=pTargCM; | |
7cd8a4ce | 553 | |
77743d11 | 554 | pProjDimu.Boost(beta); |
555 | pTargDimu.Boost(beta); | |
7cd8a4ce | 556 | |
77743d11 | 557 | yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit(); |
7cd8a4ce | 558 | xaxis=(yaxis.Cross(zaxis)).Unit(); |
559 | ||
77743d11 | 560 | pMu1Dimu=pMu1CM; |
561 | pMu2Dimu=pMu2CM; | |
562 | pMu1Dimu.Boost(beta); | |
563 | pMu2Dimu.Boost(beta); | |
7cd8a4ce | 564 | |
565 | Double_t phi; | |
77743d11 | 566 | if(charge1>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis)); |
567 | else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis)); | |
7cd8a4ce | 568 | |
569 | return phi; | |
570 | } | |
571 | ||
572 | //________________________________________________________________________ | |
573 | void AliCFMuonResTask1::Terminate(Option_t *) | |
574 | { | |
575 | // draw result of the Invariant mass MC and ESD | |
576 | ||
577 | AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2)); | |
578 | ||
579 | TH1D *kmass = cont->ShowProjection(5,0); | |
580 | TH1D *rmass = cont->ShowProjection(5,1); | |
581 | TH1D *kcostCS = cont->ShowProjection(9,0); | |
582 | TH1D *rcostCS = cont->ShowProjection(9,1); | |
583 | ||
584 | TCanvas *c1 = new TCanvas("AliCFMuonResTask1","JPSI MC & ESD",10,10,510,510); | |
585 | c1->Divide(2,2); | |
586 | c1->cd(1); | |
587 | kmass->Draw("HIST"); | |
588 | c1->cd(2); | |
589 | rmass->Draw("HIST"); | |
590 | c1->cd(3); | |
591 | kcostCS->Draw("HIST"); | |
592 | c1->cd(4); | |
593 | rcostCS->Draw("HIST"); | |
594 | } | |
595 | //________________________________________________________________________ | |
596 | ||
597 | #endif |