Fixing computing rule violations (Roberta)

diff --git a/PWG3/muon/AliAnalysisTaskMuonDistributions.cxx b/PWG3/muon/AliAnalysisTaskMuonDistributions.cxx
index 6fc32de980a2b1d3987e0011a4b02ab695909a59..f486a49772d47820674e2850a1d7d1be837143fc 100755 (executable)
--- a/PWG3/muon/AliAnalysisTaskMuonDistributions.cxx
+++ b/PWG3/muon/AliAnalysisTaskMuonDistributions.cxx
@@ -24,21 +24,15 @@
 //#ifndef ALIANALYSISTASKMUONDISTRIBUTIONS_CXX
 //#define ALIANALYSISTASKMUONDISTRIBUTIONS_CXX
 
-#include <TChain.h>
-#include <TTree.h>
 #include <TList.h>
 #include <TCanvas.h>
-#include <TClonesArray.h>
 #include <TH1.h>
-#include <TFile.h>
 #include <TStyle.h>
 #include <TLatex.h>
-#include <TPaveText.h>
 
 #include "AliAnalysisTaskMuonDistributions.h"
 #include "AliAnalysisTaskSE.h"
 #include "AliAnalysisManager.h"
-#include "AliAnalysisDataSlot.h"
 #include "AliESDEvent.h"
 #include "AliESD.h"
 #include "AliVEvent.h"
@@ -53,13 +47,11 @@
 #include "TParticle.h"
 #include "TLorentzVector.h"
 #include "AliESDMuonTrack.h"
-#include "AliESDtrack.h"
 #include "AliESDInputHandler.h"
 #include "AliAODEvent.h"
 #include "AliAODHeader.h"
 #include "AliAODHandler.h"
 #include "AliAODInputHandler.h"
-#include "AliMCEventHandler.h"
 
 ClassImp(AliAnalysisTaskMuonDistributions)
 
@@ -73,7 +65,7 @@ AliAnalysisTaskMuonDistributions::AliAnalysisTaskMuonDistributions() :
   fBckFitLimitMin(2.2),
   fBckFitLimitMax(2.85),
   fInvariantMassFit(kFALSE),
-  fAnalysisType(0x0),
+  fkAnalysisType(0x0),
   fOutput(0x0)
 {
 }
