]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/muon/AliCFMuonResTask1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixing some warnings and some coding convention rule violations (Roberta)
[u/mrichter/AliRoot.git] / PWG3 / muon / AliCFMuonResTask1.cxx
CommitLineData
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
54ClassImp(AliCFMuonResTask1)
55
56//__________________________________________________________________________
57AliCFMuonResTask1::AliCFMuonResTask1() :
58 fReadAODData(0),
59 fCFManager(0x0),
60 fQAHistList(0x0),
61 fHistEventsProcessed(0x0),
62 fNevt(0)
63{
64 //
65 //Default ctor
66 //
67}
68//___________________________________________________________________________
69AliCFMuonResTask1::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//___________________________________________________________________________
90AliCFMuonResTask1& 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//___________________________________________________________________________
107AliCFMuonResTask1::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//___________________________________________________________________________
121AliCFMuonResTask1::~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//_________________________________________________
132void AliCFMuonResTask1::UserExec(Option_t *)
133{
134 //
135 // Main loop function
136 //
137
138 Info("UserExec","") ;
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
152 fCFManager->SetEventInfo(fMCEvent);
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 349Float_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 358Float_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 372Float_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 382Double_t AliCFMuonResTask1::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
7cd8a4ce 383Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
384Double_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 432Double_t AliCFMuonResTask1::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
7cd8a4ce 433Double_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 472Double_t AliCFMuonResTask1::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
7cd8a4ce 473Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
474Double_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 523Double_t AliCFMuonResTask1::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
7cd8a4ce 524Double_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//________________________________________________________________________
573void 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