From: dainese Date: Tue, 1 Feb 2011 10:14:36 +0000 (+0000) Subject: Remove obsolte Jpsi classes X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=37d104631f847f0059f60f97f73a16f3380ca6bf;p=u%2Fmrichter%2FAliRoot.git Remove obsolte Jpsi classes --- diff --git a/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.cxx b/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.cxx deleted file mode 100644 index 4814efa964f..00000000000 --- a/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.cxx +++ /dev/null @@ -1,123 +0,0 @@ -/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * - * See cxx source for full Copyright notice */ - -//------------------------------------------------------------------------- -// Class AliAnalysisBtoJPSItoEle -// Unbinned log-likelihood fit analysis class -// -// Origin: C.Di Giglio -// Contact: Carmelo.Digiglio@ba.infn.it , giuseppe.bruno@ba.infn.it -//------------------------------------------------------------------------- -class TH1F ; -#include "TNtuple.h" -#include "TMath.h" - -#include "AliBtoJPSItoEleCDFfitFCN.h" -#include "AliBtoJPSItoEleCDFfitHandler.h" -#include "AliAnalysisBtoJPSItoEle.h" -#include "AliLog.h" - -ClassImp(AliAnalysisBtoJPSItoEle) - -//_______________________________________________________________________________ -AliAnalysisBtoJPSItoEle::AliAnalysisBtoJPSItoEle() : -fFCNfunction(0), -fPtBin(0), -fMCtemplate(0) -{ - // - // default constructor - // -} -//___________________________________________________________________________________ -AliAnalysisBtoJPSItoEle::AliAnalysisBtoJPSItoEle(const AliAnalysisBtoJPSItoEle& source) : -TNamed(source), -fFCNfunction(source.fFCNfunction), -fPtBin(source.fPtBin), -fMCtemplate(source.fMCtemplate) -{ - // - // copy constructor - // -} -//_________________________________________________________________________________________________ - -AliAnalysisBtoJPSItoEle &AliAnalysisBtoJPSItoEle::operator=(const AliAnalysisBtoJPSItoEle& source) -{ - // - // assignment operator - // - if(&source == this) return *this; - fFCNfunction = source.fFCNfunction; - fPtBin = source.fPtBin; - fMCtemplate = source.fMCtemplate; - - return *this; -} -//_________________________________________________________________________________________________ -AliAnalysisBtoJPSItoEle::~AliAnalysisBtoJPSItoEle() -{ - // - // destructor - // - delete fFCNfunction; - delete fMCtemplate; -} -//_________________________________________________________________________________________________ -Int_t AliAnalysisBtoJPSItoEle::DoMinimization() -{ - // - // performs the minimization - // - Int_t iret=fFCNfunction->DoMinimization(); - - return iret; -} -//_________________________________________________________________________________________________ -void AliAnalysisBtoJPSItoEle::ReadCandidates(TNtuple* nt, Double_t* &pseudoproper,Double_t* &invmass, Int_t& ncand) -{ - // - // Read N-tuple with X and M values - // - Float_t mJPSI = 0; Float_t x = 0; - Int_t nentries = 0; - ncand=0; - nt->SetBranchAddress("Mass",&mJPSI); - nt->SetBranchAddress("Xdecaytime",&x); - nentries = (Int_t)nt->GetEntries(); - pseudoproper = new Double_t[nentries]; - invmass = new Double_t[nentries]; - for(Int_t i = 0; i < nentries; i++) { - nt->GetEntry(i); - ncand++; - pseudoproper[i]=(Double_t)(10000*x); - invmass[i]=(Double_t)mJPSI; - } - - return; -} -//_________________________________________________________________________________________________ -void AliAnalysisBtoJPSItoEle::SetCsiMC() -{ - // - // Sets X distribution used as MC template for JPSI from B - // - fFCNfunction->LikelihoodPointer()->SetCsiMC(fMCtemplate); - - return; -} -//_________________________________________________________________________________________________ -void AliAnalysisBtoJPSItoEle::SetFitHandler(Double_t* x /*pseudoproper*/, Double_t* m /*inv mass*/, Int_t ncand /*candidates*/) -{ - // - // Create the fit handler object to play with different params of the fitting function - // - - fFCNfunction = new AliBtoJPSItoEleCDFfitHandler(x,m,ncand); - if(!fFCNfunction) { - - AliInfo("fFCNfunction not istanziated ---> nothing done"); - return; - - } -} diff --git a/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.h b/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.h deleted file mode 100644 index f432564c84b..00000000000 --- a/PWG3/vertexingHF/AliAnalysisBtoJPSItoEle.h +++ /dev/null @@ -1,50 +0,0 @@ -#ifndef ALIANALYSISBTOJPSITOELE_H -#define ALIANALYSISBTOJPSITOELE_H -/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * - * See cxx source for full Copyright notice */ - -//------------------------------------------------------------------------- -// Class AliAnalysisBtoJPSItoEle -// Unbinned log-likelihood fit analysis class -// -// Origin: C.Di Giglio -// Contact: Carmelo.Digiglio@ba.infn.it , giuseppe.bruno@ba.infn.it -//------------------------------------------------------------------------- -class TNtuple ; -//class AliBtoJPSItoEleCDFfitFCN ; -//class AliBtoJPSItoEleCDFfitHandler ; -#include "TH1F.h" -#include "AliBtoJPSItoEleCDFfitHandler.h" -#include "AliBtoJPSItoEleCDFfitFCN.h" - -//_________________________________________________________________________ -class AliAnalysisBtoJPSItoEle : public TNamed { - public: - // - AliAnalysisBtoJPSItoEle(); - AliAnalysisBtoJPSItoEle(const AliAnalysisBtoJPSItoEle& source); - AliAnalysisBtoJPSItoEle& operator=(const AliAnalysisBtoJPSItoEle& source); - virtual ~AliAnalysisBtoJPSItoEle(); - - Int_t DoMinimization(); - void ReadCandidates(TNtuple* nt, Double_t* &x, Double_t* &m, Int_t& n); // primary JPSI + secondary JPSI + bkg couples - - void SetPtBin(Int_t BinNum) { fPtBin = BinNum ; } - void SetCsiMC(); - //void SetResolutionConstants(Int_t BinNum); - void SetFitHandler(Double_t* /*pseudoproper*/, Double_t* /*inv mass*/, Int_t /*candidates*/); - void CloneMCtemplate(const TH1F* MCtemplate) {fMCtemplate = (TH1F*)MCtemplate->Clone("fMCtemplate");} - Double_t* GetResolutionConstants(Double_t* resolutionConst); - AliBtoJPSItoEleCDFfitHandler* GetCDFFitHandler() const { return fFCNfunction ; } - Int_t GetPtBin() const { return fPtBin ; } - - private: - // - AliBtoJPSItoEleCDFfitHandler* fFCNfunction; //! pointer to the interface class - Int_t fPtBin; // number of pt bin in which the analysis is performes - TH1F* fMCtemplate; //! template of the MC distribution for the x distribution of the secondary J/psi - - ClassDef(AliAnalysisBtoJPSItoEle,1); // AliAnalysisBtoJPSItoEle class -}; - -#endif diff --git a/PWG3/vertexingHF/AliAnalysisTaskSEBkgLikeSignJPSI.cxx b/PWG3/vertexingHF/AliAnalysisTaskSEBkgLikeSignJPSI.cxx deleted file mode 100644 index 9e5dccb1729..00000000000 --- a/PWG3/vertexingHF/AliAnalysisTaskSEBkgLikeSignJPSI.cxx +++ /dev/null @@ -1,372 +0,0 @@ -/************************************************************************** - * Copyright(c) 1998-2009, 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. * - **************************************************************************/ - -/////////////////////////////////////////////////////////////////////////// -// -// AliAnalysisTaskSE for reading both reconstructed JPSI -> ee candidates -// and like sign pairs and for drawing corresponding distributions -// -// Author: C.Di Giglio, carmelo.digiglio@ba.infn.it -/////////////////////////////////////////////////////////////////////////// - -#include -#include -#include -#include -#include -#include - -#include "AliAnalysisManager.h" -#include "AliAODHandler.h" -#include "AliAODEvent.h" -#include "AliAODVertex.h" -#include "AliAODTrack.h" -#include "AliAODRecoDecayHF2Prong.h" -#include "AliAnalysisVertexingHF.h" -#include "AliAnalysisTaskSE.h" -#include "AliAnalysisTaskSEBkgLikeSignJPSI.h" - -ClassImp(AliAnalysisTaskSEBkgLikeSignJPSI) - -//________________________________________________________________________ -AliAnalysisTaskSEBkgLikeSignJPSI::AliAnalysisTaskSEBkgLikeSignJPSI(): -AliAnalysisTaskSE(), -fOutput(0), -fHistMassJPSI(0), -fHistMassLS(0), -fHistCtsJPSI(0), -fHistCtsLS(0), -fHistCtsLSpos(0), -fHistCtsLSneg(0), -fHistCPtaJPSI(0), -fHistCPtaLS(0), -fHistd0d0JPSI(0), -fHistd0d0LS(0), -fHistDCAJPSI(0), -fHistDCALS(0), -fVHF(0), -fTotPosPairs(0), -fTotNegPairs(0), -fLsNormalization(1.) -{ - // - // Default constructor - // -} - -//________________________________________________________________________ -AliAnalysisTaskSEBkgLikeSignJPSI::AliAnalysisTaskSEBkgLikeSignJPSI(const char *name): -AliAnalysisTaskSE(name), -fOutput(0), -fHistMassJPSI(0), -fHistMassLS(0), -fHistCtsJPSI(0), -fHistCtsLS(0), -fHistCtsLSpos(0), -fHistCtsLSneg(0), -fHistCPtaJPSI(0), -fHistCPtaLS(0), -fHistd0d0JPSI(0), -fHistd0d0LS(0), -fHistDCAJPSI(0), -fHistDCALS(0), -fVHF(0), -fTotPosPairs(0), -fTotNegPairs(0), -fLsNormalization(1.) -{ - // - // Standard constructor - // - // Output slot #1 writes into a TList container - DefineOutput(1,TList::Class()); //My private output -} - -//________________________________________________________________________ -AliAnalysisTaskSEBkgLikeSignJPSI::~AliAnalysisTaskSEBkgLikeSignJPSI() -{ - // Destructor - if (fOutput) { - delete fOutput; - fOutput = 0; - } - - if (fVHF) { - delete fVHF; - fVHF = 0; - } - -} -//________________________________________________________________________ -void AliAnalysisTaskSEBkgLikeSignJPSI::Init() -{ - // Initialization - - if(fDebug > 1) printf("AnalysisTaskSEBkgLikeSignJPSI::Init() \n"); - - gROOT->LoadMacro("ConfigVertexingHF.C"); - - fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()"); - fVHF->PrintStatus(); - - return; -} - -//________________________________________________________________________ -void AliAnalysisTaskSEBkgLikeSignJPSI::UserCreateOutputObjects() -{ - // Create the output container - // - if(fDebug > 1) printf("AnalysisTaskSEBkgLikeSignJPSI::UserCreateOutputObjects() \n"); - - // Several histograms are more conveniently managed in a TList - fOutput = new TList(); - fOutput->SetOwner(); - - fHistMassJPSI = new TH1F("fHistMassJPSI", "J/#Psi invariant mass; M [GeV]; Entries",200,2.8,3.25); - fHistMassJPSI->Sumw2(); - fHistMassJPSI->SetMinimum(0); - fOutput->Add(fHistMassJPSI); - - fHistMassLS = new TH1F("fHistMassLS", "Like sign pairs invariant mass; M [GeV]; Entries",200,2.8,3.25); - fHistMassLS->Sumw2(); - fHistMassLS->SetMinimum(0); - fOutput->Add(fHistMassLS); - - fHistCtsJPSI = new TH1F("fHistCtsJPSI", "J/#Psi cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.); - fHistCtsJPSI->Sumw2(); - fHistCtsJPSI->SetMinimum(0); - fOutput->Add(fHistCtsJPSI); - - fHistCtsLS = new TH1F("fHistCtsLS", "Like sign pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.); - fHistCtsLS->Sumw2(); - fHistCtsLS->SetMinimum(0); - fOutput->Add(fHistCtsLS); - - fHistCtsLSpos = new TH1F("fHistCtsLSpos", "Like sign ++ pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.); - fHistCtsLSpos->Sumw2(); - fHistCtsLSpos->SetMinimum(0); - fOutput->Add(fHistCtsLSpos); - - fHistCtsLSneg = new TH1F("fHistCtsLSneg", "Like sign -- pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.); - fHistCtsLSneg->Sumw2(); - fHistCtsLSneg->SetMinimum(0); - fOutput->Add(fHistCtsLSneg); - - fHistCPtaJPSI = new TH1F("fHistCPtaJPSI", "J/#Psi cosine of pointing angle; Cos#Theta_{point}; Entries",200,-1.,1.); - fHistCPtaJPSI->Sumw2(); - fHistCPtaJPSI->SetMinimum(0); - fOutput->Add(fHistCPtaJPSI); - - fHistCPtaLS = new TH1F("fHistCPtaLS", "Like sign pairs cosine of pointing angle; Cos#Theta_{point}; Entries",200,-1.,1.); - fHistCPtaLS->Sumw2(); - fHistCPtaLS->SetMinimum(0); - fOutput->Add(fHistCPtaLS); - - fHistd0d0JPSI = new TH1F("fHistd0d0JPSI", "J/#Psi product of impact parameters; d0xd0 [#mu m^{2}]; Entries",200,-100000.,100000.); - fHistd0d0JPSI->Sumw2(); - fHistd0d0JPSI->SetMinimum(0); - fOutput->Add(fHistd0d0JPSI); - - fHistd0d0LS = new TH1F("fHistd0d0LS", "Like sign pairs product of impact parameters; d0xd0 [#mu m^{2}]; Entries",200,-100000.,100000.); - fHistd0d0LS->Sumw2(); - fHistd0d0LS->SetMinimum(0); - fOutput->Add(fHistd0d0LS); - - fHistDCAJPSI = new TH1F("fHistDCAJPSI", "J/#Psi distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.); - fHistDCAJPSI->Sumw2(); - fHistDCAJPSI->SetMinimum(0); - fOutput->Add(fHistDCAJPSI); - - fHistDCALS = new TH1F("fHistDCALS", "Like sign pairs distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.); - fHistDCALS->Sumw2(); - fHistDCALS->SetMinimum(0); - fOutput->Add(fHistDCALS); - - return; -} - -//________________________________________________________________________ -void AliAnalysisTaskSEBkgLikeSignJPSI::UserExec(Option_t */*option*/) -{ - // Execute analysis for current event: - // heavy flavor candidates association to MC truth - - AliAODEvent *aod = dynamic_cast (InputEvent()); - - TClonesArray *arrayJPSItoEle = 0; - TClonesArray *arrayLikeSign = 0; - - if(!aod && AODEvent() && IsStandardAOD()) { - // In case there is an AOD handler writing a standard AOD, use the AOD - // event in memory rather than the input (ESD) event. - aod = dynamic_cast (AODEvent()); - // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root) - // have to taken from the AOD event hold by the AliAODExtension - AliAODHandler* aodHandler = (AliAODHandler*) - ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()); - if(aodHandler->GetExtensions()) { - AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root"); - AliAODEvent *aodFromExt = ext->GetAOD(); - // load Jpsi candidates - arrayJPSItoEle=(TClonesArray*)aodFromExt->GetList()->FindObject("JPSItoEle"); - // load like sign candidates - arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign2Prong"); - } - } else { - // load Jpsi candidates - arrayJPSItoEle=(TClonesArray*)aod->GetList()->FindObject("JPSItoEle"); - // load like sign candidates - arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign2Prong"); - } - - - if(!arrayJPSItoEle) { - printf("AliAnalysisTaskSEBkgLikeSignJPSI::UserExec: JPSItoEle branch not found!\n"); - return; - } - if(!arrayLikeSign) { - printf("AliAnalysisTaskSEBkgLikeSignJPSI::UserExec: LikeSign2Prong branch not found!\n"); - return; - } - - // fix for temporary bug in ESDfilter - // the AODs with null vertex pointer didn't pass the PhysSel - if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return; - - // AOD primary vertex - AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex(); - - // make trkIDtoEntry register (temporary) - Int_t trkIDtoEntry[100000]; - for(Int_t it=0;itGetNumberOfTracks();it++) { - AliAODTrack *track = aod->GetTrack(it); - trkIDtoEntry[track->GetID()]=it; - } - - // loop over Like sign candidates - Int_t nPosPairs=0,nNegPairs=0; - Int_t nLikeSign = arrayLikeSign->GetEntriesFast(); - if(fDebug>1) printf("+++\n+++Number of like sign pairs ---> %d \n+++\n", nLikeSign); - - for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) { - AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayLikeSign->UncheckedAt(iLikeSign); - Bool_t unsetvtx=kFALSE; - if(!d->GetOwnPrimaryVtx()) { - d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables - unsetvtx=kTRUE; - } - //Int_t okBtoJPSIls=0; - //if(d->SelectBtoJPSI(fVHF->GetBtoJPSICuts(),okBtoJPSIls)) { - if(d) { - AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0); - fHistMassLS->Fill(d->InvMassJPSIee()); - fHistCPtaLS->Fill(d->CosPointingAngle()); - fHistd0d0LS->Fill(1e8*d->Prodd0d0()); - fHistDCALS->Fill(100*d->GetDCA()); - //PostData(1,fOutput); - if(!trk0) { - trk0=aod->GetTrack(trkIDtoEntry[d->GetProngID(0)]); - printf("references to standard AOD not available \n"); - } - if((trk0->Charge())==1) { - nPosPairs++; - fHistCtsLS->Fill(d->CosThetaStar(0,443,11,11)); - fHistCtsLSpos->Fill(d->CosThetaStar(0,443,11,11)); - //PostData(1,fOutput); - } else { - nNegPairs++; - fHistCtsLS->Fill(d->CosThetaStarJPSI()); - fHistCtsLSneg->Fill(d->CosThetaStarJPSI()); - //PostData(1,fOutput); - } - PostData(1,fOutput); - } - if(unsetvtx) d->UnsetOwnPrimaryVtx(); - } - - if(fDebug>1) printf("------------ N. of positive pairs in Event ----- %d \n", nPosPairs); - if(fDebug>1) printf("------------ N. of negative pairs in Event ----- %d \n", nNegPairs); - - fTotPosPairs += nPosPairs; - fTotNegPairs += nNegPairs; - - // loop over JPSI candidates - Int_t nBtoJpsiToEle = arrayJPSItoEle->GetEntriesFast(); - if(fDebug>1) printf("Number of like JPSI -> ee candidates ---> %d \n", nBtoJpsiToEle); - - for (Int_t iBtoJpsiToEle = 0; iBtoJpsiToEle < nBtoJpsiToEle; iBtoJpsiToEle++) { - AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayJPSItoEle->UncheckedAt(iBtoJpsiToEle); - Bool_t unsetvtx=kFALSE; - if(!d->GetOwnPrimaryVtx()) { - d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables - unsetvtx=kTRUE; - } - ///Int_t okBtoJPSI=0; - //if(d->SelectBtoJPSI(fVHF->GetBtoJPSICuts(),okBtoJPSI)) { - if(d) { - fHistMassJPSI->Fill(d->InvMassJPSIee()); - fHistCtsJPSI->Fill(d->CosThetaStarJPSI()); - fHistd0d0JPSI->Fill(1e8*d->Prodd0d0()); - fHistCPtaJPSI->Fill(d->CosPointingAngle()); - fHistDCAJPSI->Fill(100*d->GetDCA()); - PostData(1,fOutput); - } - if(unsetvtx) d->UnsetOwnPrimaryVtx(); - } - - return; -} - -//________________________________________________________________________ -void AliAnalysisTaskSEBkgLikeSignJPSI::Terminate(Option_t */*option*/) -{ - // Terminate analysis - // - if(fDebug > 1) printf("AnalysisTaskSEBkgLikeSignJPSI: Terminate() \n"); - - fOutput = dynamic_cast (GetOutputData(1)); - if (!fOutput) { - printf("ERROR: fOutput not available\n"); - return; - } - - fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs*fTotNegPairs); - - fHistMassJPSI = dynamic_cast(fOutput->FindObject("fHistMassJPSI")); - fHistMassLS = dynamic_cast(fOutput->FindObject("fHistMassLS")); - fHistCtsJPSI = dynamic_cast(fOutput->FindObject("fHistCtsJPSI")); - fHistCtsLS = dynamic_cast(fOutput->FindObject("fHistCtsLS")); - fHistCtsLSpos = dynamic_cast(fOutput->FindObject("fHistCtsLSpos")); - fHistCtsLSneg = dynamic_cast(fOutput->FindObject("fHistCtsLSneg")); - fHistCPtaJPSI = dynamic_cast(fOutput->FindObject("fHistCPtaJPSI")); - fHistCPtaLS = dynamic_cast(fOutput->FindObject("fHistCPtaLS")); - fHistd0d0JPSI = dynamic_cast(fOutput->FindObject("fHistd0d0JPSI")); - fHistd0d0LS = dynamic_cast(fOutput->FindObject("fHistd0d0LS")); - fHistDCAJPSI = dynamic_cast(fOutput->FindObject("fHistDCAJPSI")); - fHistDCALS = dynamic_cast(fOutput->FindObject("fHistDCALS")); - - if(fLsNormalization>0.) { - fHistMassLS->Scale((1/fLsNormalization)*fHistMassLS->GetEntries()); - fHistCtsLS->Scale((1/fLsNormalization)*fHistCtsLS->GetEntries()); - fHistCtsLSpos->Scale((1/fLsNormalization)*fHistCtsLSpos->GetEntries()); - fHistCtsLSneg->Scale((1/fLsNormalization)*fHistCtsLSneg->GetEntries()); - fHistCPtaLS->Scale((1/fLsNormalization)*fHistCPtaLS->GetEntries()); - fHistd0d0LS->Scale((1/fLsNormalization)*fHistd0d0LS->GetEntries()); - fHistDCALS->Scale((1/fLsNormalization)*fHistDCALS->GetEntries()); - } - - return; -} diff --git a/PWG3/vertexingHF/AliAnalysisTaskSEBkgLikeSignJPSI.h b/PWG3/vertexingHF/AliAnalysisTaskSEBkgLikeSignJPSI.h deleted file mode 100644 index c893c2b4c62..00000000000 --- a/PWG3/vertexingHF/AliAnalysisTaskSEBkgLikeSignJPSI.h +++ /dev/null @@ -1,67 +0,0 @@ -#ifndef AliAnalysisTaskSEBkgLikeSignJPSI_H -#define AliAnalysisTaskSEBkgLikeSignJPSI_H - -/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * - * See cxx source for full Copyright notice */ - -//************************************************************************* -// Class AliAnalysisTaskSEBkgLikeSignJPSI -// AliAnalysisTaskSE for reading both reconstructed JPSI -> ee candidates -// and like sign pairs and for drawing corresponding distributions -// Author: C.Di Giglio, carmelo.digiglio@ba.infn.it -//************************************************************************* - -#include -#include -#include -#include -#include - -#include "AliAnalysisTaskSE.h" -#include "AliAnalysisVertexingHF.h" - -class AliAnalysisTaskSEBkgLikeSignJPSI : public AliAnalysisTaskSE -{ - public: - - AliAnalysisTaskSEBkgLikeSignJPSI(); - AliAnalysisTaskSEBkgLikeSignJPSI(const char *name); - virtual ~AliAnalysisTaskSEBkgLikeSignJPSI(); - - - // Implementation of interface methods - virtual void UserCreateOutputObjects(); - virtual void Init(); - virtual void LocalInit() {Init();} - virtual void UserExec(Option_t *option); - virtual void Terminate(Option_t *option); - - private: - - AliAnalysisTaskSEBkgLikeSignJPSI(const AliAnalysisTaskSEBkgLikeSignJPSI &source); - AliAnalysisTaskSEBkgLikeSignJPSI& operator=(const AliAnalysisTaskSEBkgLikeSignJPSI& source); - - TList *fOutput; //! list send on output slot 0 - TH1F *fHistMassJPSI; //! output histograms - TH1F *fHistMassLS; //! - TH1F *fHistCtsJPSI; //! Cosine of decay angle - TH1F *fHistCtsLS; //! - TH1F *fHistCtsLSpos; //! - TH1F *fHistCtsLSneg; //! - TH1F *fHistCPtaJPSI; //! Cosine of pointing angle - TH1F *fHistCPtaLS; //! - TH1F *fHistd0d0JPSI; //! Product of impact parameters - TH1F *fHistd0d0LS; //! - TH1F *fHistDCAJPSI; //! Distance of closest approach - TH1F *fHistDCALS; //! like-sign - AliAnalysisVertexingHF *fVHF; // Vertexer heavy flavour (used to pass the cuts) - - Int_t fTotPosPairs; // - Int_t fTotNegPairs; // normalization - Double_t fLsNormalization; // - - ClassDef(AliAnalysisTaskSEBkgLikeSignJPSI,1); // comparison of unlike-sign and like-sign background for J/psi->ee -}; - -#endif - diff --git a/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.cxx b/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.cxx deleted file mode 100644 index 5cd8ce5ec2e..00000000000 --- a/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.cxx +++ /dev/null @@ -1,539 +0,0 @@ -/************************************************************************** - * Copyright(c) 1998-2009, 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. * - **************************************************************************/ -//************************************************************************* -// Class AliAnalysisTaskSEBtoJPSItoEle -// AliAnalysisTaskSE for the reconstruction -// of heavy-flavour decay candidates -// Author: C.Di Giglio, carmelo.digiglio@ba.infn.it -//************************************************************************* -class AliAnalysisTaskSE; -class AliESDEvent; -class AliVEvent; -class iostream; -class TSystem; -class TROOT; -#include -#include -#include - -#include "AliAnalysisManager.h" -#include "AliAODHandler.h" -#include "AliAODEvent.h" -#include "AliAODMCParticle.h" -#include "AliAODMCHeader.h" -#include "AliAODRecoDecayHF2Prong.h" - -#include "AliAnalysisBtoJPSItoEle.h" -#include "AliAnalysisTaskSEBtoJPSItoEle.h" - -ClassImp(AliAnalysisTaskSEBtoJPSItoEle) - -//______________________________________________________________________________ -AliAnalysisTaskSEBtoJPSItoEle::AliAnalysisTaskSEBtoJPSItoEle() : -AliAnalysisTaskSE(), -fOutput(0), -fCdfFit(0), -fNtupleJPSI(0), -fhDecayTimeMCjpsifromB(0), -fhDecayTime(0), -fhInvMass(0), -fhD0(0), -fhD0D0(0), -fhCosThetaStar(0), -fhCosThetaPointing(0), -fhCtsVsD0D0(0), -fOkAODMC(kFALSE), -fNameMCfile(""), -fOkMinimize(kFALSE) -{ - // default constructor -} -//_________________________________________________________________________________________________ -AliAnalysisTaskSEBtoJPSItoEle::AliAnalysisTaskSEBtoJPSItoEle(const char* name) : -AliAnalysisTaskSE(name), -fOutput(0), -fCdfFit(0), -fNtupleJPSI(0), -fhDecayTimeMCjpsifromB(0), -fhDecayTime(0), -fhInvMass(0), -fhD0(0), -fhD0D0(0), -fhCosThetaStar(0), -fhCosThetaPointing(0), -fhCtsVsD0D0(0), -fOkAODMC(kFALSE), -fNameMCfile(""), -fOkMinimize(kFALSE) -{ - // standard constructor - - // Output slot #1 writes into a TList container - DefineOutput(1,TList::Class()); //My private output -} -//_________________________________________________________________________________________________ -AliAnalysisTaskSEBtoJPSItoEle::~AliAnalysisTaskSEBtoJPSItoEle() -{ - // destructor - - if (fOutput) { - delete fOutput; - fOutput = 0; - } -} -//_________________________________________________________________________________________________ -void AliAnalysisTaskSEBtoJPSItoEle::Init() -{ - // Initialization - - if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle::Init() \n"); - - return; - -} -//_________________________________________________________________________________________________ -void AliAnalysisTaskSEBtoJPSItoEle::UserCreateOutputObjects() -{ - // Create the output container - - if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle::UserCreateOutputObjects() \n"); - - // Several histograms are more conveniently managed in a TList - fOutput = new TList(); - fOutput->SetOwner(); - - if(fOkAODMC){ - fhDecayTimeMCjpsifromB = new TH1F("fhDecayTimeMCjpsifromB", "Secondary J/#Psi Monte Carlo pseudo proper decay time; X [#mu m]; Entries",100,-4000.,4000.); - fhDecayTimeMCjpsifromB->Sumw2(); - fhDecayTimeMCjpsifromB->SetMinimum(0); - fOutput->Add(fhDecayTimeMCjpsifromB); - } - - fhDecayTime = new TH1F("fhDecayTime", "J/#Psi pseudo proper decay time; X [#mu m]; Entries",200,-2000.,2000.); - fhDecayTime->Sumw2(); - fhDecayTime->SetMinimum(0); - fOutput->Add(fhDecayTime); - - fhInvMass = new TH1F("fhInvMass", "J/#Psi invariant mass; M [GeV]; Entries",100,0.,3.2); - fhInvMass->Sumw2(); - fhInvMass->SetMinimum(0); - fOutput->Add(fhInvMass); - - fhD0 = new TH1F("fhD0", "Impact parameter; D0 [#mu m]; Entries",100,-5000.,5000.); - fhD0->Sumw2(); - fhD0->SetMinimum(0); - fOutput->Add(fhD0); - - fhD0D0 = new TH1F("fhD0D0", "Product of impact parameters; D0D0 [#mu m^2]; Entries",100,-100000.,100000.); - fhD0D0->Sumw2(); - fhD0D0->SetMinimum(0); - fOutput->Add(fhD0D0); - - fhCosThetaStar = new TH1F("fhCosThetaStar", "Cosine of decay angle; Cosine Theta Star; Entries",50,-1.2,1.2); - fhCosThetaStar->Sumw2(); - fhCosThetaStar->SetMinimum(0); - fOutput->Add(fhCosThetaStar); - - fhCosThetaPointing = new TH1F("fhCosThetaPointing", "Cosine of pointing angle; Cosine Pointing Angle; Entries",100,-1,1); - fhCosThetaPointing->Sumw2(); - fhCosThetaPointing->SetMinimum(0); - fOutput->Add(fhCosThetaPointing); - - fhCtsVsD0D0 = new TH2F("fhCtsVsD0D0", "Cosine theta star Vs. D0; Cosine theta star; D0D0",50,-1,1,100,-100000,100000); - fhCtsVsD0D0->Sumw2(); - fhCtsVsD0D0->SetMinimum(0); - fOutput->Add(fhCtsVsD0D0); - - fNtupleJPSI = new TNtuple("fNtupleJPSI","J/#Psi pseudo-proper decay time & invariant mass","Xdecaytime:Mass"); - fOutput->Add(fNtupleJPSI); - - return; -} -//_________________________________________________________________________________________________ -void AliAnalysisTaskSEBtoJPSItoEle::UserExec(Option_t */*option*/) -{ - // Execute analysis for current event: - - AliAODEvent *aod = dynamic_cast (InputEvent()); - - TClonesArray *inputArrayJPSItoEle = 0; - - if(!aod && AODEvent() && IsStandardAOD()) { - // In case there is an AOD handler writing a standard AOD, use the AOD - // event in memory rather than the input (ESD) event. - aod = dynamic_cast (AODEvent()); - // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root) - // have to taken from the AOD event hold by the AliAODExtension - AliAODHandler* aodHandler = (AliAODHandler*) - ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()); - if(aodHandler->GetExtensions()) { - AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root"); - AliAODEvent *aodFromExt = ext->GetAOD(); - // load Jpsi candidates - inputArrayJPSItoEle=(TClonesArray*)aodFromExt->GetList()->FindObject("JPSItoEle"); - } - } else { - // load Jpsi candidates - inputArrayJPSItoEle=(TClonesArray*)aod->GetList()->FindObject("JPSItoEle"); - } - - if(!inputArrayJPSItoEle) { - printf("AliAnalysisTaskSECompareHF::UserExec: JPSItoEle branch not found!\n"); - return; - } - - // fix for temporary bug in ESDfilter - // the AODs with null vertex pointer didn't pass the PhysSel - if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return; - - - // AOD primary vertex - AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex(); - - // load MC particles - TClonesArray* mcArray = - dynamic_cast(aod->FindListObject(AliAODMCParticle::StdBranchName())); - if (!mcArray) AliError("Could not find Monte-Carlo in AOD"); - AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast())); - - // Read AOD Monte Carlo information - if(fOkAODMC) ReadAODMCInfo(aod,inputArrayJPSItoEle); - - // loop over J/Psi->ee candidates - Int_t nInJPSItoEle = inputArrayJPSItoEle->GetEntriesFast(); - printf("Number of JPSI->e+e- candidates: %d\n",nInJPSItoEle); - - for (Int_t iJPSItoEle = 0; iJPSItoEle < nInJPSItoEle; iJPSItoEle++) { - AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayJPSItoEle->UncheckedAt(iJPSItoEle); - - Int_t mcLabel = d->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI - - //Secondary vertex - Double_t vtxSec[3] = {0.,0.,0.}; - Double_t vtxPrim[3] = {0.,0.,0.}; - d->GetSecondaryVtx(vtxSec); - Bool_t unsetvtx=kFALSE; - if(!d->GetOwnPrimaryVtx()) { - d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables - unsetvtx=kTRUE; - } - //Int_t okJPSI=0; - // here analyze only if cuts are passed - - if(d->Pt() < fPtCuts[1] && d->Pt() > fPtCuts[0]){ // apply pt cut only - if (mcLabel == -1) - { - AliDebug(2,"No MC particle found"); - } - else { - - fhInvMass->Fill(d->InvMassJPSIee()); - fhD0->Fill(10000*d->ImpParXY()); - fhD0D0->Fill(1e8*d->Prodd0d0()); - fhCosThetaStar->Fill(d->CosThetaStarJPSI()); - fhCtsVsD0D0->Fill(d->CosThetaStarJPSI(),1e8*d->Prodd0d0()); - fhCosThetaPointing->Fill(d->CosPointingAngle()); - - // compute X variable - AliAODVertex* primVertex = d->GetOwnPrimaryVtx(); - vtxPrim[0] = primVertex->GetX(); - vtxPrim[1] = primVertex->GetY(); - vtxPrim[2] = primVertex->GetZ(); - Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(d->Px()) + (vtxSec[1]-vtxPrim[1])*(d->Py()))/d->Pt(); - Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/d->Pt(); - - fhDecayTime->Fill(10000*pseudoProperDecayTime); - // Post the data already here - PostData(1,fOutput); - - fNtupleJPSI->Fill(pseudoProperDecayTime,d->InvMassJPSIee()); // fill N-tuple with X,M values - - }// - - } // end of JPSItoEle candidates selection according to cuts - - if(unsetvtx) d->UnsetOwnPrimaryVtx(); - - }// end loop on JPSI to ele candidates - -} -//_________________________________________________________________________________________________ -void AliAnalysisTaskSEBtoJPSItoEle::Terminate(Option_t */*option*/) -{ - // - // Terminate analysis - // - - if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle: Terminate() \n"); - - fOutput = dynamic_cast (GetOutputData(1)); - if (!fOutput) { - printf("ERROR: fOutput not available\n"); - return; - } - - if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast(fOutput->FindObject("fhDecayTimeMCjpsifromB")); - fhDecayTime = dynamic_cast(fOutput->FindObject("fhDecayTime")); - fhInvMass = dynamic_cast(fOutput->FindObject("fhInvMass")); - fhD0 = dynamic_cast(fOutput->FindObject("fhD0")); - fhD0D0 = dynamic_cast(fOutput->FindObject("fhD0D0")); - fhCosThetaStar = dynamic_cast(fOutput->FindObject("fhCosThetaStar")); - fhCosThetaPointing = dynamic_cast(fOutput->FindObject("fhCosThetaPointing")); - fhCtsVsD0D0 = dynamic_cast(fOutput->FindObject("fhCtsVsD0D0")); - fNtupleJPSI = dynamic_cast(fOutput->FindObject("fNtupleJPSI")); - - if(fOkMinimize){ - Double_t* pseudoProperValues=0x0; - Double_t* invMassValues=0x0; - Int_t ncand=0; - fCdfFit = new AliAnalysisBtoJPSItoEle(); - - printf("+++\n+++ AliAnalysisBtoJPSItoEle object instantiated ---> OK\n+++\n"); - fCdfFit->ReadCandidates(fNtupleJPSI,pseudoProperValues,invMassValues,ncand); - printf("+++\n+++ X and M vectors filled starting from N-tuple ---> OK\n+++\n"); - printf("+++\n+++ Number of candidates ---> %d J/psi \n+++\n",ncand); - fCdfFit->SetPtBin(0); - printf("+++\n+++ Pt bin setted n. ---> %d \n+++\n",fCdfFit->GetPtBin()); - - if(fOkAODMC) {fCdfFit->CloneMCtemplate(fhDecayTimeMCjpsifromB);} - else { - //SetPathMCfile("."); - TH1F* hLocal = OpenLocalMCFile(fNameMCfile,fCdfFit->GetPtBin()); - fCdfFit->CloneMCtemplate(hLocal); - } - - printf("+++\n+++ MC template histo copied ---> OK\n+++\n"); - //fCdfFit->DoMinimization(pseudoProperValues,invMassValues,ncand); - fCdfFit->DoMinimization(); - - if(pseudoProperValues) delete [] pseudoProperValues; - pseudoProperValues=NULL; - if(invMassValues) delete [] invMassValues; - invMassValues=NULL; - - if(fCdfFit) delete [] fCdfFit; - fCdfFit=NULL; - - } - - return; -} -//_________________________________________________________________________________________________ -void AliAnalysisTaskSEBtoJPSItoEle::SetPtCuts(const Double_t ptCuts[2]) -{ - // - // apply Pt cuts (lower cut, upper cut) - // - for(Int_t i = 0; i < 2; i++) fPtCuts[i] = ptCuts[i]; -} -//_________________________________________________________________________________________________ -void AliAnalysisTaskSEBtoJPSItoEle::SetCutsD0(const Double_t cuts[9]) -{ - // - // apply D0 and JPSI cuts - // - for(Int_t i = 0; i < 9; i++) fCuts[i] = cuts[i]; -} -//_________________________________________________________________________________________________ -void AliAnalysisTaskSEBtoJPSItoEle::ReadAODMCInfo(AliAODEvent* aodEvent, const TClonesArray* inputArray) -{ - // - // Reads MC information if needed (for example in case of parametrized PID) - // - - // load MC particles - TClonesArray* mcArray = - dynamic_cast(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); - if (!mcArray) AliError("Could not find Monte-Carlo in AOD"); - AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast())); - - // load MC header - AliAODMCHeader* mcHeader = - (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()); - if(!mcHeader){ - printf("AliAnalysisTaskSEBtoJPSItoEle::UserExec: MC AOD header branch not found!\n"); - } - - AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex(); - - //Double_t rmax = 0.000005; - Double_t mcPrimVtx[3]; - - // get MC primary Vtx - mcHeader->GetVertex(mcPrimVtx); - - Int_t nInCandidates = inputArray->GetEntriesFast(); - printf("Number of Candidates for MC analysis: %d\n",nInCandidates); - - Int_t lab0, lab1, pdgMother, labMother, pdgJPSI, pdg0, pdg1, labJPSIMother, pdgJPSIMother; - Int_t labJPSIdaugh0=0; - Int_t labJPSIdaugh1=0; - - for (Int_t iCandidate = 0; iCandidate < nInCandidates; iCandidate++) { - - AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iCandidate); - Int_t mcLabelSec = dd->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI - - Double_t vtxPrim[3] = {0.,0.,0.}; - Double_t vtxSec[3] = {0.,0.,0.}; - dd->GetSecondaryVtx(vtxSec); - Bool_t unsetvtx=kFALSE; - if(!dd->GetOwnPrimaryVtx()) { - dd->SetOwnPrimaryVtx(vtx1); // needed to compute all variables - unsetvtx=kTRUE; - } - - if(dd->Pt() < fPtCuts[1] && dd->Pt() > fPtCuts[0]){ // apply pt cut only - - AliAODTrack *trk0 = (AliAODTrack*)dd->GetDaughter(0); - AliAODTrack *trk1 = (AliAODTrack*)dd->GetDaughter(1); - lab0 = trk0->GetLabel(); - lab1 = trk1->GetLabel(); - - AliAODMCParticle* part0 = (AliAODMCParticle*)mcArray->At(lab0); - if(!part0) { - printf("no MC particle\n"); - continue; - } - - while(part0->GetMother()>=0) { - labMother=part0->GetMother(); - part0 = (AliAODMCParticle*)mcArray->At(labMother); - if(!part0) { - printf("no MC mother particle\n"); - break; - } - pdgMother = TMath::Abs(part0->GetPdgCode()); - if(pdgMother==443) {//this for JPSI - labJPSIdaugh0=labMother; - break; - } - } - - AliAODMCParticle* part1 = (AliAODMCParticle*)mcArray->At(lab1); - if(!part1) { - printf("no MC particle\n"); - continue; - } - - while(part1->GetMother()>=0) { - labMother=part1->GetMother(); - part1 = (AliAODMCParticle*)mcArray->At(labMother); - if(!part1) { - printf("no MC mother particle\n"); - break; - } - pdgMother = TMath::Abs(part1->GetPdgCode()); - if(pdgMother==443) {//this for JPSI - labJPSIdaugh1=labMother; - break; - } - } - - if (mcLabelSec == -1) - { - AliDebug(2,"No MC particle found"); - } - else { - - if(labJPSIdaugh0>=0 && labJPSIdaugh1>=0 && labJPSIdaugh0==labJPSIdaugh1) { // check if JPSI is signal - AliAODMCParticle* partJPSI = (AliAODMCParticle*)mcArray->At(labJPSIdaugh0); - pdgJPSI = partJPSI->GetPdgCode(); - Int_t pdaugh0 = partJPSI->GetDaughter(0); - Int_t pdaugh1 = partJPSI->GetDaughter(1); - AliAODMCParticle* partDaugh0 = (AliAODMCParticle*)mcArray->At(pdaugh0); - AliAODMCParticle* partDaugh1 = (AliAODMCParticle*)mcArray->At(pdaugh1); - pdg0 = partDaugh0->GetPdgCode(); - pdg1 = partDaugh1->GetPdgCode(); - if(TMath::Abs(pdg0) == 11 && TMath::Abs(pdg1) == 11){ // this is for MC JPSI -> ee - labJPSIMother = partJPSI->GetMother(); - AliAODMCParticle* partJPSIMother = (AliAODMCParticle*)mcArray->At(labJPSIMother); - pdgJPSIMother = partJPSIMother->GetPdgCode(); - - // chech if JPSI comes from B - if( pdgJPSIMother==511 || pdgJPSIMother==521 || - pdgJPSIMother==10511 || pdgJPSIMother==10521 || - pdgJPSIMother==513 || pdgJPSIMother==523 || - pdgJPSIMother==10513 || pdgJPSIMother==10523 || - pdgJPSIMother==20513 || pdgJPSIMother==20523 || - pdgJPSIMother==515 || pdgJPSIMother==525 || - pdgJPSIMother==531 || pdgJPSIMother==10531 || - pdgJPSIMother==533 || pdgJPSIMother==10533 || - pdgJPSIMother==20533 || pdgJPSIMother==535 || - pdgJPSIMother==541 || pdgJPSIMother==10541 || - pdgJPSIMother==543 || pdgJPSIMother==10543 || - pdgJPSIMother==20543 || pdgJPSIMother==545) - { //this is for MC JPSI -> ee from B-hadrons - - fhInvMass->Fill(dd->InvMassJPSIee()); - fhD0->Fill(10000*dd->ImpParXY()); - fhD0D0->Fill(1e8*dd->Prodd0d0()); - fhCosThetaStar->Fill(dd->CosThetaStarJPSI()); - fhCtsVsD0D0->Fill(dd->CosThetaStarJPSI(),1e8*dd->Prodd0d0()); - fhCosThetaPointing->Fill(dd->CosPointingAngle()); - - // compute X variable - AliAODVertex* primVertex = dd->GetOwnPrimaryVtx(); - vtxPrim[0] = primVertex->GetX(); - vtxPrim[1] = primVertex->GetY(); - vtxPrim[2] = primVertex->GetZ(); - Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dd->Px()) + (vtxSec[1]-vtxPrim[1])*(dd->Py()))/dd->Pt(); - Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dd->Pt(); - - fhDecayTime->Fill(10000*pseudoProperDecayTime); - // Post the data already here - PostData(1,fOutput); - - fNtupleJPSI->Fill(pseudoProperDecayTime,dd->InvMassJPSIee()); // fill N-tuple with X,M values - - Double_t mcLxy = ((partJPSI->Xv()-mcPrimVtx[0])*(partJPSI->Px()) + (partJPSI->Yv()-mcPrimVtx[1])*(partJPSI->Py()))/partJPSI->Pt(); - Double_t mcPseudoProperDecayTime = mcLxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/partJPSI->Pt(); - fhDecayTimeMCjpsifromB->Fill(10000*mcPseudoProperDecayTime,1); - - // Post the data already here - PostData(1,fOutput); - - } //this is for MC JPSI -> ee from B-hadrons - } //this is for MC JPSI -> ee - } - } // end of check if JPSI is signal - } - } - -} -//_________________________________________________________________________________________________ -TH1F* AliAnalysisTaskSEBtoJPSItoEle::OpenLocalMCFile(TString datadir, Int_t binNum) -{ - // - // Open a local file with MC x distribution for JPSI from B-hadron - // - - TH1F* hCsiMC = new TH1F(); - TString useFile = datadir.Data(); - useFile.Append("/CsiMCfromKine_10PtBins.root"); - TFile *f = new TFile(useFile); - Double_t scale = 0.; - char processCsiMCXHisto[1024]; - if(binNum == 0) sprintf(processCsiMCXHisto,"hPseudoProper"); - else sprintf(processCsiMCXHisto,"hPseudoProper%d",binNum); - hCsiMC = (TH1F*)f->Get(processCsiMCXHisto); - scale=1/hCsiMC->Integral(); - hCsiMC->Scale(scale); - printf ("Opening Histo with Csi_MC(x) template with n. %f entries", hCsiMC->GetEntries()); - - return hCsiMC; - -} - diff --git a/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.h b/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.h deleted file mode 100644 index 2ebb4f45708..00000000000 --- a/PWG3/vertexingHF/AliAnalysisTaskSEBtoJPSItoEle.h +++ /dev/null @@ -1,69 +0,0 @@ -#ifndef ALIANALYSISTASKSEBTOJPSITOELE_H -#define ALIANALYSISTASKSEBTOJPSITOELE_H -/* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. * - * See cxx source for full Copyright notice */ - -//************************************************************************* -// Class AliAnalysisTaskSEBtoJPSItoEle -// AliAnalysisTaskSE for the reconstruction -// of heavy-flavour decay candidates -// Author: C.Di Giglio, carmelo.digiglio@ba.infn.it -//************************************************************************* - -class TH1F; -class TList; -class AliAnalysisBtoJPSItoEle; -#include -#include -#include -#include "AliAnalysisTaskSE.h" - -class AliAnalysisTaskSEBtoJPSItoEle : public AliAnalysisTaskSE -{ - public: - - AliAnalysisTaskSEBtoJPSItoEle(); - AliAnalysisTaskSEBtoJPSItoEle(const char *name); - virtual ~AliAnalysisTaskSEBtoJPSItoEle(); - - // Implementation of interface methods - virtual void UserCreateOutputObjects(); - virtual void Init(); - virtual void LocalInit() {Init();} - virtual void UserExec(Option_t *option); - virtual void Terminate(Option_t *option); - - void SetCutsD0(const Double_t cutsD0[9]); - void SetPtCuts(const Double_t ptcuts[2]); - void SetAODMCInfo(Bool_t OkMCInfo) { fOkAODMC = OkMCInfo;} - void SetMinimize(Bool_t OkMinimize) { fOkMinimize = OkMinimize;} - void ReadAODMCInfo(AliAODEvent* aodEv, const TClonesArray* inArray); - void SetPathMCfile (TString mcfilename="CsiMCfromKine_10PtBins.root") {fNameMCfile = mcfilename;} - TH1F* OpenLocalMCFile(TString dataDir, Int_t nbin); - - private: - - AliAnalysisTaskSEBtoJPSItoEle(const AliAnalysisTaskSEBtoJPSItoEle &source); - AliAnalysisTaskSEBtoJPSItoEle& operator=(const AliAnalysisTaskSEBtoJPSItoEle& source); - // - TList *fOutput; //! list send on output slot 0 - AliAnalysisBtoJPSItoEle *fCdfFit; // Unbinned log-likelihood minimizer - TNtuple *fNtupleJPSI; //! Ntuple of pseudo-proper decay time & invariant mass values - TH1F *fhDecayTimeMCjpsifromB; //! Pseudo-proper decay time distribution used as template for JPSIs from B - TH1F *fhDecayTime; //! Pseudo-proper decay time distribution - TH1F *fhInvMass; //! Invariant mass distribution - TH1F *fhD0; //! Impact parameter distribution - TH1F *fhD0D0; //! Product of impact parameters distributions - TH1F *fhCosThetaStar; //! Cosine of decay angle distribution - TH1F *fhCosThetaPointing; //! Cosine of pointing angle distribution - TH2F *fhCtsVsD0D0; //! Cos theta star Vs. D0D0 distribution - Bool_t fOkAODMC; // Flag to read AOD monte carlo information - TString fNameMCfile; // Name of the MC file with X template - Bool_t fOkMinimize; // Flag to enable unbinned log-likelihood minimization - - Double_t fCuts[9]; // cuts for N-tuple values selection - Double_t fPtCuts[2]; // Max and min pt of the candidates - - ClassDef(AliAnalysisTaskSEBtoJPSItoEle,1); // AliAnalysisTaskSE for the reconstruction of heavy-flavour decay candidates -}; -#endif diff --git a/PWG3/vertexingHF/AliAnalysisTaskSEJPSItoEle.cxx b/PWG3/vertexingHF/AliAnalysisTaskSEJPSItoEle.cxx deleted file mode 100644 index 3af8cab88db..00000000000 --- a/PWG3/vertexingHF/AliAnalysisTaskSEJPSItoEle.cxx +++ /dev/null @@ -1,838 +0,0 @@ -/************************************************************************** - * Copyright(c) 1998-2009, 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. * - **************************************************************************/ -//************************************************************************* -// Class AliAnalysisTaskSEJPSItoEle -// AliAnalysisTaskSE class to read both J/psi -> e+e- candidates and -// like-sign pair candidates and to store them into a specific -// stand-alone AOD file for J/psi into dieletrons analysis. -// The current task: -// -// 1) produces several histograms (including invariant mass distributions) -// for both unlike sign (US) and like sign (LS) pairs. -// -// 2) selects only J/Psi to dielectrons candidates and copies it to a -// specific AOD file, namely AliAODjpsi.root. The references -// to AliAODTrack objects are copied as well. -// The specific AliAODjpsi.root is thought as the input file -// to the AliAnalysisBtoJPSItoEle.h class, which performs (locally) -// the analysis of J/Psi from beauty hadrons decays. -// -// Author: C.Di Giglio, carmelo.digiglio@ba.infn.it -//************************************************************************* -class AliAnalysisTaskSE; -class AliESDEvent; -class AliVEvent; -class iostream; -class TSystem; -class TROOT; -#include "TFile.h" -#include "TClonesArray.h" -#include "TDatabasePDG.h" -#include "TROOT.h" -#include "TCanvas.h" -#include "TBits.h" //********************** - -#include "AliAODEvent.h" -#include "AliAODMCParticle.h" -#include "AliAODMCHeader.h" -#include "AliAODRecoDecayHF2Prong.h" -#include "AliAODHandler.h" -#include "AliAnalysisManager.h" -#include "AliAnalysisTaskSEJPSItoEle.h" -#include "AliAODTrack.h" -#include "AliExternalTrackParam.h" - -ClassImp(AliAnalysisTaskSEJPSItoEle) - -//______________________________________________________________________________ -AliAnalysisTaskSEJPSItoEle::AliAnalysisTaskSEJPSItoEle() : -AliAnalysisTaskSE(), -fOutput(0), -fhDecayTimeMCjpsifromB(0), -fhDecayTime(0), -fhDecayTimeOut(0), -fhInvMass(0), -fhdEdxTPC(0), -fhD0(0), -fhD0D0(0), -fhCosThetaStar(0), -fhCosThetaPointing(0), -fhDCA(0), -fhCtsVsD0D0(0), -fHistMassLS(0), -fHistCtsLS(0), -fHistCtsLSpos(0), -fHistCtsLSneg(0), -fHistCPtaLS(0), -fHistd0d0LS(0), -fHistDCALS(0), -fVHF(0), -fTotPosPairs(0), -fTotNegPairs(0), -fLsNormalization(1.), -fOkAODMC(kFALSE), -fOkLikeSign(kFALSE), -fVerticesHFTClArr(0), -fJpsiToEleTClArr(0), -fLikeSignTClArr(0), -fTracksTClArr(0), -fOrigAOD(0), -fNewAOD(0) -{ - // default constructor -} -//_________________________________________________________________________________________________ -AliAnalysisTaskSEJPSItoEle::AliAnalysisTaskSEJPSItoEle(const char* name) : -AliAnalysisTaskSE(name), -fOutput(0), -fhDecayTimeMCjpsifromB(0), -fhDecayTime(0), -fhDecayTimeOut(0), -fhInvMass(0), -fhdEdxTPC(0), -fhD0(0), -fhD0D0(0), -fhCosThetaStar(0), -fhCosThetaPointing(0), -fhDCA(0), -fhCtsVsD0D0(0), -fHistMassLS(0), -fHistCtsLS(0), -fHistCtsLSpos(0), -fHistCtsLSneg(0), -fHistCPtaLS(0), -fHistd0d0LS(0), -fHistDCALS(0), -fVHF(0), -fTotPosPairs(0), -fTotNegPairs(0), -fLsNormalization(1.), -fOkAODMC(kFALSE), -fOkLikeSign(kFALSE), -fVerticesHFTClArr(0), -fJpsiToEleTClArr(0), -fLikeSignTClArr(0), -fTracksTClArr(0), -fOrigAOD(0), -fNewAOD(0) -{ - // standard constructor - - // Input slot #0 works with an Ntuple - DefineInput(0, TChain::Class()); - - // Output slot #0 writes into a TTree container - // Output slot #1 writes into a TList container - DefineOutput(0, TTree::Class()); - DefineOutput(1,TList::Class()); //My private output -} -//_________________________________________________________________________________________________ -AliAnalysisTaskSEJPSItoEle::~AliAnalysisTaskSEJPSItoEle() -{ - // destructor - - if (fOutput) { - delete fOutput; - fOutput = 0; - } - - if (fVHF) { - delete fVHF; - fVHF = 0; - } - -} -//_________________________________________________________________________________________________ -void AliAnalysisTaskSEJPSItoEle::Init() -{ - // Initialization - - if(fDebug > 1) printf("AliAnalysisTaskSEJPSItoEle::Init() \n"); - - gROOT->LoadMacro("ConfigVertexingHF.C"); - - fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()"); - fVHF->PrintStatus(); - - return; -} -//_________________________________________________________________________________________________ -void AliAnalysisTaskSEJPSItoEle::UserCreateOutputObjects() -{ - // Create the output container - - if(fDebug > 1) printf("AliAnalysisTaskSEJPSItoEle::UserCreateOutputObjects() \n"); - - if (!AODEvent()) { - Fatal("UserCreateOutputObjects", "This task needs an AOD handler"); - return; - } - - if(!fNewAOD) fNewAOD = new AliAODEvent(); - fNewAOD->CreateStdContent(); - - TString filename = "AliAOD.Jpsitoele.root"; - if (!IsStandardAOD()) filename = ""; - - // Array of HF vertices - fVerticesHFTClArr = new TClonesArray("AliAODVertex", 0); - fVerticesHFTClArr->SetName("VerticesHF"); - AddAODBranch("TClonesArray", &fVerticesHFTClArr); - - // Array of J/psi -> e+e- candidates - fJpsiToEleTClArr = new TClonesArray("AliAODRecoDecayHF2Prong",0); - fJpsiToEleTClArr->SetName("JpsiToEleEvents"); - AddAODBranch("TClonesArray", &fJpsiToEleTClArr, filename); - - // Array of like sign candidates - fLikeSignTClArr = new TClonesArray("AliAODRecoDecayHF2Prong",0); - fLikeSignTClArr->SetName("LikeSignEvents"); - AddAODBranch("TClonesArray", &fLikeSignTClArr, filename); - - // Array of tracks from J/psi -> e+e- candidates - fTracksTClArr = new TClonesArray("AliAODTrack", 0); - fTracksTClArr->SetName("Tracks"); - AddAODBranch("TClonesArray", &fTracksTClArr, filename); - - fOutput = new TList(); - fOutput->SetOwner(); - - if(fOkAODMC){ - fhDecayTimeMCjpsifromB = new TH1F("fhDecayTimeMCjpsifromB", "Secondary J/#Psi Monte Carlo pseudo proper decay time; X [#mu m]; Entries",100,-4000.,4000.); - fhDecayTimeMCjpsifromB->Sumw2(); - fhDecayTimeMCjpsifromB->SetMinimum(0); - fOutput->Add(fhDecayTimeMCjpsifromB); - } - - // Invariant mass - fhInvMass = new TH1F("fhInvMass", "J/#Psi invariant mass; Inv. mass [GeV]; Entries/5 MeV",100,2.75,3.25); - fhInvMass->Sumw2(); - fhInvMass->SetMinimum(0); - fOutput->Add(fhInvMass); - - fHistMassLS = new TH1F("fHistMassLS", "Like sign pairs invariant mass; Inv. mass [GeV]; Entries/5 MeV",100,2.75,3.25); - fHistMassLS->Sumw2(); - fHistMassLS->SetMinimum(0); - fOutput->Add(fHistMassLS); - - // TPC signal (only for candidate tracks) - fhdEdxTPC = new TH2F("fhdEdxTPC","TPC signal (dE/dx)",100,0.,10.,1000,0.,100.); - fhdEdxTPC->SetXTitle("p_{T} [GeV]"); - fhdEdxTPC->SetYTitle("TPC signal (arb. units)"); - fOutput->Add(fhdEdxTPC); - - // Pseudo proper decay time - fhDecayTime = new TH1F("fhDecayTime", "J/#Psi pseudo proper decay time; X [#mu m]; Entries",200,-2000.,2000.); - fhDecayTime->Sumw2(); - fhDecayTime->SetMinimum(0); - fOutput->Add(fhDecayTime); - - // Pseudo proper decay time - fhDecayTimeOut = new TH1F("fhDecayTimeOut", "J/#Psi pseudo proper decay time (output standalone AOD); X [#mu m]; Entries",200,-2000.,2000.); - fhDecayTimeOut->Sumw2(); - fhDecayTimeOut->SetMinimum(0); - fOutput->Add(fhDecayTimeOut); - - // Dictance of closest approach - fhDCA = new TH1F("fhDCA", "J/#Psi distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.); - fhDCA->Sumw2(); - fhDCA->SetMinimum(0); - fOutput->Add(fhDCA); - - fHistDCALS = new TH1F("fHistDCALS", "Like sign pairs distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.); - fHistDCALS->Sumw2(); - fHistDCALS->SetMinimum(0); - fOutput->Add(fHistDCALS); - - // Impact parameter - fhD0 = new TH1F("fhD0", "Impact parameter; d0 [#mu m]; Entries",100,-5000.,5000.); - fhD0->Sumw2(); - fhD0->SetMinimum(0); - fOutput->Add(fhD0); - - // Product of impact parameters - fhD0D0 = new TH1F("fhD0D0", "Product of impact parameters; D0D0 [#mu m^2]; Entries",100,-100000.,100000.); - fhD0D0->Sumw2(); - fhD0D0->SetMinimum(0); - fOutput->Add(fhD0D0); - - fHistd0d0LS = new TH1F("fHistd0d0LS", "Like sign pairs product of impact parameters; d0xd0 [#mu m^{2}]; Entries",200,-100000.,100000.); - fHistd0d0LS->Sumw2(); - fHistd0d0LS->SetMinimum(0); - fOutput->Add(fHistd0d0LS); - - // Cosine of the decay angle - fhCosThetaStar = new TH1F("fhCosThetaStar", "Cosine of decay angle; Cosine Theta Star; Entries",50,-1.2,1.2); - fhCosThetaStar->Sumw2(); - fhCosThetaStar->SetMinimum(0); - fOutput->Add(fhCosThetaStar); - - fHistCtsLS = new TH1F("fHistCtsLS", "Like sign pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.); - fHistCtsLS->Sumw2(); - fHistCtsLS->SetMinimum(0); - fOutput->Add(fHistCtsLS); - - fHistCtsLSpos = new TH1F("fHistCtsLSpos", "Like sign ++ pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.); - fHistCtsLSpos->Sumw2(); - fHistCtsLSpos->SetMinimum(0); - fOutput->Add(fHistCtsLSpos); - - fHistCtsLSneg = new TH1F("fHistCtsLSneg", "Like sign -- pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.); - fHistCtsLSneg->Sumw2(); - fHistCtsLSneg->SetMinimum(0); - fOutput->Add(fHistCtsLSneg); - - // Cosine of pointing angle - fhCosThetaPointing = new TH1F("fhCosThetaPointing", "Cosine of pointing angle; Cosine Pointing Angle; Entries",100,-1,1); - fhCosThetaPointing->Sumw2(); - fhCosThetaPointing->SetMinimum(0); - fOutput->Add(fhCosThetaPointing); - - fHistCPtaLS = new TH1F("fHistCPtaLS", "Like sign pairs cosine of pointing angle; Cos#Theta_{point}; Entries",200,-1.,1.); - fHistCPtaLS->Sumw2(); - fHistCPtaLS->SetMinimum(0); - fOutput->Add(fHistCPtaLS); - - // CosThetaStar vs. d0xd0 - fhCtsVsD0D0 = new TH2F("fhCtsVsD0D0", "Cosine theta star Vs. D0; Cosine theta star; D0D0",50,-1,1,100,-100000,100000); - fhCtsVsD0D0->Sumw2(); - fhCtsVsD0D0->SetMinimum(0); - fOutput->Add(fhCtsVsD0D0); - - return; -} -//_________________________________________________________________________________________________ -void AliAnalysisTaskSEJPSItoEle::UserExec(Option_t */*option*/) -{ - // Execute analysis for current event: - - AliAODEvent *aod = dynamic_cast (InputEvent()); - - TClonesArray *arrayJPSItoEle = 0; - TClonesArray *arrayLikeSign = 0; - TClonesArray *arrayTracks = 0; - - if(!aod && AODEvent() && IsStandardAOD()) { - // In case there is an AOD handler writing a standard AOD, use the AOD - // event in memory rather than the input (ESD) event. - aod = dynamic_cast (AODEvent()); - // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root) - // have to taken from the AOD event hold by the AliAODExtension - AliAODHandler* aodHandler = (AliAODHandler*) - ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()); - - if(aodHandler->GetExtensions()) { - AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root"); - AliAODEvent *aodFromExt = ext->GetAOD(); - - // load tracks from AOD event - arrayTracks=(TClonesArray*)aod->GetList()->FindObject("tracks"); - Int_t totTracks = arrayTracks->GetEntriesFast(); - if(fDebug>1) printf("Number of tracks === %d \n",totTracks); - - // load Jpsi candidates from AOD friend - arrayJPSItoEle=(TClonesArray*)aodFromExt->GetList()->FindObject("JPSItoEle"); - Int_t totJPSItoEleCand = arrayJPSItoEle->GetEntriesFast(); - if(fDebug>1) printf("Number of J/psi->ee candidates === %d \n",totJPSItoEleCand); - - // load like sign candidates from AOD friend - arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign2Prong"); - Int_t totLikeSignCand = arrayLikeSign->GetEntriesFast(); - if(fDebug>1) printf("Number of like sign candidates === %d \n",totLikeSignCand); - } - - } else { - // load tracks from AOD event - arrayTracks=(TClonesArray*)aod->GetList()->FindObject("tracks"); - Int_t totTracks = arrayTracks->GetEntriesFast(); - if(fDebug>1) printf("Number of tracks === %d \n",totTracks); - - // load Jpsi candidates - arrayJPSItoEle=(TClonesArray*)aod->GetList()->FindObject("JPSItoEle"); - Int_t totJPSItoEleCand = arrayJPSItoEle->GetEntriesFast(); - if(fDebug>1) printf("Number of J/psi->ee candidates === %d \n",totJPSItoEleCand); - - // load like sign candidates - arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign2Prong"); - Int_t totLikeSignCand = arrayLikeSign->GetEntriesFast(); - if(fDebug>1) printf("Number of like sign candidates === %d \n",totLikeSignCand); - } - - fOrigAOD = aod; // copy pointer to the current AliAODEvent in the data member fOrigAOD - if (!aod) return; - - if(!arrayTracks) { - printf("AliAnalysisTaskSEJPSItoEle::UserExec: Tracks branch not found!\n"); - return; - } - if(!arrayJPSItoEle) { - printf("AliAnalysisTaskSEJPSItoEle::UserExec: JPSItoEle branch not found!\n"); - return; - } - if(!arrayLikeSign) { - printf("AliAnalysisTaskSEJPSItoEle::UserExec: LikeSign2Prong branch not found!\n"); - return; - } - - // fix for temporary bug in ESDfilter - // the AODs with null vertex pointer didn't pass the PhysSel - if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return; - - // load MC particles and read MC info (for sim only) - TClonesArray* mcArray=0; - if(fOkAODMC){ - mcArray = dynamic_cast(aod->FindListObject(AliAODMCParticle::StdBranchName())); - if (!mcArray) AliError("Could not find Monte-Carlo in AOD"); - ReadAODMCInfo(aod,arrayJPSItoEle); - } - - // retrieve AOD primary vertex - AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex(); - - Int_t iOutVerticesHF = 0, iOutJPSItoEle = 0, iOutLikeSign = 0, iOutTracks = 0; - - fVerticesHFTClArr->Delete(); - iOutVerticesHF = fVerticesHFTClArr->GetEntriesFast(); - TClonesArray &verticesHFRef = *fVerticesHFTClArr; - - fJpsiToEleTClArr->Delete(); - iOutJPSItoEle = fJpsiToEleTClArr->GetEntriesFast(); - TClonesArray &arrayJpsiToEleRef = *fJpsiToEleTClArr; - - fLikeSignTClArr->Delete(); - iOutLikeSign = fLikeSignTClArr->GetEntriesFast(); - TClonesArray &arrayLikeSignRef = *fLikeSignTClArr; - - fTracksTClArr->Delete(); - iOutTracks = fTracksTClArr->GetEntriesFast(); - //TClonesArray &arrayTrackRef = *fTracksTClArr; - - // - // LOOP OVER LIKE SIGN CANDIDATES - // - - // Access to the new AOD container of tracks - - Int_t trkIDtoEntry[100000]; - AliAODTrack* trackin=0; - for(Int_t it=0;itGetNumberOfTracks();it++) { - trackin = aod->GetTrack(it); - trkIDtoEntry[trackin->GetID()]=it; - } -/* - for(Int_t it=0;itGetNumberOfTracks();it++) { - vtx1->AddDaughter(trackout = new(arrayTrackRef[iOutTracks++]) AliAODTrack((*(aod->GetTrack(trkIDtoEntry[it]))))); - trackin = (AliAODTrack*)arrayTracks->UncheckedAt(ntracks); - printf("trackin ==== %d \n", trackin); - AliAODTrack* trackout = new(arrayTrackRef[iOutTracks++]) AliAODTrack(*trackin); - } - const AliAODTrack* trackin; - Int_t numTracks = arrayTracks->GetEntriesFast(); - for(Int_t ntracks=0; ntracksUncheckedAt(ntracks); - printf("trackin ==== %d \n", trackin); - AliAODTrack* trackout = new(arrayTrackRef[iOutTracks++]) AliAODTrack(*trackin); - } -*/ - - Int_t nPosPairs=0,nNegPairs=0; - Int_t nLikeSign = arrayLikeSign->GetEntriesFast(); - AliAODRecoDecayHF2Prong* dlsout = 0; - //if(fDebug>1) printf("+++\n+++ Number of total like sign pairs (before cuts)---> %d \n+++\n", nLikeSign); - - for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) { - AliAODRecoDecayHF2Prong *dlsin = (AliAODRecoDecayHF2Prong*)arrayLikeSign->UncheckedAt(iLikeSign); - Bool_t unsetvtx=kFALSE; - if(!dlsin->GetOwnPrimaryVtx()) { - dlsin->SetOwnPrimaryVtx(vtx1); // needed to compute all variables - unsetvtx=kTRUE; - } - //Int_t okBtoJPSIls=0; - //if(dlsin->Pt() < fPtCuts[1] && dlsin->Pt() > fPtCuts[0]){ // apply pt cut only - //if(dlsin->SelectBtoJPSI(fVHF->GetBtoJPSICuts(),okBtoJPSIls)) { - AliAODTrack *trk0 = (AliAODTrack*)dlsin->GetDaughter(0); - fHistMassLS->Fill(dlsin->InvMassJPSIee()); - fHistCPtaLS->Fill(dlsin->CosPointingAngle()); - fHistd0d0LS->Fill(1e8*dlsin->Prodd0d0()); - fHistDCALS->Fill(100*dlsin->GetDCA()); - if(!trk0) { - trk0=aod->GetTrack(trkIDtoEntry[dlsin->GetProngID(0)]); - printf("references to standard AOD not available \n"); - } - if((trk0->Charge())==1) { - nPosPairs++; - fHistCtsLS->Fill(dlsin->CosThetaStar(0,443,11,11)); - fHistCtsLSpos->Fill(dlsin->CosThetaStar(0,443,11,11)); - } else { - nNegPairs++; - fHistCtsLS->Fill(dlsin->CosThetaStarJPSI()); - fHistCtsLSneg->Fill(dlsin->CosThetaStarJPSI()); - } - PostData(1,fOutput); - - // Clone like sign candidate for output AOD - dlsout = new(arrayLikeSignRef[iOutLikeSign++]) AliAODRecoDecayHF2Prong(*dlsin); - - //} - if(unsetvtx) dlsin->UnsetOwnPrimaryVtx(); - } - - //if(fDebug>1) printf("+++\n+++ N. of positive pairs passing cuts in Event ----- %d \n+++\n", nPosPairs); - //if(fDebug>1) printf("+++\n+++ N. of negative pairs passing cuts in Event ----- %d \n+++\n", nNegPairs); - - fTotPosPairs += nPosPairs; - fTotNegPairs += nNegPairs; - - // - // LOOP OVER J/psi->e+e- CANDIDATES - // - - Int_t nInJPSItoEle = arrayJPSItoEle->GetEntriesFast(); - //if(fDebug>1) printf("+++\n+++ Number of total like JPSI -> ee candidates (before cuts)---> %d \n+++\n", nInJPSItoEle); - - //totJPSIin += nInJPSItoEle; - Bool_t isOkEle[2] = {kFALSE,kFALSE}; - - for (Int_t iJPSItoEle = 0; iJPSItoEle < nInJPSItoEle; iJPSItoEle++) { - - AliAODRecoDecayHF2Prong *dIn = (AliAODRecoDecayHF2Prong*)arrayJPSItoEle->UncheckedAt(iJPSItoEle); - Int_t mcLabel = 0; - if(fOkAODMC) mcLabel = dIn->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI - - //Secondary vertex - Double_t vtxSec[3] = {0.,0.,0.}; - Double_t vtxPrim[3] = {0.,0.,0.}; - Double_t vtxSecOut[3] = {0.,0.,0.}; - Double_t vtxPrimOut[3] = {0.,0.,0.}; - dIn->GetSecondaryVtx(vtxSec); - Bool_t unsetvtx=kFALSE; - if(!dIn->GetOwnPrimaryVtx()) { - dIn->SetOwnPrimaryVtx(vtx1); // needed to compute all variables - unsetvtx=kTRUE; - } - - //Int_t okBtoJPSI=0; - //if(dIn->Pt() < fPtCuts[1] && dIn->Pt() > fPtCuts[0]){ // apply pt cut only - //if(dIn->SelectBtoJPSI(fVHF->GetBtoJPSICuts(),okBtoJPSI)) { - if ( fOkAODMC && mcLabel == -1){AliDebug(2,"No MC particle found");} else { - - //fhInvMass->Fill(dIn->InvMassJPSIee()); - fhD0->Fill(10000*dIn->ImpParXY()); - fhD0D0->Fill(1e8*dIn->Prodd0d0()); - fhCosThetaStar->Fill(dIn->CosThetaStarJPSI()); - fhCtsVsD0D0->Fill(dIn->CosThetaStarJPSI(),1e8*dIn->Prodd0d0()); - fhCosThetaPointing->Fill(dIn->CosPointingAngle()); - fhDCA->Fill(100*dIn->GetDCA()); - - // compute X variable - AliAODVertex* primVertex = dIn->GetOwnPrimaryVtx(); - vtxPrim[0] = primVertex->GetX(); - vtxPrim[1] = primVertex->GetY(); - vtxPrim[2] = primVertex->GetZ(); - Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dIn->Px()) + (vtxSec[1]-vtxPrim[1])*(dIn->Py()))/dIn->Pt(); - Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dIn->Pt(); - fhDecayTime->Fill(10000*pseudoProperDecayTime); - - // retrieve daughter tracks - AliAODTrack *trk0 = (AliAODTrack*)dIn->GetDaughter(0); - TBits clusterMap0 = trk0->GetTPCClusterMap(); - Int_t npoints0 = clusterMap0.CountBits(0); - AliAODPid *pid0 = trk0->GetDetPid(); - if(GetNumberOfSigmas(trk0->P(),pid0->GetTPCsignal(),npoints0,AliPID::kElectron) < 3.) isOkEle[0]=kTRUE; - - AliAODTrack *trk1 = (AliAODTrack*)dIn->GetDaughter(1); - TBits clusterMap1 = trk1->GetTPCClusterMap(); - Int_t npoints1 = clusterMap1.CountBits(0); - AliAODPid *pid1 = trk1->GetDetPid(); - if(GetNumberOfSigmas(trk1->P(),pid1->GetTPCsignal(),npoints1,AliPID::kElectron) < 3.) isOkEle[1]=kTRUE; - - if(isOkEle[0] && isOkEle[1]) { - fhInvMass->Fill(dIn->InvMassJPSIee()); - fhdEdxTPC->Fill(pid0->GetTPCmomentum(),pid0->GetTPCsignal()); - fhdEdxTPC->Fill(pid1->GetTPCmomentum(),pid1->GetTPCsignal()); - } - - //printf("TPC momentum %f\n", pid0->GetTPCmomentum()); - //printf("TPC signal %f\n", pid0->GetTPCsignal()); - - // clone candidate for output AOD - AliAODVertex *v = new(verticesHFRef[iOutVerticesHF++]) - AliAODVertex(*(dIn->GetSecondaryVtx())); - AliAODRecoDecayHF2Prong* dOut = new(arrayJpsiToEleRef[iOutJPSItoEle++]) - AliAODRecoDecayHF2Prong(*dIn); // copy constructor used - dOut->SetSecondaryVtx(v); - dOut->SetOwnPrimaryVtx((AliAODVertex*)((dIn->GetOwnPrimaryVtx())->Clone())); - v->SetParent(dOut); - - // compute X variable from the cloned candidate in the stand-alone AOD file - dOut->GetSecondaryVtx(vtxSecOut); - AliAODVertex* primVertexOut = dOut->GetOwnPrimaryVtx(); - vtxPrimOut[0] = primVertexOut->GetX(); - vtxPrimOut[1] = primVertexOut->GetY(); - vtxPrimOut[2] = primVertexOut->GetZ(); - Double_t lxyOut = ((vtxSecOut[0]-vtxPrimOut[0])*(dOut->Px()) + (vtxSecOut[1]-vtxPrimOut[1])*(dOut->Py()))/dOut->Pt(); - Double_t pseudoProperDecayTimeOut = lxyOut*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dOut->Pt(); - fhDecayTimeOut->Fill(10000*pseudoProperDecayTimeOut); - - //AliAODTrack* trk0 = (AliAODTrack*)dOut->GetDaughter(0); - //printf("+++ Pt first daughter = %d \n +++ \n", trk0->Pt()); - - } - - //} // end of JPSItoEle candidates selection according to cuts - - if(unsetvtx) dIn->UnsetOwnPrimaryVtx(); - - PostData(0,fOutput); - - }// end loop on JPSI to ele candidates - - //printf("+++\n+++ Number of selected J/psi->e+e-: %d\n+++\n",iOutJPSItoEle); - - //totJPSIout += iOutJPSItoEle; - - return; - -} -//_________________________________________________________________________________________________ -void AliAnalysisTaskSEJPSItoEle::Terminate(Option_t */*option*/) -{ - // - // Terminate analysis - // - - if(fDebug > 1) printf("AliAnalysisTaskSEJPSItoEle: Terminate() \n"); - - fOutput = dynamic_cast (GetOutputData(1)); - if (!fOutput) { - printf("ERROR: fOutput not available\n"); - return; - } - - fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs*fTotNegPairs); - - if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast(fOutput->FindObject("fhDecayTimeMCjpsifromB")); - fhDecayTime = dynamic_cast(fOutput->FindObject("fhDecayTime")); - fhDecayTimeOut = dynamic_cast(fOutput->FindObject("fhDecayTimeOut")); - fhInvMass = dynamic_cast(fOutput->FindObject("fhInvMass")); - fhdEdxTPC = dynamic_cast(fOutput->FindObject("fhdEdxTPC")); - fhD0 = dynamic_cast(fOutput->FindObject("fhD0")); - fhD0D0 = dynamic_cast(fOutput->FindObject("fhD0D0")); - fhCosThetaStar = dynamic_cast(fOutput->FindObject("fhCosThetaStar")); - fhCosThetaPointing = dynamic_cast(fOutput->FindObject("fhCosThetaPointing")); - fhCtsVsD0D0 = dynamic_cast(fOutput->FindObject("fhCtsVsD0D0")); - fHistMassLS = dynamic_cast(fOutput->FindObject("fHistMassLS")); - fHistCtsLS = dynamic_cast(fOutput->FindObject("fHistCtsLS")); - fHistCtsLSpos = dynamic_cast(fOutput->FindObject("fHistCtsLSpos")); - fHistCtsLSneg = dynamic_cast(fOutput->FindObject("fHistCtsLSneg")); - fHistCPtaLS = dynamic_cast(fOutput->FindObject("fHistCPtaLS")); - fHistd0d0LS = dynamic_cast(fOutput->FindObject("fHistd0d0LS")); - fhDCA = dynamic_cast(fOutput->FindObject("fhDCA")); - fHistDCALS = dynamic_cast(fOutput->FindObject("fHistDCALS")); - - if(fLsNormalization>0.) { - fHistMassLS->Scale((1/fLsNormalization)*fHistMassLS->GetEntries()); - fHistCtsLS->Scale((1/fLsNormalization)*fHistCtsLS->GetEntries()); - fHistCtsLSpos->Scale((1/fLsNormalization)*fHistCtsLSpos->GetEntries()); - fHistCtsLSneg->Scale((1/fLsNormalization)*fHistCtsLSneg->GetEntries()); - fHistCPtaLS->Scale((1/fLsNormalization)*fHistCPtaLS->GetEntries()); - fHistd0d0LS->Scale((1/fLsNormalization)*fHistd0d0LS->GetEntries()); - fHistDCALS->Scale((1/fLsNormalization)*fHistDCALS->GetEntries()); - } - - //printf("Tot JPSI in %d\n", totJPSIin); - //printf("Tot JPSI out %d\n", totJPSIout); - - return; -} -//_________________________________________________________________________________________________ -void AliAnalysisTaskSEJPSItoEle::SetPtCuts(const Double_t ptCuts[2]) -{ - // - // apply Pt cuts (lower cut, upper cut) - // - for(Int_t i = 0; i < 2; i++) fPtCuts[i] = ptCuts[i]; -} -//_________________________________________________________________________________________________ -void AliAnalysisTaskSEJPSItoEle::SetCutsJPSI(const Double_t cuts[9]) -{ - // - // apply JPSI cuts - // - for(Int_t i = 0; i < 9; i++) fCuts[i] = cuts[i]; -} -//_________________________________________________________________________________________________ -void AliAnalysisTaskSEJPSItoEle::ReadAODMCInfo(AliAODEvent* aodEvent, const TClonesArray* inputArray) -{ - // - // Reads MC information if needed (for example in case of parametrized PID) - // - - // load MC particles - TClonesArray* mcArray = - dynamic_cast(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); - if (!mcArray) AliError("Could not find Monte-Carlo in AOD"); - AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast())); - - // load MC header - AliAODMCHeader* mcHeader = - (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()); - if(!mcHeader){ - printf("AliAnalysisTaskSEJPSItoEle::UserExec: MC AOD header branch not found!\n"); - } - - AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex(); - - Double_t mcPrimVtx[3]; - - // get MC primary Vtx - mcHeader->GetVertex(mcPrimVtx); - - Int_t nInCandidates = inputArray->GetEntriesFast(); - printf("Number of Candidates for MC analysis: %d\n",nInCandidates); - - Int_t lab0, lab1, pdgMother, labMother, pdgJPSI, pdg0, pdg1, labJPSIMother, pdgJPSIMother; - Int_t labJPSIdaugh0=0; - Int_t labJPSIdaugh1=0; - - for (Int_t iCandidate = 0; iCandidate < nInCandidates; iCandidate++) { - - AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iCandidate); - Int_t mcLabelSec = dd->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI - - Double_t vtxPrim[3] = {0.,0.,0.}; - Double_t vtxSec[3] = {0.,0.,0.}; - dd->GetSecondaryVtx(vtxSec); - Bool_t unsetvtx=kFALSE; - if(!dd->GetOwnPrimaryVtx()) { - dd->SetOwnPrimaryVtx(vtx1); // needed to compute all variables - unsetvtx=kTRUE; - } - - if(dd->Pt() < fPtCuts[1] && dd->Pt() > fPtCuts[0]){ // apply pt cut only - - AliAODTrack *trk0 = (AliAODTrack*)dd->GetDaughter(0); - AliAODTrack *trk1 = (AliAODTrack*)dd->GetDaughter(1); - lab0 = trk0->GetLabel(); - lab1 = trk1->GetLabel(); - - AliAODMCParticle* part0 = (AliAODMCParticle*)mcArray->At(lab0); - if(!part0) { - printf("no MC particle\n"); - continue; - } - - while(part0->GetMother()>=0) { - labMother=part0->GetMother(); - part0 = (AliAODMCParticle*)mcArray->At(labMother); - if(!part0) { - printf("no MC mother particle\n"); - break; - } - pdgMother = TMath::Abs(part0->GetPdgCode()); - if(pdgMother==443) {//this for JPSI - labJPSIdaugh0=labMother; - break; - } - } - - AliAODMCParticle* part1 = (AliAODMCParticle*)mcArray->At(lab1); - if(!part1) { - printf("no MC particle\n"); - continue; - } - - while(part1->GetMother()>=0) { - labMother=part1->GetMother(); - part1 = (AliAODMCParticle*)mcArray->At(labMother); - if(!part1) { - printf("no MC mother particle\n"); - break; - } - pdgMother = TMath::Abs(part1->GetPdgCode()); - if(pdgMother==443) {//this for JPSI - labJPSIdaugh1=labMother; - break; - } - } - - if (mcLabelSec == -1) - { - AliDebug(2,"No MC particle found"); - } - else { - - if(labJPSIdaugh0>=0 && labJPSIdaugh1>=0 && labJPSIdaugh0==labJPSIdaugh1) { // check if JPSI is signal - AliAODMCParticle* partJPSI = (AliAODMCParticle*)mcArray->At(labJPSIdaugh0); - pdgJPSI = partJPSI->GetPdgCode(); - Int_t pdaugh0 = partJPSI->GetDaughter(0); - Int_t pdaugh1 = partJPSI->GetDaughter(1); - AliAODMCParticle* partDaugh0 = (AliAODMCParticle*)mcArray->At(pdaugh0); - AliAODMCParticle* partDaugh1 = (AliAODMCParticle*)mcArray->At(pdaugh1); - pdg0 = partDaugh0->GetPdgCode(); - pdg1 = partDaugh1->GetPdgCode(); - if(TMath::Abs(pdg0) == 11 && TMath::Abs(pdg1) == 11){ // this is for MC JPSI -> ee - labJPSIMother = partJPSI->GetMother(); - AliAODMCParticle* partJPSIMother = (AliAODMCParticle*)mcArray->At(labJPSIMother); - pdgJPSIMother = partJPSIMother->GetPdgCode(); - - // chech if JPSI comes from B - if( pdgJPSIMother==511 || pdgJPSIMother==521 || - pdgJPSIMother==10511 || pdgJPSIMother==10521 || - pdgJPSIMother==513 || pdgJPSIMother==523 || - pdgJPSIMother==10513 || pdgJPSIMother==10523 || - pdgJPSIMother==20513 || pdgJPSIMother==20523 || - pdgJPSIMother==515 || pdgJPSIMother==525 || - pdgJPSIMother==531 || pdgJPSIMother==10531 || - pdgJPSIMother==533 || pdgJPSIMother==10533 || - pdgJPSIMother==20533 || pdgJPSIMother==535 || - pdgJPSIMother==541 || pdgJPSIMother==10541 || - pdgJPSIMother==543 || pdgJPSIMother==10543 || - pdgJPSIMother==20543 || pdgJPSIMother==545) - { //this is for MC JPSI -> ee from B-hadrons - - fhInvMass->Fill(dd->InvMassJPSIee()); - fhD0->Fill(10000*dd->ImpParXY()); - fhD0D0->Fill(1e8*dd->Prodd0d0()); - fhCosThetaStar->Fill(dd->CosThetaStarJPSI()); - fhCtsVsD0D0->Fill(dd->CosThetaStarJPSI(),1e8*dd->Prodd0d0()); - fhCosThetaPointing->Fill(dd->CosPointingAngle()); - - // compute X variable - AliAODVertex* primVertex = dd->GetOwnPrimaryVtx(); - vtxPrim[0] = primVertex->GetX(); - vtxPrim[1] = primVertex->GetY(); - vtxPrim[2] = primVertex->GetZ(); - Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dd->Px()) + (vtxSec[1]-vtxPrim[1])*(dd->Py()))/dd->Pt(); - Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dd->Pt(); - - fhDecayTime->Fill(10000*pseudoProperDecayTime); - // Post the data already here - PostData(1,fOutput); - - Double_t mcLxy = ((partJPSI->Xv()-mcPrimVtx[0])*(partJPSI->Px()) + (partJPSI->Yv()-mcPrimVtx[1])*(partJPSI->Py()))/partJPSI->Pt(); - Double_t mcPseudoProperDecayTime = mcLxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/partJPSI->Pt(); - fhDecayTimeMCjpsifromB->Fill(10000*mcPseudoProperDecayTime,1); - - // Post the data already here - PostData(1,fOutput); - - } //this is for MC JPSI -> ee from B-hadrons - } //this is for MC JPSI -> ee - } - } // end of check if JPSI is signal - } - } -} diff --git a/PWG3/vertexingHF/AliAnalysisTaskSEJPSItoEle.h b/PWG3/vertexingHF/AliAnalysisTaskSEJPSItoEle.h deleted file mode 100644 index d7198f2d093..00000000000 --- a/PWG3/vertexingHF/AliAnalysisTaskSEJPSItoEle.h +++ /dev/null @@ -1,134 +0,0 @@ -#ifndef ALIANALYSISTASKSEJPSITOELE_H -#define ALIANALYSISTASKSEJPSITOELE_H -/* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. * - * See cxx source for full Copyright notice */ - -//************************************************************************* -// Class AliAnalysisTaskSEJPSItoEle -// AliAnalysisTaskSE class to select -// J/psi -> e+e- candidates only and save them -// into a specific stand-alone AOD file -// Author: C.Di Giglio, carmelo.digiglio@ba.infn.it -//************************************************************************* - -#include "TClonesArray.h" -#include "TH1F.h" -#include "TH2F.h" -#include "TList.h" -#include "TChain.h" - -#include "AliAnalysisTaskSE.h" -#include "AliAODEvent.h" -#include "AliAnalysisManager.h" -#include "AliAnalysisVertexingHF.h" -#include "AliPID.h" //************ -#include "AliAODPid.h" -#include "AliExternalTrackParam.h"//*************** - -class AliAnalysisTaskSEJPSItoEle : public AliAnalysisTaskSE -{ - public: - - AliAnalysisTaskSEJPSItoEle(); - AliAnalysisTaskSEJPSItoEle(const char *name); - virtual ~AliAnalysisTaskSEJPSItoEle(); - - virtual void UserCreateOutputObjects(); - virtual void Init(); - virtual void LocalInit() {Init();} - virtual void UserExec(Option_t *option); - virtual void Terminate(Option_t *option); - - void SetCutsJPSI(const Double_t cutsJPSI[9]); - void SetPtCuts(const Double_t ptcuts[2]); - void SetAODMCInfo(Bool_t OkMCInfo) { fOkAODMC = OkMCInfo;} - void ReadAODMCInfo(AliAODEvent* aodEv, const TClonesArray* inArray); - // - // - // - Double_t GetExpectedSignal(const Float_t mom, - AliPID::EParticleType n=AliPID::kKaon) const { - - Double_t mass=AliPID::ParticleMass(n); - Double_t betaGamma = mom/mass; - Double_t bb = AliExternalTrackParam::BetheBlochAleph(betaGamma,0.0283086,2.63394e+01,5.04114e-11,2.12543,4.88663); - return bb*50.; //bb*fMIP; - } - - Double_t GetExpectedSigma(const Float_t mom, const Int_t nPoints, - AliPID::EParticleType n=AliPID::kKaon) const { - if (nPoints != 0) - return GetExpectedSignal(mom,n)*0.07*sqrt(1. + 0./nPoints); - else - return GetExpectedSignal(mom,n)*0.07; - } - - Float_t GetNumberOfSigmas(const Float_t mom, const Float_t dEdx, - const Int_t nPoints, - AliPID::EParticleType n=AliPID::kKaon) const { - - Double_t bethe=GetExpectedSignal(mom,n); - Double_t sigma=GetExpectedSigma(mom,nPoints,n); - return (dEdx-bethe)/sigma; - } - // - // - // - private: - - AliAnalysisTaskSEJPSItoEle(const AliAnalysisTaskSEJPSItoEle &source); - AliAnalysisTaskSEJPSItoEle& operator=(const AliAnalysisTaskSEJPSItoEle& source); - - TList *fOutput; //! list send on output slot 0 - - //output histograms - TH1F *fhDecayTimeMCjpsifromB; //! Pseudo-proper decay time distribution used as template for JPSIs from B - TH1F *fhDecayTime; //! Pseudo-proper decay time distribution - TH1F *fhDecayTimeOut; //! Pseudo-proper decay time distribution (stand-alone AOD) - TH1F *fhInvMass; //! Invariant mass distribution - TH2F *fhdEdxTPC; - TH1F *fhD0; //! Impact parameter distribution - TH1F *fhD0D0; //! Product of impact parameters distributions - TH1F *fhCosThetaStar; //! Cosine of decay angle distribution - TH1F *fhCosThetaPointing; //! Cosine of pointing angle distribution - TH1F *fhDCA; //! Distance of closest approach - TH2F *fhCtsVsD0D0; //! Cos theta star Vs. D0D0 distribution - - //like sign pairs histograms - TH1F *fHistMassLS; //! Invariant mass distribution - TH1F *fHistCtsLS; //! Cosine of decay angle distribution - TH1F *fHistCtsLSpos; //! Cosine of decay angle distribution (++ pairs) - TH1F *fHistCtsLSneg; //! Cosine of decay angle distribution (-- pairs) - TH1F *fHistCPtaLS; //! Cosine of pointing angle distribution - TH1F *fHistd0d0LS; //! Product of impact parameters distributions - TH1F *fHistDCALS; //! Distance of closest approach - AliAnalysisVertexingHF *fVHF; //! Vertexer heavy flavour (used to pass the cuts) - - //Int_t totJPSIin; - //Int_t totJPSIout; - - //like sign spectrum normalization - Int_t fTotPosPairs; // - Int_t fTotNegPairs; // - Double_t fLsNormalization; // Like sign normalization factor - - //flags for analysis - Bool_t fOkAODMC; // Flag to read AOD monte carlo information - Bool_t fOkLikeSign; // Flag to select like sign candidates analysis - - //momentum cuts - Double_t fCuts[9]; // cuts for N-tuple values selection - Double_t fPtCuts[2]; // Max and min pt of the candidates - - TClonesArray *fVerticesHFTClArr; // Array of heavy flavor vertices to be replicated in stand-alone AOD - TClonesArray *fJpsiToEleTClArr; // Array of J/psi->e+e- candidates to be replicated in stand-alone AOD - TClonesArray *fLikeSignTClArr; // Array of like sign candidates to be replicated in stand-alone AOD - TClonesArray *fTracksTClArr; // Array of tracks belonging to J/psi->e+e- candidates to be replicated in stand-alone AOD - AliAODEvent *fOrigAOD; // original AOD event - AliAODEvent *fNewAOD; // new AOD event with only JPSItoEle candidates stored - - ClassDef(AliAnalysisTaskSEJPSItoEle,1); // AliAnalysisTaskSE for the reconstruction of heavy-flavour decay candidates -}; - - -#endif diff --git a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.cxx b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.cxx deleted file mode 100644 index 9134f468205..00000000000 --- a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.cxx +++ /dev/null @@ -1,618 +0,0 @@ -/************************************************************************** - * Copyright(c) 1998-2009, 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. * - **************************************************************************/ -#include "AliLog.h" -#include "AliBtoJPSItoEleCDFfitFCN.h" - -//_________________________________________________________________________ -// Class AliBtoJPSItoEleCDFfitFCN -// Definition of main function used in -// unbinned log-likelihood fit for -// the channel B -> JPsi + X -> e+e- + X -// -// Origin: C.Di Giglio -// Contact: Carmelo.Digiglio@ba.infn.it , Giuseppe.Bruno@ba.infn.it -//_________________________________________________________________________ - -ClassImp(AliBtoJPSItoEleCDFfitFCN) - -//_________________________________________________________________________________________________ -AliBtoJPSItoEleCDFfitFCN::AliBtoJPSItoEleCDFfitFCN() : -fFPlus(0.), -fFMinus(0.), -fFSym(0.), -fIntegral(1.), -fintxFunB(1.), -fintxDecayTimeBkgPos(1.), -fintxDecayTimeBkgNeg(1.), -fintxDecayTimeBkgSym(1.), -fintmMassSig(1.), -fintxRes(1.), -fintmMassBkg(1.), -fhCsiMC(0x0), -fMassWndHigh(0.), -fMassWndLow(0.), -fCrystalBallParam(kFALSE) -{ - // - // constructor - // - SetCrystalBallFunction(kFALSE); - SetMassWndHigh(0.2); - SetMassWndLow(0.5); - for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = 0.; - fParameters[9] = 1.;fParameters[11] = 1.;fParameters[12] = 1.; - for(Int_t index=0; index<5; index++) fResolutionConstants[index] = 0.; - AliInfo("Instance of AliBtoJPSItoEleCDFfitFCN-class created"); -} -//_________________________________________________________________________________________________ -AliBtoJPSItoEleCDFfitFCN::AliBtoJPSItoEleCDFfitFCN(const AliBtoJPSItoEleCDFfitFCN& source) : -TNamed(source), -fFPlus(source.fFPlus), -fFMinus(source.fFMinus), -fFSym(source.fFSym), -fIntegral(source.fIntegral), -fintxFunB(source.fintxFunB), -fintxDecayTimeBkgPos(source.fintxDecayTimeBkgPos), -fintxDecayTimeBkgNeg(source.fintxDecayTimeBkgNeg), -fintxDecayTimeBkgSym(source.fintxDecayTimeBkgSym), -fintmMassSig(source.fintmMassSig), -fintxRes(source.fintxRes), -fintmMassBkg(source.fintmMassBkg), -fhCsiMC(source.fhCsiMC), -fMassWndHigh(source.fMassWndHigh), -fMassWndLow(source.fMassWndLow), -fCrystalBallParam(source.fCrystalBallParam) -{ - // - // Copy constructor - // - for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = source.fParameters[iPar]; - for(Int_t index=0; index<5; index++) fResolutionConstants[index] = source.fResolutionConstants[index]; -} -//_________________________________________________________________________________________________ -AliBtoJPSItoEleCDFfitFCN& AliBtoJPSItoEleCDFfitFCN::operator=(const AliBtoJPSItoEleCDFfitFCN& source) -{ - // - // Assignment operator - // - if(&source == this) return *this; - fFPlus = source.fFPlus; - fFMinus = source.fFMinus; - fFSym = source.fFSym; - fIntegral = source.fIntegral; - fintxFunB = source.fintxFunB; - fintxDecayTimeBkgPos = source.fintxDecayTimeBkgPos; - fintxDecayTimeBkgNeg = source.fintxDecayTimeBkgNeg; - fintxDecayTimeBkgSym = source.fintxDecayTimeBkgSym; - fintmMassSig = source.fintmMassSig; - fintxRes = source.fintxRes; - fintmMassBkg = source.fintmMassBkg; - fhCsiMC = source.fhCsiMC; - fCrystalBallParam = source.fCrystalBallParam; - - for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = source.fParameters[iPar]; - for(Int_t index=0; index<5; index++) fResolutionConstants[index] = source.fResolutionConstants[index]; - - return *this; -} -//_________________________________________________________________________________________________ -AliBtoJPSItoEleCDFfitFCN::~AliBtoJPSItoEleCDFfitFCN() -{ - // - // Default destructor - // - - delete fhCsiMC; - for(Int_t iPar = 0; iPar < 16; iPar++) fParameters[iPar] = 0.; - for(Int_t index=0; index<5; index++) fResolutionConstants[index] = 0.; -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateLikelihood(const Double_t* pseudoproperdecaytime, - const Double_t* invariantmass, const Int_t ncand) -{ - // - // This function evaluates the Likelihood fnction - // It returns the -Log(of the likelihood function) - // - Double_t f = 0.; - Double_t ret = 0.; - - for(Int_t i=0; i < ncand; i++) { - f = EvaluateCDFfuncNorm(pseudoproperdecaytime[i],invariantmass[i]); - if(f < 0.) { - continue; - } - ret+=-1.*TMath::Log(f); - } - return ret; -} -//_________________________________________________________________________________________________ -void AliBtoJPSItoEleCDFfitFCN::SetAllParameters(const Double_t* parameters) -{ - // - // Sets array of FCN parameters - // - for(Int_t index = 0; index < 16; index++) fParameters[index] = parameters[index]; -} -//_________________________________________________________________________________________________ -void AliBtoJPSItoEleCDFfitFCN::ComputeIntegral() -{ - // - // this function compute the integral of the likelihood function - // (theoretical function) in order to normalize it to unity - // - Double_t np = 100.0; //number of integration steps - Double_t npres = 200.0; //number of integration steps for the resolution function - Double_t npm = 200.; - Double_t stepm;Double_t stepx;Double_t stepxres; //integration step width in variable m,x - Double_t mx=0.;Double_t xprime=0.; - Double_t xlow = -4000.; Double_t xup = 4000.; - Double_t i; Double_t j; - Double_t sumx = 0.0;Double_t intx = 0.0;Double_t intm = 0.0; - stepm = (fMassWndHigh-fMassWndLow)/npm; - stepx = (xup-xlow)/np; - stepxres = (xup-xlow)/npres; - -// compute integrals for all the terms - - Double_t iRes; - Double_t intxRes = 0.0; - Double_t sumxRes = 0.0; - for(iRes = 1.0; iRes <= npres/2.; iRes++) { - xprime = xlow + (iRes - .5)*stepxres; - sumxRes += ResolutionFunc(xprime); - xprime = xup - (iRes - .5)*stepxres; - sumxRes += ResolutionFunc(xprime); - } - intxRes = sumxRes*stepxres; - SetIntegralRes(intxRes); - -// - Double_t iFunB; - Double_t intxFunB = 0.0; - Double_t sumxFunB = 0.0; - for(iFunB = 1.0; iFunB <= np/2; iFunB++) { - xprime = xlow + (iFunB - .5)*stepx; - sumxFunB += FunB(xprime); - xprime = xup - (iFunB - .5)*stepx; - sumxFunB += FunB(xprime); - } - intxFunB = sumxFunB*stepx; - SetIntegralFunB(intxFunB); - -// - Double_t iDecayTimeBkgPos; - Double_t intxDecayTimeBkgPos = 0.0; - Double_t sumxDecayTimeBkgPos = 0.0; - for(iDecayTimeBkgPos = 1.0; iDecayTimeBkgPos <= np/2; iDecayTimeBkgPos++) { - xprime = xlow + (iDecayTimeBkgPos - .5)*stepx; - sumxDecayTimeBkgPos += FunBkgPos(xprime); - xprime = xup - (iDecayTimeBkgPos - .5)*stepx; - sumxDecayTimeBkgPos += FunBkgPos(xprime); - } - intxDecayTimeBkgPos = sumxDecayTimeBkgPos*stepx; - SetIntegralBkgPos(intxDecayTimeBkgPos); - -// - Double_t iDecayTimeBkgNeg; - Double_t intxDecayTimeBkgNeg = 0.0; - Double_t sumxDecayTimeBkgNeg = 0.0; - for(iDecayTimeBkgNeg = 1.0; iDecayTimeBkgNeg<= np/2; iDecayTimeBkgNeg++) { - xprime = xlow + (iDecayTimeBkgNeg - .5)*stepx; - sumxDecayTimeBkgNeg += FunBkgNeg(xprime); - xprime = xup - (iDecayTimeBkgNeg - .5)*stepx; - sumxDecayTimeBkgNeg += FunBkgNeg(xprime); - } - intxDecayTimeBkgNeg = sumxDecayTimeBkgNeg*stepx; - SetIntegralBkgNeg(intxDecayTimeBkgNeg); -// - Double_t iDecayTimeBkgSym; - Double_t intxDecayTimeBkgSym = 0.0; - Double_t sumxDecayTimeBkgSym = 0.0; - for(iDecayTimeBkgSym = 1.0; intxDecayTimeBkgSym <= np/2; intxDecayTimeBkgSym++) { - xprime = xlow + (intxDecayTimeBkgSym - .5)*stepx; - sumxDecayTimeBkgSym += FunBkgSym(xprime); - xprime = xup - (intxDecayTimeBkgSym - .5)*stepx; - sumxDecayTimeBkgSym += FunBkgSym(xprime); - } - intxDecayTimeBkgSym = sumxDecayTimeBkgSym*stepx; - SetIntegralBkgSym(intxDecayTimeBkgSym); - -// - Double_t iMassSig; - Double_t intmMassSig = 0.0; - Double_t summMassSig = 0.0; - for(iMassSig = 1.0; iMassSig<= npm/2.; iMassSig++) { - mx = fMassWndLow + (iMassSig - .5)*stepm; - summMassSig += EvaluateCDFInvMassSigDistr(mx); - mx = fMassWndHigh - (iMassSig - .5)*stepm; - summMassSig += EvaluateCDFInvMassSigDistr(mx); - } - intmMassSig = summMassSig*stepm; - SetIntegralMassSig(intmMassSig); - -// - Double_t iMassBkg; - Double_t intmMassBkg = 0.0; - Double_t summMassBkg = 0.0; - for(iMassBkg = 1.0; iMassBkg <= npm/2.; iMassBkg++) { - mx = fMassWndLow + (iMassBkg - .5)*stepm; - summMassBkg += EvaluateCDFInvMassBkgDistr(mx); - mx = fMassWndHigh - (iMassBkg - .5)*stepm; - summMassBkg += EvaluateCDFInvMassBkgDistr(mx); - } - intmMassBkg = summMassBkg*stepm; - SetIntegralMassBkg(intmMassBkg); - -// -// Compute integral of the whole distribution function -// - for(i = 1.0; i <= np; i++) { - Double_t summ = 0.0; - xprime = xlow + (i - .5)*stepx; - for(j = 1.0; j <= npm/2; j++) { - mx = fMassWndLow + (j - .5)*stepm; - summ += EvaluateCDFfunc(xprime,mx); - mx = fMassWndHigh - (j - .5)*stepm; - summ += EvaluateCDFfunc(xprime,mx); - } - intm = summ*stepm; - sumx += intm; - } - intx = sumx*stepx; - SetIntegral(intx); - -} -//_________________________________________________________________________________________________ -void AliBtoJPSItoEleCDFfitFCN::PrintStatus() -{ - // - // Print the parameters of the fits - // - printf("\n"); - printf("actual value of fRadius---------------------------------------->> | %f \n", GetRadius()); - printf("actual value of fTheta ---------------------------------------->> | %f \n", GetTheta()); - printf("actual value of fPhi ------------------------------------------>> | %f \n", GetPhi()); - printf("actual value of fFPlus ---------------------------------------->> | %f \n", GetFPlus()); - printf("actual value of fFMinus --------------------------------------->> | %f \n", GetFMinus()); - printf("actual value of fFSym ----------------------------------------->> | %f \n", GetFSym()); - printf("actual value of fOneOvLamPlus --------------------------------->> | %f \n", GetLamPlus()); - printf("actual value of fOneOvLamMinus -------------------------------->> | %f \n", GetLamMinus()); - printf("actual value of fOneOvLamSym ---------------------------------->> | %f \n", GetLamSym()); - printf("actual value of fMassBkgSlope --------------------------------->> | %f \n", GetMassSlope()); - printf("actual value of fFractionJpsiFromBeauty ----------------------->> | %f \n", GetFractionJpsiFromBeauty()); - printf("actual value of fFsig ----------------------------------------->> | %f \n", GetFsig()); - if(fCrystalBallParam){ - printf("actual value of fCrystalBallMmean ----------------------------->> | %f \n", GetCrystalBallMmean()); - printf("actual value of fCrystalBallNexp ------------------------------>> | %f \n", GetCrystalBallNexp()); - printf("actual value of fCrystalBallSigma ----------------------------->> | %f \n", GetCrystalBallSigma()); - printf("actual value of fCrystalBallAlpha ----------------------------->> | %f \n", GetCrystalBallAlpha()); - printf("actual value of fCrystalBallNorm ----------------------------->> | %f \n", GetCrystalBallNorm()); - }else{ - printf("actual value of fMpv ------------------------------------------>> | %f \n", GetCrystalBallMmean()); - printf("actual value of fConstRovL ------------------------------------>> | %f \n", GetCrystalBallNexp()); - printf("actual value of fSigmaL --------------------------------------->> | %f \n", GetCrystalBallSigma()); - printf("actual value of fSigmaR --------------------------------------->> | %f \n", GetCrystalBallAlpha()); - } - printf("actual value of fSigmaResol ----------------------------------->> | %f \n", GetSigmaResol()); - printf("actual value of fNResol --------------------------------------->> | %f \n", GetNResol()); - printf("\n"); - printf("Actual value of normalization integral for FunB ------------------->> | %f \n", GetIntegralFunB()); - printf("Actual value of normalization integral for BkgPos ----------------->> | %f \n", GetIntegralBkgPos()); - printf("Actual value of normalization integral for BkgNeg ----------------->> | %f \n", GetIntegralBkgNeg()); - printf("Actual value of normalization integral for BkgSym ----------------->> | %f \n", GetIntegralBkgSym()); - printf("Actual value of normalization integral for MassSig ---------------->> | %f \n", GetIntegralMassSig()); - printf("Actual value of normalization integral for MassBkg ---------------->> | %f \n", GetIntegralMassBkg()); - printf("Actual value of normalization integral for Resolution ------------->> | %f \n", GetIntegralRes()); - printf("Actual value of normalization integral for FCN -------------------->> | %f \n", GetIntegral()); - - printf("\n"); -} -//_________________________________________________________________________________________________ -void AliBtoJPSItoEleCDFfitFCN::SetResolutionConstants() -{ - // - // This method must be update: - // for the time beeing the values are hard-wired. - // Implementations have to be done to set the values from outside - // (e.g. from a ConfigHF file) starting from an indipendent fit - // of primary JPSI distribution. - // - - fResolutionConstants[0] = 8.; // mean sigma2/sigma1 - fResolutionConstants[1] = 0.1675; // mean Integral2/Integral1 - fResolutionConstants[2] = 1374.; // sigma2 - fResolutionConstants[3] = 0.001022; // N2 - fResolutionConstants[4] = 686.6; // mu2 -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfunc(Double_t x, Double_t m) const -{ - return fParameters[8]*EvaluateCDFfuncSignalPart(x,m) + (1. - fParameters[8])*EvaluateCDFfuncBkgPart(x,m); -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncNorm(Double_t x, Double_t m) const -{ - return EvaluateCDFfunc(x,m)/fIntegral; -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncSignalPart(Double_t x, Double_t m) const -{ - return EvaluateCDFDecayTimeSigDistr(x)*(EvaluateCDFInvMassSigDistr(m)/fintmMassSig); -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeSigDistr(Double_t x) const -{ - // - // Implementation of the Background part of the Likelyhood function - // - - Double_t retvalue = 0.; - Double_t FunBnorm = FunB(x)/fintxFunB; - Double_t FunPnorm = ResolutionFunc(x)/fintxRes; - retvalue = fParameters[7]*FunBnorm + (1. - fParameters[7])*FunPnorm; - return retvalue; -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassSigDistr(Double_t m) const -{ - // - // Parametrization of signal part invariant mass distribution - // It can be either Crystal Ball function or sum of two Landau - // - - Double_t fitval = 0.; - - if(fCrystalBallParam){ - Double_t t = (m-fParameters[9])/fParameters[11]; ; - if (fParameters[12] < 0) t = -t; - - Double_t absAlpha = TMath::Abs((Double_t)fParameters[12]); - - if (t >= -absAlpha) { - return fParameters[13]*TMath::Exp(-0.5*t*t); - } - else { - Double_t a = TMath::Power(fParameters[10]/absAlpha,fParameters[10])* TMath::Exp(-0.5*absAlpha*absAlpha); - Double_t b= fParameters[10]/absAlpha - absAlpha; - fitval = (fParameters[13]*a/TMath::Power(b - t, fParameters[10])); - return fitval; - } - }else{ - Double_t t=-1*m; - Double_t tmpv=-1*fParameters[9]; - fitval=TMath::Sqrt(TMath::Landau(t,tmpv,fParameters[11])); - fitval += fParameters[10]*(TMath::Landau(m,fParameters[9],fParameters[12])); - return fitval; - } -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::FunB(Double_t x) const -{ - // - // Parameterisation of the fit function for the x part of the Background - // - - Double_t np = 100.0; - Double_t sc = 10.; - Double_t sigma3 = fResolutionConstants[2]; - Double_t xprime; - Double_t sum = 0.0; - Double_t xlow,xupp; - Double_t step; - Double_t i; - xlow = x - sc * sigma3 ; - xupp = x + sc * sigma3 ; - step = (xupp-xlow) / np; - Double_t CsiMCxprime = 0.; - Double_t Resolutionxdiff = 0.; - Double_t xdiff = 0.; - - for(i=1.0; i<=np; i++) { - - xprime = xlow + (i-.5) * step; - CsiMCxprime = CsiMC(xprime); - xdiff = xprime - x; - Resolutionxdiff = ResolutionFunc(xdiff)/fintxRes; // normalized value - sum += CsiMCxprime * Resolutionxdiff; - - } - - return step * sum ; -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::FunP(Double_t x) const -{ - // - // Parameterisation of the Prompt part for the x distribution - // - - return ResolutionFunc(x)/fintxRes; // normalized value -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::CsiMC(Double_t x) const -{ - // - // Distribution (template) of the x distribution for the x variable - // for the J/psi coming from Beauty hadrons - // - - Double_t returnvalue = 0.; - returnvalue = fhCsiMC->GetBinContent(fhCsiMC->FindBin(x)); - return returnvalue; -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFfuncBkgPart(Double_t x,Double_t m) const -{ - // - // Return the part of the likelihood function for the background hypothesis - // - - return EvaluateCDFDecayTimeBkgDistr(x)*(EvaluateCDFInvMassBkgDistr(m)/fintmMassBkg); -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFDecayTimeBkgDistr(Double_t x) const -{ - // - // it returns the value of the probability to have a given x for the background - // - - Double_t ret = (1. - TMath::Power(fParameters[0],2.))*(ResolutionFunc(x)/fintxRes) - + TMath::Power(fParameters[0]*TMath::Cos(fParameters[1]),2.)* - (FunBkgPos(x)/fintxDecayTimeBkgPos) - + TMath::Power(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2]),2.)* - (FunBkgNeg(x)/fintxDecayTimeBkgNeg) - + TMath::Power(fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2]),2.)* - (FunBkgSym(x)/fintxDecayTimeBkgSym); - return ret; -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::EvaluateCDFInvMassBkgDistr(Double_t m) const -{ - // - // it returns the value of the probability to have a given mass for the background - // - - return 1/(fMassWndHigh-fMassWndLow) + - fParameters[6] * m - - fParameters[6] * ((fMassWndHigh+fMassWndLow)/2); -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgPos(Double_t x) const -{ - // - // exponential with positive slopes for the background part (x) - // - - Double_t np = 100.0; - Double_t sc = 10.; - Double_t sigma3 = fResolutionConstants[2]; - Double_t xprime; - Double_t sum = 0.0; - Double_t xlow,xupp; - Double_t step; - Double_t i; - xlow = x - sc * sigma3 ; - xupp = x + sc * sigma3 ; - step = (xupp-xlow) / np; - - for(i=1.0; i<=np/2; i++) { - - xprime = xlow + (i-.5) * step; - if (xprime > 0) {sum += fParameters[3] * TMath::Exp(-1*xprime*fParameters[3]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;} - - xprime = xupp - (i-.5) * step; - if (xprime > 0) {sum += fParameters[3] * TMath::Exp(-1*xprime*fParameters[3]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;} - - } - - return step * sum ; -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgNeg(Double_t x) const -{ - // - // exponential with negative slopes for the background part (x) - // - - Double_t np = 100.0; - Double_t sc = 10.; - Double_t sigma3 = fResolutionConstants[2]; - Double_t xprime; - Double_t sum = 0.0; - Double_t xlow,xupp; - Double_t step; - Double_t i; - xlow = x - sc * sigma3 ; - xupp = x + sc * sigma3 ; - step = (xupp-xlow) / np; - - for(i=1.0; i<=np/2; i++) { - - xprime = xlow + (i-.5) * step; - if (xprime < 0) {sum += fParameters[4] * TMath::Exp(xprime*fParameters[4]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;} - - xprime = xupp - (i-.5) * step; - if (xprime < 0) {sum += fParameters[4] * TMath::Exp(xprime*fParameters[4]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum = 0.;} - } - - return step * sum ; -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::FunBkgSym(Double_t x) const -{ - // - // exponential with both positive and negative slopes for the background part (x) - // - - Double_t np = 100.0; - Double_t sc = 10.; - Double_t sigma3 = fResolutionConstants[2]; - Double_t xprime; - Double_t sum1 = 0.0; - Double_t sum2 = 0.0; - Double_t xlow,xupp; - Double_t step; - Double_t i; - xlow = x - sc * sigma3 ; - xupp = x + sc * sigma3 ; - step = (xupp-xlow) / np; - - for(i=1.0; i<=np/2; i++) { - - xprime = xlow + (i-.5) * step; - if (xprime > 0) {sum1 += 0.5 * fParameters[5]*TMath::Exp(-1*xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum1 = 0.;} - - xprime = xupp - (i-.5) * step; - if (xprime > 0) {sum1 += 0.5 * fParameters[5]*TMath::Exp(-1*xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum1 = 0.;} - } - - for(i=1.0; i<=np/2; i++) { - - xprime = xlow + (i-.5) * step; - if (xprime < 0) {sum2 += 0.5 * fParameters[5]*TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum2 = 0.;} - - xprime = xupp - (i-.5) * step; - if (xprime < 0) {sum2 += 0.5 * fParameters[5]*TMath::Exp(xprime*fParameters[5]) * (ResolutionFunc(xprime-x)/fintxRes);} else {sum2 = 0.;} - } - - return step*(sum1 + sum2) ; -} -//_________________________________________________________________________________________________ -Double_t AliBtoJPSItoEleCDFfitFCN::ResolutionFunc(Double_t x) const -{ - // - // parametrization with 2 gaus - // - - Double_t ret = 0.; - Double_t x1 = x; - Double_t x2 = x; - //Double_t mean1 = 0.; - Double_t mean2 = fResolutionConstants[4]; - Double_t sigma1 = fParameters[14]; - Double_t sigma2 = fResolutionConstants[2]; - Double_t n1 = fParameters[15]; - Double_t n2 = fResolutionConstants[3]; - Double_t arg1 = x1/sigma1; - Double_t arg2 = (x2-mean2)/sigma2; - Double_t sqrt2Pi = TMath::Sqrt(2*TMath::Pi()); - - ret = n2*((n1/n2)*TMath::Exp(-0.5*arg1*arg1) + TMath::Exp(-0.5*arg2*arg2)); - - return ret/(sqrt2Pi*(n1*sigma1+n2*sigma2));//return value is normalized - -} - diff --git a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.h b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.h deleted file mode 100644 index 47daa98860d..00000000000 --- a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitFCN.h +++ /dev/null @@ -1,197 +0,0 @@ -#ifndef ALIBTOJPSITOELECDFFITFCN_H -#define ALIBTOJPSITOELECDFFITFCN_H -/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * - * See cxx source for full Copyright notice */ - -//_________________________________________________________________________ -// Class AliBtoJPSItoEleCDFfitFCN -// Definition of main function used in -// unbinned log-likelihood fit for -// the channel B -> JPsi + X -> e+e- + X -// -// Origin: C.Di Giglio -// Contact: Carmelo.Digiglio@ba.infn.it , Giuseppe.Bruno@ba.infn.it -//_________________________________________________________________________ - -#include -#include -#include "TH1F.h" -#include "TMath.h" -#include "TRandom3.h" - -enum IntegralType {kBkg, - kBkgNorm, - kSig, - kSigNorm}; - -enum PtBins {kallpt, - kptbin1,kptbin2,kptbin3, - kptbin4,kptbin5,kptbin6, - kptbin7,kptbin8,kptbin9}; -//_________________________________________________________________________________________________ -class AliBtoJPSItoEleCDFfitFCN : public TNamed { - public: - // - AliBtoJPSItoEleCDFfitFCN(); - AliBtoJPSItoEleCDFfitFCN(const AliBtoJPSItoEleCDFfitFCN& source); - AliBtoJPSItoEleCDFfitFCN& operator=(const AliBtoJPSItoEleCDFfitFCN& source); - virtual ~AliBtoJPSItoEleCDFfitFCN(); - - Double_t EvaluateLikelihood(const Double_t* pseudoproperdecaytime, - const Double_t* invariantmass, const Int_t ncand); - - Double_t GetFPlus() const { return fFPlus; } - Double_t GetFMinus() const { return fFMinus; } - Double_t GetFSym() const { return fFSym; } - Double_t GetRadius() const { return fParameters[0]; } - Double_t GetTheta() const { return fParameters[1]; } - Double_t GetPhi() const { return fParameters[2]; } - Double_t GetLamPlus() const { return fParameters[3]; } - Double_t GetLamMinus() const { return fParameters[4]; } - Double_t GetLamSym() const { return fParameters[5]; } - Double_t GetMassSlope() const { return fParameters[6]; } - Double_t GetFractionJpsiFromBeauty() const { return fParameters[7]; } - Double_t GetFsig() const { return fParameters[8]; } - Double_t GetCrystalBallMmean() const { return fParameters[9]; } - Double_t GetCrystalBallNexp() const { return fParameters[10]; } - Double_t GetCrystalBallSigma() const { return fParameters[11]; } - Double_t GetCrystalBallAlpha() const { return fParameters[12]; } - Double_t GetCrystalBallNorm() const { return fParameters[13]; } - Double_t GetSigmaResol() const { return fParameters[14]; } - Double_t GetNResol() const { return fParameters[15]; } - Double_t GetIntegral() const { return fIntegral; } - Double_t GetIntegralFunB() const { return fintxFunB; } - Double_t GetIntegralBkgPos() const { return fintxDecayTimeBkgPos; } - Double_t GetIntegralBkgNeg() const { return fintxDecayTimeBkgNeg; } - Double_t GetIntegralBkgSym() const { return fintxDecayTimeBkgSym; } - Double_t GetIntegralMassSig() const { return fintmMassSig; } - Double_t GetIntegralRes() const { return fintxRes; } - Double_t GetIntegralMassBkg() const { return fintmMassBkg; } - Bool_t GetCrystalBallParam() const { return fCrystalBallParam; } - - void SetFPlus(Double_t plus) {fFPlus = plus;} - void SetFMinus(Double_t minus) {fFMinus = minus;} - void SetFSym(Double_t sym) {fFSym = sym;} - void SetRadius(Double_t radius) {fParameters[0] = radius;} - void SetTheta(Double_t theta) {fParameters[1] = theta;} - void SetPhi(Double_t phi) {fParameters[2] = phi;} - void SetLamPlus(Double_t lamplus) {fParameters[3] = lamplus;} - void SetLamMinus(Double_t lamminus) {fParameters[4] = lamminus;} - void SetLamSym(Double_t lamsym) {fParameters[5] = lamsym;} - void SetMassSlope(Double_t slope) {fParameters[6] = slope;} - void SetFractionJpsiFromBeauty(Double_t B) {fParameters[7] = B;} - void SetFsig(Double_t Fsig) {fParameters[8] = Fsig;} - void SetCrystalBallMmean(Double_t CrystalBallMmean) {fParameters[9] = CrystalBallMmean;} - void SetCrystalBallNexp(Double_t CrystalBallNexp) {fParameters[10] = CrystalBallNexp;} - void SetCrystalBallSigma(Double_t CrystalBallSigma) {fParameters[11] = CrystalBallSigma;} - void SetCrystalBallAlpha(Double_t CrystalBallAlpha) {fParameters[12] = CrystalBallAlpha;} - void SetCrystalBallNorm(Double_t CrystalBallNorm) {fParameters[13] = CrystalBallNorm;} - void SetSigmaResol(Double_t SigmaResol) {fParameters[14] = SigmaResol;} - void SetNResol(Double_t NResol) {fParameters[15] = NResol;} - - void SetAllParameters(const Double_t* parameters); - void SetIntegral(Double_t integral) { fIntegral = integral; } - void SetIntegralFunB(Double_t integral) { fintxFunB = integral; } - void SetIntegralBkgPos(Double_t integral) { fintxDecayTimeBkgPos = integral; } - void SetIntegralBkgNeg(Double_t integral) { fintxDecayTimeBkgNeg = integral; } - void SetIntegralBkgSym(Double_t integral) { fintxDecayTimeBkgSym = integral; } - void SetIntegralMassSig(Double_t integral) { fintmMassSig = integral; } - void SetIntegralRes(Double_t integral) { fintxRes = integral; } - void SetIntegralMassBkg(Double_t integral) { fintmMassBkg = integral; } - - void SetCsiMC(const TH1F* MCtemplate) {fhCsiMC = (TH1F*)MCtemplate->Clone("fhCsiMC");} - - void SetResolutionConstants(); - void SetMassWndHigh(Double_t limit) { fMassWndHigh = TDatabasePDG::Instance()->GetParticle(443)->Mass() + limit ;} - void SetMassWndLow(Double_t limit) { fMassWndLow = TDatabasePDG::Instance()->GetParticle(443)->Mass() - limit ;} - void SetCrystalBallFunction(Bool_t okCB) {fCrystalBallParam = okCB;} - - void ConvertFromSpherical() { fFPlus = TMath::Power((fParameters[0]*TMath::Cos(fParameters[1])),2.); - fFMinus = TMath::Power((fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2])),2.); - fFSym = TMath::Power((fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2])),2.);} - - void ComputeIntegral(); - - void ReadMCtemplates(Int_t BinNum); - - void PrintStatus(); - - private: - // - Double_t fParameters[16]; /* par[0] = fRadius; - par[1] = fTheta; - par[2] = fPhi; - par[3] = fOneOvLamPlus; - par[4] = fOneOvLamMinus; - par[5] = fOneOvLamSym; - par[6] = fMassBkgSlope; - par[7] = fFractionJpsiFromBeauty; - par[8] = fFsig; - par[9] = fCrystalBallMmean; - par[10] = fCrystalBallNexp; - par[11] = fCrystalBallSigma; - par[12] = fCrystalBallAlpha; - par[13] = fCrystalBallNorm; - par[14] = fSigmaResol; - par[15] = fNResol; */ - - - Double_t fFPlus; // parameters of the log-likelihood function - Double_t fFMinus; // Slopes of the x distributions of the background - Double_t fFSym; // functions - - Double_t fIntegral; // integral values of log-likelihood function terms - Double_t fintxFunB; - Double_t fintxDecayTimeBkgPos; - Double_t fintxDecayTimeBkgNeg; - Double_t fintxDecayTimeBkgSym; - Double_t fintmMassSig; - Double_t fintxRes; - Double_t fintmMassBkg; - TH1F *fhCsiMC; // X distribution used as MC template for JPSI from B - Double_t fResolutionConstants[5]; // constants for the parametrized resolution function R(X) - Double_t fMassWndHigh; // JPSI Mass window higher limit - Double_t fMassWndLow; // JPSI Mass window lower limit - Bool_t fCrystalBallParam; // Boolean to switch to Crystall Ball parameterisation - - //// - - Double_t EvaluateCDFfunc(Double_t x, Double_t m) const ; - - Double_t EvaluateCDFfuncNorm(Double_t x, Double_t m) const ; - - //// - - Double_t EvaluateCDFfuncSignalPart(Double_t x, Double_t m) const ; // Signal part - - Double_t EvaluateCDFDecayTimeSigDistr(Double_t x) const ; - - Double_t EvaluateCDFInvMassSigDistr(Double_t m) const ; - - Double_t EvaluateCDFfuncBkgPart(Double_t x,Double_t m) const ; // Background part - - Double_t EvaluateCDFDecayTimeBkgDistr(Double_t x) const ; - - Double_t EvaluateCDFInvMassBkgDistr(Double_t m) const ; - - //// - - Double_t FunB(Double_t x) const ; - - Double_t FunP(Double_t x) const ; - - Double_t CsiMC(Double_t x) const ; - - Double_t FunBkgPos(Double_t x) const ; - - Double_t FunBkgNeg(Double_t x) const ; - - Double_t FunBkgSym(Double_t x) const ; - - Double_t ResolutionFunc(Double_t x) const ; - // - ClassDef (AliBtoJPSItoEleCDFfitFCN,1); // Unbinned log-likelihood fit - -}; - -#endif diff --git a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.cxx b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.cxx deleted file mode 100644 index 6c7851444c5..00000000000 --- a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.cxx +++ /dev/null @@ -1,238 +0,0 @@ -/************************************************************************** - * Copyright(c) 1998-2009, 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. * - **************************************************************************/ -#include -#include -#include "AliLog.h" -#include "AliBtoJPSItoEleCDFfitHandler.h" -#include "AliBtoJPSItoEleCDFfitFCN.h" - -//------------------------------------------------------------------------- -// Class AliBtoJPSItoEleCDFfitHandler -// Class to perform unbinned log-likelihood fit -// -// Origin: C. Di Giglio -// Contact: carmelo.digiglio@ba.infn.it , giuseppe.bruno@ba.infn.it -//------------------------------------------------------------------------- - -void CDFFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag); - -ClassImp(AliBtoJPSItoEleCDFfitHandler) - -//______________________________________________________________________________ -void CDFFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) -{ - // This function is called by minuit - // The corresponding member method is called - // using SetObjectFit/GetObjectFit methods of TMinuit - AliBtoJPSItoEleCDFfitHandler* dummy = (AliBtoJPSItoEleCDFfitHandler *)TVirtualFitter::GetFitter()->GetObjectFit(); - dummy->CdfFCN(npar, gin, f, par, iflag); -} - - -//_________________________________________________________________________________________________ -AliBtoJPSItoEleCDFfitHandler::AliBtoJPSItoEleCDFfitHandler(): -fIsParamFixed(16), -fPrintStatus(kFALSE), -fUp(0), -fX(0x0), -fM(0x0), -fLikely(0x0), -fNcand(0) -{ - // - // default constructor - // -} -//_________________________________________________________________________________________________ -AliBtoJPSItoEleCDFfitHandler::AliBtoJPSItoEleCDFfitHandler(Double_t* decaytime, - Double_t* invariantmass, Int_t ncand) : -fIsParamFixed(16), -fPrintStatus(kFALSE), -fUp(0), -fX(decaytime), -fM(invariantmass), -fLikely(0x0), -fNcand(ncand) -{ - // - // constructor - // - AliInfo("\n+++\n+++ Minimization object AliBtoJPSItoEleCDFfitHandler created\n+++\n"); - fLikely = new AliBtoJPSItoEleCDFfitFCN(); - AliInfo("\n+++\n+++ CDF fit function object AliBtoJPSItoEleCDFfitFCN created\n+++\n"); - AliInfo("Parameter 0 ----> fRadius"); - AliInfo("Parameter 1 ----> fTheta"); - AliInfo("Parameter 2 ----> fPhi"); - AliInfo("Parameter 3 ----> fOneOvLamPlus"); - AliInfo("Parameter 4 ----> fOneOvLamMinus"); - AliInfo("Parameter 5 ----> fOneOvLamSym"); - AliInfo("Parameter 6 ----> fMSlope"); - AliInfo("Parameter 7 ----> fB"); - AliInfo("Parameter 8 ----> fFsig"); - AliInfo("Parameter 9 ----> fMmean"); - AliInfo("Parameter 10 ----> fNexp"); - AliInfo("Parameter 11 ----> fSigma"); - AliInfo("Parameter 12 ----> fAlpha"); - AliInfo("Parameter 13 ----> fNorm"); - AliInfo("Parameter 14 ----> fSigmaResol"); - AliInfo("Parameter 15 ----> fNResol"); - AliInfo(Form("\n+++\n+++ Number of candidates ---> %d\n+++\n ", ncand)); -} -//___________________________________________________________________________ -AliBtoJPSItoEleCDFfitHandler& AliBtoJPSItoEleCDFfitHandler::operator=(const AliBtoJPSItoEleCDFfitHandler& c) -{ - // - // Assignment operator - // - if (this!=&c) { - fIsParamFixed = c.fIsParamFixed; - fPrintStatus = c.fPrintStatus; - fUp = c.fUp; - fX = c.fX; - fM = c.fM; - fLikely = c.fLikely; - fNcand = c.fNcand; - } - return *this; -} - -//_______________________________________________________________________________________ -AliBtoJPSItoEleCDFfitHandler::AliBtoJPSItoEleCDFfitHandler(const AliBtoJPSItoEleCDFfitHandler& c) : -TNamed(c), -fIsParamFixed(c.fIsParamFixed), -fPrintStatus(c.fPrintStatus), -fUp(c.fUp), -fX(c.fX), -fM(c.fM), -fLikely(c.fLikely), -fNcand(c.fNcand) -{ - // - // Copy Constructor - // -} -//_______________________________________________________________________________________ -AliBtoJPSItoEleCDFfitHandler::~AliBtoJPSItoEleCDFfitHandler() -{ - // - //destructor - // - delete fLikely; -} -//_______________________________________________________________________________________ -Int_t AliBtoJPSItoEleCDFfitHandler::DoMinimization() -{ - // - // performs the minimization - // - static TVirtualFitter *fitter = TVirtualFitter::Fitter(this,16); - fitter->SetFCN(CDFFunction); - - fitter->SetParameter(0,"fRadius",fParamStartValues[0], 1.e-06, 0., 1.); - fitter->SetParameter(1,"fTheta",fParamStartValues[1], 1.e-06, 0.,2*TMath::Pi()); - fitter->SetParameter(2,"fPhi",fParamStartValues[2], 1.e-06, 0.,2*TMath::Pi()); -// fitter->SetParameter(3,"fOneOvLamPlus",fParamStartValues[3], 1.e-10, 0., 5.e+01); -// fitter->SetParameter(4,"fOneOvLamMinus",fParamStartValues[4], 1.e-10, 0., 5.e+01); -// fitter->SetParameter(5,"fOneOvLamSym",fParamStartValues[5], 1.e-10, 0., 5.e+01); - fitter->SetParameter(3,"fOneOvLamPlus",fParamStartValues[3], 1.e-10, 0.0000001, 5.e+01); - fitter->SetParameter(4,"fOneOvLamMinus",fParamStartValues[4], 1.e-10, 0.00000001, 5.e+01); - fitter->SetParameter(5,"fOneOvLamSym",fParamStartValues[5], 1.e-10, 0.00000001, 5.e+01); - fitter->SetParameter(6,"fMSlope",fParamStartValues[6], 1.e-04, -2.5, 2.5); - fitter->SetParameter(7,"fB",fParamStartValues[7], 1.e-08, 0., 1.); - fitter->SetParameter(8,"fFsig",fParamStartValues[8], 1.e-08, 0., 1.); - fitter->SetParameter(9,"fMmean",fParamStartValues[9], 1.e-08, 0., 1.e+04); - fitter->SetParameter(10,"fNexp",fParamStartValues[10], 1.e-08, 0., 1.e+02); - fitter->SetParameter(11,"fSigma",fParamStartValues[11], 1.e-08, 0., 1.e+04); - fitter->SetParameter(12,"fAlpha",fParamStartValues[12], 1.e-08, 0., 1.e+04); - fitter->SetParameter(13,"fNorm",fParamStartValues[13], 1.e-08, 0., 1.e+01); - fitter->SetParameter(14,"fSigmaResol",fParamStartValues[14], 1.e-08, 0., 1.e+04); - fitter->SetParameter(15,"fNResol",fParamStartValues[15], 1.e-08, 0., 1.e+05); - - for(UInt_t indexparam = 0; indexparam < 16; indexparam++){ - if(IsParamFixed(indexparam))fitter->FixParameter((Int_t)indexparam); - } - - Double_t arglist[2]={10000,1.0}; - Int_t iret=fitter->ExecuteCommand("MIGRAD", arglist ,2); - fitter->PrintResults(4,0); - - AliInfo("Minimization procedure finished\n"); - - return iret; -} -//_______________________________________________________________________________________ -void AliBtoJPSItoEleCDFfitHandler::CdfFCN(Int_t & /* npar */, - Double_t * /* gin */,Double_t &f,Double_t *par,Int_t /* iflag */) -{ -// -// Definition of the FCN to be used by minuit -// - fLikely->SetAllParameters(par); - fLikely->ConvertFromSpherical(); - fLikely->ComputeIntegral(); - if(fPrintStatus)fLikely->PrintStatus(); - - TStopwatch t; - t.Start(); - - f = fLikely->EvaluateLikelihood(fX,fM,fNcand); - - t.Stop(); - AliDebug(2,Form("Real time spent to calculate function == %f \n", t.RealTime())); - AliDebug(2,Form("CPU time spent to calculate function == %f \n", t.CpuTime())); - AliDebug(2,Form("Actual value of the AliBtoJPSItoEleCDFfitFCN function == %f \n", f)); - - return; -} -//_______________________________________________________________________________________ -void AliBtoJPSItoEleCDFfitHandler::SetParamStartValues(Double_t inputparamvalues[16]) -{ - for(Int_t index=0; index < 16; index++) fParamStartValues[index] = inputparamvalues[index]; -} -//_______________________________________________________________________________________ -void AliBtoJPSItoEleCDFfitHandler::SetResolutionConstants() -{ - // - // Sets constants for the resolution function - // - fLikely->SetResolutionConstants(); - -} -//_______________________________________________________________________________________ -void AliBtoJPSItoEleCDFfitHandler::SetCrystalBallFunction(Bool_t okCB) -{ - // - // Sets the CB as the parametrization for the signal invariant mass spectrum - // (otherwise Landau is chosen) - // - fLikely->SetCrystalBallFunction(okCB); -} -//_______________________________________________________________________________________ -void AliBtoJPSItoEleCDFfitHandler::SetMassWndHigh(Double_t limit) -{ - // - // Sets upper limit for the invariant mass window (under J/PSI mass peak) - // - fLikely->SetMassWndHigh(limit); -} -//_______________________________________________________________________________________ -void AliBtoJPSItoEleCDFfitHandler::SetMassWndLow(Double_t limit) -{ - // - // Sets lower limit for the invariant mass window (under J/PSI mass peak) - // - fLikely->SetMassWndLow(limit); -} - diff --git a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.h b/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.h deleted file mode 100644 index 9c55f255e4a..00000000000 --- a/PWG3/vertexingHF/AliBtoJPSItoEleCDFfitHandler.h +++ /dev/null @@ -1,68 +0,0 @@ -#ifndef ALIBTOJPSITOELECDFFITHANDLER_H -#define ALIBTOJPSITOELECDFFITHANDLER_H -/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * - * See cxx source for full Copyright notice */ - -//___________________________________________________________________________ -// Class AliBtoJPSItoEleCDFfitHandler -// Class to perform unbinned log-likelihood fit -// -// Origin: C. Di Giglio -// Contact: Carmelo.Digiglio@ba.infn.it; Giuseppe.Bruno@ba.infn.it -//____________________________________________________________________________ - -#include -#include -#include "TMath.h" -#include "AliBtoJPSItoEleCDFfitFCN.h" -#include "AliLog.h" - -class TBits; - -class AliBtoJPSItoEleCDFfitHandler : public TNamed { - public: - // - AliBtoJPSItoEleCDFfitHandler(); - AliBtoJPSItoEleCDFfitHandler& operator= (const AliBtoJPSItoEleCDFfitHandler& c); - AliBtoJPSItoEleCDFfitHandler(const AliBtoJPSItoEleCDFfitHandler& c); - AliBtoJPSItoEleCDFfitHandler(Double_t* decaytime, Double_t* invariantmass, Int_t ncand); - ~AliBtoJPSItoEleCDFfitHandler(); - //Double_t Up() const { return fUp*TMath::Sqrt(fLikely->GetIntegral()); } - Double_t Up() const { return fUp; } - void SetErrorDef(Double_t up) {fUp = up;} - void SetPrintStatus(Bool_t okPrintStatus) { fPrintStatus = okPrintStatus; } - Bool_t GetPrintStatus() { return fPrintStatus ; } - void SetParamStartValues (Double_t*); - Double_t* GetStartParamValues() { return fParamStartValues; } - TBits GetFixedParamList() { return fIsParamFixed; } - void FixParam(UInt_t param, Bool_t value) { fIsParamFixed.SetBitNumber(param,value); } - void FixAllParam(Bool_t value) { for(UInt_t par=0;par<16;par++) fIsParamFixed.SetBitNumber(par,value); } - Bool_t IsParamFixed(UInt_t param) { return fIsParamFixed.TestBitNumber(param); } - void SetResolutionConstants(); - void SetCrystalBallFunction(Bool_t okCB); - void SetMassWndHigh(Double_t limit); - void SetMassWndLow(Double_t limit); - - Double_t operator()(const Double_t* par) const ; - void CdfFCN(Int_t & /* npar */, Double_t * /* gin */, Double_t &f, Double_t *par, Int_t /* iflag */); - - Double_t* Decaytime() const { return fX; } - Double_t* InvariantMass() const { return fM; } - AliBtoJPSItoEleCDFfitFCN* LikelihoodPointer() const { return fLikely; } - Int_t DoMinimization(); - - private: - // - TBits fIsParamFixed; //array of bits: 0 = param free; 1 = param fixed; - Bool_t fPrintStatus; //flag to enable the prit out of the algorithm at each step - Double_t fParamStartValues[16]; //array of parameters input value - Double_t fUp; //error definition - Double_t* fX; //pseudo-proper decay time X - Double_t* fM; //invariant mass M - AliBtoJPSItoEleCDFfitFCN* fLikely; //Log likelihood function - Int_t fNcand; //number of candidates - // - ClassDef(AliBtoJPSItoEleCDFfitHandler,1); - -}; -#endif diff --git a/PWG3/vertexingHF/macros/AddTaskBkgLikeSignJPSI.C b/PWG3/vertexingHF/macros/AddTaskBkgLikeSignJPSI.C deleted file mode 100644 index 4891eed20f3..00000000000 --- a/PWG3/vertexingHF/macros/AddTaskBkgLikeSignJPSI.C +++ /dev/null @@ -1,40 +0,0 @@ -AliAnalysisTaskSEBkgLikeSignJPSI *AddTaskBkgLikeSignJPSI() -{ - // - // Test macro for the AliAnalysisTaskSEBkgLikeSignJPSI - // starting from AliAOD.root file with HF + Like Sign candidates. - // C.Di Giglio, carmelo.digiglio@ba.infn.it - // - - - // Get the pointer to the existing analysis manager via the static access method. - //============================================================================== - AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); - if (!mgr) { - ::Error("AddTaskBkgLikeSignJPSI", "No analysis manager to connect to."); - return NULL; - } - - - // Like-sign background analysis task - AliAnalysisTaskSEBkgLikeSignJPSI *lsTask = new AliAnalysisTaskSEBkgLikeSignJPSI("CmpLikeSignAnalysis"); - lsTask->SetDebugLevel(0); - - mgr->AddTask(lsTask); - - // - // Create containers for input/output - AliAnalysisDataContainer *cinputLS = mgr->CreateContainer("cinputLikeSignJPSI",TChain::Class(), - AliAnalysisManager::kInputContainer); - TString outputfile = AliAnalysisManager::GetCommonFileName(); - outputfile += ":PWG3_D2H_CmpLikesignJPSI"; - AliAnalysisDataContainer *coutputLS = mgr->CreateContainer("coutputLikeSignJPSI",TList::Class(), - AliAnalysisManager::kOutputContainer, - outputfile.Data()); - - mgr->ConnectInput(lsTask,0,mgr->GetCommonInputContainer()); - mgr->ConnectOutput(lsTask,1,coutputLS); - - - return lsTask; -} diff --git a/PWG3/vertexingHF/macros/AddTaskBtoJPSItoEle.C b/PWG3/vertexingHF/macros/AddTaskBtoJPSItoEle.C deleted file mode 100644 index d475897d658..00000000000 --- a/PWG3/vertexingHF/macros/AddTaskBtoJPSItoEle.C +++ /dev/null @@ -1,33 +0,0 @@ -AliAnalysisTaskSEBtoJPSItoEle *AddTaskBtoJPSItoEle() -{ - // - // Test macro for the AliAnalysisTaskSEBtoJPSItoEle - // starting from AliAOD.root file with HF + Like Sign candidates. - // C.Di Giglio, carmelo.digiglio@ba.infn.it - // - - // Get the pointer to the existing analysis manager via the static access method. - //============================================================================== - AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); - if (!mgr) { - ::Error("AddTaskBtoJPSItoEle", "No analysis manager to connect to."); - return NULL; - } - - // Cdf unbinned log-likelihood fit analysis task - AliAnalysisTaskSEBtoJPSItoEle *hfTask = new AliAnalysisTaskSEBtoJPSItoEle("CdfFitAnalysis"); - hfTask->SetDebugLevel(0); - - mgr->AddTask(hfTask); - - // - // Create containers for input/output - mgr->ConnectInput(hfTask,0,mgr->GetCommonInputContainer()); - - AliAnalysisDataContainer *coutput = mgr->CreateContainer("coutputCdfFit",TList::Class(), - AliAnalysisManager::kOutputContainer, - "CdfFit.root"); - mgr->ConnectOutput(hfTask,1,coutput); - - return hfTask; -} diff --git a/PWG3/vertexingHF/macros/AddTaskJPSItoEle.C b/PWG3/vertexingHF/macros/AddTaskJPSItoEle.C deleted file mode 100644 index 0205cd3b574..00000000000 --- a/PWG3/vertexingHF/macros/AddTaskJPSItoEle.C +++ /dev/null @@ -1,47 +0,0 @@ -AliAnalysisTaskSEJPSItoEle *AddTaskJPSItoEle(char* fileout = "AliAOD.Jpsitoele.root") -{ - //*********************************************************************************************** - // Test macro for the AliAnalysisTaskSEJPSItoEle - // Starting from AliAOD.root + AliAOD.VertexingHF.root with HF + like sign candidates, - // it produces a specific AliAODjpsi.root file with a replica of J/psi->e+e- candidates only - // (and references to the corresponding decay tracks) + several histograms for - // both unlike sign and like sign candidates. - // - // C.Di Giglio, carmelo.digiglio@ba.infn.it - //*********************************************************************************************** - - // Get the pointer to the existing analysis manager via the static access method. - //============================================================================== - AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); - if (!mgr) { - ::Error("AddTaskBtoJPSItoEle", "No analysis manager to connect to."); - return NULL; - } - - AliAODHandler* aodHandler = new AliAODHandler(); - aodHandler->SetOutputFileName(fileout); - - mgr->SetOutputEventHandler(aodHandler); - mgr-> SetDebugLevel(10); - - AliAnalysisTaskSEJPSItoEle *hJPSItoEleTask = new AliAnalysisTaskSEJPSItoEle("AOD_JPSItoEle_filter"); - hJPSItoEleTask->SetDebugLevel(2); - - Double_t ptCuts[2] = {0.,500.}; // the cut is on the J/psi pT (--> change this) - hJPSItoEleTask->SetPtCuts(ptCuts); - hJPSItoEleTask->SetAODMCInfo(kFALSE); // only for sim - - mgr->AddTask(hJPSItoEleTask); - - mgr->ConnectInput(hJPSItoEleTask,0,mgr->GetCommonInputContainer()); - - AliAnalysisDataContainer *coutput0 = mgr->CreateContainer("tree",TTree::Class(), - AliAnalysisManager::kOutputContainer, - "default"); - AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("histos",TList::Class(), - AliAnalysisManager::kOutputContainer, - "JPSItoEleHistos.root"); - mgr->ConnectOutput(hJPSItoEleTask,0,coutput0); - mgr->ConnectOutput(hJPSItoEleTask,1,coutput1); - return hJPSItoEleTask; -} diff --git a/PWG3/vertexingHF/macros/ReadAODVertexingHF.C b/PWG3/vertexingHF/macros/ReadAODVertexingHF.C index 219e15762b8..be79b266698 100644 --- a/PWG3/vertexingHF/macros/ReadAODVertexingHF.C +++ b/PWG3/vertexingHF/macros/ReadAODVertexingHF.C @@ -130,6 +130,15 @@ void ReadAODVertexingHF(const char *aodFileName="AliAOD.root", } + // Fix references to daughter tracks + //AliAnalysisVertexingHF *fixer = new AliAnalysisVertexingHF(); + //fixer->FixReferences(aod); + //delete fixer; + // + //AliRDHFCutsD0toKpi *mycuts=new AliRDHFCutsD0toKpi(); + //mycuts->SetFixRefs(kTRUE); + //mycuts->IsEventSelected(aod); + // loop over D0->Kpi candidates Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast(); nTotD0toKpi += nD0toKpi; @@ -138,42 +147,11 @@ void ReadAODVertexingHF(const char *aodFileName="AliAOD.root", for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) { AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->UncheckedAt(iD0toKpi); - d->GetPrimaryVtx()->Print(); - - Bool_t unsetvtx=kFALSE; - //if(!d->GetOwnPrimaryVtx()) { - //d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables - //unsetvtx=kTRUE; - //} - Int_t okD0=0,okD0bar=0; - if(d->SelectD0(cutsD0,okD0,okD0bar)) { - //cout<<1e8*d->Prodd0d0()<Fill(d->InvMassD0(),0.5); - hMass->Fill(d->InvMassD0bar(),0.5); - hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle()); - hSecVtxZ->Fill(d->GetSecVtxZ()); - //cout<GetSecVtxX() <GetDaughter(0); AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1); - if(!trk0 || !trk1) { - trk0=aod->GetTrack(trkIDtoEntry[d->GetProngID(0)]); - trk1=aod->GetTrack(trkIDtoEntry[d->GetProngID(1)]); - cout<<"references to standard AOD not available"<Pt()<Pt()<<" (AliAODRecoDecay); "<GetMagneticField(),1000.,dz,covdz); - //cout<<"D0 impact parameter rphi: "<UnsetOwnPrimaryVtx(); + + if(trk0->GetStatus()) printf("ok %d\n",iD0toKpi); } @@ -182,27 +160,42 @@ void ReadAODVertexingHF(const char *aodFileName="AliAOD.root", nTotDstar += nDstar; cout<<"Number of D*->D0pi: "<SetRequireSigmaToVertex(kFALSE); + //default + esdTrackCuts->SetRequireTPCRefit(kTRUE); + //esdTrackCuts->SetRequireITSRefit(kTRUE); + //esdTrackCuts->SetMinNClustersTPC(70); + esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, + AliESDtrackCuts::kAny); + // default is kBoth, otherwise kAny + esdTrackCuts->SetMinDCAToVertexXY(0.); + esdTrackCuts->SetPtRange(0.3,1.e10); + + // soft pion pre-selections + AliESDtrackCuts* esdSoftPicuts=new AliESDtrackCuts(); + //esdSoftPicuts->SetRequireSigmaToVertex(kFALSE); + //default + //esdSoftPicuts->SetRequireTPCRefit(kFALSE); + //esdSoftPicuts->SetRequireITSRefit(kFALSE); + //esdSoftPicuts->SetMinNClustersITS(4); // default is 5 + //esdSoftPicuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, + // AliESDtrackCuts::kAny); //test d0 asimmetry + //esdSoftPicuts->SetPtRange(0.0,15.0); + + cuts->AddTrackCuts(esdTrackCuts); + cuts->AddTrackCutsSoftPi(esdSoftPicuts); + Float_t cutsArrayDStartoKpipi[14]={0.15,0.07.,0.85,0.8,0.8,0.06.,0.06.,0.001,0.6,0.15, 0.03, 0.2, 5, 0.5}; // first 9 for D0 from D*, last 5 for D* + cuts->SetCuts(14,cutsArrayDStartoKpipi); + for (Int_t iDstar = 0; iDstar < nDstar; iDstar++) { AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)arrayDstar->UncheckedAt(iDstar); - Bool_t unsetvtx=kFALSE; - //if(!c->GetOwnPrimaryVtx()) { - //c->SetOwnPrimaryVtx(vtx1); // needed to compute all variables - //c->Get2Prong()->SetOwnPrimaryVtx(vtx1); - //unsetvtx=kTRUE; - //} - if(c->SelectDstar(cutsDstar,cutsD0)) { - hDeltaMassDstar->Fill(c->DeltaInvMass()); - // get daughters - AliAODTrack *trk = (AliAODTrack*)c->GetBachelor(); - AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)c->Get2Prong(); - //cout<<"pt of soft pi: "<Pt()<Pt()<UnsetOwnPrimaryVtx(); - c->Get2Prong()->UnsetOwnPrimaryVtx(); - } + if(cuts->IsSelected(c,AliRDHFCuts::kTracks)) printf("passed\n"); + } @@ -215,7 +208,7 @@ void ReadAODVertexingHF(const char *aodFileName="AliAOD.root", Int_t nVtxsHF = arrayVerticesHF->GetEntriesFast(); nTotHF += nVtxsHF; cout<<"Number of heavy-flavour vertices: "<UncheckedAt(iVtx); // print info //cout << iVtx << ": vertex z position: " << vtxHF->GetZ() << endl;