PWG/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 /* $Id$ */
  17
  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
 140   Info("UserExec"," ") ;
 141   if (!fMCEvent) {
 142     Error("UserExec","NO MC EVENT FOUND!");
 143     return;
 144   }
 145
 146   fNevt++;
 147   Double_t containerInput[13] ;
 148   Double_t beamEnergy=7000;
 149
 150 ////////
 151 //// MC
 152 ////////
 153
 154   fCFManager->SetMCEventInfo(fMCEvent);
 155   AliStack* stack = fMCEvent->Stack();
 156
 157   // loop on the MC event
 158   for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
 159     AliMCParticle *mcPart  = (AliMCParticle*) fMCEvent->GetTrack(ipart);
 160
 161     TParticle *part = 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     TParticle *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     TParticle *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   // CHECK
 277   if ( ! esdH ) {
 278     AliError("Cannot get input event handler");
 279     return;
 280   }
 281
 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
 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);
 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);
 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);
 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 //________________________________________________________________________
 355 Float_t AliCFMuonResTask1::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
 356                                    Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const
 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 //________________________________________________________________________
 364 Float_t AliCFMuonResTask1::Rap(Float_t e, Float_t pz) const
 365 {
 366 // calculate rapidity
 367     Float_t rap;
 368     if(TMath::Abs(e-pz)>1e-7){
 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 //________________________________________________________________________
 378 Float_t AliCFMuonResTask1::Phideg(Float_t phi) const
 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 //________________________________________________________________________
 388 Double_t AliCFMuonResTask1::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 389 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
 390 Double_t Energy)
 391 {
 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
 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   //
 400   pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 401   pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 402   //
 403   // --- Get the muons parameters in the CM frame
 404   //
 405   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 406   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 407   //
 408   // --- Obtain the dimuon parameters in the CM frame
 409   //
 410   pDimuCM=pMu1CM+pMu2CM;
 411   //
 412   // --- Translate the dimuon parameters in the dimuon rest frame
 413   //
 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);
 423   //
 424   // --- Determine the z axis for the CS angle
 425   //
 426   zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 427   //
 428   // --- Determine the CS angle (angle between mu+ and the z axis defined above)
 429   Double_t cost;
 430
 431   if(charge1>0) cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());
 432   else cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());
 433   return cost;
 434 }
 435 //________________________________________________________________________
 436
 437 //________________________________________________________________________
 438 Double_t AliCFMuonResTask1::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 439 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 440 {
 441 // calculate costHE
 442   TLorentzVector pMu1CM, pMu2CM, pDimuCM; // In the CM frame
 443   TLorentzVector pMu1Dimu, pMu2Dimu; // In the dimuon rest frame
 444   TVector3 beta,zaxisCS;
 445   //
 446   // --- Get the muons parameters in the CM frame
 447   //
 448   pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 449   pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 450   //
 451   // --- Obtain the dimuon parameters in the CM frame
 452   //
 453   pDimuCM=pMu1CM+pMu2CM;
 454   //
 455   // --- Translate the muon parameters in the dimuon rest frame
 456   //
 457   beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 458   pMu1Dimu=pMu1CM;
 459   pMu2Dimu=pMu2CM;
 460   pMu1Dimu.Boost(beta);
 461   pMu2Dimu.Boost(beta);
 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;
 466   zaxis=(pDimuCM.Vect()).Unit();
 467
 468   // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
 469   Double_t cost;
 470   if(charge1>0) cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());
 471   else cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());
 472
 473   return cost;
 474 }
 475 //________________________________________________________________________
 476
 477 //________________________________________________________________________
 478 Double_t AliCFMuonResTask1::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 479 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
 480 Double_t Energy)
 481 {
 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
 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    //
 490    pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 491    pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 492    //
 493    // --- Get the muons parameters in the CM frame
 494    //
 495    pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 496    pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 497    //
 498    // --- Obtain the dimuon parameters in the CM frame
 499    //
 500    pDimuCM=pMu1CM+pMu2CM;
 501    //
 502    // --- Translate the dimuon parameters in the dimuon rest frame
 503    //
 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);
 513    //
 514    // --- Determine the z axis for the CS angle
 515    //
 516    zaxisCS=(((pProjDimu.Vect()).Unit())-((pTargDimu.Vect()).Unit())).Unit();
 517    yaxisCS=(((pProjDimu.Vect()).Unit()).Cross((pTargDimu.Vect()).Unit())).Unit();
 518    xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
 519
 520    Double_t phi;
 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)));
 523
 524    return phi;
 525 }
 526 //________________________________________________________________________
 527
 528 //________________________________________________________________________
 529 Double_t AliCFMuonResTask1::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
 530 Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t Energy)
 531 {
 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
 535    TVector3 beta,xaxis,yaxis,zaxis;
 536    Double_t mp=0.93827231;
 537
 538    //
 539    // --- Get the muons parameters in the CM frame
 540    //
 541    pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
 542    pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
 543    //
 544    // --- Obtain the dimuon parameters in the CM frame
 545    //
 546    pDimuCM=pMu1CM+pMu2CM;
 547    //
 548    // --- Translate the muon parameters in the dimuon rest frame
 549    //
 550    zaxis=(pDimuCM.Vect()).Unit();
 551
 552    beta=(-1./pDimuCM.E())*pDimuCM.Vect();
 553
 554    pProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 555    pTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp));
 556
 557    pProjDimu=pProjCM;
 558    pTargDimu=pTargCM;
 559
 560    pProjDimu.Boost(beta);
 561    pTargDimu.Boost(beta);
 562
 563    yaxis=((pProjDimu.Vect()).Cross(pTargDimu.Vect())).Unit();
 564    xaxis=(yaxis.Cross(zaxis)).Unit();
 565
 566    pMu1Dimu=pMu1CM;
 567    pMu2Dimu=pMu2CM;
 568    pMu1Dimu.Boost(beta);
 569    pMu2Dimu.Boost(beta);
 570
 571    Double_t phi;
 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));
 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));
 584     if ( ! cont ) {
 585       AliError("Cannot find container in file");
 586       return;
 587     }
 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