[u/mrichter/AliRoot.git] / PWG3 / muon / AliCFMuonResTask1.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//-----------------------------------------------------------------------
// Example of task (running locally, on AliEn and CAF),
// which provides standard way of calculating acceptance and efficiency
// between different steps of the procedure.
// The ouptut of the task is a AliCFContainer from which the efficiencies
// can be calculated
//-----------------------------------------------------------------------
// Author : R. Vernet, Consorzio Cometa - Catania (it)
//-----------------------------------------------------------------------
// Modification done by X. Lopez - LPC Clermont (fr)
//-----------------------------------------------------------------------


#ifndef ALICFMUONRESTASK1_CXX
#define ALICFMUONRESTASK1_CXX

#include "AliCFMuonResTask1.h"
#include "AliHeader.h"
#include "AliESDHeader.h"
#include "AliStack.h"
#include "TParticle.h"
#include "TLorentzVector.h"
#include "TH1I.h"
#include "AliMCEvent.h"
#include "AliAnalysisManager.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliCFManager.h"
#include "AliCFCutBase.h"
#include "AliCFContainer.h"
#include "TChain.h"
#include "AliESDtrack.h"
#include "AliLog.h"
#include "AliESDMuonTrack.h"
#include "AliESDtrack.h"
#include "AliESDInputHandler.h"
#include "TCanvas.h"

ClassImp(AliCFMuonResTask1)

//__________________________________________________________________________
AliCFMuonResTask1::AliCFMuonResTask1() :
  fReadAODData(0),
  fCFManager(0x0),
  fQAHistList(0x0),
  fHistEventsProcessed(0x0),
  fNevt(0)
{
  //
  //Default ctor
  //
}
//___________________________________________________________________________
AliCFMuonResTask1::AliCFMuonResTask1(const Char_t* name) :
  AliAnalysisTaskSE(name),
  fReadAODData(0),
  fCFManager(0x0),
  fQAHistList(0x0),
  fHistEventsProcessed(0x0),
  fNevt(0)
{
  //
  // Constructor. Initialization of Inputs and Outputs
  //
  Info("AliCFMuonResTask1","Calling Constructor");

  fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;

  DefineOutput(1,TH1I::Class());
  DefineOutput(2,AliCFContainer::Class());

}

//___________________________________________________________________________
AliCFMuonResTask1& AliCFMuonResTask1::operator=(const AliCFMuonResTask1& c) 
{
  //
  // Assignment operator
  //
  if (this!=&c) {
    AliAnalysisTaskSE::operator=(c) ;
    fReadAODData = c.fReadAODData ;
    fCFManager  = c.fCFManager;
    fQAHistList = c.fQAHistList ;
    fHistEventsProcessed = c.fHistEventsProcessed;
    fNevt = c.fNevt ;
  }
  return *this;
}

//___________________________________________________________________________
AliCFMuonResTask1::AliCFMuonResTask1(const AliCFMuonResTask1& c) :
  AliAnalysisTaskSE(c),
  fReadAODData(c.fReadAODData),
  fCFManager(c.fCFManager),
  fQAHistList(c.fQAHistList),
  fHistEventsProcessed(c.fHistEventsProcessed),
  fNevt(c.fNevt) 
{
  //
  // Copy Constructor
  //
}

//___________________________________________________________________________
AliCFMuonResTask1::~AliCFMuonResTask1() {
  //
  //destructor
  //
  Info("~AliCFMuonResTask1","Calling Destructor");
  if (fCFManager)           delete fCFManager ;
  if (fHistEventsProcessed) delete fHistEventsProcessed ;
  if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}
}

