first template task for muons
authorrvernet <rvernet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 11 Feb 2009 10:21:35 +0000 (10:21 +0000)
committerrvernet <rvernet@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 11 Feb 2009 10:21:35 +0000 (10:21 +0000)
CORRFW/test/muons/AliCFMuonResTask1.C [new file with mode: 0644]
CORRFW/test/muons/AliCFMuonResTask1.cxx [new file with mode: 0644]
CORRFW/test/muons/AliCFMuonResTask1.h [new file with mode: 0644]

diff --git a/CORRFW/test/muons/AliCFMuonResTask1.C b/CORRFW/test/muons/AliCFMuonResTask1.C
new file mode 100644 (file)
index 0000000..8cd40dc
--- /dev/null
@@ -0,0 +1,270 @@
+// DEFINITION OF A FEW CONSTANTS
+
+const Double_t nevtmin= 1;
+const Double_t nevtmax = 15000;
+
+// Muons 
+const Double_t ymin  = -4.0 ;
+const Double_t ymax  =  -2.5 ;
+
+const Double_t phimin = -180;
+const Double_t phimax = 180;
+
+// Resonance
+const Int_t    PDG = 443;
+
+const Double_t ptmin =  0.0 ;
+const Double_t ptmax =  30 ;
+const Double_t pmin =  0.0 ;
+const Double_t pmax =  700 ;
+const Int_t    charge  = 0 ;
+const Double_t mmin =  0.1 ;
+const Double_t mmax =  6 ;
+const Double_t mymin =  -5 ;
+const Double_t mymax =  -1.5 ;
+
+//----------------------------------------------------
+
+Bool_t AliCFMuonResTask1(
+                           const Bool_t useGrid = 0,
+                           const Bool_t readAOD = 0,
+                           const char * kTagXMLFile="wn.xml" // XML file containing tags
+                           )
+{
+  
+  TBenchmark benchmark;
+  benchmark.Start("AliMuonResTask1");
+
+  AliLog::SetGlobalDebugLevel(0);
+
+  Load() ; // load the required libraries
+
+  TChain * analysisChain ;
+
+///// INPUT
+
+  if (useGrid) { // data located on AliEn
+    TGrid::Connect("alien://") ;    //  Create an AliRunTagCuts and an AliEventTagCuts Object 
+                                    //  and impose some selection criteria
+    AliRunTagCuts      *runCuts   = new AliRunTagCuts(); 
+    AliEventTagCuts    *eventCuts = new AliEventTagCuts(); 
+    AliLHCTagCuts      *lhcCuts   = new AliLHCTagCuts(); 
+    AliDetectorTagCuts *detCuts   = new AliDetectorTagCuts(); 
+    eventCuts->SetMultiplicityRange(0,2000);
+
+    //  Create an AliTagAnalysis Object and chain the tags
+    AliTagAnalysis   *tagAna = new AliTagAnalysis(); 
+    if (readAOD) tagAna->SetType("AOD");  // for aliroot > v4-05
+    else         tagAna->SetType("ESD");  // for aliroot > v4-05
+    TAlienCollection *coll   = TAlienCollection::Open(kTagXMLFile); 
+    TGridResult      *tagResult = coll->GetGridResult("",0,0);
+    tagResult->Print();
+    tagAna->ChainGridTags(tagResult);
+
+    // Create a new esd chain and assign the chain that is returned by querying the tags
+    analysisChain = tagAna->QueryTags(runCuts,lhcCuts,detCuts,eventCuts); 
+  }
+
+  else {// local data
+    // here put your input data path
+    printf("\n\nRunning on local file, please check the path\n\n");
+
+    if (readAOD) {
+      analysisChain = new TChain("aodTree");
+      analysisChain->Add("AliAOD.root");
+    }
+    else {
+      analysisChain = new TChain("esdTree");
+      analysisChain->Add("/scratch/lopez/PDC08jpsi/run2-300/AliESDs.root");
+   }
+  }
+  
+///// END INPUT
+
+
+  Info("AliCFMuonResTask1",Form("CHAIN HAS %d ENTRIES",(Int_t)analysisChain->GetEntries()));
+
+  // CONTAINER DEFINITION
+  Info("AliCFMuonResTask1","SETUP CONTAINER");
+  
+  // The sensitive variables (9 in this example), their indices
+  UInt_t nevt  = 0;
+  UInt_t y1  = 1;
+  UInt_t phi1  = 2;
+  UInt_t y2  = 3;
+  UInt_t phi2  = 4;
+  UInt_t imass  = 5;
+  UInt_t y  = 6;
+  UInt_t pt = 7;
+  UInt_t p = 8;
+
+
+  // Setting up the container grid
+  UInt_t nstep = 2 ; //number of selection steps : MC and ESD 
+  const Int_t nvar   = 9 ;     //number of variables on the grid
+  const Int_t nbin1  = nevtmax ;  
+  const Int_t nbin2  = 100 ;  
+  const Int_t nbin3  = 360 ;  
+  const Int_t nbin4  = 100 ;  
+  const Int_t nbin5  = 360 ;  
+  const Int_t nbin6  = 100 ;  
+  const Int_t nbin7  = 100 ;  
+  const Int_t nbin8  = 100 ;  
+  const Int_t nbin9  = 100 ;  
+
+  // arrays for the number of bins in each dimension
+  Int_t iBin[nvar];
+  iBin[0]=nbin1;
+  iBin[1]=nbin2;
+  iBin[2]=nbin3;
+  iBin[3]=nbin4;
+  iBin[4]=nbin5;
+  iBin[5]=nbin6;
+  iBin[6]=nbin7;
+  iBin[7]=nbin8;
+  iBin[8]=nbin9;
+
+  // arrays for lower bounds :
+  Double_t *binLim1=new Double_t[nbin1+1];
+  Double_t *binLim2=new Double_t[nbin2+1];
+  Double_t *binLim3=new Double_t[nbin3+1];
+  Double_t *binLim4=new Double_t[nbin4+1];
+  Double_t *binLim5=new Double_t[nbin5+1];
+  Double_t *binLim6=new Double_t[nbin6+1];
+  Double_t *binLim7=new Double_t[nbin7+1];
+  Double_t *binLim8=new Double_t[nbin8+1];
+  Double_t *binLim9=new Double_t[nbin9+1];
+
+  // values for bin lower bounds
+  for(Int_t i=0; i<=nbin1; i++) binLim1[i]=(Double_t)nevtmin  + (nevtmax-nevtmin)  /nbin1*(Double_t)i ;
+  for(Int_t i=0; i<=nbin2; i++) binLim2[i]=(Double_t)ymin  + (ymax-ymin)  /nbin2*(Double_t)i ;
+  for(Int_t i=0; i<=nbin3; i++) binLim3[i]=(Double_t)phimin  + (phimax-phimin)  /nbin3*(Double_t)i ;
+  for(Int_t i=0; i<=nbin4; i++) binLim4[i]=(Double_t)ymin  + (ymax-ymin)  /nbin4*(Double_t)i ;
+  for(Int_t i=0; i<=nbin5; i++) binLim5[i]=(Double_t)phimin  + (phimax-phimin)  /nbin5*(Double_t)i ;
+  for(Int_t i=0; i<=nbin6; i++) binLim6[i]=(Double_t)mmin  + (mmax-mmin)  /nbin6*(Double_t)i ;
+  for(Int_t i=0; i<=nbin7; i++) binLim7[i]=(Double_t)mymin  + (mymax-mymin)  /nbin7*(Double_t)i ;
+  for(Int_t i=0; i<=nbin8; i++) binLim8[i]=(Double_t)ptmin + (ptmax-ptmin)/nbin8*(Double_t)i ; 
+  for(Int_t i=0; i<=nbin9; i++) binLim9[i]=(Double_t)pmin + (pmax-pmin)/nbin9*(Double_t)i ; 
+
+  // one container  of 2 steps (MC and ESD) with 9 variables
+  AliCFContainer* container = new AliCFContainer("container","container for tracks",nstep,nvar,iBin);
+  // setting the bin limits
+  container -> SetBinLimits(nevt,binLim1);
+  container -> SetBinLimits(y1,binLim2);
+  container -> SetBinLimits(phi1,binLim3);
+  container -> SetBinLimits(y2,binLim4);
+  container -> SetBinLimits(phi2,binLim5);
+  container -> SetBinLimits(imass,binLim6);
+  container -> SetBinLimits(y,binLim7);
+  container -> SetBinLimits(pt,binLim8);
+  container -> SetBinLimits(p,binLim9);
+
+  // Set list
+  TList* qaList = new TList();
+
+  //CREATE THE CUTS
+  // Choice of the Resonance
+  AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts("mcGenCuts","MC particle generation cuts");
+  mcGenCuts->SetRequirePdgCode(PDG);
+  mcGenCuts->SetQAOn(qaList);
+
+  // Set a pt range of the resonance
+  AliCFTrackKineCuts *mcKineCuts = new AliCFTrackKineCuts("mcKineCuts","MC-level kinematic cuts");
+  mcKineCuts->SetChargeMC(charge);
+  mcKineCuts->SetPtRange(ptmin,ptmax);
+  mcKineCuts->SetQAOn(qaList);
+
+  // Create and fill the list associated 
+  TObjArray* mcList = new TObjArray(0) ;
+  mcList->AddLast(mcKineCuts);
+  mcList->AddLast(mcGenCuts);
+
+  // kinematic cuts on muons rapidity 
+  AliCFTrackKineCuts *recKineCuts = new AliCFTrackKineCuts("recKineCuts","rec-level kine cuts");
+  recKineCuts->SetRapidityRange(ymin,ymax);
+  recKineCuts->SetQAOn(qaList);
+  TObjArray* recList = new TObjArray(0) ;
+  recList->AddLast(recKineCuts);
+
+  // CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK
+  printf("CREATE INTERFACE AND CUTS\n");
+  AliCFManager* man = new AliCFManager() ;
+  man->SetParticleContainer     (container);
+  man->SetParticleCutsList(AliCFManager::kPartGenCuts,mcList);
+  man->SetParticleCutsList(AliCFManager::kPartRecCuts,recList);
+
+
+  //CREATE THE TASK
+  printf("CREATE TASK\n");
+  // create the task
+  AliCFMuonResTask1 *task = new AliCFMuonResTask1("AliMuonResTask1");
+  task->SetCFManager(man); //here is set the CF manager
+  task->SetQAList(qaList);
+  if (readAOD)       task->SetReadAODData() ;
+
+  //SETUP THE ANALYSIS MANAGER TO READ INPUT CHAIN AND WRITE DESIRED OUTPUTS
+  printf("CREATE ANALYSIS MANAGER\n");
+  // Make the analysis manager
+  AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
+
+  if (useGrid) mgr->SetAnalysisType(AliAnalysisManager::kGridAnalysis);
+  else mgr->SetAnalysisType(AliAnalysisManager::kLocalAnalysis);
+
+
+  AliMCEventHandler*  mcHandler = new AliMCEventHandler();
+  mgr->SetMCtruthEventHandler(mcHandler);
+  AliInputEventHandler* dataHandler ;
+  
+  if   (readAOD) dataHandler = new AliAODInputHandler();
+  else           dataHandler = new AliESDInputHandler();
+  mgr->SetInputEventHandler(dataHandler);
+
+  // Create and connect containers for input/output
+
+  // input data 
+  AliAnalysisDataContainer *cinput0  = mgr->CreateContainer("cchain0",TChain::Class(),AliAnalysisManager::kInputContainer);
+
+  // output data
+  Char_t file[256];
+  sprintf(file,"CFMuonResTask1.root");
+  printf("Analysis output in %s \n",file);
+
+  // output TH1I for event counting
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("chist0", TH1I::Class(),AliAnalysisManager::kOutputContainer,file);
+  // output Correction Framework Container (for acceptance & efficiency calculations)
+  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("ccontainer0", AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,file);
+
+  cinput0->SetData(analysisChain);
+  mgr->AddTask(task);
+  mgr->ConnectInput(task,0,cinput0);
+  mgr->ConnectOutput(task,1,coutput1);
+  mgr->ConnectOutput(task,2,coutput2);
+
+  printf("READY TO RUN\n");
+  //RUN !!!
+  if (mgr->InitAnalysis()) {
+    mgr->PrintStatus();
+    mgr->StartAnalysis("local",analysisChain);
+  }
+
+  benchmark.Stop("AliMuonResTask1");
+  benchmark.Show("AliMuonResTask1");
+
+  return kTRUE ;
+}
+
+void Load() {
+
+  //load the required aliroot libraries
+  gSystem->Load("libANALYSIS") ;
+
+  gSystem->Load("libANALYSISalice") ;
+
+  //gSystem->Load("libCORRFW.so") ;
+  gSystem->Load("$ALICE_ROOT/lib/tgt_linux/libCORRFW.so") ;
+
+  //compile online the task class
+  gSystem->SetIncludePath("-I. -I$ALICE_ROOT/include -I$ROOTSYS/include");
+  gROOT->LoadMacro("./AliCFMuonResTask1.cxx+");
+}
diff --git a/CORRFW/test/muons/AliCFMuonResTask1.cxx b/CORRFW/test/muons/AliCFMuonResTask1.cxx
new file mode 100644 (file)
index 0000000..cd5bbca
--- /dev/null
@@ -0,0 +1,349 @@
+/**************************************************************************
+ * 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
diff --git a/CORRFW/test/muons/AliCFMuonResTask1.h b/CORRFW/test/muons/AliCFMuonResTask1.h
new file mode 100644 (file)
index 0000000..357ca63
--- /dev/null
@@ -0,0 +1,73 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------------------------
+// Author : R. Vernet, Consorzio Cometa - Catania (it)
+//-----------------------------------------------------------------------
+// Modification done by X. Lopez - LPC Clermont (fr)
+//-----------------------------------------------------------------------
+
+#ifndef ALICFMUONRESTASK1_H
+#define ALICFMUONRESTASK1_H
+
+#include "AliAnalysisTaskSE.h"
+
+class TH1I;
+class TParticle ;
+class TFile ;
+class AliStack ;
+class AliCFManager;
+class AliESDtrack;
+class AliVParticle;
+
+class AliCFMuonResTask1 : public AliAnalysisTaskSE {
+  public:
+
+  AliCFMuonResTask1();
+  AliCFMuonResTask1(const Char_t* name);
+  AliCFMuonResTask1& operator= (const AliCFMuonResTask1& c);
+  AliCFMuonResTask1(const AliCFMuonResTask1& c);
+  virtual ~AliCFMuonResTask1();
+
+  // ANALYSIS FRAMEWORK STUFF to loop on data and fill output objects
+  void     UserExec(Option_t *option);
+  void     Terminate(Option_t *);
+   
+  // CORRECTION FRAMEWORK RELATED FUNCTIONS
+  void           SetCFManager(AliCFManager* const io) {fCFManager = io;}   // global correction manager
+  AliCFManager * GetCFManager() const {return fCFManager;}           // get corr manager
+  void           SetQAList(TList* const list) {fQAHistList = list;}
+
+  // Data types
+  Bool_t IsReadAODData()   const {return fReadAODData;}
+  void   SetReadAODData   (Bool_t flag=kTRUE) {fReadAODData=flag;}
+
+ protected:
+  
+  Bool_t          fReadAODData ;   // flag for AOD/ESD input files
+  AliCFManager   *fCFManager   ;   // pointer to the CF manager
+  TList          *fQAHistList  ;   // list of QA histograms
+  TH1I           *fHistEventsProcessed; // simple histo for monitoring the number of events processed
+  Int_t           fNevt        ;   // event countor
+  
+  const Float_t Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1,
+               Float_t e2, Float_t px2, Float_t py2, Float_t p2);
+  const  Float_t Rap(Float_t e, Float_t pz);
+  const Float_t Phideg(Float_t phi);
+  
+  ClassDef(AliCFMuonResTask1,1);
+};
+
+#endif