QA mods to fix warnings
[u/mrichter/AliRoot.git] / PWG3 / muon / AliCFMuonResTask1.cxx
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   
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++) { 
157     AliMCParticle *mcPart  = (AliMCParticle*) fMCEvent->GetTrack(ipart);
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 //________________________________________________________________________
349 Float_t AliCFMuonResTask1::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
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 //________________________________________________________________________
358 Float_t AliCFMuonResTask1::Rap(Float_t e, Float_t pz) 
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 //________________________________________________________________________
372 Float_t AliCFMuonResTask1::Phideg(Float_t phi) 
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 //________________________________________________________________________
382 Double_t AliCFMuonResTask1::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
383 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
384 Double_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 //________________________________________________________________________
431 Double_t AliCFMuonResTask1::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
432 Double_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 //________________________________________________________________________
470 Double_t AliCFMuonResTask1::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
471 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
472 Double_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 //________________________________________________________________________
520 Double_t AliCFMuonResTask1::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
521 Double_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 //________________________________________________________________________
569 void 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