]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/muon/AliCFMuonResTask1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
o small fix
[u/mrichter/AliRoot.git] / PWG / 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
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
56ClassImp(AliCFMuonResTask1)
57
58//__________________________________________________________________________
59AliCFMuonResTask1::AliCFMuonResTask1() :
60 fReadAODData(0),
61 fCFManager(0x0),
62 fQAHistList(0x0),
63 fHistEventsProcessed(0x0),
64 fNevt(0)
65{
66 //
67 //Default ctor
68 //
69}
70//___________________________________________________________________________
71AliCFMuonResTask1::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//___________________________________________________________________________
92AliCFMuonResTask1& 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//___________________________________________________________________________
109AliCFMuonResTask1::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//___________________________________________________________________________
123AliCFMuonResTask1::~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//_________________________________________________
134void 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 355Float_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 364Float_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 378Float_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 388Double_t AliCFMuonResTask1::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
7cd8a4ce 389Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
390Double_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 438Double_t AliCFMuonResTask1::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
7cd8a4ce 439Double_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 478Double_t AliCFMuonResTask1::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
7cd8a4ce 479Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
480Double_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 529Double_t AliCFMuonResTask1::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
7cd8a4ce 530Double_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//________________________________________________________________________
579void 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