//_________________________________________________
void AliCFMuonResTask1::UserExec(Option_t *)
{
  //
  // Main loop function
  // 
  
  Info("UserExec","") ;
  if (!fMCEvent) {
    Error("UserExec","NO MC EVENT FOUND!");
    return;
  }

  fNevt++; 
  Double_t containerInput[13] ;
  Double_t BeamEnergy=7000;
 
////////
//// MC
//////// 

  fCFManager->SetEventInfo(fMCEvent);
  AliStack* stack = fMCEvent->Stack();

  // loop on the MC event
  for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) { 
    AliMCParticle *mcPart  = fMCEvent->GetTrack(ipart);
 
    TParticle *part = mcPart->Particle(); 
    TParticle *part0 = mcPart->Particle();
    TParticle *part1 = mcPart->Particle();
 
    // Selection of the resonance
    if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;

    // Mother kinematics
    Float_t e = part->Energy();
    Float_t pz = part->Pz();           
    //Float_t py = part->Py();
    //Float_t px = part->Px();
    Float_t rapmc = Rap(e,pz);
    //Float_t mass = part->GetCalcMass();

    // Decays kinematics

    Int_t p0 = part->GetDaughter(0);
    part0 = stack->Particle(p0); 
   // selection of the rapidity for first muon
    AliMCParticle *mcpart0 = new AliMCParticle(part0);
    if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart0)) continue;
    Int_t pdg0 = part0->GetPdgCode();

    Int_t p1 = part->GetDaughter(1);
    part1 = stack->Particle(p1);
   // selection of the rapidity for second muon
    AliMCParticle *mcpart1 = new AliMCParticle(part1);
    if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart1)) continue;
    Int_t pdg1 = part1->GetPdgCode();
 
   // 0 mu- 1 mu+
    Float_t e0=-999, pz0=-999, py0=-999, px0=-999, phi0=-999, rapmc0=-999;
    Float_t e1=-999, pz1=-999, py1=-999, px1=-999, phi1=-999, rapmc1=-999;
    Double_t charge0=-999, charge1=-999;

    // ordering sign: first = mu-
    if(pdg0==13){
	e0 = part0->Energy();
	pz0 = part0->Pz();
	py0 = part0->Py();
	px0 = part0->Px();
	phi0 = part0->Phi(); 
	phi0 = Phideg(phi0);    
	rapmc0=Rap(e0,pz0);
	charge0 = part0->GetPDG()->Charge()/3;

	e1 = part1->Energy();
	pz1 = part1->Pz();
	py1 = part1->Py();
	px1 = part1->Px();
	phi1 = part1->Phi();
	phi1 = Phideg(phi1);
	rapmc1=Rap(e1,pz1);
	charge1 = part1->GetPDG()->Charge()/3;
    }
    else{
	if(pdg0==-13){
	    e1 = part0->Energy();
	    pz1 = part0->Pz();
	    py1 = part0->Py();
	    px1 = part0->Px();
	    phi1 = part0->Phi();
	    phi1 = Phideg(phi1);    
	    rapmc1=Rap(e1,pz1);
	    charge1 = part0->GetPDG()->Charge()/3;
	    
	    e0 = part1->Energy();
	    pz0 = part1->Pz();
	    py0 = part1->Py();
	    px0 = part1->Px();
	    phi0 = part1->Phi();
	    phi0 = Phideg(phi0);
	    rapmc0=Rap(e0,pz0); 
	    charge0 = part1->GetPDG()->Charge()/3;
	}
    }
    
    if(pdg0==13 || pdg1==13) { 

	Float_t pmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1)+
				   (pz0+pz1)*(pz0+pz1));
	Float_t ptmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1));
	Float_t imassmc = Imass(e0,px0,py0,pz0,e1,px1,py1,pz1);
	//Float_t rapmc_check=Rap((e0+e1),(pz0+pz1));
	
	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);
	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);
	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);
	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);

	containerInput[0] = fNevt+0.5 ;   
	containerInput[1] = rapmc0 ;   
	containerInput[2] = phi0 ;   
	containerInput[3] = rapmc1 ;   
	containerInput[4] = phi1 ;   
	containerInput[5] = imassmc ;   
	containerInput[6] = rapmc ;   
	containerInput[7] = ptmc;
	containerInput[8] = pmc;	
	containerInput[9] = costCS;	
	containerInput[10] = phiCS;	
	containerInput[11] = costHE;	
	containerInput[12] = phiHE;	

	// fill the container at the first step
	fCFManager->GetParticleContainer()->Fill(containerInput,0);
    }
  }    

////////
//// ESD
////////
  
  AliESDEvent *fESD; 
  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>
      (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  fESD = esdH->GetEvent();
  Int_t mult1 = fESD->GetNumberOfMuonTracks() ;

 
    for (Int_t j = 0; j<mult1; j++) { 
	AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));
	Float_t zr1 = mu1->Charge();
// Select mu-
	if(zr1<0){
	    Float_t pxr1 = mu1->Px();
	    Float_t pyr1 = mu1->Py();
	    Float_t pzr1 = mu1->Pz();
	    Float_t phir1 = mu1->Phi(); 
	    phir1=Phideg(phir1);
 	    Float_t er1 = mu1->E();
	    Float_t rap1=Rap(er1,pzr1);
	    // selection of the rapidity for first muon
	    if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mu1)) continue;

	    for (Int_t jj = 0; jj<mult1; jj++) {
		AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));
		Float_t zr2 = mu2->Charge();
