update on muon resonance task
authorrvernet <rvernet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Feb 2009 19:17:27 +0000 (19:17 +0000)
committerrvernet <rvernet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Feb 2009 19:17:27 +0000 (19:17 +0000)
CORRFW/test/muons/AliCFMuonResTask1.cxx

index cd5bbca..9a539eb 100644 (file)
-/**************************************************************************
- * 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 "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[9] ;
-////////
-//// 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);
-    Int_t pdg0 = part0->GetPdgCode();
-    Float_t e0 = part0->Energy();
-    Float_t pz0 = part0->Pz();
-    Float_t py0 = part0->Py();
-    Float_t px0 = part0->Px();
-    Float_t phi0 = part0->Phi(); // Warning in TParticle Phi = pi + ATan2(Py,Px) = [0,2pi] 
-    phi0 = Phideg(phi0);    
-    Float_t rapmc0=Rap(e0,pz0);
-    AliMCParticle *mcpart0 = new AliMCParticle(part0);
-
-    // selection of the rapidity for first muon
-    if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,mcpart0)) continue;                
-
-    Int_t p1 = part->GetDaughter(1);
-    part1 = stack->Particle(p1);
-    Int_t pdg1 = part1->GetPdgCode();
-    Float_t e1 = part1->Energy();
-    Float_t pz1 = part1->Pz();
-    Float_t py1 = part1->Py();
-    Float_t px1 = part1->Px();
-    Float_t phi1 = part1->Phi();
-    phi1 = Phideg(phi1);
-    Float_t rapmc1=Rap(e1,pz1);
-    AliMCParticle *mcpart1 = new AliMCParticle(part1);
-
-    // selection of the rapidity for second muon
-    if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,mcpart1)) continue;
-
-    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=Rap((e0+e1),(pz0+pz1));
-
-       containerInput[0] = fNevt ;   
-       containerInput[1] = rapmc0 ;   
-       containerInput[2] = phi0 ;   
-       containerInput[3] = rapmc1 ;   
-       containerInput[4] = phi1 ;   
-       containerInput[5] = imassmc ;   
-       containerInput[6] = rapmc ;   
-       containerInput[7] = ptmc;
-       containerInput[8] = pmc;        
-
-       // fill the container at the first step
-       fCFManager->GetParticleContainer()->Fill(containerInput,0);
-
-    } // one muon is positive
-  }    
-
-////////
-//// 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(); // Warning in AliESDMuonTrack Phi = ATan2(Py,Px) = [-pi,pi]
-           phir1 = phir1 * 57.296;
-           Float_t er1 = mu1->E();
-           Float_t rap1=Rap(er1,pzr1);
-           // selection of the rapidity for first muon
-           if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,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 = phir2 * 57.296;
-                   Float_t er2 = mu2->E();
-                   Float_t rap2=Rap(er2,pzr2);
-                   // selection of the rapidity for second muon
-                   if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,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);
-
-                   containerInput[0] = fNevt ;   
-                   containerInput[1] = rap1 ;   
-                   containerInput[2] = phir1 ;   
-                   containerInput[3] = rap2 ;
-                   containerInput[4] = phir2 ;   
-                   containerInput[5] = imassrec ;   
-                   containerInput[6] = raprec ;   
-                   containerInput[7] = ptrec;
-                   containerInput[8] = prec;
-                   
-                   // 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 from TParticle in range [-180,180] 
-    Float_t phideg;
-    
-       phideg = phi-TMath::Pi();
-       phideg = phideg*57.296;
-       return phideg;
-}
-//________________________________________________________________________
-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);
-
-    TCanvas *c1 = new TCanvas("AliCFMuonResTask1","JPSI MC & ESD",10,10,510,510);
-    c1->Divide(1,2);
-    c1->cd(1);
-    kmass->Draw("HIST");
-    c1->cd(2);
-    rmass->Draw("HIST");
-}
-//________________________________________________________________________
-
-#endif
+/**************************************************************************\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ *                                                                        *\r
+ * Author: The ALICE Off-line Project.                                    *\r
+ * Contributors are mentioned in the code where appropriate.              *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes is hereby granted   *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+\r
+//-----------------------------------------------------------------------\r
+// Example of task (running locally, on AliEn and CAF),\r
+// which provides standard way of calculating acceptance and efficiency\r
+// between different steps of the procedure.\r
+// The ouptut of the task is a AliCFContainer from which the efficiencies\r
+// can be calculated\r
+//-----------------------------------------------------------------------\r
+// Author : R. Vernet, Consorzio Cometa - Catania (it)\r
+//-----------------------------------------------------------------------\r
+// Modification done by X. Lopez - LPC Clermont (fr)\r
+//-----------------------------------------------------------------------\r
+\r
+\r
+#ifndef ALICFMUONRESTASK1_CXX\r
+#define ALICFMUONRESTASK1_CXX\r
+\r
+#include "AliCFMuonResTask1.h"\r
+#include "AliHeader.h"\r
+#include "AliESDHeader.h"\r
+#include "AliStack.h"\r
+#include "TParticle.h"\r
+#include "TH1I.h"\r
+#include "AliMCEvent.h"\r
+#include "AliAnalysisManager.h"\r
+#include "AliESDEvent.h"\r
+#include "AliAODEvent.h"\r
+#include "AliCFManager.h"\r
+#include "AliCFCutBase.h"\r
+#include "AliCFContainer.h"\r
+#include "TChain.h"\r
+#include "AliESDtrack.h"\r
+#include "AliLog.h"\r
+#include "AliESDMuonTrack.h"\r
+#include "AliESDtrack.h"\r
+#include "AliESDInputHandler.h"\r
+#include "TCanvas.h"\r
+\r
+ClassImp(AliCFMuonResTask1)\r
+\r
+//__________________________________________________________________________\r
+AliCFMuonResTask1::AliCFMuonResTask1() :\r
+  fReadAODData(0),\r
+  fCFManager(0x0),\r
+  fQAHistList(0x0),\r
+  fHistEventsProcessed(0x0),\r
+  fNevt(0)\r
+{\r
+  //\r
+  //Default ctor\r
+  //\r
+}\r
+//___________________________________________________________________________\r
+AliCFMuonResTask1::AliCFMuonResTask1(const Char_t* name) :\r
+  AliAnalysisTaskSE(name),\r
+  fReadAODData(0),\r
+  fCFManager(0x0),\r
+  fQAHistList(0x0),\r
+  fHistEventsProcessed(0x0),\r
+  fNevt(0)\r
+{\r
+  //\r
+  // Constructor. Initialization of Inputs and Outputs\r
+  //\r
+  Info("AliCFMuonResTask1","Calling Constructor");\r
+\r
+  fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;\r
+\r
+  DefineOutput(1,TH1I::Class());\r
+  DefineOutput(2,AliCFContainer::Class());\r
+\r
+}\r
+\r
+//___________________________________________________________________________\r
+AliCFMuonResTask1& AliCFMuonResTask1::operator=(const AliCFMuonResTask1& c) \r
+{\r
+  //\r
+  // Assignment operator\r
+  //\r
+  if (this!=&c) {\r
+    AliAnalysisTaskSE::operator=(c) ;\r
+    fReadAODData = c.fReadAODData ;\r
+    fCFManager  = c.fCFManager;\r
+    fQAHistList = c.fQAHistList ;\r
+    fHistEventsProcessed = c.fHistEventsProcessed;\r
+    fNevt = c.fNevt ;\r
+  }\r
+  return *this;\r
+}\r
+\r
+//___________________________________________________________________________\r
+AliCFMuonResTask1::AliCFMuonResTask1(const AliCFMuonResTask1& c) :\r
+  AliAnalysisTaskSE(c),\r
+  fReadAODData(c.fReadAODData),\r
+  fCFManager(c.fCFManager),\r
+  fQAHistList(c.fQAHistList),\r
+  fHistEventsProcessed(c.fHistEventsProcessed),\r
+  fNevt(c.fNevt) \r
+{\r
+  //\r
+  // Copy Constructor\r
+  //\r
+}\r
+\r
+//___________________________________________________________________________\r
+AliCFMuonResTask1::~AliCFMuonResTask1() {\r
+  //\r
+  //destructor\r
+  //\r
+  Info("~AliCFMuonResTask1","Calling Destructor");\r
+  if (fCFManager)           delete fCFManager ;\r
+  if (fHistEventsProcessed) delete fHistEventsProcessed ;\r
+  if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}\r
+}\r
+\r
+//_________________________________________________\r
+void AliCFMuonResTask1::UserExec(Option_t *)\r
+{\r
+  //\r
+  // Main loop function\r
+  // \r
+  \r
+  Info("UserExec","") ;\r
+  if (!fMCEvent) {\r
+    Error("UserExec","NO MC EVENT FOUND!");\r
+    return;\r
+  }\r
+\r
+  fNevt++; \r
+  Double_t containerInput[9] ;\r
\r
+////////\r
+//// MC\r
+//////// \r
+\r
+  fCFManager->SetEventInfo(fMCEvent);\r
+  AliStack* stack = fMCEvent->Stack();\r
+\r
+  // loop on the MC event\r
+  for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) { \r
+    AliMCParticle *mcPart  = fMCEvent->GetTrack(ipart);\r
\r
+    TParticle *part = mcPart->Particle(); \r
+    TParticle *part0 = mcPart->Particle();\r
+    TParticle *part1 = mcPart->Particle();\r
\r
+    // Selection of the resonance\r
+    if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue;\r
+\r
+    // Mother kinematics\r
+    Float_t e = part->Energy();\r
+    Float_t pz = part->Pz();           \r
+    Float_t py = part->Py();\r
+    Float_t px = part->Px();\r
+    Float_t rapmc = Rap(e,pz);\r
+    Float_t mass = part->GetCalcMass();\r
+\r
+    // Decays kinematics\r
+\r
+    Int_t p0 = part->GetDaughter(0);\r
+    part0 = stack->Particle(p0); \r
+   // selection of the rapidity for first muon\r
+    AliMCParticle *mcpart0 = new AliMCParticle(part0);\r
+    if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,mcpart0)) continue;\r
+    Int_t pdg0 = part0->GetPdgCode();\r
+\r
+    Int_t p1 = part->GetDaughter(1);\r
+    part1 = stack->Particle(p1);\r
+   // selection of the rapidity for second muon\r
+    AliMCParticle *mcpart1 = new AliMCParticle(part1);\r
+    if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,mcpart1)) continue;\r
+    Int_t pdg1 = part1->GetPdgCode();\r
\r
+   // 0 mu- 1 mu+\r
+    Float_t e0, pz0, py0, px0, phi0, rapmc0;\r
+    Float_t e1, pz1, py1, px1, phi1, rapmc1;\r
+\r
+    // ordering sign: first = mu-\r
+    if(pdg0==13){\r
+       e0 = part0->Energy();\r
+       pz0 = part0->Pz();\r
+       py0 = part0->Py();\r
+       px0 = part0->Px();\r
+       phi0 = part0->Phi(); // Warning in TParticle Phi = pi + ATan2(Py,Px) = [0,2pi] \r
+       phi0 = Phideg(phi0);    \r
+       rapmc0=Rap(e0,pz0);\r
+\r
+       e1 = part1->Energy();\r
+       pz1 = part1->Pz();\r
+       py1 = part1->Py();\r
+       px1 = part1->Px();\r
+       phi1 = part1->Phi();\r
+       phi1 = Phideg(phi1);\r
+       rapmc1=Rap(e1,pz1);\r
+    }\r
+    else{\r
+       e1 = part0->Energy();\r
+       pz1 = part0->Pz();\r
+       py1 = part0->Py();\r
+       px1 = part0->Px();\r
+       phi1 = part0->Phi();\r
+       phi1 = Phideg(phi1);    \r
+       rapmc1=Rap(e1,pz1);\r
+\r
+       e0 = part1->Energy();\r
+       pz0 = part1->Pz();\r
+       py0 = part1->Py();\r
+       px0 = part1->Px();\r
+       phi0 = part1->Phi();\r
+       phi0 = Phideg(phi0);\r
+       rapmc0=Rap(e0,pz0); \r
+    }\r
+\r
+       Float_t pmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1)+\r
+                                  (pz0+pz1)*(pz0+pz1));\r
+       Float_t ptmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1));\r
+       Float_t imassmc = Imass(e0,px0,py0,pz0,e1,px1,py1,pz1);\r
+       Float_t rapmc_check=Rap((e0+e1),(pz0+pz1));\r
+\r
+       containerInput[0] = fNevt ;   \r
+       containerInput[1] = rapmc0 ;   \r
+       containerInput[2] = phi0 ;   \r
+       containerInput[3] = rapmc1 ;   \r
+       containerInput[4] = phi1 ;   \r
+       containerInput[5] = imassmc ;   \r
+       containerInput[6] = rapmc ;   \r
+       containerInput[7] = ptmc;\r
+       containerInput[8] = pmc;        \r
+\r
+       // fill the container at the first step\r
+       fCFManager->GetParticleContainer()->Fill(containerInput,0);\r
\r
+  }    \r
+\r
+////////\r
+//// ESD\r
+////////\r
+  \r
+  AliESDEvent *fESD; \r
+  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>\r
+      (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
+  fESD = esdH->GetEvent();\r
+\r
+  Int_t mult1 = fESD->GetNumberOfMuonTracks() ;\r
+\r
+    for (Int_t j = 0; j<mult1; j++) { \r
+       AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));\r
+       Float_t zr1 = mu1->Charge();\r
+// Select mu-\r
+       if(zr1<0){\r
+           Float_t pxr1 = mu1->Px();\r
+           Float_t pyr1 = mu1->Py();\r
+           Float_t pzr1 = mu1->Pz();\r
+           Float_t phir1 = mu1->Phi(); // Warning in AliESDMuonTrack Phi = ATan2(Py,Px) = [-pi,pi]\r
+           phir1 = phir1 * 57.296;\r
+           Float_t er1 = mu1->E();\r
+           Float_t rap1=Rap(er1,pzr1);\r
+           // selection of the rapidity for first muon\r
+           if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,mu1)) continue;\r
+\r
+           for (Int_t jj = 0; jj<mult1; jj++) {\r
+               AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));\r
+               Float_t zr2 = mu2->Charge();\r
+// Select mu+\r
+               if(zr2>0 && jj !=j){\r
+                   Float_t pxr2 = mu2->Px();\r
+                   Float_t pyr2 = mu2->Py();\r
+                   Float_t pzr2 = mu2->Pz();\r
+                   Float_t phir2 = mu2->Phi();\r
+                   phir2 = phir2 * 57.296;\r
+                   Float_t er2 = mu2->E();\r
+                   Float_t rap2=Rap(er2,pzr2);\r
+                   // selection of the rapidity for second muon\r
+                   if (!fCFManager->CheckParticleCuts(AliCFManager::kPartRecCuts,mu2)) continue;\r
+\r
+                   Float_t prec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2)+\r
+                                               (pzr1+pzr2)*(pzr1+pzr2));\r
+                   Float_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));\r
+                   Float_t raprec=Rap((er1+er2),(pzr1+pzr2));\r
+                   Float_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);\r
+\r
+                   containerInput[0] = fNevt ;   \r
+                   containerInput[1] = rap1 ;   \r
+                   containerInput[2] = phir1 ;   \r
+                   containerInput[3] = rap2 ;\r
+                   containerInput[4] = phir2 ;   \r
+                   containerInput[5] = imassrec ;   \r
+                   containerInput[6] = raprec ;   \r
+                   containerInput[7] = ptrec;\r
+                   containerInput[8] = prec;\r
+                   \r
+                   // fill the container at the second step\r
+                   fCFManager->GetParticleContainer()->Fill(containerInput,1);\r
+\r
+               }  // mu+ Selection\r
+           }      // second mu Loop\r
+       }          // mu- Selection\r
+    }        \r
+\r
+//  ----------\r
+  fHistEventsProcessed->Fill(0);\r
+  PostData(1,fHistEventsProcessed) ;\r
+  PostData(2,fCFManager->GetParticleContainer()) ;\r
+}\r
+//________________________________________________________________________\r
+const Float_t AliCFMuonResTask1::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,\r
+                                  Float_t e2, Float_t px2, Float_t py2, Float_t pz2) \r
+{\r
+// invariant mass calculation\r
+    Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+\r
+                                    (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));\r
+    return imassrec;\r
+}\r
+//________________________________________________________________________\r
+const Float_t AliCFMuonResTask1::Rap(Float_t e, Float_t pz) \r
+{\r
+// calculate rapidity\r
+    Float_t rap;\r
+    if(e!=pz){\r
+       rap = 0.5*TMath::Log((e+pz)/(e-pz));\r
+       return rap;\r
+    }\r
+    else{\r
+       rap = -200;\r
+       return rap;\r
+    }\r
+}\r
+//________________________________________________________________________\r
+const Float_t AliCFMuonResTask1::Phideg(Float_t phi) \r
+{\r
+// calculate Phi from TParticle in range [-180,180] \r
+    Float_t phideg;\r
+    \r
+       phideg = phi-TMath::Pi();\r
+       phideg = phideg*57.296;\r
+       return phideg;\r
+}\r
+//________________________________________________________________________\r
+void AliCFMuonResTask1::Terminate(Option_t *) \r
+{\r
+  // draw result of the Invariant mass MC and ESD\r
+\r
+    AliCFContainer *cont = dynamic_cast<AliCFContainer*> (GetOutputData(2));   \r
+\r
+    TH1D *kmass = cont->ShowProjection(5,0);\r
+    TH1D *rmass = cont->ShowProjection(5,1);\r
+\r
+    TCanvas *c1 = new TCanvas("AliCFMuonResTask1","JPSI MC & ESD",10,10,510,510);\r
+    c1->Divide(1,2);\r
+    c1->cd(1);\r
+    kmass->Draw("HIST");\r
+    c1->cd(2);\r
+    rmass->Draw("HIST");\r
+}\r
+//________________________________________________________________________\r
+\r
+#endif\r