]>
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 | ||
27de2dfb | 16 | /* $Id$ */ |
17 | ||
7cd8a4ce | 18 | //----------------------------------------------------------------------- |
19 | // Example of task (running locally, on AliEn and CAF), | |
20 | // which provides standard way of calculating acceptance and efficiency | |
21 | // between different steps of the procedure. | |
22 | // The ouptut of the task is a AliCFContainer from which the efficiencies | |
23 | // can be calculated | |
24 | //----------------------------------------------------------------------- | |
25 | // Author : R. Vernet, Consorzio Cometa - Catania (it) | |
26 | //----------------------------------------------------------------------- | |
27 | // Modification done by X. Lopez - LPC Clermont (fr) | |
28 | //----------------------------------------------------------------------- | |
29 | ||
30 | ||
31 | #ifndef ALICFMUONRESTASK1_CXX | |
32 | #define ALICFMUONRESTASK1_CXX | |
33 | ||
34 | #include "AliCFMuonResTask1.h" | |
35 | #include "AliHeader.h" | |
36 | #include "AliESDHeader.h" | |
37 | #include "AliStack.h" | |
38 | #include "TParticle.h" | |
39 | #include "TLorentzVector.h" | |
40 | #include "TH1I.h" | |
41 | #include "AliMCEvent.h" | |
42 | #include "AliAnalysisManager.h" | |
43 | #include "AliESDEvent.h" | |
44 | #include "AliAODEvent.h" | |
45 | #include "AliCFManager.h" | |
46 | #include "AliCFCutBase.h" | |
47 | #include "AliCFContainer.h" | |
48 | #include "TChain.h" | |
49 | #include "AliESDtrack.h" | |
50 | #include "AliLog.h" | |
51 | #include "AliESDMuonTrack.h" | |
52 | #include "AliESDtrack.h" | |
53 | #include "AliESDInputHandler.h" | |
54 | #include "TCanvas.h" | |
55 | ||
56 | ClassImp(AliCFMuonResTask1) | |
57 | ||
58 | //__________________________________________________________________________ | |
59 | AliCFMuonResTask1::AliCFMuonResTask1() : | |
60 | fReadAODData(0), | |
61 | fCFManager(0x0), | |
62 | fQAHistList(0x0), | |
63 | fHistEventsProcessed(0x0), | |
64 | fNevt(0) | |
65 | { | |
66 | // | |
67 | //Default ctor | |
68 | // | |
69 | } | |
70 | //___________________________________________________________________________ | |
71 | AliCFMuonResTask1::AliCFMuonResTask1(const Char_t* name) : | |
72 | AliAnalysisTaskSE(name), | |
73 | fReadAODData(0), | |
74 | fCFManager(0x0), | |
75 | fQAHistList(0x0), | |
76 | fHistEventsProcessed(0x0), | |
77 | fNevt(0) | |
78 | { | |
79 | // | |
80 | // Constructor. Initialization of Inputs and Outputs | |
81 | // | |
82 | Info("AliCFMuonResTask1","Calling Constructor"); | |
83 | ||
84 | fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ; | |
85 | ||
86 | DefineOutput(1,TH1I::Class()); | |
87 | DefineOutput(2,AliCFContainer::Class()); | |
88 | ||
89 | } | |
90 | ||
91 | //___________________________________________________________________________ | |
92 | AliCFMuonResTask1& AliCFMuonResTask1::operator=(const AliCFMuonResTask1& c) | |
93 | { | |
94 | // | |
95 | // Assignment operator | |
96 | // | |
97 | if (this!=&c) { | |
98 | AliAnalysisTaskSE::operator=(c) ; | |
99 | fReadAODData = c.fReadAODData ; | |
100 | fCFManager = c.fCFManager; | |
101 | fQAHistList = c.fQAHistList ; | |
102 | fHistEventsProcessed = c.fHistEventsProcessed; | |
103 | fNevt = c.fNevt ; | |
104 | } | |
105 | return *this; | |
106 | } | |
107 | ||
108 | //___________________________________________________________________________ | |
109 | AliCFMuonResTask1::AliCFMuonResTask1(const AliCFMuonResTask1& c) : | |
110 | AliAnalysisTaskSE(c), | |
111 | fReadAODData(c.fReadAODData), | |
112 | fCFManager(c.fCFManager), | |
113 | fQAHistList(c.fQAHistList), | |
114 | fHistEventsProcessed(c.fHistEventsProcessed), | |
115 | fNevt(c.fNevt) | |
116 | { | |
117 | // | |
118 | // Copy Constructor | |
119 | // | |
120 | } | |
121 | ||
122 | //___________________________________________________________________________ | |
123 | AliCFMuonResTask1::~AliCFMuonResTask1() { | |
124 | // | |
125 | //destructor | |
126 | // | |
127 | Info("~AliCFMuonResTask1","Calling Destructor"); | |
128 | if (fCFManager) delete fCFManager ; | |
129 | if (fHistEventsProcessed) delete fHistEventsProcessed ; | |
130 | if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;} | |
131 | } | |
132 | ||
133 | //_________________________________________________ | |
134 | void AliCFMuonResTask1::UserExec(Option_t *) | |
135 | { | |
136 | // | |
137 | // Main loop function | |
138 | // | |
139 | ||
dfbceb5b | 140 | Info("UserExec"," ") ; |
7cd8a4ce | 141 | if (!fMCEvent) { |
142 | Error("UserExec","NO MC EVENT FOUND!"); | |
143 | return; | |
144 | } | |
145 | ||
146 | fNevt++; | |
147 | Double_t containerInput[13] ; | |
77743d11 | 148 | Double_t beamEnergy=7000; |
7cd8a4ce | 149 | |
150 | //////// | |
151 | //// MC | |
152 | //////// | |
153 | ||
297f59ea | 154 | fCFManager->SetMCEventInfo(fMCEvent); |
7cd8a4ce | 155 | AliStack* stack = fMCEvent->Stack(); |
156 | ||
157 | // loop on the MC event | |
158 | for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) { | |
7aad0c47 | 159 | AliMCParticle *mcPart = (AliMCParticle*) fMCEvent->GetTrack(ipart); |
7cd8a4ce | 160 | |
161 | TParticle *part = mcPart->Particle(); | |
7cd8a4ce | 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); | |
3c9a0892 | 177 | TParticle *part0 = stack->Particle(p0); |
7cd8a4ce | 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); | |
3c9a0892 | 184 | TParticle *part1 = stack->Particle(p1); |
7cd8a4ce | 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()); | |
197fefa8 | 276 | // CHECK |
277 | if ( ! esdH ) { | |
278 | AliError("Cannot get input event handler"); | |
279 | return; | |
280 | } | |
281 | ||
7cd8a4ce | 282 | fESD = esdH->GetEvent(); |
283 | Int_t mult1 = fESD->GetNumberOfMuonTracks() ; | |
284 | ||
285 | ||
286 | for (Int_t j = 0; j<mult1; j++) { | |
287 | AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j))); | |
288 | Float_t zr1 = mu1->Charge(); | |
289 | // Select mu- | |
290 | if(zr1<0){ | |
291 | Float_t pxr1 = mu1->Px(); | |
292 | Float_t pyr1 = mu1->Py(); | |
293 | Float_t pzr1 = mu1->Pz(); | |
294 | Float_t phir1 = mu1->Phi(); | |
295 | phir1=Phideg(phir1); | |
296 | Float_t er1 = mu1->E(); | |
297 | Float_t rap1=Rap(er1,pzr1); | |
298 | // selection of the rapidity for first muon | |
299 | if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mu1)) continue; | |
300 | ||
301 | for (Int_t jj = 0; jj<mult1; jj++) { | |
302 | AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj))); | |
303 | Float_t zr2 = mu2->Charge(); | |
304 | // Select mu+ | |
305 | if(zr2>0 && jj !=j){ | |
306 | Float_t pxr2 = mu2->Px(); | |
307 | Float_t pyr2 = mu2->Py(); | |
308 | Float_t pzr2 = mu2->Pz(); | |
309 | Float_t phir2 = mu2->Phi(); | |
310 | phir2 = Phideg(phir2); | |
311 | Float_t er2 = mu2->E(); | |
312 | Float_t rap2=Rap(er2,pzr2); | |
313 | // selection of the rapidity for second muon | |
314 | if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mu2)) continue; | |
315 | ||
316 | Float_t prec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2)+ | |
317 | (pzr1+pzr2)*(pzr1+pzr2)); | |
318 | Float_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2)); | |
319 | Float_t raprec=Rap((er1+er2),(pzr1+pzr2)); | |
320 | Float_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2); | |
321 | ||
77743d11 | 322 | 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 | 323 | 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 | 324 | 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); |
325 | 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 | 326 | |
327 | containerInput[0] = fNevt+0.5 ; | |
328 | containerInput[1] = rap1 ; | |
329 | containerInput[2] = phir1 ; | |
330 | containerInput[3] = rap2 ; | |
331 | containerInput[4] = phir2 ; | |
332 | containerInput[5] = imassrec ; | |
333 | containerInput[6] = raprec ; | |
334 | containerInput[7] = ptrec; | |
335 | containerInput[8] = prec; | |
336 | containerInput[9] = costCSrec; | |
337 | containerInput[10] = phiCSrec; | |
338 | containerInput[11] = costHErec; | |
339 | containerInput[12] = phiHErec; | |
340 | ||
341 | // fill the container at the second step | |
342 | fCFManager->GetParticleContainer()->Fill(containerInput,1); | |
343 | ||
344 | } // mu+ Selection | |
345 | } // second mu Loop | |
346 | } // mu- Selection | |
347 | } | |
348 | ||
349 | // ---------- | |
350 | fHistEventsProcessed->Fill(0); | |
351 | PostData(1,fHistEventsProcessed) ; | |
352 | PostData(2,fCFManager->GetParticleContainer()) ; | |
353 | } | |
354 | //________________________________________________________________________ | |
eedcd663 | 355 | Float_t AliCFMuonResTask1::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1, |
77743d11 | 356 | Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const |
7cd8a4ce | 357 | { |
358 | // invariant mass calculation | |
359 | Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+ | |
360 | (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2))); | |
361 | return imassrec; | |
362 | } | |
363 | //________________________________________________________________________ | |
77743d11 | 364 | Float_t AliCFMuonResTask1::Rap(Float_t e, Float_t pz) const |
7cd8a4ce | 365 | { |
366 | // calculate rapidity | |
367 | Float_t rap; | |
77743d11 | 368 | if(TMath::Abs(e-pz)>1e-7){ |
7cd8a4ce | 369 | rap = 0.5*TMath::Log((e+pz)/(e-pz)); |
370 | return rap; | |
371 | } | |
372 | else{ | |
373 | rap = -200; | |
374 | return rap; | |
375 | } | |
376 | } | |
377 | //________________________________________________________________________ | |
77743d11 | 378 | Float_t AliCFMuonResTask1::Phideg(Float_t phi) const |
7cd8a4ce | 379 | { |
380 | // calculate Phi in range [-180,180] | |
381 | Float_t phideg; | |
382 | ||
383 | phideg = phi-TMath::Pi(); | |
384 | phideg = phideg*57.296; | |
385 | return phideg; | |
386 | } | |
387 | //________________________________________________________________________ | |
eedcd663 | 388 | Double_t AliCFMuonResTask1::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, |
7cd8a4ce | 389 | Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, |
390 | Double_t Energy) | |
391 | { | |
77743d11 | 392 | // calculate costCS |
393 | TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM. frame | |
394 | TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame | |
7cd8a4ce | 395 | TVector3 beta,zaxisCS; |
396 | Double_t mp=0.93827231; | |
397 | // | |
398 | // --- Fill the Lorentz vector for projectile and target in the CM frame | |
399 | // | |
77743d11 | 400 | pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp)); |
401 | pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp)); | |
7cd8a4ce | 402 | // |
403 | // --- Get the muons parameters in the CM frame | |
404 | // | |
77743d11 | 405 | pMu1CM.SetPxPyPzE(px1,py1,pz1,e1); |
406 | pMu2CM.SetPxPyPzE(px2,py2,pz2,e2); | |
7cd8a4ce | 407 | // |
408 | // --- Obtain the dimuon parameters in the CM frame | |
409 | // | |
77743d11 | 410 | pDimuCM=pMu1CM+pMu2CM; |
7cd8a4ce | 411 | // |
412 | // --- Translate the dimuon parameters in the dimuon rest frame | |
413 | // | |
77743d11 | 414 | beta=(-1./pDimuCM.E())*pDimuCM.Vect(); |
415 | pMu1Dimu=pMu1CM; | |
416 | pMu2Dimu=pMu2CM; | |
417 | pProjDimu=pProjCM; | |
418 | pTargDimu=pTargCM; | |
419 | pMu1Dimu.Boost(beta); | |
420 | pMu2Dimu.Boost(beta); | |
421 | pProjDimu.Boost(beta); | |
422 | pTargDimu.Boost(beta); | |
7cd8a4ce | 423 | // |
424 | // --- Determine the z axis for the CS angle | |
425 | // | |
77743d11 | 426 | zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit(); |
7cd8a4ce | 427 | // |
428 | // --- Determine the CS angle (angle between mu+ and the z axis defined above) | |
429 | Double_t cost; | |
430 | ||
77743d11 | 431 | if(charge1>0) cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit()); |
432 | else cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit()); | |
7cd8a4ce | 433 | return cost; |
434 | } | |
435 | //________________________________________________________________________ | |
436 | ||
437 | //________________________________________________________________________ | |
eedcd663 | 438 | Double_t AliCFMuonResTask1::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, |
7cd8a4ce | 439 | Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2) |
440 | { | |
77743d11 | 441 | // calculate costHE |
442 | TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame | |
443 | TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame | |
7cd8a4ce | 444 | TVector3 beta,zaxisCS; |
445 | // | |
446 | // --- Get the muons parameters in the CM frame | |
447 | // | |
77743d11 | 448 | pMu1CM.SetPxPyPzE(px1,py1,pz1,e1); |
449 | pMu2CM.SetPxPyPzE(px2,py2,pz2,e2); | |
7cd8a4ce | 450 | // |
451 | // --- Obtain the dimuon parameters in the CM frame | |
452 | // | |
77743d11 | 453 | pDimuCM=pMu1CM+pMu2CM; |
7cd8a4ce | 454 | // |
455 | // --- Translate the muon parameters in the dimuon rest frame | |
456 | // | |
77743d11 | 457 | beta=(-1./pDimuCM.E())*pDimuCM.Vect(); |
458 | pMu1Dimu=pMu1CM; | |
459 | pMu2Dimu=pMu2CM; | |
460 | pMu1Dimu.Boost(beta); | |
461 | pMu2Dimu.Boost(beta); | |
7cd8a4ce | 462 | // |
463 | // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system) | |
464 | // | |
465 | TVector3 zaxis; | |
77743d11 | 466 | zaxis=(pDimuCM.Vect()).Unit(); |
7cd8a4ce | 467 | |
468 | // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above) | |
469 | Double_t cost; | |
77743d11 | 470 | if(charge1>0) cost = zaxis.Dot((pMu1Dimu.Vect()).Unit()); |
471 | else cost = zaxis.Dot((pMu2Dimu.Vect()).Unit()); | |
7cd8a4ce | 472 | |
473 | return cost; | |
474 | } | |
475 | //________________________________________________________________________ | |
476 | ||
477 | //________________________________________________________________________ | |
eedcd663 | 478 | Double_t AliCFMuonResTask1::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, |
7cd8a4ce | 479 | Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, |
480 | Double_t Energy) | |
481 | { | |
77743d11 | 482 | // calculate phiCS |
483 | TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame | |
484 | TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame | |
7cd8a4ce | 485 | TVector3 beta,yaxisCS, xaxisCS, zaxisCS; |
486 | Double_t mp=0.93827231; | |
487 | // | |
488 | // --- Fill the Lorentz vector for projectile and target in the CM frame | |
489 | // | |
77743d11 | 490 | pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp)); |
491 | pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp)); | |
7cd8a4ce | 492 | // |
493 | // --- Get the muons parameters in the CM frame | |
494 | // | |
77743d11 | 495 | pMu1CM.SetPxPyPzE(px1,py1,pz1,e1); |
496 | pMu2CM.SetPxPyPzE(px2,py2,pz2,e2); | |
7cd8a4ce | 497 | // |
498 | // --- Obtain the dimuon parameters in the CM frame | |
499 | // | |
77743d11 | 500 | pDimuCM=pMu1CM+pMu2CM; |
7cd8a4ce | 501 | // |
502 | // --- Translate the dimuon parameters in the dimuon rest frame | |
503 | // | |
77743d11 | 504 | beta=(-1./pDimuCM.E())*pDimuCM.Vect(); |
505 | pMu1Dimu=pMu1CM; | |
506 | pMu2Dimu=pMu2CM; | |
507 | pProjDimu=pProjCM; | |
508 | pTargDimu=pTargCM; | |
509 | pMu1Dimu.Boost(beta); | |
510 | pMu2Dimu.Boost(beta); | |
511 | pProjDimu.Boost(beta); | |
512 | pTargDimu.Boost(beta); | |
7cd8a4ce | 513 | // |
514 | // --- Determine the z axis for the CS angle | |
515 | // | |
77743d11 | 516 | zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit(); |
517 | yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit(); | |
7cd8a4ce | 518 | xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit(); |
519 | ||
520 | Double_t phi; | |
77743d11 | 521 | if(charge1>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxisCS),((pMu1Dimu.Vect()).Dot(xaxisCS))); |
522 | else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxisCS),((pMu2Dimu.Vect()).Dot(xaxisCS))); | |
7cd8a4ce | 523 | |
524 | return phi; | |
525 | } | |
526 | //________________________________________________________________________ | |
527 | ||
528 | //________________________________________________________________________ | |
eedcd663 | 529 | Double_t AliCFMuonResTask1::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, |
7cd8a4ce | 530 | Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t Energy) |
531 | { | |
77743d11 | 532 | // calculate phiHE |
533 | TLorentzVector pMu1CM, pMu2CM, pProjCM, pTargCM, pDimuCM; // In the CM frame | |
534 | TLorentzVector pMu1Dimu, pMu2Dimu, pProjDimu, pTargDimu; // In the dimuon rest frame | |
7cd8a4ce | 535 | TVector3 beta,xaxis,yaxis,zaxis; |
536 | Double_t mp=0.93827231; | |
537 | ||
538 | // | |
539 | // --- Get the muons parameters in the CM frame | |
540 | // | |
77743d11 | 541 | pMu1CM.SetPxPyPzE(px1,py1,pz1,e1); |
542 | pMu2CM.SetPxPyPzE(px2,py2,pz2,e2); | |
7cd8a4ce | 543 | // |
544 | // --- Obtain the dimuon parameters in the CM frame | |
545 | // | |
77743d11 | 546 | pDimuCM=pMu1CM+pMu2CM; |
7cd8a4ce | 547 | // |
548 | // --- Translate the muon parameters in the dimuon rest frame | |
549 | // | |
77743d11 | 550 | zaxis=(pDimuCM.Vect()).Unit(); |
7cd8a4ce | 551 | |
77743d11 | 552 | beta=(-1./pDimuCM.E())*pDimuCM.Vect(); |
7cd8a4ce | 553 | |
77743d11 | 554 | pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp)); |
555 | pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp)); | |
7cd8a4ce | 556 | |
77743d11 | 557 | pProjDimu=pProjCM; |
558 | pTargDimu=pTargCM; | |
7cd8a4ce | 559 | |
77743d11 | 560 | pProjDimu.Boost(beta); |
561 | pTargDimu.Boost(beta); | |
7cd8a4ce | 562 | |
77743d11 | 563 | yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit(); |
7cd8a4ce | 564 | xaxis=(yaxis.Cross(zaxis)).Unit(); |
565 | ||
77743d11 | 566 | pMu1Dimu=pMu1CM; |
567 | pMu2Dimu=pMu2CM; | |
568 | pMu1Dimu.Boost(beta); | |
569 | pMu2Dimu.Boost(beta); | |
7cd8a4ce | 570 | |
571 | Double_t phi; | |
77743d11 | 572 | if(charge1>0) phi = TMath::ATan2((pMu1Dimu.Vect()).Dot(yaxis),(pMu1Dimu.Vect()).Dot(xaxis)); |
573 | else phi = TMath::ATan2((pMu2Dimu.Vect()).Dot(yaxis),(pMu2Dimu.Vect()).Dot(xaxis)); | |
7cd8a4ce | 574 | |
575 | return phi; | |
576 | } | |
577 | ||
578 | //________________________________________________________________________ | |
579 | void AliCFMuonResTask1::Terminate(Option_t *) | |
580 | { | |
581 | // draw result of the Invariant mass MC and ESD | |
582 | ||
583 | AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2)); | |
197fefa8 | 584 | if ( ! cont ) { |
585 | AliError("Cannot find container in file"); | |
586 | return; | |
587 | } | |
7cd8a4ce | 588 | |
589 | TH1D *kmass = cont->ShowProjection(5,0); | |
590 | TH1D *rmass = cont->ShowProjection(5,1); | |
591 | TH1D *kcostCS = cont->ShowProjection(9,0); | |
592 | TH1D *rcostCS = cont->ShowProjection(9,1); | |
593 | ||
594 | TCanvas *c1 = new TCanvas("AliCFMuonResTask1","JPSI MC & ESD",10,10,510,510); | |
595 | c1->Divide(2,2); | |
596 | c1->cd(1); | |
597 | kmass->Draw("HIST"); | |
598 | c1->cd(2); | |
599 | rmass->Draw("HIST"); | |
600 | c1->cd(3); | |
601 | kcostCS->Draw("HIST"); | |
602 | c1->cd(4); | |
603 | rcostCS->Draw("HIST"); | |
604 | } | |
605 | //________________________________________________________________________ | |
606 | ||
607 | #endif |