// Select mu+
		if(zr2>0 && jj !=j){
		    Float_t pxr2 = mu2->Px();
		    Float_t pyr2 = mu2->Py();
		    Float_t pzr2 = mu2->Pz();
		    Float_t phir2 = mu2->Phi();
		    phir2 = Phideg(phir2);
		    Float_t er2 = mu2->E();
		    Float_t rap2=Rap(er2,pzr2);
		    // selection of the rapidity for second muon
		    if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mu2)) continue;

		    Float_t prec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2)+
						(pzr1+pzr2)*(pzr1+pzr2));
		    Float_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
		    Float_t raprec=Rap((er1+er2),(pzr1+pzr2));
		    Float_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
		    
		    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);
		    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);
		    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);
		    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);

		    containerInput[0] = fNevt+0.5 ;   
		    containerInput[1] = rap1 ;   
		    containerInput[2] = phir1 ;   
		    containerInput[3] = rap2 ;
		    containerInput[4] = phir2 ;   
		    containerInput[5] = imassrec ;   
		    containerInput[6] = raprec ;   
		    containerInput[7] = ptrec;
		    containerInput[8] = prec;
		    containerInput[9] = costCSrec;	
		    containerInput[10] = phiCSrec;	
		    containerInput[11] = costHErec;	
		    containerInput[12] = phiHErec;	
		    
		    // fill the container at the second step
		    fCFManager->GetParticleContainer()->Fill(containerInput,1);

		}  // mu+ Selection
	    }      // second mu Loop
	}          // mu- Selection
    }        
 
//  ----------
  fHistEventsProcessed->Fill(0);
  PostData(1,fHistEventsProcessed) ;
  PostData(2,fCFManager->GetParticleContainer()) ;
}
//________________________________________________________________________
const Float_t AliCFMuonResTask1::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
				   Float_t e2, Float_t px2, Float_t py2, Float_t pz2) 
{
// invariant mass calculation
    Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+
                                    (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
    return imassrec;
}
//________________________________________________________________________
const Float_t AliCFMuonResTask1::Rap(Float_t e, Float_t pz) 
{
// calculate rapidity
    Float_t rap;
    if(e!=pz){
	rap = 0.5*TMath::Log((e+pz)/(e-pz));
	return rap;
    }
    else{
	rap = -200;
	return rap;
    }
}
//________________________________________________________________________
const Float_t AliCFMuonResTask1::Phideg(Float_t phi) 
{
// calculate Phi in range [-180,180] 
    Float_t phideg;
    
	phideg = phi-TMath::Pi();
	phideg = phideg*57.296;
	return phideg;
}
//________________________________________________________________________
const Double_t AliCFMuonResTask1::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
Double_t Energy)
{
  TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM. frame
  TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
  TVector3 beta,zaxisCS;
  Double_t mp=0.93827231;
  //
  // --- Fill the Lorentz vector for projectile and target in the CM frame
  //
  PProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp)); 
  PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp)); 
  //
  // --- Get the muons parameters in the CM frame 
  //
  PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
  PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
  //
  // --- Obtain the dimuon parameters in the CM frame
  //
  PDimuCM=PMu1CM+PMu2CM;
  //
  // --- Translate the dimuon parameters in the dimuon rest frame
  //
  beta=(-1./PDimuCM.E())*PDimuCM.Vect();
  PMu1Dimu=PMu1CM;
  PMu2Dimu=PMu2CM;
  PProjDimu=PProjCM;
  PTargDimu=PTargCM;
  PMu1Dimu.Boost(beta);
  PMu2Dimu.Boost(beta);
  PProjDimu.Boost(beta);
  PTargDimu.Boost(beta);
  //
  // --- Determine the z axis for the CS angle 
  //
  zaxisCS=(((PProjDimu.Vect()).Unit())-((PTargDimu.Vect()).Unit())).Unit();
  //				     
  // --- Determine the CS angle (angle between mu+ and the z axis defined above)
  Double_t cost;
  
  if(charge1>0) cost = zaxisCS.Dot((PMu1Dimu.Vect()).Unit());
  else cost = zaxisCS.Dot((PMu2Dimu.Vect()).Unit());
  return cost;
}
//________________________________________________________________________

//________________________________________________________________________
const Double_t AliCFMuonResTask1::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
{
  TLorentzVector PMu1CM, PMu2CM, PDimuCM; // In the CM frame 
  TLorentzVector PMu1Dimu, PMu2Dimu; // In the dimuon rest frame
  TVector3 beta,zaxisCS;
  //
  // --- Get the muons parameters in the CM frame
  //
  PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
  PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
  //
  // --- Obtain the dimuon parameters in the CM frame
  //
  PDimuCM=PMu1CM+PMu2CM;
  //
  // --- Translate the muon parameters in the dimuon rest frame
  //
  beta=(-1./PDimuCM.E())*PDimuCM.Vect();
  PMu1Dimu=PMu1CM;
  PMu2Dimu=PMu2CM;
  PMu1Dimu.Boost(beta);
  PMu2Dimu.Boost(beta);
  //
  // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system)
  //
  TVector3 zaxis;
  zaxis=(PDimuCM.Vect()).Unit();
  
  // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above)
  Double_t cost;
  if(charge1>0) cost = zaxis.Dot((PMu1Dimu.Vect()).Unit()); 
  else cost = zaxis.Dot((PMu2Dimu.Vect()).Unit()); 
  
  return cost;
}
//________________________________________________________________________

