New task and container for dimuon pairs (Livio, Enrico, Roberta)
authormartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 28 Feb 2010 18:49:59 +0000 (18:49 +0000)
committermartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 28 Feb 2010 18:49:59 +0000 (18:49 +0000)
PWG3/PWG3muonLinkDef.h
PWG3/libPWG3muon.pkg
PWG3/muon/AddPWG3MuonTrain.C
PWG3/muon/AddTaskDimuonCFContainerBuilder.C [new file with mode: 0644]
PWG3/muon/AliAnalysisTaskDimuonCFContainerBuilder.cxx [new file with mode: 0644]
PWG3/muon/AliAnalysisTaskDimuonCFContainerBuilder.h [new file with mode: 0644]

index 793deed..2dbac69 100644 (file)
@@ -26,5 +26,6 @@
 #pragma link C++ class AliAODMuonTrack+;
 #pragma link C++ class AliAODMuonPair+;
 #pragma link C++ class AliAnalysisTaskSEMuonsHF+;
+#pragma link C++ class AliAnalysisTaskDimuonCFContainerBuilder+;
 #endif
 
index 6754e4b..639a22e 100644 (file)
@@ -21,7 +21,8 @@ SRCS:=         muon/AliAnalysisTaskESDMuonFilter.cxx \
                 muon/AliMCMuonPair.cxx \
                 muon/AliAODMuonTrack.cxx \
                 muon/AliAODMuonPair.cxx \
-                muon/AliAnalysisTaskSEMuonsHF.cxx
+                muon/AliAnalysisTaskSEMuonsHF.cxx \
+                muon/AliAnalysisTaskDimuonCFContainerBuilder.cxx
 
 
 HDRS:= $(SRCS:.cxx=.h)