@@ -88,9 +80,10 @@ AliAnalysisTaskMuonDistributions::AliAnalysisTaskMuonDistributions(const Char_t*
   fBckFitLimitMin(2.2),
   fBckFitLimitMax(2.85),
   fInvariantMassFit(kFALSE),
-  fAnalysisType(0x0),
+  fkAnalysisType(0x0),
   fOutput(0x0)
 {
+  //
   // Constructor. Initialization of Inputs and Outputs
   //
   Info("AliAnalysisTaskMuonDistributions","Calling Constructor");
@@ -122,7 +115,7 @@ AliAnalysisTaskMuonDistributions::AliAnalysisTaskMuonDistributions(const AliAnal
   fBckFitLimitMin(c.fBckFitLimitMin),
   fBckFitLimitMax(c.fBckFitLimitMax),
   fInvariantMassFit(c.fInvariantMassFit),
-  fAnalysisType(c.fAnalysisType),
+  fkAnalysisType(c.fkAnalysisType),
   fOutput(c.fOutput)
  {
   //
@@ -140,7 +133,9 @@ AliAnalysisTaskMuonDistributions::~AliAnalysisTaskMuonDistributions() {
 
 //___________________________________________________________________________
 void AliAnalysisTaskMuonDistributions::UserCreateOutputObjects(){
-        
+ //
+ // output objects creation
+ //     
  fOutput = new TList();
  fOutput->SetOwner(); 
  //
@@ -150,11 +145,11 @@ void AliAnalysisTaskMuonDistributions::UserCreateOutputObjects(){
  //
  // dimuon histos
  //
- TH1D *hMass_Dimu   = new TH1D("hMass_Dimu","hMass_Dimu;M_{#mu#mu} (GeV/c^{2})",180,0,9.);     
- TH1D *hPt_Dimu  = new TH1D("hPt_Dimu","hPt_Dimu;p_{T} (GeV/c)",100,0,20);     
- TH1D *hRapidity_Dimu  = new TH1D("hRapidity_Dimu","hRapidity_Dimu;y",100,-5,-2);      
- TH1D *hCosThetaCS_Dimu  = new TH1D("hCosThetaCS_Dimu","hCosThetaCS_Dimu;cos#theta_{CS}",100,-1.,1.);  
- TH1D *hCosThetaHE_Dimu  = new TH1D("hCosThetaHE_Dimu","hCosThetaHE_Dimu;cos#theta_{HE}",100,-1.,1.);  
+ TH1D *hMassDimu   = new TH1D("hMassDimu","hMassDimu;M_{#mu#mu} (GeV/c^{2})",180,0,9.);        
+ TH1D *hPtDimu  = new TH1D("hPtDimu","hPtDimu;p_{T} (GeV/c)",100,0,20);        
+ TH1D *hRapidityDimu  = new TH1D("hRapidityDimu","hRapidityDimu;y",100,-5,-2); 
+ TH1D *hCosThetaCSDimu  = new TH1D("hCosThetaCSDimu","hCosThetaCSDimu;cos#theta_{CS}",100,-1.,1.);     
+ TH1D *hCosThetaHEDimu  = new TH1D("hCosThetaHEDimu","hCosThetaHEDimu;cos#theta_{HE}",100,-1.,1.);     
  //
  // muon histos
  //
@@ -163,11 +158,11 @@ void AliAnalysisTaskMuonDistributions::UserCreateOutputObjects(){
  TH1D *hRapidity  = new TH1D("hRapidity","hRapidity;y",100,-5,-2);     
        
  fOutput->Add(hNumberMuonTracks);      
- fOutput->Add(hMass_Dimu);     
- fOutput->Add(hPt_Dimu);       
- fOutput->Add(hRapidity_Dimu);         
- fOutput->Add(hCosThetaCS_Dimu);       
- fOutput->Add(hCosThetaHE_Dimu);       
+ fOutput->Add(hMassDimu);      
+ fOutput->Add(hPtDimu);        
+ fOutput->Add(hRapidityDimu);  
+ fOutput->Add(hCosThetaCSDimu);        
+ fOutput->Add(hCosThetaHEDimu);        
  fOutput->Add(hP);     
  fOutput->Add(hPt);    
  fOutput->Add(hRapidity);      
@@ -177,84 +172,87 @@ void AliAnalysisTaskMuonDistributions::UserCreateOutputObjects(){
 //_________________________________________________
 void AliAnalysisTaskMuonDistributions::UserExec(Option_t *)
 {
+//
+// Execute analysis for current event
+//
   AliESDEvent *esd=0x0;
   AliAODEvent *aod=0x0;
   
-  if(strcmp(fAnalysisType,"ESD")==0){
+  if(strcmp(fkAnalysisType,"ESD")==0){
     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>
         (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
     esd = esdH->GetEvent();
-  } else if(strcmp(fAnalysisType,"AOD")==0){
+  } else if(strcmp(fkAnalysisType,"AOD")==0){
     aod = dynamic_cast<AliAODEvent*> (InputEvent());
   }
   
   Int_t ntracks=-999;
-  if(strcmp(fAnalysisType,"ESD")==0) ntracks=esd->GetNumberOfMuonTracks();
-  else if(strcmp(fAnalysisType,"AOD")==0) ntracks=aod->GetNumberOfTracks();
+  if(strcmp(fkAnalysisType,"ESD")==0) ntracks=esd->GetNumberOfMuonTracks();
+  else if(strcmp(fkAnalysisType,"AOD")==0) ntracks=aod->GetNumberOfTracks();
   Int_t nmuontracks=0;
   
   for (Int_t j = 0; j<ntracks; j++) {
-    Float_t px_mu1=-999; Float_t py_mu1=-999; Float_t pz_mu1=-999; Float_t pt_mu1=-999; Float_t p_mu1=-999;
-    Float_t e_mu1=-999; Float_t rapidity_mu1=-999;  Float_t charge_mu1=-999; 
-    if(strcmp(fAnalysisType,"ESD")==0){ 
+    Float_t pxmu1=-999; Float_t pymu1=-999; Float_t pzmu1=-999; Float_t ptmu1=-999; Float_t pmu1=-999;
+    Float_t emu1=-999; Float_t rapiditymu1=-999;  Float_t chargemu1=-999; 
+    if(strcmp(fkAnalysisType,"ESD")==0){ 
       AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(esd->GetMuonTrack(j)));
       if (!mu1->ContainTrackerData()) continue;
-      charge_mu1 = mu1->Charge();
-      px_mu1 = mu1->Px();
-      py_mu1 = mu1->Py();
-      pz_mu1 = mu1->Pz();
-      e_mu1 = mu1->E();
-      p_mu1 = mu1->P();
-      pt_mu1 = mu1->Pt();
-      rapidity_mu1 = Rapidity(e_mu1,pz_mu1);
-    } else if(strcmp(fAnalysisType,"AOD")==0){
+      chargemu1 = mu1->Charge();
+      pxmu1 = mu1->Px();
+      pymu1 = mu1->Py();
+      pzmu1 = mu1->Pz();
+      emu1 = mu1->E();
+      pmu1 = mu1->P();
+      ptmu1 = mu1->Pt();
+      rapiditymu1 = Rapidity(emu1,pzmu1);
+    } else if(strcmp(fkAnalysisType,"AOD")==0){
       AliAODTrack *mu1 = aod->GetTrack(j);
       if(!mu1->IsMuonTrack()) continue;
-      charge_mu1 = mu1->Charge();
-      px_mu1 = mu1->Px();
-      py_mu1 = mu1->Py();
-      pz_mu1 = mu1->Pz();
-      e_mu1 = mu1->E();
-      p_mu1 = mu1->P();
-      pt_mu1 = mu1->Pt();
-      rapidity_mu1 = Rapidity(e_mu1,pz_mu1);
+      chargemu1 = mu1->Charge();
+      pxmu1 = mu1->Px();
+      pymu1 = mu1->Py();
+      pzmu1 = mu1->Pz();
+      emu1 = mu1->E();
+      pmu1 = mu1->P();
+      ptmu1 = mu1->Pt();
+      rapiditymu1 = Rapidity(emu1,pzmu1);
     }
-    ((TH1D*)(fOutput->FindObject("hP")))->Fill(p_mu1);    
-    ((TH1D*)(fOutput->FindObject("hPt")))->Fill(pt_mu1);    
-    ((TH1D*)(fOutput->FindObject("hRapidity")))->Fill(rapidity_mu1);     
+    ((TH1D*)(fOutput->FindObject("hP")))->Fill(pmu1);    
+    ((TH1D*)(fOutput->FindObject("hPt")))->Fill(ptmu1);    
+    ((TH1D*)(fOutput->FindObject("hRapidity")))->Fill(rapiditymu1);      
     nmuontracks++;
-    if(charge_mu1<0){
+    if(chargemu1<0){
       for (Int_t jj = 0; jj<ntracks; jj++) {
-        Float_t px_mu2=-999; Float_t py_mu2=-999; Float_t pz_mu2=-999;
-        Float_t e_mu2=-999;Float_t charge_mu2=-999; 
-        if(strcmp(fAnalysisType,"ESD")==0){ 
+        Float_t pxmu2=-999; Float_t pymu2=-999; Float_t pzmu2=-999;
+        Float_t emu2=-999;Float_t chargemu2=-999; 
+        if(strcmp(fkAnalysisType,"ESD")==0){ 
           AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(esd->GetMuonTrack(jj)));
           if (!mu2->ContainTrackerData()) continue;
-         charge_mu2 = mu2->Charge();
-          px_mu2 = mu2->Px();
-          py_mu2 = mu2->Py();
-          pz_mu2 = mu2->Pz();
-         e_mu2 = mu2->E();
-        } else if(strcmp(fAnalysisType,"AOD")==0){
+         chargemu2 = mu2->Charge();
+          pxmu2 = mu2->Px();
+          pymu2 = mu2->Py();
+          pzmu2 = mu2->Pz();
+         emu2 = mu2->E();
+        } else if(strcmp(fkAnalysisType,"AOD")==0){
           AliAODTrack *mu2 = aod->GetTrack(jj);
           if(!mu2->IsMuonTrack()) continue; 
-         charge_mu2 = mu2->Charge();
-          px_mu2 = mu2->Px();
-          py_mu2 = mu2->Py();
-          pz_mu2 = mu2->Pz();
-         e_mu2 = mu2->E();
+         chargemu2 = mu2->Charge();
+          pxmu2 = mu2->Px();
+          pymu2 = mu2->Py();
+          pzmu2 = mu2->Pz();
+         emu2 = mu2->E();
         }
-        if(charge_mu2>0){
-         Float_t pt_dimu = TMath::Sqrt((px_mu1+px_mu2)*(px_mu1+px_mu2)+(py_mu1+py_mu2)*(py_mu1+py_mu2));
-         Float_t mass_dimu = InvMass(e_mu1,px_mu1,py_mu1,pz_mu1,e_mu2,px_mu2,py_mu2,pz_mu2);
-         Float_t rapidity_dimu = Rapidity((e_mu1+e_mu2),(pz_mu1+pz_mu2));
-         Double_t costhetaCS_dimu = CostCS((Double_t) px_mu1,(Double_t) py_mu1,(Double_t)pz_mu1,(Double_t) e_mu1,(Double_t)charge_mu1,(Double_t) px_mu2,(Double_t) py_mu2,(Double_t)pz_mu2,(Double_t) e_mu2);
-         Double_t costhetaHE_dimu = CostHE((Double_t) px_mu1,(Double_t) py_mu1,(Double_t)pz_mu1,(Double_t) e_mu1,(Double_t)charge_mu1,(Double_t) px_mu2,(Double_t) py_mu2,(Double_t)pz_mu2,(Double_t) e_mu2);
-         ((TH1D*)(fOutput->FindObject("hMass_Dimu")))->Fill(mass_dimu);
-         ((TH1D*)(fOutput->FindObject("hPt_Dimu")))->Fill(pt_dimu);    
-         ((TH1D*)(fOutput->FindObject("hRapidity_Dimu")))->Fill(rapidity_dimu);        
-         ((TH1D*)(fOutput->FindObject("hCosThetaCS_Dimu")))->Fill(costhetaCS_dimu);    
-         ((TH1D*)(fOutput->FindObject("hCosThetaHE_Dimu")))->Fill(costhetaHE_dimu);    
+        if(chargemu2>0){
+         Float_t ptdimu = TMath::Sqrt((pxmu1+pxmu2)*(pxmu1+pxmu2)+(pymu1+pymu2)*(pymu1+pymu2));
+         Float_t massdimu = InvMass(emu1,pxmu1,pymu1,pzmu1,emu2,pxmu2,pymu2,pzmu2);
+         Float_t rapiditydimu = Rapidity((emu1+emu2),(pzmu1+pzmu2));
+         Double_t costhetaCSdimu = CostCS((Double_t) pxmu1,(Double_t) pymu1,(Double_t)pzmu1,(Double_t) emu1,(Double_t)chargemu1,(Double_t) pxmu2,(Double_t) pymu2,(Double_t)pzmu2,(Double_t) emu2);
+         Double_t costhetaHEdimu = CostHE((Double_t) pxmu1,(Double_t) pymu1,(Double_t)pzmu1,(Double_t) emu1,(Double_t)chargemu1,(Double_t) pxmu2,(Double_t) pymu2,(Double_t)pzmu2,(Double_t) emu2);
+         ((TH1D*)(fOutput->FindObject("hMassDimu")))->Fill(massdimu);
+         ((TH1D*)(fOutput->FindObject("hPtDimu")))->Fill(ptdimu);      
+         ((TH1D*)(fOutput->FindObject("hRapidityDimu")))->Fill(rapiditydimu);  
+         ((TH1D*)(fOutput->FindObject("hCosThetaCSDimu")))->Fill(costhetaCSdimu);      
+         ((TH1D*)(fOutput->FindObject("hCosThetaHEDimu")))->Fill(costhetaHEdimu);      
         }
         //delete mu2;
       }      // second mu Loop
@@ -270,6 +268,9 @@ void AliAnalysisTaskMuonDistributions::UserExec(Option_t *)
 //________________________________________________________________________
 void AliAnalysisTaskMuonDistributions::Terminate(Option_t *) 
 {
+//
+// Draw histos
+//
   gStyle->SetCanvasColor(10);
   gStyle->SetFrameFillColor(10);
   Int_t xmin=20; 
@@ -277,70 +278,74 @@ void AliAnalysisTaskMuonDistributions::Terminate(Option_t *)
   
   printf("Using beam Energy=%f \n",fBeamEnergy);
 
-  TH1D *h_NumberMuonTracks = dynamic_cast<TH1D*> (fOutput->FindObject("hNumberMuonTracks"));  
-  TH1D *h_Mass_Dimu = dynamic_cast<TH1D*> (fOutput->FindObject("hMass_Dimu"));  
-  TH1D *h_Pt_Dimu = dynamic_cast<TH1D*> (fOutput->FindObject("hPt_Dimu"));  
-  TH1D *h_Rapidity_Dimu = dynamic_cast<TH1D*> (fOutput->FindObject("hRapidity_Dimu"));  
-  TH1D *h_CostCS_Dimu = dynamic_cast<TH1D*> (fOutput->FindObject("hCosThetaCS_Dimu"));  
-  TH1D *h_CostHE_Dimu = dynamic_cast<TH1D*> (fOutput->FindObject("hCosThetaHE_Dimu"));  
-  TH1D *h_P = dynamic_cast<TH1D*> (fOutput->FindObject("hP"));  
-  TH1D *h_Pt = dynamic_cast<TH1D*> (fOutput->FindObject("hPt"));  
-  TH1D *h_Rapidity = dynamic_cast<TH1D*> (fOutput->FindObject("hRapidity"));  
+  TH1D *hNumberMuonTracks = dynamic_cast<TH1D*> (fOutput->FindObject("hNumberMuonTracks"));  
+  TH1D *hMassDimu = dynamic_cast<TH1D*> (fOutput->FindObject("hMassDimu"));  
+  TH1D *hPtDimu = dynamic_cast<TH1D*> (fOutput->FindObject("hPtDimu"));  
+  TH1D *hRapidityDimu = dynamic_cast<TH1D*> (fOutput->FindObject("hRapidityDimu"));  
+  TH1D *hCostCSDimu = dynamic_cast<TH1D*> (fOutput->FindObject("hCosThetaCSDimu"));  
+  TH1D *hCostHEDimu = dynamic_cast<TH1D*> (fOutput->FindObject("hCosThetaHEDimu"));  
+  TH1D *hP = dynamic_cast<TH1D*> (fOutput->FindObject("hP"));  
+  TH1D *hPt = dynamic_cast<TH1D*> (fOutput->FindObject("hPt"));  
+  TH1D *hRapidity = dynamic_cast<TH1D*> (fOutput->FindObject("hRapidity"));  
 
   TCanvas *c0 = new TCanvas("c0","Plots",xmin,ymin,600,600);
   c0->Divide(2,2);
   c0->cd(1);
-  h_NumberMuonTracks->Draw();
+  hNumberMuonTracks->Draw();
   
   xmin+=20; ymin+=20;
   TCanvas *c1 = new TCanvas("c1","Muon kinematic distributions Plots",xmin,ymin,600,600);
   c1->Divide(2,2);  
   c1->cd(1);
   gPad->SetLogy(1);
-  h_P->Draw();
+  hP->Draw();
   c1->cd(2);
   gPad->SetLogy(1);
-  h_Pt->Draw();
+  hPt->Draw();
   c1->cd(3);
-  h_Rapidity->Draw();
+  hRapidity->Draw();
   
   xmin+=20; ymin+=20;
   TCanvas *c2 = new TCanvas("c2","Dimuon kinematic distributions Plots",xmin,ymin,600,600);
   c2->Divide(2,2);  
   c2->cd(1);
   gPad->SetLogy(1);
-  h_Pt_Dimu->Draw();
+  hPtDimu->Draw();
   c2->cd(2);
-  h_Rapidity_Dimu->Draw();
+  hRapidityDimu->Draw();
   c2->cd(3);
-  h_CostCS_Dimu->Draw();
+  hCostCSDimu->Draw();
   c2->cd(4);
-  h_CostHE_Dimu->Draw();
+  hCostHEDimu->Draw();
   
   xmin+=20; ymin+=20;
   TCanvas *c3 = new TCanvas("c3","Invariant Mass Plots",xmin,ymin,600,600);
   gPad->SetLogy(1);
-  h_Mass_Dimu->Draw();  
-  if(fInvariantMassFit) FitInvMass(h_Mass_Dimu);
+  hMassDimu->Draw();  
+  if(fInvariantMassFit) FitInvMass(hMassDimu);
   c3->Update();
 
 }
 
 //________________________________________________________________________
 Float_t AliAnalysisTaskMuonDistributions::InvMass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
-                                  Float_t e2, Float_t px2, Float_t py2, Float_t pz2) 
+                                  Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const 
 {
+//
 // 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;
 }
 //________________________________________________________________________
-Float_t AliAnalysisTaskMuonDistributions::Rapidity(Float_t e, Float_t pz) 
+Float_t AliAnalysisTaskMuonDistributions::Rapidity(Float_t e, Float_t pz) const 
 {
+//
 // calculate rapidity
+//
     Float_t rap;
-    if(e!=pz){
+    if(TMath::Abs(e-pz)>1e-7){
       rap = 0.5*TMath::Log((e+pz)/(e-pz));
       return rap;
     } else {
@@ -352,46 +357,49 @@ Float_t AliAnalysisTaskMuonDistributions::Rapidity(Float_t e, Float_t pz)
 Double_t AliAnalysisTaskMuonDistributions::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)
 {
-  TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM. frame
-  TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame
+//TMath
+// costCS calculation
+//
+  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.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); 
-  PTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); 
+  pProjCM.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); 
+  pTargCM.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); 
   //
   // --- Get the muons parameters in the CM frame 
   //
-  PMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
-  PMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
+  pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
+  pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
   //
   // --- Obtain the dimuon parameters in the CM frame
   //
-  PDimuCM=PMu1CM+PMu2CM;
+  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);
+  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();
+  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());}
+  if(charge1>0) {cost = zaxisCS.Dot((pMu1Dimu.Vect()).Unit());}
+  else {cost = zaxisCS.Dot((pMu2Dimu.Vect()).Unit());}
   return cost;
 }
 
@@ -399,48 +407,53 @@ Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2)
 Double_t AliAnalysisTaskMuonDistributions::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
+//
+// costHE calculation
+//
+  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);
+  pMu1CM.SetPxPyPzE(px1,py1,pz1,e1);
+  pMu2CM.SetPxPyPzE(px2,py2,pz2,e2);
   //
   // --- Obtain the dimuon parameters in the CM frame
   //
-  PDimuCM=PMu1CM+PMu2CM;
+  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);
+  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();
+  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());} 
+  if(charge1>0) {cost = zaxis.Dot((pMu1Dimu.Vect()).Unit());} 
+  else {cost = zaxis.Dot((pMu2Dimu.Vect()).Unit());} 
   return cost;
 }  
 
 //________________________________________________________________________
 void AliAnalysisTaskMuonDistributions::FitInvMass(TH1D *histo)
 {
+//
+// Fit to the Invariant Mass Spectrum
+//
   TF1 *gau = new TF1("gau","gaus",fPsiFitLimitMin,fPsiFitLimitMax);
   TF1 *ex = new TF1("ex","expo",fBckFitLimitMin,fBckFitLimitMax);    
   TF1 *tot = new TF1("mtot","gaus(0)+expo(3)",fInvMassFitLimitMin,fInvMassFitLimitMax);
   Double_t par[5];
-  Double_t BinWidth= histo->GetBinWidth(1);
-  printf("BW=%f\n",BinWidth);
+  Double_t binWidth= histo->GetBinWidth(1);
   gau->SetLineColor(3);
   gau->SetLineWidth(2);
   histo->Fit(gau,"RlQ"); 
@@ -454,52 +467,52 @@ void AliAnalysisTaskMuonDistributions::FitInvMass(TH1D *histo)
   tot->SetLineWidth(2);
   histo->Fit(tot,"Rl+");
   histo->Draw("e");
-  Double_t Chi2 = tot->GetChisquare();
-  Double_t NDF = tot->GetNDF();
-  Float_t MeanPsi= tot->GetParameter(1);
-  Float_t SigPsi= tot->GetParameter(2)*1000.;
-  Double_t NPsiFit = TMath::Sqrt(2*3.1415)*tot->GetParameter(0)*tot->GetParameter(2)/BinWidth;
+  Double_t chi2 = tot->GetChisquare();
+  Double_t ndf = tot->GetNDF();
+  Float_t meanPsi= tot->GetParameter(1);
+  Float_t sigPsi= tot->GetParameter(2)*1000.;
+  Double_t nPsiFit = TMath::Sqrt(2*3.1415)*tot->GetParameter(0)*tot->GetParameter(2)/binWidth;
   TF1 *exfix = new TF1("exfix","expo",2.,5.);  
   exfix->SetParameter(0,tot->GetParameter(3));  
   exfix->SetParameter(1,tot->GetParameter(4));  
-  Double_t NBck = exfix->Integral(2.9,3.3)/BinWidth;
+  Double_t nBck = exfix->Integral(2.9,3.3)/binWidth;
   
   printf("\n\n****************************************************************************\n");
-  char psi_text[100];
-  sprintf(psi_text,"N. J/#psi (2.9-3.3)=%10.0f",NPsiFit);
-  printf("\nN. J/psi (2.9-3.3)=%10.0f\n",NPsiFit);
-  TLatex *psi_latex = new TLatex(4.5,0.85*histo->GetMaximum(),psi_text);
-  psi_latex->SetTextColor(2);
-  psi_latex->SetTextSize(0.03);
-  psi_latex->SetTextAlign(2);
-  psi_latex->Draw();
+  char psitext[100];
+  sprintf(psitext,"N. J/#psi (2.9-3.3)=%10.0f",nPsiFit);
+  printf("\nN. J/psi (2.9-3.3)=%10.0f\n",nPsiFit);
+  TLatex *psilatex = new TLatex(4.5,0.85*histo->GetMaximum(),psitext);
+  psilatex->SetTextColor(2);
+  psilatex->SetTextSize(0.03);
+  psilatex->SetTextAlign(2);
+  psilatex->Draw();
   
-  char psi2_text[100];
-  sprintf(psi2_text,"J/#psi m=%4.3f GeV #sigma=%4.2f MeV",MeanPsi,SigPsi);
-  printf("J/psi m=%4.3f GeV sigma=%4.2f MeV\n",MeanPsi,SigPsi);
-  TLatex *psi2_latex = new TLatex(4.5,0.425*histo->GetMaximum(),psi2_text);
-  psi2_latex->SetTextColor(2);
-  psi2_latex->SetTextSize(0.03);
-  psi2_latex->SetTextAlign(2);
-  psi2_latex->Draw();
+  char psi2text[100];
+  sprintf(psi2text,"J/#psi m=%4.3f GeV #sigma=%4.2f MeV",meanPsi,sigPsi);
+  printf("J/psi m=%4.3f GeV sigma=%4.2f MeV\n",meanPsi,sigPsi);
+  TLatex *psi2latex = new TLatex(4.5,0.425*histo->GetMaximum(),psi2text);
+  psi2latex->SetTextColor(2);
+  psi2latex->SetTextSize(0.03);
+  psi2latex->SetTextAlign(2);
+  psi2latex->Draw();
   
-  char sb_text[100];
-  sprintf(sb_text,"S/B (2.9-3.3)=%4.2f ",NPsiFit/NBck);
-  printf("S/B (2.9-3.3)=%4.2f\n",NPsiFit/NBck);
-  TLatex *sb_latex = new TLatex(4.5,0.212*histo->GetMaximum(),sb_text);
-  sb_latex->SetTextColor(2);
-  sb_latex->SetTextSize(0.03);
-  sb_latex->SetTextAlign(2);
-  sb_latex->Draw();
+  char sbtext[100];
+  sprintf(sbtext,"S/B (2.9-3.3)=%4.2f ",nPsiFit/nBck);
+  printf("S/B (2.9-3.3)=%4.2f\n",nPsiFit/nBck);
+  TLatex *sblatex = new TLatex(4.5,0.212*histo->GetMaximum(),sbtext);
+  sblatex->SetTextColor(2);
+  sblatex->SetTextSize(0.03);
+  sblatex->SetTextAlign(2);
+  sblatex->Draw();
   
-  char chi2_text[100];
-  sprintf(chi2_text,"#chi^2/ndf =%4.2f ",Chi2/NDF);
-  printf("chi^2/ndf =%4.2f\n",Chi2/NDF);
-  TLatex *chi2_latex = new TLatex(4.5,0.106*histo->GetMaximum(),chi2_text);
-  chi2_latex->SetTextColor(2);
-  chi2_latex->SetTextSize(0.03);
-  chi2_latex->SetTextAlign(2);
-  chi2_latex->Draw();
+  char chi2text[100];
+  sprintf(chi2text,"#chi^2/ndf =%4.2f ",chi2/ndf);
+  printf("chi^2/ndf =%4.2f\n",chi2/ndf);
+  TLatex *chi2latex = new TLatex(4.5,0.106*histo->GetMaximum(),chi2text);
+  chi2latex->SetTextColor(2);
+  chi2latex->SetTextSize(0.03);
+  chi2latex->SetTextAlign(2);
+  chi2latex->Draw();
   printf("\n****************************************************************************\n");
   
 }

diff --git a/PWG3/muon/AliAnalysisTaskMuonDistributions.h b/PWG3/muon/AliAnalysisTaskMuonDistributions.h
index c485e12b7546f235e4cd5dabbd6daba8138809b9..a6164903d9852af597a1f222799390ba5a474e14 100755 (executable)
--- a/PWG3/muon/AliAnalysisTaskMuonDistributions.h
+++ b/PWG3/muon/AliAnalysisTaskMuonDistributions.h
@@ -1,8 +1,12 @@
 #ifndef ALIANALYSISTASKMUONDISTRIBUTIONS_H
 #define ALIANALYSISTASKMUONDISTRIBUTIONS_H
 
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* Analysis Task for Muon/Dimuon distributions */
+
 #include "AliAnalysisTaskSE.h"
-#include "TMath.h"
 
 class TH1D;
 class TParticle ;
@@ -27,7 +31,7 @@ class AliAnalysisTaskMuonDistributions : public AliAnalysisTaskSE {
   void UserCreateOutputObjects();
   
   void SetBeamEnergy(Double_t en) {fBeamEnergy=en;}
-  void SetAnalysisType(const char* type) {fAnalysisType=type;}
+  void SetAnalysisType(const char* type) {fkAnalysisType=type;}
   void SetInvMassFitLimits(Double_t xmin, Double_t xmax) {fInvMassFitLimitMin=xmin; fInvMassFitLimitMax=xmax;}
   void SetPsiFitLimits(Double_t xmin, Double_t xmax) {fPsiFitLimitMin=xmin; fPsiFitLimitMax=xmax;}
   void SetBckFitLimits(Double_t xmin, Double_t xmax) {fBckFitLimitMin=xmin; fBckFitLimitMax=xmax;}
@@ -36,22 +40,22 @@ class AliAnalysisTaskMuonDistributions : public AliAnalysisTaskSE {
  protected:
   
   Double_t fBeamEnergy;   // Energy of the beam (required for the CS angle)
-  Double_t fInvMassFitLimitMin;  // invariant mass spectrum fit limits 
-  Double_t fInvMassFitLimitMax;
-  Double_t fPsiFitLimitMin;  // psi fit limits 
-  Double_t fPsiFitLimitMax;
-  Double_t fBckFitLimitMin;  // bck fit limits 
-  Double_t fBckFitLimitMax;
+  Double_t fInvMassFitLimitMin;  // invariant mass spectrum fit lower limit 
+  Double_t fInvMassFitLimitMax;  // invariant mass spectrum fit upper limit 
+  Double_t fPsiFitLimitMin;  // psi fit lower limits 
+  Double_t fPsiFitLimitMax;  // psi fit upper limits 
+  Double_t fBckFitLimitMin;  // bck fit lower limits 
+  Double_t fBckFitLimitMax;  // bck fit upper limits
   Bool_t fInvariantMassFit; // flag to perform or not inv. mass fit
     
-  const char* fAnalysisType; //ESD or AOD based analysis
-  TList *fOutput;
+  const char* fkAnalysisType; //ESD or AOD based analysis
+  TList *fOutput;  // output file
       
-  Float_t InvMass (Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, Float_t);
-  Float_t Rapidity (Float_t, Float_t);
-  Double_t CostCS (Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t);
-  Double_t CostHE (Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t);
-  void FitInvMass(TH1D* );
+  Float_t InvMass (Float_t e1, Float_t px1, Float_t py1, Float_t pz1, Float_t e2, Float_t px2, Float_t py2, Float_t pz2) const;
+  Float_t Rapidity (Float_t e, Float_t pz) const;
+  Double_t 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 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);
+  void FitInvMass(TH1D* histo);
   
   ClassDef(AliAnalysisTaskMuonDistributions,1);
 };