QA mods to fix warnings
[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] ;
146 Double_t BeamEnergy=7000;
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
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);
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);
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);
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
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);
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);
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);
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,
7cd8a4ce 350 Float_t e2, Float_t px2, Float_t py2, Float_t pz2)
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//________________________________________________________________________
eedcd663 358Float_t AliCFMuonResTask1::Rap(Float_t e, Float_t pz)
7cd8a4ce 359{
360// calculate rapidity
361 Float_t rap;
362 if(e!=pz){
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//________________________________________________________________________
eedcd663 372Float_t AliCFMuonResTask1::Phideg(Float_t phi)
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{
386 TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM. frame
387 TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
388 TVector3 beta,zaxisCS;
389 Double_t mp=0.93827231;
390 //
391 // --- Fill the Lorentz vector for projectile and target in the CM frame
392 //
393 PProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
394 PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
395 //
396 // --- Get the muons parameters in the CM frame
397 //
398 PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
399 PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
400 //
401 // --- Obtain the dimuon parameters in the CM frame
402 //
403 PDimuCM=PMu1CM+PMu2CM;
404 //
405 // --- Translate the dimuon parameters in the dimuon rest frame
406 //
407 beta=(-1./PDimuCM.E())*PDimuCM.Vect();
408 PMu1Dimu=PMu1CM;
409 PMu2Dimu=PMu2CM;
410 PProjDimu=PProjCM;
411 PTargDimu=PTargCM;
412 PMu1Dimu.Boost(beta);
413 PMu2Dimu.Boost(beta);
414 PProjDimu.Boost(beta);
415 PTargDimu.Boost(beta);
416 //
417 // --- Determine the z axis for the CS angle
418 //
419 zaxisCS=(((PProjDimu.Vect()).Unit())-((PTargDimu.Vect()).Unit())).Unit();
420 //
421 // --- Determine the CS angle (angle between mu+ and the z axis defined above)
422 Double_t cost;
423
424 if(charge1>0) cost = zaxisCS.Dot((PMu1Dimu.Vect()).Unit());
425 else cost = zaxisCS.Dot((PMu2Dimu.Vect()).Unit());
426 return cost;
427}
428//________________________________________________________________________
429
430//________________________________________________________________________
eedcd663 431Double_t AliCFMuonResTask1::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
7cd8a4ce 432Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
433{
434 TLorentzVector PMu1CM, PMu2CM, PDimuCM; // In the CM frame
435 TLorentzVector PMu1Dimu, PMu2Dimu; // In the dimuon rest frame
436 TVector3 beta,zaxisCS;
437 //
438 // --- Get the muons parameters in the CM frame
439 //
440 PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
441 PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
442 //
443 // --- Obtain the dimuon parameters in the CM frame
444 //
445 PDimuCM=PMu1CM+PMu2CM;
446 //
447 // --- Translate the muon parameters in the dimuon rest frame
448 //
449 beta=(-1./PDimuCM.E())*PDimuCM.Vect();
450 PMu1Dimu=PMu1CM;
451 PMu2Dimu=PMu2CM;
452 PMu1Dimu.Boost(beta);
453 PMu2Dimu.Boost(beta);
454 //
455 // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
456 //
457 TVector3 zaxis;
458 zaxis=(PDimuCM.Vect()).Unit();
459
460 // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
461 Double_t cost;
462 if(charge1>0) cost = zaxis.Dot((PMu1Dimu.Vect()).Unit());
463 else cost = zaxis.Dot((PMu2Dimu.Vect()).Unit());
464
465 return cost;
466}
467//________________________________________________________________________
468
469//________________________________________________________________________
eedcd663 470Double_t AliCFMuonResTask1::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
7cd8a4ce 471Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
472Double_t Energy)
473{
474 TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM frame
475 TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
476 TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
477 Double_t mp=0.93827231;
478 //
479 // --- Fill the Lorentz vector for projectile and target in the CM frame
480 //
481 PProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
482 PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
483 //
484 // --- Get the muons parameters in the CM frame
485 //
486 PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
487 PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
488 //
489 // --- Obtain the dimuon parameters in the CM frame
490 //
491 PDimuCM=PMu1CM+PMu2CM;
492 //
493 // --- Translate the dimuon parameters in the dimuon rest frame
494 //
495 beta=(-1./PDimuCM.E())*PDimuCM.Vect();
496 PMu1Dimu=PMu1CM;
497 PMu2Dimu=PMu2CM;
498 PProjDimu=PProjCM;
499 PTargDimu=PTargCM;
500 PMu1Dimu.Boost(beta);
501 PMu2Dimu.Boost(beta);
502 PProjDimu.Boost(beta);
503 PTargDimu.Boost(beta);
504 //
505 // --- Determine the z axis for the CS angle
506 //
507 zaxisCS=(((PProjDimu.Vect()).Unit())-((PTargDimu.Vect()).Unit())).Unit();
508 yaxisCS=(((PProjDimu.Vect()).Unit()).Cross((PTargDimu.Vect()).Unit())).Unit();
509 xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
510
511 Double_t phi;
512 if(charge1>0) phi = TMath::ATan2((PMu1Dimu.Vect()).Dot(yaxisCS),((PMu1Dimu.Vect()).Dot(xaxisCS)));
513 else phi = TMath::ATan2((PMu2Dimu.Vect()).Dot(yaxisCS),((PMu2Dimu.Vect()).Dot(xaxisCS)));
514
515 return phi;
516}
517//________________________________________________________________________
518
519//________________________________________________________________________
eedcd663 520Double_t AliCFMuonResTask1::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
7cd8a4ce 521Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t Energy)
522{
523 TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM frame
524 TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
525 TVector3 beta,xaxis,yaxis,zaxis;
526 Double_t mp=0.93827231;
527
528 //
529 // --- Get the muons parameters in the CM frame
530 //
531 PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
532 PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
533 //
534 // --- Obtain the dimuon parameters in the CM frame
535 //
536 PDimuCM=PMu1CM+PMu2CM;
537 //
538 // --- Translate the muon parameters in the dimuon rest frame
539 //
540 zaxis=(PDimuCM.Vect()).Unit();
541
542 beta=(-1./PDimuCM.E())*PDimuCM.Vect();
543
544 PProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
545 PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
546
547 PProjDimu=PProjCM;
548 PTargDimu=PTargCM;
549
550 PProjDimu.Boost(beta);
551 PTargDimu.Boost(beta);
552
553 yaxis=((PProjDimu.Vect()).Cross(PTargDimu.Vect())).Unit();
554 xaxis=(yaxis.Cross(zaxis)).Unit();
555
556 PMu1Dimu=PMu1CM;
557 PMu2Dimu=PMu2CM;
558 PMu1Dimu.Boost(beta);
559 PMu2Dimu.Boost(beta);
560
561 Double_t phi;
562 if(charge1>0) phi = TMath::ATan2((PMu1Dimu.Vect()).Dot(yaxis),(PMu1Dimu.Vect()).Dot(xaxis));
563 else phi = TMath::ATan2((PMu2Dimu.Vect()).Dot(yaxis),(PMu2Dimu.Vect()).Dot(xaxis));
564
565 return phi;
566}
567
568//________________________________________________________________________
569void AliCFMuonResTask1::Terminate(Option_t *)
570{
571 // draw result of the Invariant mass MC and ESD
572
573 AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2));
574
575 TH1D *kmass = cont->ShowProjection(5,0);
576 TH1D *rmass = cont->ShowProjection(5,1);
577 TH1D *kcostCS = cont->ShowProjection(9,0);
578 TH1D *rcostCS = cont->ShowProjection(9,1);
579
580 TCanvas *c1 = new TCanvas("AliCFMuonResTask1","JPSI MC & ESD",10,10,510,510);
581 c1->Divide(2,2);
582 c1->cd(1);
583 kmass->Draw("HIST");
584 c1->cd(2);
585 rmass->Draw("HIST");
586 c1->cd(3);
587 kcostCS->Draw("HIST");
588 c1->cd(4);
589 rcostCS->Draw("HIST");
590}
591//________________________________________________________________________
592
593#endif