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 const 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 const 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 const 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 const 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 const 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 const 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 const 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