//________________________________________________________________________
const Double_t AliCFMuonResTask1::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2,
Double_t Energy)
{
   TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM frame
   TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
   TVector3 beta,yaxisCS, xaxisCS, zaxisCS;
   Double_t mp=0.93827231;
   //
   // --- Fill the Lorentz vector for projectile and target in the CM frame
   //
   PProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp)); 
   PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp)); 
   //
   // --- Get the muons parameters in the CM frame 
   //
   PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
   PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
   //
   // --- Obtain the dimuon parameters in the CM frame
   //
   PDimuCM=PMu1CM+PMu2CM;
   //
   // --- Translate the dimuon parameters in the dimuon rest frame
   //
   beta=(-1./PDimuCM.E())*PDimuCM.Vect();
   PMu1Dimu=PMu1CM;
   PMu2Dimu=PMu2CM;
   PProjDimu=PProjCM;
   PTargDimu=PTargCM;
   PMu1Dimu.Boost(beta);
   PMu2Dimu.Boost(beta);
   PProjDimu.Boost(beta);
   PTargDimu.Boost(beta);
   //
   // --- Determine the z axis for the CS angle 
   //
   zaxisCS=(((PProjDimu.Vect()).Unit())-((PTargDimu.Vect()).Unit())).Unit();
   yaxisCS=(((PProjDimu.Vect()).Unit()).Cross((PTargDimu.Vect()).Unit())).Unit();
   xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit();
 
   Double_t phi;
   if(charge1>0) phi = TMath::ATan2((PMu1Dimu.Vect()).Dot(yaxisCS),((PMu1Dimu.Vect()).Dot(xaxisCS)));
   else phi = TMath::ATan2((PMu2Dimu.Vect()).Dot(yaxisCS),((PMu2Dimu.Vect()).Dot(xaxisCS)));
     
   return phi;
}
//________________________________________________________________________

//________________________________________________________________________
const Double_t AliCFMuonResTask1::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1,
Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t Energy)
{
   TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM frame 
   TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
   TVector3 beta,xaxis,yaxis,zaxis;
   Double_t mp=0.93827231;
 
   //
   // --- Get the muons parameters in the CM frame
   //
   PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
   PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
   //
   // --- Obtain the dimuon parameters in the CM frame
   //
   PDimuCM=PMu1CM+PMu2CM;
   //
   // --- Translate the muon parameters in the dimuon rest frame
   // 
   zaxis=(PDimuCM.Vect()).Unit();
 
   beta=(-1./PDimuCM.E())*PDimuCM.Vect();
 
   PProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp));
   PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp)); 
 
   PProjDimu=PProjCM;
   PTargDimu=PTargCM;
 
   PProjDimu.Boost(beta);
   PTargDimu.Boost(beta);
   
   yaxis=((PProjDimu.Vect()).Cross(PTargDimu.Vect())).Unit();
   xaxis=(yaxis.Cross(zaxis)).Unit();
   
   PMu1Dimu=PMu1CM;
   PMu2Dimu=PMu2CM;
   PMu1Dimu.Boost(beta);
   PMu2Dimu.Boost(beta);
 
   Double_t phi;
   if(charge1>0) phi = TMath::ATan2((PMu1Dimu.Vect()).Dot(yaxis),(PMu1Dimu.Vect()).Dot(xaxis));
   else phi = TMath::ATan2((PMu2Dimu.Vect()).Dot(yaxis),(PMu2Dimu.Vect()).Dot(xaxis));
   
   return phi;
}

//________________________________________________________________________
void AliCFMuonResTask1::Terminate(Option_t *) 
{
  // draw result of the Invariant mass MC and ESD

    AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2));   

    TH1D *kmass = cont->ShowProjection(5,0);
    TH1D *rmass = cont->ShowProjection(5,1);
    TH1D *kcostCS = cont->ShowProjection(9,0);
    TH1D *rcostCS = cont->ShowProjection(9,1);

    TCanvas *c1 = new TCanvas("AliCFMuonResTask1","JPSI MC & ESD",10,10,510,510);
    c1->Divide(2,2);
    c1->cd(1);
    kmass->Draw("HIST");
    c1->cd(2);
    rmass->Draw("HIST");
    c1->cd(3);
    kcostCS->Draw("HIST");
    c1->cd(4);
    rcostCS->Draw("HIST");
}
//________________________________________________________________________

#endif
Commit	Line	Data
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
	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 = 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(EnergyEnergy+mpmp));
394	PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(EnergyEnergy+mpmp));
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(EnergyEnergy+mpmp));
482	PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(EnergyEnergy+mpmp));
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(EnergyEnergy+mpmp));
545	PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(EnergyEnergy+mpmp));
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