From: martinez Date: Thu, 2 Jul 2009 14:05:32 +0000 (+0000) Subject: Classes for efficiency corrections (Xaver, Roberta) X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=7cd8a4ce948ff356f7f0c34075d8e8f005b768db Classes for efficiency corrections (Xaver, Roberta) --- diff --git a/PWG3/CMake_libPWG3muon.txt b/PWG3/CMake_libPWG3muon.txt index 331c65adc32..82aa2f08514 100644 --- a/PWG3/CMake_libPWG3muon.txt +++ b/PWG3/CMake_libPWG3muon.txt @@ -12,6 +12,7 @@ set(SRCS muon/AliAODEventInfo.cxx muon/AliESDMuonTrackCuts.cxx muon/AliAnalysisTaskSingleMuESD.cxx + muon/AliCFMuonResTask1.cxx ) # fill list of header files from list of source files diff --git a/PWG3/PWG3muonLinkDef.h b/PWG3/PWG3muonLinkDef.h index ef76b856bc7..2fa8994a07a 100644 --- a/PWG3/PWG3muonLinkDef.h +++ b/PWG3/PWG3muonLinkDef.h @@ -15,6 +15,7 @@ #pragma link C++ class AliESDMuonTrackCuts+; #pragma link C++ class AliAnalysisTaskSingleMuESD+; #pragma link C++ class AliAnalysisTaskLinkToMC+; +#pragma link C++ class AliCFMuonResTask1+; #endif diff --git a/PWG3/libPWG3muon.pkg b/PWG3/libPWG3muon.pkg index 44fd78ec84a..5b13ef15abf 100644 --- a/PWG3/libPWG3muon.pkg +++ b/PWG3/libPWG3muon.pkg @@ -10,7 +10,8 @@ SRCS:= muon/AliAnalysisTaskESDMuonFilter.cxx \ muon/AliAODDimuon.cxx \ muon/AliAODEventInfo.cxx \ muon/AliESDMuonTrackCuts.cxx \ - muon/AliAnalysisTaskSingleMuESD.cxx + muon/AliAnalysisTaskSingleMuESD.cxx \ + muon/AliCFMuonResTask1.cxx HDRS:= $(SRCS:.cxx=.h) diff --git a/PWG3/muon/AliCFMuonResTask1.C b/PWG3/muon/AliCFMuonResTask1.C new file mode 100644 index 00000000000..682b0276072 --- /dev/null +++ b/PWG3/muon/AliCFMuonResTask1.C @@ -0,0 +1,307 @@ +// 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; //JPsi + +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 = 12 ; +const Double_t mymin = -4 ; +const Double_t mymax = -2.5 ; +const Double_t costCSmin = -1.; +const Double_t costCSmax = 1.; +const Double_t costHEmin = -1.; +const Double_t costHEmax = 1.; +const Double_t phiCSmin = -TMath::Pi(); +const Double_t phiCSmax = TMath::Pi(); +const Double_t phiHEmin = -TMath::Pi(); +const Double_t phiHEmax = TMath::Pi(); + +//---------------------------------------------------- + +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("/alidata/alice/arnaldi/CORRFW/polar-1.1000.1000/AliESDs.root"); + analysisChain->Add("/alidata/alice/arnaldi/CORRFW/polar-1.1001.1000/AliESDs.root"); + analysisChain->Add("/alidata/alice/arnaldi/CORRFW/polar-1.1002.1000/AliESDs.root"); + } + } + +///// END INPUT + + + Info("AliCFMuonResTask1",Form("CHAIN HAS %d ENTRIES",(Int_t)analysisChain->GetEntries())); + + // CONTAINER DEFINITION + Info("AliCFMuonResTask1","SETUP CONTAINER"); + + // The sensitive variables (13 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; + UInt_t costCS = 9; + UInt_t phiCS = 10; + UInt_t costHE = 11; + UInt_t phiHE = 12; + + + // Setting up the container grid + UInt_t nstep = 2 ; //number of selection steps : MC and ESD + const Int_t nvar = 13 ; //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 = 50 ; + const Int_t nbin8 = 50 ; + const Int_t nbin9 = 50 ; + const Int_t nbin10 = 20 ; + const Int_t nbin11 = 20 ; + const Int_t nbin12 = 20 ; + const Int_t nbin13 = 20 ; + + // 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; + iBin[9]=nbin10; + iBin[10]=nbin11; + iBin[11]=nbin12; + iBin[12]=nbin13; + + // 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]; + 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]; + + // 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 ; + for(Int_t i=0; i<=nbin10; i++) binLim10[i]=(Double_t)costCSmin + (costCSmax-costCSmin)/nbin10*(Double_t)i ; + for(Int_t i=0; i<=nbin11; i++) binLim11[i]=(Double_t)phiCSmin + (phiCSmax-phiCSmin)/nbin11*(Double_t)i ; + for(Int_t i=0; i<=nbin12; i++) binLim12[i]=(Double_t)costHEmin + (costHEmax-costHEmin)/nbin12*(Double_t)i ; + for(Int_t i=0; i<=nbin13; i++) binLim13[i]=(Double_t)phiHEmin + (phiHEmax-phiHEmin)/nbin13*(Double_t)i ; + + // one container of 2 steps (MC and ESD) with 12 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); + container -> SetBinLimits(costCS,binLim10); + container -> SetBinLimits(phiCS,binLim11); + container -> SetBinLimits(costHE,binLim12); + container -> SetBinLimits(phiHE,binLim13); + + // 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->SetRapidityRange(-4,-2.5); + 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::kPartAccCuts,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,mgr->GetCommonInputContainer()); + + 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/PWG3/muon/AliCFMuonResTask1.cxx b/PWG3/muon/AliCFMuonResTask1.cxx new file mode 100644 index 00000000000..356a6b3066d --- /dev/null +++ b/PWG3/muon/AliCFMuonResTask1.cxx @@ -0,0 +1,593 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +//----------------------------------------------------------------------- +// Example of task (running locally, on AliEn and CAF), +// which provides standard way of calculating acceptance and efficiency +// between different steps of the procedure. +// The ouptut of the task is a AliCFContainer from which the efficiencies +// can be calculated +//----------------------------------------------------------------------- +// Author : R. Vernet, Consorzio Cometa - Catania (it) +//----------------------------------------------------------------------- +// Modification done by X. Lopez - LPC Clermont (fr) +//----------------------------------------------------------------------- + + +#ifndef ALICFMUONRESTASK1_CXX +#define ALICFMUONRESTASK1_CXX + +#include "AliCFMuonResTask1.h" +#include "AliHeader.h" +#include "AliESDHeader.h" +#include "AliStack.h" +#include "TParticle.h" +#include "TLorentzVector.h" +#include "TH1I.h" +#include "AliMCEvent.h" +#include "AliAnalysisManager.h" +#include "AliESDEvent.h" +#include "AliAODEvent.h" +#include "AliCFManager.h" +#include "AliCFCutBase.h" +#include "AliCFContainer.h" +#include "TChain.h" +#include "AliESDtrack.h" +#include "AliLog.h" +#include "AliESDMuonTrack.h" +#include "AliESDtrack.h" +#include "AliESDInputHandler.h" +#include "TCanvas.h" + +ClassImp(AliCFMuonResTask1) + +//__________________________________________________________________________ +AliCFMuonResTask1::AliCFMuonResTask1() : + fReadAODData(0), + fCFManager(0x0), + fQAHistList(0x0), + fHistEventsProcessed(0x0), + fNevt(0) +{ + // + //Default ctor + // +} +//___________________________________________________________________________ +AliCFMuonResTask1::AliCFMuonResTask1(const Char_t* name) : + AliAnalysisTaskSE(name), + fReadAODData(0), + fCFManager(0x0), + fQAHistList(0x0), + fHistEventsProcessed(0x0), + fNevt(0) +{ + // + // Constructor. Initialization of Inputs and Outputs + // + Info("AliCFMuonResTask1","Calling Constructor"); + + fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ; + + DefineOutput(1,TH1I::Class()); + DefineOutput(2,AliCFContainer::Class()); + +} + +//___________________________________________________________________________ +AliCFMuonResTask1& AliCFMuonResTask1::operator=(const AliCFMuonResTask1& c) +{ + // + // Assignment operator + // + if (this!=&c) { + AliAnalysisTaskSE::operator=(c) ; + fReadAODData = c.fReadAODData ; + fCFManager = c.fCFManager; + fQAHistList = c.fQAHistList ; + fHistEventsProcessed = c.fHistEventsProcessed; + fNevt = c.fNevt ; + } + return *this; +} + +//___________________________________________________________________________ +AliCFMuonResTask1::AliCFMuonResTask1(const AliCFMuonResTask1& c) : + AliAnalysisTaskSE(c), + fReadAODData(c.fReadAODData), + fCFManager(c.fCFManager), + fQAHistList(c.fQAHistList), + fHistEventsProcessed(c.fHistEventsProcessed), + fNevt(c.fNevt) +{ + // + // Copy Constructor + // +} + +//___________________________________________________________________________ +AliCFMuonResTask1::~AliCFMuonResTask1() { + // + //destructor + // + Info("~AliCFMuonResTask1","Calling Destructor"); + if (fCFManager) delete fCFManager ; + if (fHistEventsProcessed) delete fHistEventsProcessed ; + if (fQAHistList) {fQAHistList->Clear(); delete fQAHistList;} +} + +//_________________________________________________ +void AliCFMuonResTask1::UserExec(Option_t *) +{ + // + // Main loop function + // + + Info("UserExec","") ; + if (!fMCEvent) { + Error("UserExec","NO MC EVENT FOUND!"); + return; + } + + fNevt++; + Double_t containerInput[13] ; + Double_t BeamEnergy=7000; + +//////// +//// MC +//////// + + fCFManager->SetEventInfo(fMCEvent); + AliStack* stack = fMCEvent->Stack(); + + // loop on the MC event + for (Int_t ipart=0; ipartGetNumberOfTracks(); ipart++) { + AliMCParticle *mcPart = fMCEvent->GetTrack(ipart); + + TParticle *part = mcPart->Particle(); + TParticle *part0 = mcPart->Particle(); + TParticle *part1 = mcPart->Particle(); + + // Selection of the resonance + if (!fCFManager->CheckParticleCuts(AliCFManager::kPartGenCuts,mcPart)) continue; + + // Mother kinematics + Float_t e = part->Energy(); + Float_t pz = part->Pz(); + //Float_t py = part->Py(); + //Float_t px = part->Px(); + Float_t rapmc = Rap(e,pz); + //Float_t mass = part->GetCalcMass(); + + // Decays kinematics + + Int_t p0 = part->GetDaughter(0); + part0 = stack->Particle(p0); + // selection of the rapidity for first muon + AliMCParticle *mcpart0 = new AliMCParticle(part0); + if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart0)) continue; + Int_t pdg0 = part0->GetPdgCode(); + + Int_t p1 = part->GetDaughter(1); + part1 = stack->Particle(p1); + // selection of the rapidity for second muon + AliMCParticle *mcpart1 = new AliMCParticle(part1); + if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mcpart1)) continue; + Int_t pdg1 = part1->GetPdgCode(); + + // 0 mu- 1 mu+ + Float_t e0=-999, pz0=-999, py0=-999, px0=-999, phi0=-999, rapmc0=-999; + Float_t e1=-999, pz1=-999, py1=-999, px1=-999, phi1=-999, rapmc1=-999; + Double_t charge0=-999, charge1=-999; + + // ordering sign: first = mu- + if(pdg0==13){ + e0 = part0->Energy(); + pz0 = part0->Pz(); + py0 = part0->Py(); + px0 = part0->Px(); + phi0 = part0->Phi(); + phi0 = Phideg(phi0); + rapmc0=Rap(e0,pz0); + charge0 = part0->GetPDG()->Charge()/3; + + e1 = part1->Energy(); + pz1 = part1->Pz(); + py1 = part1->Py(); + px1 = part1->Px(); + phi1 = part1->Phi(); + phi1 = Phideg(phi1); + rapmc1=Rap(e1,pz1); + charge1 = part1->GetPDG()->Charge()/3; + } + else{ + if(pdg0==-13){ + e1 = part0->Energy(); + pz1 = part0->Pz(); + py1 = part0->Py(); + px1 = part0->Px(); + phi1 = part0->Phi(); + phi1 = Phideg(phi1); + rapmc1=Rap(e1,pz1); + charge1 = part0->GetPDG()->Charge()/3; + + e0 = part1->Energy(); + pz0 = part1->Pz(); + py0 = part1->Py(); + px0 = part1->Px(); + phi0 = part1->Phi(); + phi0 = Phideg(phi0); + rapmc0=Rap(e0,pz0); + charge0 = part1->GetPDG()->Charge()/3; + } + } + + if(pdg0==13 || pdg1==13) { + + Float_t pmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1)+ + (pz0+pz1)*(pz0+pz1)); + Float_t ptmc = TMath::Sqrt((px0+px1)*(px0+px1)+(py0+py1)*(py0+py1)); + Float_t imassmc = Imass(e0,px0,py0,pz0,e1,px1,py1,pz1); + //Float_t rapmc_check=Rap((e0+e1),(pz0+pz1)); + + Double_t costCS = CostCS((Double_t) px0,(Double_t) py0,(Double_t) pz0,(Double_t) e0, (Double_t)charge0,(Double_t) px1,(Double_t) py1,(Double_t) pz1,(Double_t) e1, BeamEnergy); + Double_t costHE = CostHE((Double_t) px0,(Double_t) py0,(Double_t) pz0,(Double_t) e0, (Double_t)charge0,(Double_t) px1,(Double_t) py1,(Double_t) pz1,(Double_t) e1); + Double_t phiCS = PhiCS((Double_t) px0,(Double_t) py0,(Double_t) pz0,(Double_t) e0, (Double_t)charge0,(Double_t) px1,(Double_t) py1,(Double_t) pz1,(Double_t) e1, BeamEnergy); + Double_t phiHE = PhiHE((Double_t) px0,(Double_t) py0,(Double_t) pz0,(Double_t) e0, (Double_t)charge0,(Double_t) px1,(Double_t) py1,(Double_t) pz1,(Double_t) e1, BeamEnergy); + + containerInput[0] = fNevt+0.5 ; + containerInput[1] = rapmc0 ; + containerInput[2] = phi0 ; + containerInput[3] = rapmc1 ; + containerInput[4] = phi1 ; + containerInput[5] = imassmc ; + containerInput[6] = rapmc ; + containerInput[7] = ptmc; + containerInput[8] = pmc; + containerInput[9] = costCS; + containerInput[10] = phiCS; + containerInput[11] = costHE; + containerInput[12] = phiHE; + + // fill the container at the first step + fCFManager->GetParticleContainer()->Fill(containerInput,0); + } + } + +//////// +//// ESD +//////// + + AliESDEvent *fESD; + AliESDInputHandler *esdH = dynamic_cast + (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); + fESD = esdH->GetEvent(); + Int_t mult1 = fESD->GetNumberOfMuonTracks() ; + + + for (Int_t j = 0; jGetMuonTrack(j))); + Float_t zr1 = mu1->Charge(); +// Select mu- + if(zr1<0){ + Float_t pxr1 = mu1->Px(); + Float_t pyr1 = mu1->Py(); + Float_t pzr1 = mu1->Pz(); + Float_t phir1 = mu1->Phi(); + phir1=Phideg(phir1); + Float_t er1 = mu1->E(); + Float_t rap1=Rap(er1,pzr1); + // selection of the rapidity for first muon + if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mu1)) continue; + + for (Int_t jj = 0; jjGetMuonTrack(jj))); + Float_t zr2 = mu2->Charge(); +// Select mu+ + if(zr2>0 && jj !=j){ + Float_t pxr2 = mu2->Px(); + Float_t pyr2 = mu2->Py(); + Float_t pzr2 = mu2->Pz(); + Float_t phir2 = mu2->Phi(); + phir2 = Phideg(phir2); + Float_t er2 = mu2->E(); + Float_t rap2=Rap(er2,pzr2); + // selection of the rapidity for second muon + if (!fCFManager->CheckParticleCuts(AliCFManager::kPartAccCuts,mu2)) continue; + + Float_t prec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2)+ + (pzr1+pzr2)*(pzr1+pzr2)); + Float_t ptrec = TMath::Sqrt((pxr1+pxr2)*(pxr1+pxr2)+(pyr1+pyr2)*(pyr1+pyr2)); + Float_t raprec=Rap((er1+er2),(pzr1+pzr2)); + Float_t imassrec = Imass(er1,pxr1,pyr1,pzr1,er2,pxr2,pyr2,pzr2); + + Double_t costCSrec = CostCS((Double_t) pxr1,(Double_t) pyr1,(Double_t)pzr1,(Double_t) er1, (Double_t)zr1,(Double_t) pxr2,(Double_t) pyr2,(Double_t)pzr2,(Double_t) er2, BeamEnergy); + Double_t costHErec = CostHE((Double_t) pxr1,(Double_t) pyr1,(Double_t)pzr1,(Double_t) er1, (Double_t)zr1,(Double_t) pxr2,(Double_t) pyr2,(Double_t)pzr2,(Double_t) er2); + Double_t phiCSrec = PhiCS((Double_t) pxr1,(Double_t) pyr1,(Double_t)pzr1,(Double_t) er1, (Double_t)zr1,(Double_t) pxr2,(Double_t) pyr2,(Double_t)pzr2,(Double_t) er2, BeamEnergy); + Double_t phiHErec = PhiHE((Double_t) pxr1,(Double_t) pyr1,(Double_t)pzr1,(Double_t) er1, (Double_t)zr1,(Double_t) pxr2,(Double_t) pyr2,(Double_t)pzr2,(Double_t) er2, BeamEnergy); + + containerInput[0] = fNevt+0.5 ; + containerInput[1] = rap1 ; + containerInput[2] = phir1 ; + containerInput[3] = rap2 ; + containerInput[4] = phir2 ; + containerInput[5] = imassrec ; + containerInput[6] = raprec ; + containerInput[7] = ptrec; + containerInput[8] = prec; + containerInput[9] = costCSrec; + containerInput[10] = phiCSrec; + containerInput[11] = costHErec; + containerInput[12] = phiHErec; + + // fill the container at the second step + fCFManager->GetParticleContainer()->Fill(containerInput,1); + + } // mu+ Selection + } // second mu Loop + } // mu- Selection + } + +// ---------- + fHistEventsProcessed->Fill(0); + PostData(1,fHistEventsProcessed) ; + PostData(2,fCFManager->GetParticleContainer()) ; +} +//________________________________________________________________________ +const Float_t AliCFMuonResTask1::Imass(Float_t e1, Float_t px1, Float_t py1, Float_t pz1, + Float_t e2, Float_t px2, Float_t py2, Float_t pz2) +{ +// invariant mass calculation + Float_t imassrec = TMath::Sqrt((e1+e2)*(e1+e2)-((px1+px2)*(px1+px2)+ + (py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2))); + return imassrec; +} +//________________________________________________________________________ +const Float_t AliCFMuonResTask1::Rap(Float_t e, Float_t pz) +{ +// calculate rapidity + Float_t rap; + if(e!=pz){ + rap = 0.5*TMath::Log((e+pz)/(e-pz)); + return rap; + } + else{ + rap = -200; + return rap; + } +} +//________________________________________________________________________ +const Float_t AliCFMuonResTask1::Phideg(Float_t phi) +{ +// calculate Phi in range [-180,180] + Float_t phideg; + + phideg = phi-TMath::Pi(); + phideg = phideg*57.296; + return phideg; +} +//________________________________________________________________________ +const Double_t AliCFMuonResTask1::CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, +Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, +Double_t Energy) +{ + TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM. frame + TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame + TVector3 beta,zaxisCS; + Double_t mp=0.93827231; + // + // --- Fill the Lorentz vector for projectile and target in the CM frame + // + PProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp)); + PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp)); + // + // --- Get the muons parameters in the CM frame + // + PMu1CM.SetPxPyPzE(px1,py1,pz1,e1); + PMu2CM.SetPxPyPzE(px2,py2,pz2,e2); + // + // --- Obtain the dimuon parameters in the CM frame + // + PDimuCM=PMu1CM+PMu2CM; + // + // --- Translate the dimuon parameters in the dimuon rest frame + // + beta=(-1./PDimuCM.E())*PDimuCM.Vect(); + PMu1Dimu=PMu1CM; + PMu2Dimu=PMu2CM; + PProjDimu=PProjCM; + PTargDimu=PTargCM; + PMu1Dimu.Boost(beta); + PMu2Dimu.Boost(beta); + PProjDimu.Boost(beta); + PTargDimu.Boost(beta); + // + // --- Determine the z axis for the CS angle + // + zaxisCS=(((PProjDimu.Vect()).Unit())-((PTargDimu.Vect()).Unit())).Unit(); + // + // --- Determine the CS angle (angle between mu+ and the z axis defined above) + Double_t cost; + + if(charge1>0) cost = zaxisCS.Dot((PMu1Dimu.Vect()).Unit()); + else cost = zaxisCS.Dot((PMu2Dimu.Vect()).Unit()); + return cost; +} +//________________________________________________________________________ + +//________________________________________________________________________ +const Double_t AliCFMuonResTask1::CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, +Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2) +{ + TLorentzVector PMu1CM, PMu2CM, PDimuCM; // In the CM frame + TLorentzVector PMu1Dimu, PMu2Dimu; // In the dimuon rest frame + TVector3 beta,zaxisCS; + // + // --- Get the muons parameters in the CM frame + // + PMu1CM.SetPxPyPzE(px1,py1,pz1,e1); + PMu2CM.SetPxPyPzE(px2,py2,pz2,e2); + // + // --- Obtain the dimuon parameters in the CM frame + // + PDimuCM=PMu1CM+PMu2CM; + // + // --- Translate the muon parameters in the dimuon rest frame + // + beta=(-1./PDimuCM.E())*PDimuCM.Vect(); + PMu1Dimu=PMu1CM; + PMu2Dimu=PMu2CM; + PMu1Dimu.Boost(beta); + PMu2Dimu.Boost(beta); + // + // --- Determine the z axis for the calculation of the polarization angle (i.e. the direction of the dimuon in the CM system) + // + TVector3 zaxis; + zaxis=(PDimuCM.Vect()).Unit(); + + // --- Calculation of the polarization angle (angle between mu+ and the z axis defined above) + Double_t cost; + if(charge1>0) cost = zaxis.Dot((PMu1Dimu.Vect()).Unit()); + else cost = zaxis.Dot((PMu2Dimu.Vect()).Unit()); + + return cost; +} +//________________________________________________________________________ + +//________________________________________________________________________ +const Double_t AliCFMuonResTask1::PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, +Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, +Double_t Energy) +{ + TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM frame + TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame + TVector3 beta,yaxisCS, xaxisCS, zaxisCS; + Double_t mp=0.93827231; + // + // --- Fill the Lorentz vector for projectile and target in the CM frame + // + PProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp)); + PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp)); + // + // --- Get the muons parameters in the CM frame + // + PMu1CM.SetPxPyPzE(px1,py1,pz1,e1); + PMu2CM.SetPxPyPzE(px2,py2,pz2,e2); + // + // --- Obtain the dimuon parameters in the CM frame + // + PDimuCM=PMu1CM+PMu2CM; + // + // --- Translate the dimuon parameters in the dimuon rest frame + // + beta=(-1./PDimuCM.E())*PDimuCM.Vect(); + PMu1Dimu=PMu1CM; + PMu2Dimu=PMu2CM; + PProjDimu=PProjCM; + PTargDimu=PTargCM; + PMu1Dimu.Boost(beta); + PMu2Dimu.Boost(beta); + PProjDimu.Boost(beta); + PTargDimu.Boost(beta); + // + // --- Determine the z axis for the CS angle + // + zaxisCS=(((PProjDimu.Vect()).Unit())-((PTargDimu.Vect()).Unit())).Unit(); + yaxisCS=(((PProjDimu.Vect()).Unit()).Cross((PTargDimu.Vect()).Unit())).Unit(); + xaxisCS=(yaxisCS.Cross(zaxisCS)).Unit(); + + Double_t phi; + if(charge1>0) phi = TMath::ATan2((PMu1Dimu.Vect()).Dot(yaxisCS),((PMu1Dimu.Vect()).Dot(xaxisCS))); + else phi = TMath::ATan2((PMu2Dimu.Vect()).Dot(yaxisCS),((PMu2Dimu.Vect()).Dot(xaxisCS))); + + return phi; +} +//________________________________________________________________________ + +//________________________________________________________________________ +const Double_t AliCFMuonResTask1::PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, +Double_t charge1, Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t Energy) +{ + TLorentzVector PMu1CM, PMu2CM, PProjCM, PTargCM, PDimuCM; // In the CM frame + TLorentzVector PMu1Dimu, PMu2Dimu, PProjDimu, PTargDimu; // In the dimuon rest frame + TVector3 beta,xaxis,yaxis,zaxis; + Double_t mp=0.93827231; + + // + // --- Get the muons parameters in the CM frame + // + PMu1CM.SetPxPyPzE(px1,py1,pz1,e1); + PMu2CM.SetPxPyPzE(px2,py2,pz2,e2); + // + // --- Obtain the dimuon parameters in the CM frame + // + PDimuCM=PMu1CM+PMu2CM; + // + // --- Translate the muon parameters in the dimuon rest frame + // + zaxis=(PDimuCM.Vect()).Unit(); + + beta=(-1./PDimuCM.E())*PDimuCM.Vect(); + + PProjCM.SetPxPyPzE(0.,0.,-Energy,TMath::Sqrt(Energy*Energy+mp*mp)); + PTargCM.SetPxPyPzE(0.,0.,Energy,TMath::Sqrt(Energy*Energy+mp*mp)); + + PProjDimu=PProjCM; + PTargDimu=PTargCM; + + PProjDimu.Boost(beta); + PTargDimu.Boost(beta); + + yaxis=((PProjDimu.Vect()).Cross(PTargDimu.Vect())).Unit(); + xaxis=(yaxis.Cross(zaxis)).Unit(); + + PMu1Dimu=PMu1CM; + PMu2Dimu=PMu2CM; + PMu1Dimu.Boost(beta); + PMu2Dimu.Boost(beta); + + Double_t phi; + if(charge1>0) phi = TMath::ATan2((PMu1Dimu.Vect()).Dot(yaxis),(PMu1Dimu.Vect()).Dot(xaxis)); + else phi = TMath::ATan2((PMu2Dimu.Vect()).Dot(yaxis),(PMu2Dimu.Vect()).Dot(xaxis)); + + return phi; +} + +//________________________________________________________________________ +void AliCFMuonResTask1::Terminate(Option_t *) +{ + // draw result of the Invariant mass MC and ESD + + AliCFContainer *cont = dynamic_cast (GetOutputData(2)); + + TH1D *kmass = cont->ShowProjection(5,0); + TH1D *rmass = cont->ShowProjection(5,1); + TH1D *kcostCS = cont->ShowProjection(9,0); + TH1D *rcostCS = cont->ShowProjection(9,1); + + TCanvas *c1 = new TCanvas("AliCFMuonResTask1","JPSI MC & ESD",10,10,510,510); + c1->Divide(2,2); + c1->cd(1); + kmass->Draw("HIST"); + c1->cd(2); + rmass->Draw("HIST"); + c1->cd(3); + kcostCS->Draw("HIST"); + c1->cd(4); + rcostCS->Draw("HIST"); +} +//________________________________________________________________________ + +#endif diff --git a/PWG3/muon/AliCFMuonResTask1.h b/PWG3/muon/AliCFMuonResTask1.h new file mode 100644 index 00000000000..1a930470633 --- /dev/null +++ b/PWG3/muon/AliCFMuonResTask1.h @@ -0,0 +1,83 @@ +/************************************************************************** + * 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 TLorentzVector ; +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 + Double_t fNevt ; // event counter + + 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); + + const Double_t CostCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, Double_t charge1, + Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t energy); + const Double_t CostHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, Double_t charge1, + Double_t px2, Double_t py2, Double_t pz2, Double_t e2); + const Double_t PhiCS(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, Double_t charge1, + Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t energy); + const Double_t PhiHE(Double_t px1, Double_t py1, Double_t pz1, Double_t e1, Double_t charge1, + Double_t px2, Double_t py2, Double_t pz2, Double_t e2, Double_t energy); + + ClassDef(AliCFMuonResTask1,1); +}; + +#endif