index 531cab8..b3a2baf 100644 (file)
@@ -1,8 +1,9 @@
-Int_t AddPWG3MuonTrain(Int_t iESDAnalysis,
-                       Int_t iAODAnalyis,
-                       Int_t addMuonDistributions,
-                      Int_t addSingleMuonAnalysis,
-                      Int_t addMuonHFAnalysis) {
+Int_t AddPWG3MuonTrain(Int_t iESDAnalysis=1,
+                       Int_t iAODAnalyis=0,
+                       Int_t addMuonDistributions=0,
+                      Int_t addSingleMuonAnalysis=0,
+                      Int_t addMuonHFAnalysis=0,
+                      Int_t addDimuonCFContainer=0) {
 
   // Analysis wagons for PWG3Muon (Roberta)
 
@@ -41,5 +42,20 @@ Int_t AddPWG3MuonTrain(Int_t iESDAnalysis,
     ntasks++;
   }
 
+  if(addDimuonCFContainer) {
+    taskName="AddTaskDimuonCFContainerBuilder.C"; taskName.Prepend(loadMacroPath.Data());
+    gROOT->LoadMacro(taskName.Data());
+    Int_t runMode = 0;
+    Bool_t isAOD;
+    if(iAODAnalysis==1) isAOD=kTRUE;
+    else isAOD=kFALSE;
+    Bool_t isAcceptance = kFALSE;
+    Bool_t readMC   = kFALSE;
+    Double_t ebeam = 3500.;
+    AliAnalysisTaskDimuonCFContainerBuilder *dimuonCFtask =
+    AddTaskDimuonCFContainerBuilder(isAOD,readMC,isAcceptance,ebeam);
+    ntasks++;
+  }
+
   return ntasks;
 }
diff --git a/PWG3/muon/AddTaskDimuonCFContainerBuilder.C b/PWG3/muon/AddTaskDimuonCFContainerBuilder.C
new file mode 100644 (file)
index 0000000..f6e5096
--- /dev/null
@@ -0,0 +1,214 @@
+// VARIABLES RANGES
+const Double_t ymin          =  -4.0       ;           
+const Double_t ymax          =  -2.5       ;
+const Double_t ptmin         =   0.0       ;
+const Double_t ptmax         =  20.        ;
+const Double_t cCSmin        =  -1.        ;
+const Double_t cCSmax        =   1.        ;
+const Double_t cHEmin        =  -1.         ;
+const Double_t cHEmax        =   1.         ;
+const Double_t pCSmin        =   0.         ;           
+const Double_t pCSmax        =   TMath::Pi();  
+const Double_t pHEmin        =   0.         ;           
+const Double_t pHEmax        =   TMath::Pi(); 
+const Double_t massmin       =   0.         ; 
+const Double_t massmax       =  12.         ;          
+const Double_t trigmin       =   0.        ;
+const Double_t trigmax       =   4.         ;
+const Double_t ptmuminMIN    =   0.        ;
+const Double_t ptmuminMAX    =   100.      ;
+const Double_t ptmumaxMIN    =   0.         ; 
+const Double_t ptmumaxMAX    =   100.       ; 
+const Double_t thetamuminMIN =   0.        ;
+const Double_t thetamuminMAX =   180.      ;
+const Double_t thetamumaxMIN =   0.         ;
+const Double_t thetamumaxMAX =   180.       ;
+const Double_t pmuminMIN     =   0.        ;
+const Double_t pmuminMAX     =   100.      ;
+const Double_t pmumaxMIN     =   0.         ;
+const Double_t pmumaxMAX     =   100.       ;
+const Double_t trigsideMIN   =   0         ;
+const Double_t trigsideMAX   =   4         ;
+
+
+AliAnalysisTaskDimuonCFContainerBuilder *AddTaskDimuonCFContainerBuilder(Bool_t readAOD=kTRUE, Bool_t readMC=kTRUE, 
+                                               Bool_t isaccept = kTRUE, Double_t beamEn=3500)
+{
+
+   // Check and Info printings
+   //==============================================================================
+    if(!readMC && isaccept) {
+       printf("ERROR: incompatible choice readMC-isaccept. If isaccept you must readMC!\n"); 
+       return NULL;
+    } else if (readMC && isaccept) printf("Creating task for filling a CFcontainer with acceptance data.\n");
+    else if (readMC && !isaccept) printf("Creating task for filling a CFcontainer with simulated data.\n");
+    else printf("Creating task for filling a CFcontainer with real data.\n");
+       
+
+   // Get the pointer to the existing analysis manager via the static access method.
+   //==============================================================================
+    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+    if (!mgr) {
+      ::Error("AddDataTaskMuonPolarCF", "No analysis manager to connect to.");
+      return NULL;
+    }   
+   
+
+   // MC handler if needed
+   //==============================================================================
+    if(!readAOD && readMC){
+     AliMCEventHandler *mcH = (AliMCEventHandler*)mgr->GetMCtruthEventHandler();
+     if (!mcH) {
+       ::Error("AddDataTaskMuonPolarCF", "No MC handler connected");
+       return NULL;
+     } 
+    }
+
+
+   // DEFINING CONTAINER
+   //==============================================================================
+
+     UInt_t y                = 0;              // Association of variables with int numbers
+     UInt_t pt        = 1;
+     UInt_t costHE    = 2;
+     UInt_t phiHE     = 3;
+     UInt_t costCS    = 4;
+     UInt_t phiCS     = 5;
+     UInt_t mass      = 6;
+     UInt_t trig      = 7;
+     UInt_t ptmumin   = 8;
+     UInt_t ptmumax   = 9;
+     UInt_t thetamumin=10;
+     UInt_t thetamumax=11;
+     UInt_t pmumin    =12;
+     UInt_t pmumax    =13;
+     UInt_t trigside  =14;
+
+     UInt_t nstep = 8;                 // Number of layers (always 8 - not always filled). 4 with CINT and 4 with CMU
+
+     const Int_t nvar   = 15 ;         // Number of variables of the grid
+
+     const Int_t nbin1  = 10 ;         // Number of bins for each variable
+     const Int_t nbin2  = 10 ;  
+     const Int_t nbin3  = 20 ;  
+     const Int_t nbin4  = 20 ;  
+     const Int_t nbin5  = 20 ;  
+     const Int_t nbin6  = 20 ;  
+     const Int_t nbin7  = 240;  
+     const Int_t nbin8  = 40 ;  
+     const Int_t nbin9  = 100;  
+     const Int_t nbin10 = 100;  
+     const Int_t nbin11 = 180;  
+     const Int_t nbin12 = 180;  
+     const Int_t nbin13 = 100;  
+     const Int_t nbin14 = 100;  
+     const Int_t nbin15 = 4  ;  
+
+     Int_t iBin[nvar];                 // Array containing the number of bins for each variable
+     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;
+     iBin[9] =nbin10;
+     iBin[10]=nbin11;
+     iBin[11]=nbin12;
+     iBin[12]=nbin13;
+     iBin[13]=nbin14;
+     iBin[14]=nbin15;
+  
+     Double_t *binLim1 = new Double_t[nbin1+1];                // Arrays for lower bounds 
+     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];
+     Double_t *binLim10= new Double_t[nbin10+1];
+     Double_t *binLim11= new Double_t[nbin11+1];
+     Double_t *binLim12= new Double_t[nbin12+1];
+     Double_t *binLim13= new Double_t[nbin13+1];
+     Double_t *binLim14= new Double_t[nbin14+1];
+     Double_t *binLim15= new Double_t[nbin15+1];
+     for(Int_t i=0; i<=nbin1; i++) binLim1[i] =(Double_t)ymin+(ymax-ymin)/nbin1*(Double_t)i;
+     for(Int_t i=0; i<=nbin2; i++) binLim2[i] =(Double_t)ptmin+(ptmax-ptmin)/nbin2*(Double_t)i; 
+     for(Int_t i=0; i<=nbin3; i++) binLim3[i] =(Double_t)cHEmin+(cHEmax-cHEmin)/nbin3*(Double_t)i ; 
+     for(Int_t i=0; i<=nbin4; i++) binLim4[i] =(Double_t)pHEmin+(pHEmax-pHEmin)/nbin4*(Double_t)i ; 
+     for(Int_t i=0; i<=nbin5; i++) binLim5[i] =(Double_t)cCSmin+(cCSmax-cCSmin)/nbin5*(Double_t)i ; 
+     for(Int_t i=0; i<=nbin6; i++) binLim6[i] =(Double_t)pCSmin+(pCSmax-pCSmin)/nbin6*(Double_t)i ; 
+     for(Int_t i=0; i<=nbin7; i++) binLim7[i] =(Double_t)massmin+(massmax-massmin)/nbin7*(Double_t)i ; 
+     for(Int_t i=0; i<=nbin8; i++) binLim8[i] =(Double_t)trigmin+(trigmax-trigmin)/nbin8*(Double_t)i ; 
+     for(Int_t i=0; i<=nbin9; i++) binLim9[i] =(Double_t)ptmuminMIN+(ptmuminMAX-ptmuminMIN)/nbin9*(Double_t)i ; 
+     for(Int_t i=0; i<=nbin10;i++) binLim10[i]=(Double_t)ptmumaxMIN+(ptmumaxMAX-ptmumaxMIN)/nbin10*(Double_t)i ; 
+     for(Int_t i=0; i<=nbin11;i++) binLim11[i]=(Double_t)thetamuminMIN+(thetamuminMAX-thetamuminMIN)/nbin11*(Double_t)i ; 
+     for(Int_t i=0; i<=nbin12;i++) binLim12[i]=(Double_t)thetamumaxMIN+(thetamumaxMAX-thetamumaxMIN)/nbin12*(Double_t)i ; 
+     for(Int_t i=0; i<=nbin13;i++) binLim13[i]=(Double_t)pmuminMIN+(pmuminMAX-pmuminMIN)/nbin13*(Double_t)i ; 
+     for(Int_t i=0; i<=nbin14;i++) binLim14[i]=(Double_t)pmumaxMIN+(pmumaxMAX-pmumaxMIN)/nbin14*(Double_t)i ; 
+     for(Int_t i=0; i<=nbin15;i++) binLim15[i]=(Double_t)trigsideMIN+(trigsideMAX-trigsideMIN)/nbin15*(Double_t)i ; 
+
+     AliCFContainer* container = new AliCFContainer("container","Container for Dimuons",nstep,nvar,iBin);
+  
+     container -> SetBinLimits(y,binLim1);             // setting the bin limits
+     container -> SetBinLimits(pt,binLim2);
+     container -> SetBinLimits(costHE,binLim3);
+     container -> SetBinLimits(phiHE,binLim4);
+     container -> SetBinLimits(costCS,binLim5);
+     container -> SetBinLimits(phiCS,binLim6);
+     container -> SetBinLimits(mass,binLim7);
+     container -> SetBinLimits(trig,binLim8);
+     container -> SetBinLimits(ptmumin,binLim9);
+     container -> SetBinLimits(ptmumax,binLim10);
+     container -> SetBinLimits(thetamumin,binLim11);
+     container -> SetBinLimits(thetamumax,binLim12);
+     container -> SetBinLimits(pmumin,binLim13);
+     container -> SetBinLimits(pmumax,binLim14);
+     container -> SetBinLimits(trigside,binLim15);
+
+
+   // CF Manager
+   //==============================================================================
+     AliCFManager* man = new AliCFManager() ;
+     man->SetParticleContainer(container);
+
+  
+   // Outputs: list of histograms + CFContainer
+   //==============================================================================
+     TString outputfile = AliAnalysisManager::GetCommonFileName();
+     outputfile += ":PWG3Muon_DimuonCFContainer";
+
+     AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("Histos",TList::Class(),AliAnalysisManager::kOutputContainer,outputfile);
+     AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("DimuonCFContainer",AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,outputfile);
+
+  // The task with the associtated CF manager
+   //==============================================================================
+     AliAnalysisTaskDimuonCFContainerBuilder *task = new AliAnalysisTaskDimuonCFContainerBuilder("AliAnalysisTaskDimuonCFContainerBuilder",readAOD,readMC,isaccept,beamEn);
+     task->SetCFManager(man);
+
+
+   // Additional settings for the task (including some cuts)
+   //==============================================================================
+    //task->SetDistinguishTrigClass(kTRUE);
+    //task->SetReadMCinfo(kTRUE);
+    //Double_t ptlimits[2]={1.,1000.};
+    //task->SetPtSingMuLimits(ptlimits);
+    //task->SetCutonZvtxSPD(kTRUE);
+    //Double_t vtxlimits[2]={-10.,10.};
+    //task->SetZprimVertLimits(vtxlimits);
+
+
+   // Adding the task to the analysis manager and connecting inputs and outputs
+   //==============================================================================
+    mgr->AddTask(task);
+    mgr->ConnectInput(task,  0, mgr->GetCommonInputContainer());
+    mgr->ConnectOutput(task,1,coutput1);
+    mgr->ConnectOutput(task,2,coutput2);
+
+   return task;
+}
diff --git a/PWG3/muon/AliAnalysisTaskDimuonCFContainerBuilder.cxx b/PWG3/muon/AliAnalysisTaskDimuonCFContainerBuilder.cxx
new file mode 100644 (file)
index 0000000..e0b51bc
--- /dev/null
@@ -0,0 +1,1027 @@
+#ifndef ALIANALYSISTASKDIMUONCFCONTAINERBUILDER_CXX
+#define ALIANALYSISTASKDIMUONCFCONTAINERBUILDER_CXX
+
+#include "AliAnalysisTaskDimuonCFContainerBuilder.h"
+#include "AliHeader.h"
+#include "AliESDHeader.h"
+#include "AliStack.h"
+#include "TParticle.h"
+#include "TString.h"
+#include "TLorentzVector.h"
+#include "TFile.h"
+#include "TH1I.h"
+#include "TH2D.h"
+#include "AliMCEvent.h"
+#include "AliAnalysisManager.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliCFManager.h"
+#include "AliCFCutBase.h"
+#include "AliCFContainer.h"
+#include "AliCFEffGrid.h"
+#include "AliCFDataGrid.h"
+#include "TChain.h"
+#include "AliESDtrack.h"
+#include "AliLog.h"
+#include "AliESDMuonTrack.h"
+#include "AliESDtrack.h"
+#include "AliESDInputHandler.h"
+#include "TCanvas.h"
+#include "AliAODInputHandler.h"
+#include "AliAODMCParticle.h"
+
+//     Analysis task for the building of a dimuon CF container
+//     Also some single-muon variables are stored
+//     L. Bianchi - Universita' & INFN Torino
+
+
+ClassImp(AliAnalysisTaskDimuonCFContainerBuilder)
+
+//__________________________________________________________________________
+AliAnalysisTaskDimuonCFContainerBuilder::AliAnalysisTaskDimuonCFContainerBuilder() :
+  fReadAODData(0),
+  fReadMCInfo(kTRUE),
+  fIsAccProduction(kTRUE),
+  fCFManager(0x0),
+  fQAHistList(0x0),
+  fNevt(0),
+  fBeamEnergy(3500.),
+  fOutput(0x0),
+  fCutOnzVtxSPD(kFALSE),
+  fCutOnNContributors(kFALSE),
+  fTrigClassMuon(""),
+  fTrigClassInteraction(""),
+  fDistinguishTrigClass(kFALSE)
+{
+  
+  //Default constructor
+  
+  Double_t chilims[2]={0.,10000.};
+  Double_t ptlims[2]={0.,100.};
+  Double_t thetalims[2]={0.,180.};
+  Double_t vtxlims[2]={-1000.,1000.};
+  SetChi2Limits(chilims);
+  SetChi2MatchLimits(chilims);
+  SetPtSingMuLimits(ptlims);
+  SetThetaSingMuLimits(thetalims);
+  SetZprimVertLimits(vtxlims);
+  SetTrigClassMuonName();
+  SetTrigClassInteracName();
+  TString namemuonside[4]={"CMUS1A-","CMUS1B-","CMUS1C-","CMUS1-E-"};
+  TString nameinteractionside[4]={"CINT1A-","CINT1B-","CINT1C-","CINT1-E"};
+  SetTrigClassMuonSideName(namemuonside);
+  SetTrigClassInteracSideName(nameinteractionside);
+}
+//___________________________________________________________________________
+AliAnalysisTaskDimuonCFContainerBuilder::AliAnalysisTaskDimuonCFContainerBuilder(const Char_t* name, Bool_t readaod, Bool_t readMC, Bool_t isaccept, Double_t beamEn) :
+  AliAnalysisTaskSE(name),
+  fReadAODData(0),
+  fReadMCInfo(kTRUE),
+  fIsAccProduction(kTRUE),
+  fCFManager(0x0),
+  fQAHistList(0x0),
+  fNevt(0),
+  fBeamEnergy(3500.),
+  fOutput(0x0),
+  fCutOnzVtxSPD(kFALSE),
+  fCutOnNContributors(kFALSE),
+  fTrigClassMuon(""),
+  fTrigClassInteraction(""),
+  fDistinguishTrigClass(kFALSE)
+{
+  //
+  // Constructor. Initialization of Inputs and Outputs
+  //
+  Info("AliAnalysisTaskDimuonCFContainerBuilder","Calling Constructor");
+
+  SetReadAODData(readaod);
+  SetReadMCinfo(readMC);
+  SetIsAccProd(isaccept);
+  SetBeamEnergy(beamEn);
+
+  Double_t chilims[2]={0.,10000.};
+  Double_t ptlims[2]={0.,100.};
+  Double_t thetalims[2]={0.,180.};
+  Double_t vtxlims[2]={-1000.,1000.};
+  SetChi2Limits(chilims);
+  SetChi2MatchLimits(chilims);
+  SetPtSingMuLimits(ptlims);
+  SetThetaSingMuLimits(thetalims);
+  SetZprimVertLimits(vtxlims);
+  SetTrigClassMuonName();
+  SetTrigClassInteracName();
+  TString namemuonside[4]={"CMUS1A-","CMUS1B-","CMUS1C-","CMUS1-E-"};
+  TString nameinteractionside[4]={"CINT1A-","CINT1B-","CINT1C-","CINT1-E"};
+  SetTrigClassMuonSideName(namemuonside);
+  SetTrigClassInteracSideName(nameinteractionside);
+  
+  DefineOutput(1,TList::Class());
+  DefineOutput(2,AliCFContainer::Class());
+
+}
+
+//___________________________________________________________________________
+AliAnalysisTaskDimuonCFContainerBuilder& AliAnalysisTaskDimuonCFContainerBuilder::operator=(const AliAnalysisTaskDimuonCFContainerBuilder& c) 
+{
+  //
+  // Assignment operator
+  //
+  if (this!=&c) {
+    AliAnalysisTaskSE::operator=(c) ;
+    fReadAODData = c.fReadAODData ;
+    fCFManager  = c.fCFManager;
+    fQAHistList = c.fQAHistList ;
+    fNevt = c.fNevt ;
+  }
+  return *this;
+}
+
+//___________________________________________________________________________
+AliAnalysisTaskDimuonCFContainerBuilder::AliAnalysisTaskDimuonCFContainerBuilder(const AliAnalysisTaskDimuonCFContainerBuilder& c) :
+  AliAnalysisTaskSE(c),
+  fReadAODData(c.fReadAODData),
+  fReadMCInfo(c.fReadMCInfo),
+  fIsAccProduction(c.fIsAccProduction),
+  fCFManager(c.fCFManager),
+  fQAHistList(c.fQAHistList),
+  fNevt(c.fNevt),
+  fBeamEnergy(c.fBeamEnergy),
+  fOutput(c.fOutput),
+  fCutOnzVtxSPD(c.fCutOnzVtxSPD),
+  fCutOnNContributors(c.fCutOnNContributors),
+  fTrigClassMuon(c.fTrigClassMuon),
+  fTrigClassInteraction(c.fTrigClassInteraction),
+  fDistinguishTrigClass(c.fDistinguishTrigClass)
+{
+  
+  // Copy Constructor
+  
+}
+
+//___________________________________________________________________________
+AliAnalysisTaskDimuonCFContainerBuilder::~AliAnalysisTaskDimuonCFContainerBuilder() {
+  //
+  //destructor
+  //
+  Info("~AliAnalysisTaskDimuonCFContainerBuilder","Calling Destructor");
+  if (fCFManager)           delete fCFManager ;
+  if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;}
+}
+
+//___________________________________________________________________________
+void AliAnalysisTaskDimuonCFContainerBuilder::UserCreateOutputObjects(){
+       
+ fOutput = new TList();
+ fOutput->SetOwner(); 
+ TH1D *hnevts  = new TH1D("hnevts","hnevts",1,0,1);                                    // Stat check histos
+ TH1D *JpsiMult = new TH1D("JpsiMult","JpsiMult",20,0,20);     
+
+ TH1D *ptdimuREC    = new TH1D("ptdimuREC","ptdimuREC",200,0,20);                      // Dimu check histos
+ TH1D *ydimuREC     = new TH1D("ydimuREC","ydimuREC",200,-10.,10.);    
+ TH1D *ydimuFail     = new TH1D("ydimuFail","ydimuFail",10,660,670);
+ TH1D *costHEdimuREC= new TH1D("costHEdimuREC","costHEdimuREC",200,-1.,1.);    
+ TH1D *costCSdimuREC= new TH1D("costCSdimuREC","costCSdimuREC",200,-1.,1.);    
+ TH1D *costHEdimuFail= new TH1D("costHEdimuFail","costHEdimuFail",10,660.,670.);       
+ TH1D *costCSdimuFail= new TH1D("costCSdimuFail","costCSdimuFail",10,660.,670.);       
+ TH1D *phiHEdimuREC = new TH1D("phiHEdimuREC","phiHEdimuREC",100,0.,TMath::Pi());      
+ TH1D *phiCSdimuREC = new TH1D("phiCSdimuREC","phiCSdimuREC",100,0.,TMath::Pi());      
+ TH1D *phiHEdimuFail = new TH1D("phiHEdimuFail","phiHEdimuFail",10,660.,670.); 
+ TH1D *phiCSdimuFail = new TH1D("phiCSdimuFail","phiCSdimuFail",10,660.,670.); 
+ TH1D *imassTot     = new TH1D("imassTot","imassTot",100,0,8); 
+ TH1D *trigCond     = new TH1D("trigCond","trigCond",40,0,4);  
+
+ TH1D *EmuonREC        = new TH1D("EmuonREC","EmuonREC",200,0.,20.);                           // Mu check histos
+ TH1D *ptmuonREC= new TH1D("ptmuonREC","ptmuonREC",200,0.,20.);
+ TH1D *ymuonREC        = new TH1D("ymuonREC","ymuonREC",200,-10.,10.); 
+ TH1D *hdca    = new TH1D("hdca","hdca",20,0,200);
+ TH2D *hdcay   = new TH2D("hdcay","hdcay",200,0,200,20,-5,0);
+
+ TH1D *zvSPD   = new TH1D("zvSPD","zvSPD",100,-50,50.);                                // Event check histos
+ TH1D *zvSPDcut        = new TH1D("zvSPDcut","zvSPDcut",100,-50,50.);
+ TH1D *NContrib        = new TH1D("NContrib","NContrib",100,0,100);
+
+
+ fOutput->Add(hnevts);
+ fOutput->Add(JpsiMult); 
+
+ fOutput->Add(ptdimuREC); 
+ fOutput->Add(ydimuREC); 
+ fOutput->Add(ydimuFail); 
+ fOutput->Add(costHEdimuREC);
+ fOutput->Add(costCSdimuREC);
+ fOutput->Add(costHEdimuFail);
+ fOutput->Add(costCSdimuFail);
+ fOutput->Add(phiHEdimuREC);
+ fOutput->Add(phiCSdimuREC);
+ fOutput->Add(phiHEdimuFail);
+ fOutput->Add(phiCSdimuFail);
+ fOutput->Add(imassTot); 
+ fOutput->Add(trigCond); 
+
+ fOutput->Add(EmuonREC); 
+ fOutput->Add(ptmuonREC);
+ fOutput->Add(ymuonREC);
+ fOutput->Add(hdca);
+ fOutput->Add(hdcay);
+ fOutput->Add(zvSPD); 
+ fOutput->Add(zvSPDcut); 
+ fOutput->Add(NContrib); 
+
+       
+ fOutput->ls();
+
+} 
+
+
+//_________________________________________________
+void AliAnalysisTaskDimuonCFContainerBuilder::UserExec(Option_t *)
+{
+  
+  
+  //Info("UserExec","");
+  fNevt++;
+  ((TH1D*)(fOutput->FindObject("hnevts")))->Fill(0.5);
+
+  Double_t containerInput[15]={666,666,666,666,666,666,666,666,666,666,666,666,666,666,666};   //y, pT, costHE, phiHE, costCS, phiCS, mass, TrigCond
+  Int_t cutAccept =1;
+  Int_t NumbJpsis =0;
+   
+  if (!fReadAODData){     // ESD-based ANALYSIS
+
+    if (fReadMCInfo){
+
+      if (!fMCEvent) {
+       Error("UserExec","NO MC EVENT FOUND!");
+       //return;
+      }
+      // MC part  ---------------------------------------------------------------------------------------------------
+
+      //fCFManager->SetEventInfo(fMCEvent);                                    CHANGES IN NEW VERSIONS - MANUAL CUT ON PDG
+      AliStack* stack = fMCEvent->Stack();
+      // loop on the MC event
+      for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) { 
+
+       AliMCParticle *mcPart  = (AliMCParticle*) fMCEvent->GetTrack(ipart);
+  
+       TParticle *part = mcPart->Particle(); 
+       TParticle *part0 = mcPart->Particle();
+       TParticle *part1 = mcPart->Particle();
+
+       // Mother kinematics
+       Double_t e = part->Energy();
+       Double_t pz = part->Pz();                    
+       Double_t rapmc = Rap(e,pz);
+
+       // Selection of the resonance
+       //if (!fCFManager->CheckParticleCuts(0,mcPart)) continue;
+       if (part->GetPdgCode()!=443) continue;                                  // MANUAL CUT ON PDG
+       NumbJpsis++;
+       
+       // Decays kinematics
+       Int_t p0 = part->GetDaughter(0);
+       part0 = stack->Particle(p0); 
+       Int_t pdg0 = part0->GetPdgCode();
+       Int_t p1 = part->GetDaughter(1);
+       part1 = stack->Particle(p1);
+       Int_t pdg1 = part1->GetPdgCode();
+  
+       Double_t e0 = part0->Energy();
+       Double_t pz0 = part0->Pz();
+       Double_t py0 = part0->Py();
+       Double_t px0 = part0->Px();
+       Double_t charge0 = part0->GetPDG()->Charge()/3;
+       Double_t theta0 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px0*px0+py0*py0),pz0);
+       Double_t pt0 = TMath::Sqrt(px0*px0+py0*py0);
+       Double_t Mom0 = part0->P();
+
+       Double_t e1 = part1->Energy();
+       Double_t pz1 = part1->Pz();
+       Double_t py1 = part1->Py();
+       Double_t px1 = part1->Px();
+       //Double_t charge1 = part1->GetPDG()->Charge()/3;
+       Double_t theta1 = (180./TMath::Pi())*TMath::ATan2(TMath::Sqrt(px1*px1+py1*py1),pz1);
+       Double_t pt1 = TMath::Sqrt(px1*px1+py1*py1);
+       Double_t Mom1 = part1->P();
+
+     
+       if(pdg0==13 || pdg1==13) { 
+         Double_t ptmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1));
+         Double_t imassmc = Imass(e0,px0,py0,pz0,e1,px1,py1,pz1);
+      
+         Double_t costCS = CostCS(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
+         Double_t costHE = CostHE(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
+         Double_t phiCS  = PhiCS(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
+         Double_t phiHE  = PhiHE(px0,py0,pz0,e0,charge0,px1,py1,pz1,e1);
+         containerInput[0]  = rapmc ;   
+         containerInput[1]  = ptmc;
+         containerInput[2]  = costHE;    
+         containerInput[3]  = TMath::Abs(phiHE);     
+         containerInput[4]  = costCS;  
+         containerInput[5]  = TMath::Abs(phiCS);
+         containerInput[6]  = imassmc;      
+         containerInput[7]  = 1.;              //for generated no trigger condition
+         if (pt0<pt1) {
+           containerInput[8]=pt0; 
+           containerInput[9]=pt1;
+         } else {
+           containerInput[8]=pt1; 
+           containerInput[9]=pt0;
+         }
+         if (theta0<theta1) {
+           containerInput[10]=theta0; 
+           containerInput[11]=theta1;
+         } else {
+           containerInput[10]=theta1; 
+           containerInput[11]=theta0;
+         }
+         if (Mom0<Mom1) {
+           containerInput[12]=Mom0; 
+           containerInput[13]=Mom1;
+         } else {
+           containerInput[12]=Mom1; 
+           containerInput[13]=Mom0;
+         }
+         containerInput[14]=1.;
+         // fill the container at the first step
+         fCFManager->GetParticleContainer()->Fill(containerInput,0);
+               
+         if(fIsAccProduction){
+           // acceptance cuts on single mu
+           if(theta0<fThetaSingMuCut[0]||theta0>fThetaSingMuCut[1]||theta1<fThetaSingMuCut[0]||theta1>fThetaSingMuCut[1]) cutAccept=0;
+           if(pt0<fPtSingMuCut[0] || pt0>fPtSingMuCut[1] || pt1<fPtSingMuCut[0] || pt1>fPtSingMuCut[1]) cutAccept=0;
+           if (cutAccept ==1) fCFManager->GetParticleContainer()->Fill(containerInput,2);
+         }
+       }
+      } 
+    } //fReadMCInfo
+
+
+      // ESD part  ---------------------------------------------------------------------------------------------------
+  
+      AliESDEvent *fESD; 
+      AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+      fESD = esdH->GetEvent();
+      Int_t mult1 = fESD->GetNumberOfMuonTracks() ;
+
+      Int_t trigfired=-1;
+      Int_t trigside=-1;
+      if(fDistinguishTrigClass){
+       //aod->GetHeader()->SetFiredTriggerClasses("CMU");
+       TString trigclass = fESD->GetFiredTriggerClasses();
+       printf("TrigClass: %s\n",trigclass.Data());
+        if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
+        else if (trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
+       printf("Muon container: %d\n",trigfired);
+       for(Int_t i=0;i<4;i++){
+         if(trigfired==1 && trigclass.Contains(fTrigClassMuonSide[i])) {trigside=i;printf("Muon Side: %d\n",trigside);}
+         if(trigfired==0 && trigclass.Contains(fTrigClassInteractionSide[i])) {trigside=i;printf("Interaction Side: %d\n",trigside);}
+       }
+      }
+
+      ((TH1D*)(fOutput->FindObject("zvSPD")))->Fill(fESD->GetPrimaryVertexSPD()->GetZ());
+      ((TH1D*)(fOutput->FindObject("NContrib")))->Fill(fESD->GetPrimaryVertexSPD()->GetNContributors());
+      if (!fCutOnzVtxSPD || (fCutOnzVtxSPD && (fESD->GetPrimaryVertexSPD()->GetZ()>fzPrimVertexSPD[0] && fESD->GetPrimaryVertexSPD()->GetZ()<fzPrimVertexSPD[1]))){
+      ((TH1D*)(fOutput->FindObject("zvSPDcut")))->Fill(fESD->GetPrimaryVertexSPD()->GetZ());
+      if (!fCutOnNContributors || (fCutOnNContributors && (fESD->GetPrimaryVertexSPD()->GetNContributors()>fNContributors[0] && fESD->GetPrimaryVertexSPD()->GetNContributors()<fNContributors[1]))){
+      for (Int_t j = 0; j<mult1; j++) { 
+
+       AliESDMuonTrack* mu1 = new AliESDMuonTrack(*(fESD->GetMuonTrack(j)));
+       if (!mu1->ContainTrackerData()) continue;
+       if (mu1->GetChi2()<fChi2Track[0] || mu1->GetChi2()>fChi2Track[1]) continue;
+       if (mu1->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu1->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
+       Double_t zr1 = mu1->Charge();
+       Double_t pxr1  = mu1->Px();
+       Double_t pyr1  = mu1->Py();
+       Double_t pzr1  = mu1->Pz();
+        Double_t ptmu1 = TMath::Sqrt(pxr1*pxr1+pyr1*pyr1);
+       Double_t er1 = mu1->E();
+       Double_t Mom1 = mu1->P();
+       Double_t theta1 = (180./TMath::Pi())*mu1->Theta();
+       Double_t rapiditymu1 = Rap(er1,pzr1); 
+       ((TH1D*)(fOutput->FindObject("EmuonREC")))->Fill(er1);
+       ((TH1D*)(fOutput->FindObject("ptmuonREC")))->Fill(ptmu1);
+       ((TH1D*)(fOutput->FindObject("ymuonREC")))->Fill(rapiditymu1);
+
+       if(zr1<0){
+
+           for (Int_t jj = 0; jj<mult1; jj++) {
+
+             AliESDMuonTrack* mu2 = new AliESDMuonTrack(*(fESD->GetMuonTrack(jj)));
+             if (!mu2->ContainTrackerData()) continue;
+             if (mu2->GetChi2()<fChi2Track[0] || mu2->GetChi2()>fChi2Track[1]) continue;
+             if (mu2->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu2->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
+             Double_t zr2 = mu2->Charge();
+
+             if(zr2>0){
+               Double_t TrigCondition=0;
+               if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) TrigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
+               else TrigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();
+               ((TH1D*)(fOutput->FindObject("trigCond")))->Fill(TrigCondition);
+               Double_t pxr2 = mu2->Px();
+               Double_t pyr2 = mu2->Py();
+               Double_t pzr2 = mu2->Pz();
+                Double_t ptmu2 = TMath::Sqrt(pxr2*pxr2+pyr2*pyr2);
+               Double_t er2 = mu2->E();
+               Double_t Mom2 = mu2->P();
+               Double_t theta2 = (180./TMath::Pi())*mu2->Theta();
+
+               Double_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2));
+               ((TH1D*)(fOutput->FindObject("ptdimuREC")))->Fill(ptrec);
+               Double_t raprec= Rap((er1+er2),(pzr1+pzr2));
+               ((TH1D*)(fOutput->FindObject("ydimuREC")))->Fill(raprec);
+               Double_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2);
+               ((TH1D*)(fOutput->FindObject("imassTot")))->Fill(imassrec);
+
+               if (imassrec>12.) continue;
+                   
+               Double_t costCSrec = CostCS(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
+               ((TH1D*)(fOutput->FindObject("costCSdimuREC")))->Fill(costCSrec);
+               if(costCSrec==666) ((TH1D*)(fOutput->FindObject("costCSdimuFail")))->Fill(costCSrec);
+               Double_t costHErec = CostHE(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
+               ((TH1D*)(fOutput->FindObject("costHEdimuREC")))->Fill(costHErec);
+               if(costHErec==666) ((TH1D*)(fOutput->FindObject("costHEdimuFail")))->Fill(costHErec);
+               Double_t phiCSrec  = PhiCS(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
+               ((TH1D*)(fOutput->FindObject("phiCSdimuREC")))->Fill(phiCSrec);
+               if(phiCSrec==666) ((TH1D*)(fOutput->FindObject("phiCSdimuFail")))->Fill(phiCSrec);
+               Double_t phiHErec  = PhiHE(pxr1,pyr1,pzr1,er1,zr1,pxr2,pyr2,pzr2,er2);
+               ((TH1D*)(fOutput->FindObject("phiHEdimuREC")))->Fill(phiHErec);
+               if(phiHErec==666) ((TH1D*)(fOutput->FindObject("phiHEdimuFail")))->Fill(phiHErec);
+
+               containerInput[0] = raprec ;   
+               containerInput[1] = ptrec;
+               containerInput[2] = costHErec;  
+               containerInput[3] = TMath::Abs(phiHErec);       
+               containerInput[4] = costCSrec;  
+               containerInput[5] = TMath::Abs(phiCSrec);
+               containerInput[6] = imassrec;  
+               containerInput[7] = TrigCondition+0.05;
+               if (ptmu1<ptmu2) {
+                 containerInput[8]=ptmu1; 
+                 containerInput[9]=ptmu2;
+               } else {
+                 containerInput[8]=ptmu2; 
+                 containerInput[9]=ptmu1;
+               }
+               if (theta1<theta2) {
+                 containerInput[10]=theta1; 
+                 containerInput[11]=theta2;
+               } else {
+                 containerInput[10]=theta2; 
+                 containerInput[11]=theta1;
+               }
+               if (Mom1<Mom2) {
+                 containerInput[12]=Mom1; 
+                 containerInput[13]=Mom2;
+               } else {
+                 containerInput[12]=Mom2; 
+                 containerInput[13]=Mom1;
+               }
+               if (fDistinguishTrigClass && trigside==0) containerInput[14]=0.;
+               else if (fDistinguishTrigClass && trigside==1) containerInput[14]=1.;
+               else if (fDistinguishTrigClass && trigside==2) containerInput[14]=2.;
+               else if (fDistinguishTrigClass && trigside==3) containerInput[14]=3.;
+               else containerInput[14]=0.;
+                   
+               if(fDistinguishTrigClass){    
+                 if (trigfired==1) fCFManager->GetParticleContainer()->Fill(containerInput,5);
+                 else if (trigfired==0) fCFManager->GetParticleContainer()->Fill(containerInput,1);
+               } else {
+                 fCFManager->GetParticleContainer()->Fill(containerInput,1);
+                 if(fIsAccProduction){
+                   if(cutAccept==1) fCFManager->GetParticleContainer()->Fill(containerInput,3);
+                 }
+               }
+             }  // mu+ Selection
+
+           }      // second mu Loop
+        }          // mu- Selection
+     }  
+       }
+      }
+  } else {     // AOD-based ANALYSIS
+
+    AliAODEvent *aod;
+    AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+    aod = aodH->GetEvent();
+    Int_t ntracks=aod->GetNumberOfTracks(); 
+
+    // MC part  ---------------------------------------------------------------------------------------------------
+
+    if (fReadMCInfo){
+      TClonesArray *mcarray = dynamic_cast<TClonesArray*> (aod->FindListObject(AliAODMCParticle::StdBranchName()));  //array of MC particles in this event
+      for(int ii=0;ii<mcarray->GetEntries();ii++){
+       AliAODMCParticle *mctrack = (AliAODMCParticle*) mcarray->At(ii);
+       if(mctrack->GetPdgCode()!=13) continue;
+       Int_t NumbMother = mctrack->GetMother();
+       if (NumbMother==-1) continue;
+       AliAODMCParticle *mother = (AliAODMCParticle*) mcarray->At(NumbMother);
+       if (mother->GetPdgCode()!=443) continue;
+       NumbJpsis++;
+       Int_t daught0 = mother->GetDaughter(0);
+       Int_t daught1 = mother->GetDaughter(1);
+       AliAODMCParticle *mcDaughter0 = (AliAODMCParticle*) mcarray->At(daught0);
+       Double_t pxmc0 = mcDaughter0->Px();
+       Double_t pymc0 = mcDaughter0->Py();
+       Double_t pzmc0 = mcDaughter0->Pz();
+       Double_t ptmc0 = mcDaughter0->Pt();
+       Double_t Mommc0 = mcDaughter0->P();
+       Double_t Emc0  = mcDaughter0->E();
+       Double_t thetamc0 = (180./TMath::Pi())*mcDaughter0->Theta();
+       Int_t charge0 = (Int_t) mcDaughter0->Charge()/3;
+       AliAODMCParticle *mcDaughter1 = (AliAODMCParticle*) mcarray->At(daught1);
+       Double_t pxmc1 = mcDaughter1->Px();
+       Double_t pymc1 = mcDaughter1->Py();
+       Double_t pzmc1 = mcDaughter1->Pz();
+       Double_t ptmc1 = mcDaughter1->Pt();
+       Double_t Mommc1 = mcDaughter1->P();
+       Double_t Emc1  = mcDaughter1->E();
+       Double_t thetamc1 = (180./TMath::Pi())*mcDaughter1->Theta();
+       Int_t charge1 = (Int_t) mcDaughter1->Charge()/3;
+      
+       if (charge0==charge1 || TMath::Abs(mcDaughter0->GetPdgCode())!=13 || TMath::Abs(mcDaughter1->GetPdgCode())!=13) continue;
+       //Double_t rapJpsi = mother->Y();               //PROBLEM!!!
+       Double_t rapJpsi = Rap(mother->E(),mother->Pz());
+       Double_t ptJpsi = mother->Pt();
+       Double_t costHE = CostHE(pxmc0,pymc0,pzmc0,Emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,Emc1);
+       Double_t phiHE  = PhiHE(pxmc0,pymc0,pzmc0,Emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,Emc1);
+       Double_t costCS = CostCS(pxmc0,pymc0,pzmc0,Emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,Emc1);
+       Double_t phiCS  = PhiCS(pxmc0,pymc0,pzmc0,Emc0,(Double_t)charge0,pxmc1,pymc1,pzmc1,Emc1);
+       Double_t massJpsi = mother->M();
+       containerInput[0] = rapJpsi;
+       containerInput[1] = ptJpsi;
+       containerInput[2] = costHE;
+       containerInput[3] = TMath::Abs(phiHE);
+       containerInput[4] = costCS;
+       containerInput[5] = TMath::Abs(phiCS);
+       containerInput[6] = massJpsi;
+       containerInput[7] = 1.;         //for generated no trigger condition
+       if (ptmc0<ptmc1) {
+         containerInput[8]=ptmc0; 
+         containerInput[9]=ptmc1;
+       } else {
+         containerInput[8]=ptmc1; 
+         containerInput[9]=ptmc0;
+       }
+       if (thetamc0<thetamc1) {
+         containerInput[10]=thetamc0; 
+         containerInput[11]=thetamc1;
+       } else {
+         containerInput[10]=thetamc1; 
+         containerInput[11]=thetamc0;
+       }
+       if (Mommc0<Mommc1) {
+         containerInput[12]=Mommc0; 
+         containerInput[13]=Mommc1;
+       } else {
+         containerInput[12]=Mommc1; 
+         containerInput[13]=Mommc0;
+       }
+       containerInput[14]=1.;
+      
+       if(rapJpsi!=666 && costHE!=666 && phiHE!=666 && costCS!=666 && phiCS!=666){
+         fCFManager->GetParticleContainer()->Fill(containerInput,0);
+         if(fIsAccProduction){
+           // acceptance cuts on single mu
+           if(thetamc0<fThetaSingMuCut[0]||thetamc0>fThetaSingMuCut[1]||thetamc1<fThetaSingMuCut[0]||thetamc1>fThetaSingMuCut[1]) cutAccept=0;
+           if(ptmc0<fPtSingMuCut[0] || ptmc0>fPtSingMuCut[1] || ptmc1<fPtSingMuCut[0] || ptmc1>fPtSingMuCut[1]) cutAccept=0;
+           if (cutAccept ==1) fCFManager->GetParticleContainer()->Fill(containerInput,2);
+         }
+         
+       }
+
+       }
+      }
+      
+      
+      // AOD part  ---------------------------------------------------------------------------------------------------
+      
+      Int_t trigfired=-1;
+      Int_t trigside=-1;
+      if(fDistinguishTrigClass){
+       TString trigclass = aod->GetFiredTriggerClasses();
+       printf("TrigClass: %s\n",trigclass.Data());
+        if(trigclass.Contains(fTrigClassMuon)) trigfired = 1;
+        else if (trigclass.Contains(fTrigClassInteraction)) trigfired = 0;
+       printf("Muon container: %d\n",trigfired);
+       for(Int_t i=0;i<4;i++){
+         if(trigfired==1 && trigclass.Contains(fTrigClassMuonSide[i])) {trigside=i;printf("Muon Side: %d\n",trigside);}
+         if(trigfired==0 && trigclass.Contains(fTrigClassInteractionSide[i])) {trigside=i;printf("Interaction Side: %d\n",trigside);}
+       }
+      }
+      
+      ((TH1D*)(fOutput->FindObject("zvSPD")))->Fill(aod->GetPrimaryVertex()->GetZ());
+      Int_t ncontr = aod->GetPrimaryVertex()->GetNContributors();
+      ((TH1D*)(fOutput->FindObject("NContrib")))->Fill(ncontr);
+      if (!fCutOnzVtxSPD || (fCutOnzVtxSPD && (aod->GetPrimaryVertex()->GetZ()>fzPrimVertexSPD[0] && aod->GetPrimaryVertex()->GetZ()<fzPrimVertexSPD[1]))){    //NOT REALLY SPD VERTEX
+      ((TH1D*)(fOutput->FindObject("zvSPDcut")))->Fill(aod->GetPrimaryVertex()->GetZ());
+      if (!fCutOnNContributors || (fCutOnNContributors && (aod->GetPrimaryVertex()->GetNContributors()>fNContributors[0] && aod->GetPrimaryVertex()->GetNContributors()<fNContributors[1]))){
+      for (Int_t j = 0; j<ntracks; j++) {
+       AliAODTrack *mu1 = aod->GetTrack(j);
+       if(!mu1->IsMuonTrack()) continue;
+       if (mu1->Chi2perNDF()<fChi2Track[0] || mu1->Chi2perNDF()>fChi2Track[1]) continue;
+       if (mu1->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu1->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
+       Double_t chargemu1 = mu1->Charge();
+       Double_t pxmu1 = mu1->Px();
+       Double_t pymu1 = mu1->Py();
+       Double_t pzmu1 = mu1->Pz();
+       Double_t emu1 = mu1->E();
+       Double_t pmu1 = mu1->P();
+       Double_t ptmu1 = mu1->Pt();
+       Double_t rapiditymu1 = Rap(emu1,pzmu1); 
+       Double_t thetamu1 = (180./TMath::Pi())*mu1->Theta();
+       Double_t Rdcamu1 = mu1->DCA();
+       ((TH1D*)(fOutput->FindObject("EmuonREC")))->Fill(emu1);
+       ((TH1D*)(fOutput->FindObject("ptmuonREC")))->Fill(ptmu1);
+       ((TH1D*)(fOutput->FindObject("ymuonREC")))->Fill(rapiditymu1);
+       if(chargemu1<0){
+         for (Int_t jj = 0; jj<ntracks; jj++) { 
+           AliAODTrack *mu2 = aod->GetTrack(jj);
+           if(!mu2->IsMuonTrack()) continue;
+           if (mu2->Chi2perNDF()<fChi2Track[0] || mu2->Chi2perNDF()>fChi2Track[1]) continue;
+           if (mu2->GetChi2MatchTrigger()<fChi2MatchTrig[0] || mu2->GetChi2MatchTrigger()>fChi2MatchTrig[1]) continue;
+           Double_t chargemu2 = mu2->Charge();
+           Double_t pxmu2 = mu2->Px();
+           Double_t pymu2 = mu2->Py();
+           Double_t pzmu2 = mu2->Pz();
+           Double_t emu2 = mu2->E();
+           Double_t pmu2 = mu2->P();
+           Double_t ptmu2 = mu2->Pt();
+           //Double_t rapiditymu2 = Rap(emu2,pzmu2); 
+           Double_t thetamu2 = (180./TMath::Pi())*mu2->Theta();
+           Double_t Rdcamu2 = mu2->DCA();
+           if(chargemu2>0){
+           if (mu1->GetMatchTrigger()>0) printf("Mu1: charge: %f, match: %d\n",chargemu1,1);
+           else  printf("Mu1: charge: %f, match: %d\n",chargemu1,0);
+           if (mu2->GetMatchTrigger()>0) printf("Mu2: charge: %f, match: %d\n",chargemu2,1);
+           else printf("Mu2: charge: %f, match: %d\n",chargemu2,0);
+             ((TH1D*)(fOutput->FindObject("hdca")))->Fill(Rdcamu2);
+             ((TH1D*)(fOutput->FindObject("hdca")))->Fill(Rdcamu1);
+             Double_t TrigCondition=0;
+             if (mu1->GetMatchTrigger()>=mu2->GetMatchTrigger()) TrigCondition = mu1->GetMatchTrigger()+0.1*mu2->GetMatchTrigger();
+               else TrigCondition = mu2->GetMatchTrigger()+0.1*mu1->GetMatchTrigger();     
+             containerInput[0] = (Double_t) Rap((emu1+emu2),(pzmu1+pzmu2));  
+             ((TH1D*)(fOutput->FindObject("ydimuREC")))->Fill(containerInput[0]);
+             ((TH2D*)(fOutput->FindObject("hdcay")))->Fill(TMath::Max(Rdcamu1,Rdcamu2),containerInput[0]);
+             //((TH2D*)(fOutput->FindObject("hdcay")))->Fill(Rdcamu2,containerInput[0]);
+             containerInput[1] = (Double_t) TMath::Sqrt((pxmu1+pxmu2)*(pxmu1+pxmu2)+(pymu1+pymu2)*(pymu1+pymu2));
+             ((TH1D*)(fOutput->FindObject("ptdimuREC")))->Fill(containerInput[1]);
+             containerInput[2] = CostHE(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2);
+             ((TH1D*)(fOutput->FindObject("costHEdimuREC")))->Fill(containerInput[2]);
+             if(containerInput[2]==666) ((TH1D*)(fOutput->FindObject("costHEdimuFail")))->Fill(containerInput[2]);
+             containerInput[3] = TMath::Abs(PhiHE(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2));
+             ((TH1D*)(fOutput->FindObject("phiHEdimuREC")))->Fill(containerInput[3]);
+             if(containerInput[3]==666) ((TH1D*)(fOutput->FindObject("phiHEdimuFail")))->Fill(containerInput[3]);
+             containerInput[4] = CostCS(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2);
+             ((TH1D*)(fOutput->FindObject("costCSdimuREC")))->Fill(containerInput[4]);
+             if(containerInput[4]==666) ((TH1D*)(fOutput->FindObject("costCSdimuFail")))->Fill(containerInput[4]);
+             containerInput[5] = TMath::Abs(PhiCS(pxmu1,pymu1,pzmu1,emu1,chargemu1,pxmu2,pymu2,pzmu2,emu2));
+             ((TH1D*)(fOutput->FindObject("phiCSdimuREC")))->Fill(containerInput[5]);
+             if(containerInput[5]==666) ((TH1D*)(fOutput->FindObject("phiCSdimuFail")))->Fill(containerInput[5]);
+             containerInput[6] = Imass(emu1,pxmu1,pymu1,pzmu1,emu2,pxmu2,pymu2,pzmu2);
+             ((TH1D*)(fOutput->FindObject("imassTot")))->Fill(containerInput[6]);
+             containerInput[7] = TrigCondition+0.05;
+             ((TH1D*)(fOutput->FindObject("trigCond")))->Fill(containerInput[7]);
+             if (ptmu1<ptmu2) {
+               containerInput[8]=ptmu1; 
+               containerInput[9]=ptmu2;
+             } else {
+               containerInput[8]=ptmu2; 
+               containerInput[9]=ptmu1;
+             }
+             if (thetamu1<thetamu2) {
+               containerInput[10]=thetamu1; 
+               containerInput[11]=thetamu2;
+             } else {
+               containerInput[10]=thetamu2; 
+               containerInput[11]=thetamu1;
+             }
+             if (pmu1<pmu2) {
+               containerInput[12]=pmu1; 
+               containerInput[13]=pmu2;
+             } else {
+               containerInput[12]=pmu2; 
+               containerInput[13]=pmu1;
+             }
+             if (fDistinguishTrigClass && trigside==0) containerInput[14]=0.;
+             else if (fDistinguishTrigClass && trigside==1) containerInput[14]=1.;
+             else if (fDistinguishTrigClass && trigside==2) containerInput[14]=2.;
+             else if (fDistinguishTrigClass && trigside==3) containerInput[14]=3.;
+             else containerInput[14]=0.;
+             
+             if(containerInput[2]!=666 && containerInput[3]!=666 && containerInput[4]!=666 && containerInput[5]!=666){
+               if(fDistinguishTrigClass){
+                 if (trigfired==1) fCFManager->GetParticleContainer()->Fill(containerInput,5);
+                 else if (trigfired==0) fCFManager->GetParticleContainer()->Fill(containerInput,1);
+               } else {
+                 fCFManager->GetParticleContainer()->Fill(containerInput,1);
+                 if(fIsAccProduction){
+                   if (cutAccept ==1) fCFManager->GetParticleContainer()->Fill(containerInput,3);
+                 }
+               }
+             } else if (containerInput[0]==666) ((TH1D*)(fOutput->FindObject("ydimuFail")))->Fill(containerInput[0]);
+           
+         }
+
+       }
+
+      }
+
+    } 
+    }
+    }
+   
+
+  }
+
+
+//  ----------
+  if (fReadMCInfo) ((TH1D*)(fOutput->FindObject("JpsiMult")))->Fill(NumbJpsis);
+
+  PostData(1,fOutput) ;
+  PostData(2,fCFManager->GetParticleContainer()) ;
+  
+  
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskDimuonCFContainerBuilder::Imass(Double_t e1, Double_t px1, Double_t py1, Double_t pz1,
+                                  Double_t e2, Double_t px2, Double_t py2, Double_t pz2) 
+{
+// invariant mass calculation
+    Double_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+
+                                    (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2)));
+    return imassrec;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskDimuonCFContainerBuilder::Rap(Double_t e, Double_t pz) 
+{
+// calculate rapidity
+    Double_t rap;
+    if(e>TMath::Abs(pz)){                         // in order to avoid problems with AODs
+       rap = 0.5*TMath::Log((e+pz)/(e-pz));
+       return rap;
+    }
+    else{
+       rap = 666.;
+       return rap;
+    }
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskDimuonCFContainerBuilder::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)
+{
+// Cosine of the theta decay angle (mu+) in the Collins-Soper frame
+
+  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)); 
+  //
+  // --- 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();
+  if(beta.Mag()>=1) return 666.;
+  PMu1Dimu=PMu1CM;
+  PMu2Dimu=PMu2CM;
+  PProjDimu=PProjCM;
+  PTargDimu=PTargCM;
+  PMu1Dimu.Boost(beta);
+  PMu2Dimu.Boost(beta);
+  PProjDimu.Boost(beta);
+  PTargDimu.Boost(beta);
+  
+  //Debugging part -------------------------------------
+  Double_t debugProj[4]={0.,0.,0.,0.};
+  Double_t debugTarg[4]={0.,0.,0.,0.};
+  Double_t debugMu1[4]={0.,0.,0.,0.};
+  Double_t debugMu2[4]={0.,0.,0.,0.};
+  PMu1Dimu.GetXYZT(debugMu1);
+  PMu2Dimu.GetXYZT(debugMu2);
+  PProjDimu.GetXYZT(debugProj);
+  PTargDimu.GetXYZT(debugTarg);
+  if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666; 
+  if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666; 
+  if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666; 
+  if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666; 
+  //----------------------------------------------------
+
+  // --- 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;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskDimuonCFContainerBuilder::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)
+{
+// Cosine of the theta decay angle (mu+) in the Helicity frame
+  
+  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();
+  if(beta.Mag()>=1) return 666.;
+  PMu1Dimu=PMu1CM;
+  PMu2Dimu=PMu2CM;
+  PMu1Dimu.Boost(beta);
+  PMu2Dimu.Boost(beta);
+  
+  //Debugging part -------------------------------------
+  Double_t debugMu1[4]={0.,0.,0.,0.};
+  Double_t debugMu2[4]={0.,0.,0.,0.};
+  PMu1Dimu.GetXYZT(debugMu1);
+  PMu2Dimu.GetXYZT(debugMu2);
+  if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666; 
+  if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666; 
+  //----------------------------------------------------
+  // --- 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;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskDimuonCFContainerBuilder::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)
+{
+// Phi decay angle (mu+) in the Collins-Soper frame
+
+   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.,-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);
+   
+   // --- 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();
+   if(beta.Mag()>=1) return 666.;
+   PMu1Dimu=PMu1CM;
+   PMu2Dimu=PMu2CM;
+   PProjDimu=PProjCM;
+   PTargDimu=PTargCM;
+   PMu1Dimu.Boost(beta);
+   PMu2Dimu.Boost(beta);
+   PProjDimu.Boost(beta);
+   PTargDimu.Boost(beta);
+
+   //Debugging part -------------------------------------
+   Double_t debugProj[4]={0.,0.,0.,0.};
+   Double_t debugTarg[4]={0.,0.,0.,0.};
+   Double_t debugMu1[4]={0.,0.,0.,0.};
+   Double_t debugMu2[4]={0.,0.,0.,0.};
+   PMu1Dimu.GetXYZT(debugMu1);
+   PMu2Dimu.GetXYZT(debugMu2);
+   PProjDimu.GetXYZT(debugProj);
+   PTargDimu.GetXYZT(debugTarg);
+   if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666; 
+   if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666; 
+   if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666; 
+   if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666; 
+   //----------------------------------------------------
+
+   // --- 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=0.;
+   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)));
+   }
+   if (phi>TMath::Pi()) phi=phi-TMath::Pi();
+   
+   return phi;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskDimuonCFContainerBuilder::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)
+{
+// Phi decay angle (mu+) in the Helicity frame
+  TLorentzVector PMu1Lab, PMu2Lab, PProjLab, PTargLab, PDimuLab; // In the lab. 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 LAB frame
+  PMu1Lab.SetPxPyPzE(px1,py1,pz1,e1);
+  PMu2Lab.SetPxPyPzE(px2,py2,pz2,e2);
+  
+  // --- Obtain the dimuon parameters in the LAB frame
+  PDimuLab=PMu1Lab+PMu2Lab;
+  zaxis=(PDimuLab.Vect()).Unit();
+  
+  // --- Translate the muon parameters in the dimuon rest frame
+  beta=(-1./PDimuLab.E())*PDimuLab.Vect();
+  if(beta.Mag()>=1.) return 666.;
+
+  PProjLab.SetPxPyPzE(0.,0.,-fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // proiettile
+  PTargLab.SetPxPyPzE(0.,0.,fBeamEnergy,TMath::Sqrt(fBeamEnergy*fBeamEnergy+mp*mp)); // bersaglio
+
+  PProjDimu=PProjLab;
+  PTargDimu=PTargLab;
+
+  PProjDimu.Boost(beta);
+  PTargDimu.Boost(beta);
+  
+  yaxis=((PProjDimu.Vect()).Cross(PTargDimu.Vect())).Unit();
+  xaxis=(yaxis.Cross(zaxis)).Unit();
+  
+  PMu1Dimu=PMu1Lab;
+  PMu2Dimu=PMu2Lab;
+  PMu1Dimu.Boost(beta);
+  PMu2Dimu.Boost(beta);
+  
+  //Debugging part -------------------------------------
+  Double_t debugProj[4]={0.,0.,0.,0.};
+  Double_t debugTarg[4]={0.,0.,0.,0.};
+  Double_t debugMu1[4]={0.,0.,0.,0.};
+  Double_t debugMu2[4]={0.,0.,0.,0.};
+  PMu1Dimu.GetXYZT(debugMu1);
+  PMu2Dimu.GetXYZT(debugMu2);
+  PProjDimu.GetXYZT(debugProj);
+  PTargDimu.GetXYZT(debugTarg);
+  if (debugProj[0]!=debugProj[0] ||debugProj[1]!=debugProj[1] || debugProj[2]!=debugProj[2] ||debugProj[3]!=debugProj[3]) return 666; 
+  if (debugTarg[0]!=debugTarg[0] ||debugTarg[1]!=debugTarg[1] || debugTarg[2]!=debugTarg[2] ||debugTarg[3]!=debugTarg[3]) return 666; 
+  if (debugMu1[0]!=debugMu1[0] ||debugMu1[1]!=debugMu1[1] || debugMu1[2]!=debugMu1[2] ||debugMu1[3]!=debugMu1[3]) return 666; 
+  if (debugMu2[0]!=debugMu2[0] ||debugMu2[1]!=debugMu2[1] || debugMu2[2]!=debugMu2[2] ||debugMu2[3]!=debugMu2[3]) return 666; 
+  //----------------------------------------------------
+  
+  Double_t phi=0.;
+   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 AliAnalysisTaskDimuonCFContainerBuilder::Terminate(Option_t *) 
+{
+}
+
+#endif
diff --git a/PWG3/muon/AliAnalysisTaskDimuonCFContainerBuilder.h b/PWG3/muon/AliAnalysisTaskDimuonCFContainerBuilder.h
new file mode 100644 (file)
index 0000000..2f2a8ac
--- /dev/null
@@ -0,0 +1,98 @@
+#ifndef ALIANALYSISTASKDIMUONCFCONTAINERBUILDER_H
+#define ALIANALYSISTASKDIMUONCFCONTAINERBUILDER_H
+
+#include "AliAnalysisTaskSE.h"
+#include "AliCFContainer.h"
+#include "TMath.h"
+#include "TString.h"
+#include "AliCFManager.h"
+
+//     Analysis task for the building of a dimuon CF container
+//     Also some single-muon variables are stored
+//     L. Bianchi - Universita' & INFN Torino
+
+
+
+class AliAnalysisTaskDimuonCFContainerBuilder : public AliAnalysisTaskSE {
+  public:
+
+  AliAnalysisTaskDimuonCFContainerBuilder();
+  AliAnalysisTaskDimuonCFContainerBuilder(const Char_t* name, Bool_t readaod, Bool_t readMC, Bool_t isaccept, Double_t beamEn);
+  AliAnalysisTaskDimuonCFContainerBuilder& operator= (const AliAnalysisTaskDimuonCFContainerBuilder& c);
+  AliAnalysisTaskDimuonCFContainerBuilder(const AliAnalysisTaskDimuonCFContainerBuilder& c);
+  virtual ~AliAnalysisTaskDimuonCFContainerBuilder();
+
+  // ANALYSIS FRAMEWORK STUFF to loop on data and fill output objects
+  void     UserExec(Option_t *option);
+  void     Terminate(Option_t *);
+  void     UserCreateOutputObjects();
+  
+  // 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;}
+
+  // Setters and Getters
+  Bool_t IsReadAODData   ()                    const {return fReadAODData;}
+  void   SetReadAODData        (Bool_t flag=kTRUE)       {fReadAODData=flag;}
+  void  SetReadMCinfo          (Bool_t flag=kTRUE)       {fReadMCInfo=flag;}
+  void   SetIsAccProd          (Bool_t flag=kTRUE)       {fIsAccProduction=flag;}
+  void  SetBeamEnergy          (Double_t en)             {fBeamEnergy=en;}
+  void   SetChi2Limits         (Double_t chi2track[])    {fChi2Track[0]=chi2track[0];fChi2Track[1]=chi2track[1];}
+  void   SetChi2MatchLimits    (Double_t chi2match[])    {fChi2MatchTrig[0]=chi2match[0];fChi2MatchTrig[1]=chi2match[1];}
+  void   SetPtSingMuLimits     (Double_t PtSingle[])     {fPtSingMuCut[0]=PtSingle[0];fPtSingMuCut[1]=PtSingle[1];}
+  void   SetThetaSingMuLimits  (Double_t ThetaSingle[])  {fThetaSingMuCut[0]=ThetaSingle[0];fThetaSingMuCut[1]=ThetaSingle[1];}
+  void   SetZprimVertLimits    (Double_t Zprimvtx[])     {fzPrimVertexSPD[0]=Zprimvtx[0];fzPrimVertexSPD[1]=Zprimvtx[1];}
+  void  SetCutonZvtxSPD        (Bool_t   cut=kFALSE)     {fCutOnzVtxSPD=cut;}
+  void   SetNContributorsLimits        (Double_t NContr[])       {fNContributors[0]=NContr[0];fNContributors[1]=NContr[1];}
+  void  SetCutonNContributors  (Bool_t   cut=kFALSE)     {fCutOnNContributors=cut;}
+  void  SetDistinguishTrigClass(Bool_t   dist=kFALSE)    {fDistinguishTrigClass=dist;}
+  void   SetTrigClassMuonName  (TString name = "CMU")    {fTrigClassMuon=name;}
+  void   SetTrigClassInteracName(TString name = "CINT")          {fTrigClassInteraction=name;}
+  void   SetTrigClassMuonSideName(TString name[4])       {for(Int_t i=0;i<4;i++) fTrigClassMuonSide[i]=name[i];}
+  void   SetTrigClassInteracSideName(TString name[4])    {for(Int_t i=0;i<4;i++) fTrigClassInteractionSide[i]=name[i];}
+ protected:
+  
+  Bool_t               fReadAODData            ;    // flag for AOD/ESD input files
+  Bool_t               fReadMCInfo             ;    // flag for reading MC info (ESD->Kinematics, AOD->MCbranch)
+  Bool_t               fIsAccProduction        ;    // flag to activate in case of acceptance MC production (in this case fReadMCInfo==kTRUE)
+  AliCFManager         *fCFManager             ;    // pointer to the CF manager
+  TList                *fQAHistList            ;    // list of QA histograms
+  Double_t             fNevt                   ;    // event counter
+  Double_t             fBeamEnergy             ;    // Energy of the beam (required for the CS angle)
+  TList                *fOutput                ;    // list of TH in output
+  
+                                                    // CUTS ON TRACKS
+  Double_t             fChi2Track[2]           ;    // Cut on chi2 of the tracks ([0]==chi2min, [1]==chi2max)
+  Double_t             fChi2MatchTrig[2]       ;    // Cut on chi2matchtrigger of the tracks ([0]==chi2Matchmin, [1]==chi2Matchmax)
+  Double_t             fPtSingMuCut[2]         ;    // Cut on pt of single-mu tracks ([0]==ptmin, [1]==ptmax)
+  Double_t             fThetaSingMuCut[2]      ;    // Cut on polar angle (wrt beam axis) of single-mu tracks ([0]==thetamin, [1]==thetamax)
+
+                                                    // CUTS ON EVENT
+  Double_t             fzPrimVertexSPD[2]      ;    // Cut on the z coordinate of the primary vertex in SPD (full ITS for AODs)
+  Bool_t               fCutOnzVtxSPD           ;    // flag to activate the cut on the z of the primary vertex
+  Int_t                        fNContributors[2]       ;    // Cut on NContributors in SPD
+  Bool_t               fCutOnNContributors     ;    // flag to activate the cut on NContributors in SPD
+
+                                                    // CUTS ON THE FIRED TRIGGER CLASS
+  TString              fTrigClassMuon          ;    // name of the muon trigger class (CMU by default)
+  TString              fTrigClassInteraction   ;    // name of the interaction trigger class (CINT by default)
+  TString              fTrigClassMuonSide[4]   ;    // name of the muon trigger classes containing the side
+  TString              fTrigClassInteractionSide[4];// name of the interaction trigger classes containing the side
+  Bool_t               fDistinguishTrigClass   ;    // flag to activate the cut on the fired trigger class
+  
+  
+  
+  Double_t Imass  (Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t);
+  Double_t Rap   (Double_t, Double_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);
+  Double_t PhiCS  (Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t);
+  Double_t PhiHE  (Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t);
+  
+  ClassDef(AliAnalysisTaskDimuonCFContainerBuilder,1);
+};
+
+#endif