+++ /dev/null
-/* 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;
-
- }
-}
+++ /dev/null
-#ifndef ALIANALYSISBTOJPSITOELE_H\r
-#define ALIANALYSISBTOJPSITOELE_H\r
-/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
- * See cxx source for full Copyright notice */\r
-\r
-//-------------------------------------------------------------------------\r
-// Class AliAnalysisBtoJPSItoEle\r
-// Unbinned log-likelihood fit analysis class\r
-//\r
-// Origin: C.Di Giglio\r
-// Contact: Carmelo.Digiglio@ba.infn.it , giuseppe.bruno@ba.infn.it\r
-//-------------------------------------------------------------------------\r
-class TNtuple ;\r
-//class AliBtoJPSItoEleCDFfitFCN ;\r
-//class AliBtoJPSItoEleCDFfitHandler ; \r
-#include "TH1F.h"\r
-#include "AliBtoJPSItoEleCDFfitHandler.h"\r
-#include "AliBtoJPSItoEleCDFfitFCN.h"\r
-\r
-//_________________________________________________________________________ \r
-class AliAnalysisBtoJPSItoEle : public TNamed {\r
- public:\r
- //\r
- AliAnalysisBtoJPSItoEle();\r
- AliAnalysisBtoJPSItoEle(const AliAnalysisBtoJPSItoEle& source);\r
- AliAnalysisBtoJPSItoEle& operator=(const AliAnalysisBtoJPSItoEle& source);\r
- virtual ~AliAnalysisBtoJPSItoEle();\r
- \r
- Int_t DoMinimization();\r
- void ReadCandidates(TNtuple* nt, Double_t* &x, Double_t* &m, Int_t& n); // primary JPSI + secondary JPSI + bkg couples\r
-\r
- void SetPtBin(Int_t BinNum) { fPtBin = BinNum ; }\r
- void SetCsiMC();\r
- //void SetResolutionConstants(Int_t BinNum);\r
- void SetFitHandler(Double_t* /*pseudoproper*/, Double_t* /*inv mass*/, Int_t /*candidates*/); \r
- void CloneMCtemplate(const TH1F* MCtemplate) {fMCtemplate = (TH1F*)MCtemplate->Clone("fMCtemplate");}\r
- Double_t* GetResolutionConstants(Double_t* resolutionConst);\r
- AliBtoJPSItoEleCDFfitHandler* GetCDFFitHandler() const { return fFCNfunction ; }\r
- Int_t GetPtBin() const { return fPtBin ; }\r
-\r
- private:\r
- //\r
- AliBtoJPSItoEleCDFfitHandler* fFCNfunction; //! pointer to the interface class\r
- Int_t fPtBin; // number of pt bin in which the analysis is performes\r
- TH1F* fMCtemplate; //! template of the MC distribution for the x distribution of the secondary J/psi\r
-\r
- ClassDef(AliAnalysisBtoJPSItoEle,1); // AliAnalysisBtoJPSItoEle class\r
-};\r
-\r
-#endif\r
+++ /dev/null
-/**************************************************************************
- * 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 <TSystem.h>
-#include <TROOT.h>
-#include <TClonesArray.h>
-#include <TNtuple.h>
-#include <TList.h>
-#include <TH1F.h>
-
-#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<AliAODEvent*> (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<AliAODEvent*> (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;it<aod->GetNumberOfTracks();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<TList*> (GetOutputData(1));
- if (!fOutput) {
- printf("ERROR: fOutput not available\n");
- return;
- }
-
- fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs*fTotNegPairs);
-
- fHistMassJPSI = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMassJPSI"));
- fHistMassLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMassLS"));
- fHistCtsJPSI = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsJPSI"));
- fHistCtsLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLS"));
- fHistCtsLSpos = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSpos"));
- fHistCtsLSneg = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSneg"));
- fHistCPtaJPSI = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCPtaJPSI"));
- fHistCPtaLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCPtaLS"));
- fHistd0d0JPSI = dynamic_cast<TH1F*>(fOutput->FindObject("fHistd0d0JPSI"));
- fHistd0d0LS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistd0d0LS"));
- fHistDCAJPSI = dynamic_cast<TH1F*>(fOutput->FindObject("fHistDCAJPSI"));
- fHistDCALS = dynamic_cast<TH1F*>(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;
-}
+++ /dev/null
-#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 <TROOT.h>
-#include <TSystem.h>
-#include <TList.h>
-#include <TH1F.h>
-#include <TH2F.h>
-
-#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
-
+++ /dev/null
-/**************************************************************************
- * 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 <TFile.h>
-#include <TClonesArray.h>
-#include <TDatabasePDG.h>
-
-#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<AliAODEvent*> (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<AliAODEvent*> (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<TClonesArray*>(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<TList*> (GetOutputData(1));
- if (!fOutput) {
- printf("ERROR: fOutput not available\n");
- return;
- }
-
- if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeMCjpsifromB"));
- fhDecayTime = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTime"));
- fhInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMass"));
- fhD0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0"));
- fhD0D0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0D0"));
- fhCosThetaStar = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaStar"));
- fhCosThetaPointing = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaPointing"));
- fhCtsVsD0D0 = dynamic_cast<TH2F*>(fOutput->FindObject("fhCtsVsD0D0"));
- fNtupleJPSI = dynamic_cast<TNtuple*>(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<TClonesArray*>(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;
-
-}
-
+++ /dev/null
-#ifndef ALIANALYSISTASKSEBTOJPSITOELE_H\r
-#define ALIANALYSISTASKSEBTOJPSITOELE_H\r
-/* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *\r
- * See cxx source for full Copyright notice */\r
-\r
-//*************************************************************************\r
-// Class AliAnalysisTaskSEBtoJPSItoEle\r
-// AliAnalysisTaskSE for the reconstruction \r
-// of heavy-flavour decay candidates\r
-// Author: C.Di Giglio, carmelo.digiglio@ba.infn.it\r
-//*************************************************************************\r
-\r
-class TH1F;\r
-class TList;\r
-class AliAnalysisBtoJPSItoEle;\r
-#include <TClonesArray.h>\r
-#include <TNtuple.h>\r
-#include <TH2F.h> \r
-#include "AliAnalysisTaskSE.h"\r
-\r
-class AliAnalysisTaskSEBtoJPSItoEle : public AliAnalysisTaskSE\r
-{\r
- public:\r
-\r
- AliAnalysisTaskSEBtoJPSItoEle();\r
- AliAnalysisTaskSEBtoJPSItoEle(const char *name);\r
- virtual ~AliAnalysisTaskSEBtoJPSItoEle();\r
-\r
- // Implementation of interface methods\r
- virtual void UserCreateOutputObjects();\r
- virtual void Init();\r
- virtual void LocalInit() {Init();}\r
- virtual void UserExec(Option_t *option);\r
- virtual void Terminate(Option_t *option);\r
-\r
- void SetCutsD0(const Double_t cutsD0[9]);\r
- void SetPtCuts(const Double_t ptcuts[2]);\r
- void SetAODMCInfo(Bool_t OkMCInfo) { fOkAODMC = OkMCInfo;}\r
- void SetMinimize(Bool_t OkMinimize) { fOkMinimize = OkMinimize;}\r
- void ReadAODMCInfo(AliAODEvent* aodEv, const TClonesArray* inArray);\r
- void SetPathMCfile (TString mcfilename="CsiMCfromKine_10PtBins.root") {fNameMCfile = mcfilename;}\r
- TH1F* OpenLocalMCFile(TString dataDir, Int_t nbin); \r
-\r
- private:\r
-\r
- AliAnalysisTaskSEBtoJPSItoEle(const AliAnalysisTaskSEBtoJPSItoEle &source);\r
- AliAnalysisTaskSEBtoJPSItoEle& operator=(const AliAnalysisTaskSEBtoJPSItoEle& source); \r
- //\r
- TList *fOutput; //! list send on output slot 0\r
- AliAnalysisBtoJPSItoEle *fCdfFit; // Unbinned log-likelihood minimizer \r
- TNtuple *fNtupleJPSI; //! Ntuple of pseudo-proper decay time & invariant mass values\r
- TH1F *fhDecayTimeMCjpsifromB; //! Pseudo-proper decay time distribution used as template for JPSIs from B\r
- TH1F *fhDecayTime; //! Pseudo-proper decay time distribution\r
- TH1F *fhInvMass; //! Invariant mass distribution\r
- TH1F *fhD0; //! Impact parameter distribution\r
- TH1F *fhD0D0; //! Product of impact parameters distributions\r
- TH1F *fhCosThetaStar; //! Cosine of decay angle distribution\r
- TH1F *fhCosThetaPointing; //! Cosine of pointing angle distribution\r
- TH2F *fhCtsVsD0D0; //! Cos theta star Vs. D0D0 distribution\r
- Bool_t fOkAODMC; // Flag to read AOD monte carlo information\r
- TString fNameMCfile; // Name of the MC file with X template\r
- Bool_t fOkMinimize; // Flag to enable unbinned log-likelihood minimization\r
- \r
- Double_t fCuts[9]; // cuts for N-tuple values selection\r
- Double_t fPtCuts[2]; // Max and min pt of the candidates\r
-\r
- ClassDef(AliAnalysisTaskSEBtoJPSItoEle,1); // AliAnalysisTaskSE for the reconstruction of heavy-flavour decay candidates\r
-};\r
-#endif\r
+++ /dev/null
-/**************************************************************************
- * 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<AliAODEvent*> (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<AliAODEvent*> (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<TClonesArray*>(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;it<aod->GetNumberOfTracks();it++) {
- trackin = aod->GetTrack(it);
- trkIDtoEntry[trackin->GetID()]=it;
- }
-/*
- for(Int_t it=0;it<aod->GetNumberOfTracks();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; ntracks<numTracks; ntracks++){
- trackin = (AliAODTrack*)arrayTracks->UncheckedAt(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<TList*> (GetOutputData(1));
- if (!fOutput) {
- printf("ERROR: fOutput not available\n");
- return;
- }
-
- fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs*fTotNegPairs);
-
- if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeMCjpsifromB"));
- fhDecayTime = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTime"));
- fhDecayTimeOut = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeOut"));
- fhInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMass"));
- fhdEdxTPC = dynamic_cast<TH2F*>(fOutput->FindObject("fhdEdxTPC"));
- fhD0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0"));
- fhD0D0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0D0"));
- fhCosThetaStar = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaStar"));
- fhCosThetaPointing = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaPointing"));
- fhCtsVsD0D0 = dynamic_cast<TH2F*>(fOutput->FindObject("fhCtsVsD0D0"));
- fHistMassLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMassLS"));
- fHistCtsLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLS"));
- fHistCtsLSpos = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSpos"));
- fHistCtsLSneg = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSneg"));
- fHistCPtaLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCPtaLS"));
- fHistd0d0LS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistd0d0LS"));
- fhDCA = dynamic_cast<TH1F*>(fOutput->FindObject("fhDCA"));
- fHistDCALS = dynamic_cast<TH1F*>(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<TClonesArray*>(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
- }
- }
-}
+++ /dev/null
-#ifndef ALIANALYSISTASKSEJPSITOELE_H\r
-#define ALIANALYSISTASKSEJPSITOELE_H\r
-/* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *\r
- * See cxx source for full Copyright notice */\r
-\r
-//*************************************************************************\r
-// Class AliAnalysisTaskSEJPSItoEle\r
-// AliAnalysisTaskSE class to select \r
-// J/psi -> e+e- candidates only and save them \r
-// into a specific stand-alone AOD file \r
-// Author: C.Di Giglio, carmelo.digiglio@ba.infn.it\r
-//*************************************************************************\r
-\r
-#include "TClonesArray.h"\r
-#include "TH1F.h" \r
-#include "TH2F.h"\r
-#include "TList.h"\r
-#include "TChain.h"\r
-\r
-#include "AliAnalysisTaskSE.h"\r
-#include "AliAODEvent.h"\r
-#include "AliAnalysisManager.h"\r
-#include "AliAnalysisVertexingHF.h"\r
-#include "AliPID.h" //************\r
-#include "AliAODPid.h"\r
-#include "AliExternalTrackParam.h"//***************\r
-\r
-class AliAnalysisTaskSEJPSItoEle : public AliAnalysisTaskSE\r
-{\r
- public:\r
-\r
- AliAnalysisTaskSEJPSItoEle();\r
- AliAnalysisTaskSEJPSItoEle(const char *name);\r
- virtual ~AliAnalysisTaskSEJPSItoEle();\r
-\r
- virtual void UserCreateOutputObjects();\r
- virtual void Init();\r
- virtual void LocalInit() {Init();}\r
- virtual void UserExec(Option_t *option);\r
- virtual void Terminate(Option_t *option);\r
-\r
- void SetCutsJPSI(const Double_t cutsJPSI[9]);\r
- void SetPtCuts(const Double_t ptcuts[2]);\r
- void SetAODMCInfo(Bool_t OkMCInfo) { fOkAODMC = OkMCInfo;}\r
- void ReadAODMCInfo(AliAODEvent* aodEv, const TClonesArray* inArray);\r
- //\r
- //\r
- //\r
- Double_t GetExpectedSignal(const Float_t mom,\r
- AliPID::EParticleType n=AliPID::kKaon) const {\r
-\r
- Double_t mass=AliPID::ParticleMass(n);\r
- Double_t betaGamma = mom/mass;\r
- Double_t bb = AliExternalTrackParam::BetheBlochAleph(betaGamma,0.0283086,2.63394e+01,5.04114e-11,2.12543,4.88663);\r
- return bb*50.; //bb*fMIP;\r
- }\r
-\r
- Double_t GetExpectedSigma(const Float_t mom, const Int_t nPoints,\r
- AliPID::EParticleType n=AliPID::kKaon) const {\r
- if (nPoints != 0) \r
- return GetExpectedSignal(mom,n)*0.07*sqrt(1. + 0./nPoints);\r
- else\r
- return GetExpectedSignal(mom,n)*0.07;\r
- }\r
-\r
- Float_t GetNumberOfSigmas(const Float_t mom, const Float_t dEdx, \r
- const Int_t nPoints,\r
- AliPID::EParticleType n=AliPID::kKaon) const {\r
-\r
- Double_t bethe=GetExpectedSignal(mom,n);\r
- Double_t sigma=GetExpectedSigma(mom,nPoints,n);\r
- return (dEdx-bethe)/sigma;\r
- }\r
- //\r
- //\r
- //\r
- private:\r
-\r
- AliAnalysisTaskSEJPSItoEle(const AliAnalysisTaskSEJPSItoEle &source);\r
- AliAnalysisTaskSEJPSItoEle& operator=(const AliAnalysisTaskSEJPSItoEle& source); \r
- \r
- TList *fOutput; //! list send on output slot 0\r
-\r
- //output histograms\r
- TH1F *fhDecayTimeMCjpsifromB; //! Pseudo-proper decay time distribution used as template for JPSIs from B\r
- TH1F *fhDecayTime; //! Pseudo-proper decay time distribution\r
- TH1F *fhDecayTimeOut; //! Pseudo-proper decay time distribution (stand-alone AOD)\r
- TH1F *fhInvMass; //! Invariant mass distribution\r
- TH2F *fhdEdxTPC;\r
- TH1F *fhD0; //! Impact parameter distribution\r
- TH1F *fhD0D0; //! Product of impact parameters distributions\r
- TH1F *fhCosThetaStar; //! Cosine of decay angle distribution\r
- TH1F *fhCosThetaPointing; //! Cosine of pointing angle distribution\r
- TH1F *fhDCA; //! Distance of closest approach\r
- TH2F *fhCtsVsD0D0; //! Cos theta star Vs. D0D0 distribution\r
-\r
- //like sign pairs histograms \r
- TH1F *fHistMassLS; //! Invariant mass distribution\r
- TH1F *fHistCtsLS; //! Cosine of decay angle distribution\r
- TH1F *fHistCtsLSpos; //! Cosine of decay angle distribution (++ pairs)\r
- TH1F *fHistCtsLSneg; //! Cosine of decay angle distribution (-- pairs)\r
- TH1F *fHistCPtaLS; //! Cosine of pointing angle distribution\r
- TH1F *fHistd0d0LS; //! Product of impact parameters distributions\r
- TH1F *fHistDCALS; //! Distance of closest approach\r
- AliAnalysisVertexingHF *fVHF; //! Vertexer heavy flavour (used to pass the cuts)\r
-\r
- //Int_t totJPSIin;\r
- //Int_t totJPSIout;\r
-\r
- //like sign spectrum normalization\r
- Int_t fTotPosPairs; //\r
- Int_t fTotNegPairs; //\r
- Double_t fLsNormalization; // Like sign normalization factor\r
-\r
- //flags for analysis\r
- Bool_t fOkAODMC; // Flag to read AOD monte carlo information\r
- Bool_t fOkLikeSign; // Flag to select like sign candidates analysis\r
-\r
- //momentum cuts \r
- Double_t fCuts[9]; // cuts for N-tuple values selection\r
- Double_t fPtCuts[2]; // Max and min pt of the candidates\r
-\r
- TClonesArray *fVerticesHFTClArr; // Array of heavy flavor vertices to be replicated in stand-alone AOD\r
- TClonesArray *fJpsiToEleTClArr; // Array of J/psi->e+e- candidates to be replicated in stand-alone AOD\r
- TClonesArray *fLikeSignTClArr; // Array of like sign candidates to be replicated in stand-alone AOD\r
- TClonesArray *fTracksTClArr; // Array of tracks belonging to J/psi->e+e- candidates to be replicated in stand-alone AOD \r
- AliAODEvent *fOrigAOD; // original AOD event\r
- AliAODEvent *fNewAOD; // new AOD event with only JPSItoEle candidates stored\r
-\r
- ClassDef(AliAnalysisTaskSEJPSItoEle,1); // AliAnalysisTaskSE for the reconstruction of heavy-flavour decay candidates\r
-};\r
-\r
-\r
-#endif\r
+++ /dev/null
-/**************************************************************************
- * 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
-
-}
-
+++ /dev/null
-#ifndef ALIBTOJPSITOELECDFFITFCN_H\r
-#define ALIBTOJPSITOELECDFFITFCN_H\r
-/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
- * See cxx source for full Copyright notice */\r
-\r
-//_________________________________________________________________________\r
-// Class AliBtoJPSItoEleCDFfitFCN\r
-// Definition of main function used in \r
-// unbinned log-likelihood fit for\r
-// the channel B -> JPsi + X -> e+e- + X\r
-// \r
-// Origin: C.Di Giglio\r
-// Contact: Carmelo.Digiglio@ba.infn.it , Giuseppe.Bruno@ba.infn.it\r
-//_________________________________________________________________________\r
-\r
-#include <TNamed.h>\r
-#include <TDatabasePDG.h>\r
-#include "TH1F.h"\r
-#include "TMath.h"\r
-#include "TRandom3.h"\r
-\r
-enum IntegralType {kBkg, \r
- kBkgNorm, \r
- kSig, \r
- kSigNorm};\r
-\r
-enum PtBins {kallpt, \r
- kptbin1,kptbin2,kptbin3,\r
- kptbin4,kptbin5,kptbin6,\r
- kptbin7,kptbin8,kptbin9};\r
-//_________________________________________________________________________________________________\r
-class AliBtoJPSItoEleCDFfitFCN : public TNamed {\r
- public:\r
- //\r
- AliBtoJPSItoEleCDFfitFCN();\r
- AliBtoJPSItoEleCDFfitFCN(const AliBtoJPSItoEleCDFfitFCN& source); \r
- AliBtoJPSItoEleCDFfitFCN& operator=(const AliBtoJPSItoEleCDFfitFCN& source); \r
- virtual ~AliBtoJPSItoEleCDFfitFCN();\r
-\r
- Double_t EvaluateLikelihood(const Double_t* pseudoproperdecaytime,\r
- const Double_t* invariantmass, const Int_t ncand);\r
- \r
- Double_t GetFPlus() const { return fFPlus; }\r
- Double_t GetFMinus() const { return fFMinus; }\r
- Double_t GetFSym() const { return fFSym; }\r
- Double_t GetRadius() const { return fParameters[0]; }\r
- Double_t GetTheta() const { return fParameters[1]; }\r
- Double_t GetPhi() const { return fParameters[2]; }\r
- Double_t GetLamPlus() const { return fParameters[3]; }\r
- Double_t GetLamMinus() const { return fParameters[4]; }\r
- Double_t GetLamSym() const { return fParameters[5]; }\r
- Double_t GetMassSlope() const { return fParameters[6]; }\r
- Double_t GetFractionJpsiFromBeauty() const { return fParameters[7]; }\r
- Double_t GetFsig() const { return fParameters[8]; }\r
- Double_t GetCrystalBallMmean() const { return fParameters[9]; }\r
- Double_t GetCrystalBallNexp() const { return fParameters[10]; }\r
- Double_t GetCrystalBallSigma() const { return fParameters[11]; }\r
- Double_t GetCrystalBallAlpha() const { return fParameters[12]; }\r
- Double_t GetCrystalBallNorm() const { return fParameters[13]; }\r
- Double_t GetSigmaResol() const { return fParameters[14]; }\r
- Double_t GetNResol() const { return fParameters[15]; }\r
- Double_t GetIntegral() const { return fIntegral; }\r
- Double_t GetIntegralFunB() const { return fintxFunB; }\r
- Double_t GetIntegralBkgPos() const { return fintxDecayTimeBkgPos; }\r
- Double_t GetIntegralBkgNeg() const { return fintxDecayTimeBkgNeg; }\r
- Double_t GetIntegralBkgSym() const { return fintxDecayTimeBkgSym; }\r
- Double_t GetIntegralMassSig() const { return fintmMassSig; }\r
- Double_t GetIntegralRes() const { return fintxRes; }\r
- Double_t GetIntegralMassBkg() const { return fintmMassBkg; }\r
- Bool_t GetCrystalBallParam() const { return fCrystalBallParam; }\r
-\r
- void SetFPlus(Double_t plus) {fFPlus = plus;}\r
- void SetFMinus(Double_t minus) {fFMinus = minus;}\r
- void SetFSym(Double_t sym) {fFSym = sym;}\r
- void SetRadius(Double_t radius) {fParameters[0] = radius;}\r
- void SetTheta(Double_t theta) {fParameters[1] = theta;}\r
- void SetPhi(Double_t phi) {fParameters[2] = phi;}\r
- void SetLamPlus(Double_t lamplus) {fParameters[3] = lamplus;}\r
- void SetLamMinus(Double_t lamminus) {fParameters[4] = lamminus;}\r
- void SetLamSym(Double_t lamsym) {fParameters[5] = lamsym;}\r
- void SetMassSlope(Double_t slope) {fParameters[6] = slope;}\r
- void SetFractionJpsiFromBeauty(Double_t B) {fParameters[7] = B;}\r
- void SetFsig(Double_t Fsig) {fParameters[8] = Fsig;}\r
- void SetCrystalBallMmean(Double_t CrystalBallMmean) {fParameters[9] = CrystalBallMmean;}\r
- void SetCrystalBallNexp(Double_t CrystalBallNexp) {fParameters[10] = CrystalBallNexp;}\r
- void SetCrystalBallSigma(Double_t CrystalBallSigma) {fParameters[11] = CrystalBallSigma;}\r
- void SetCrystalBallAlpha(Double_t CrystalBallAlpha) {fParameters[12] = CrystalBallAlpha;}\r
- void SetCrystalBallNorm(Double_t CrystalBallNorm) {fParameters[13] = CrystalBallNorm;}\r
- void SetSigmaResol(Double_t SigmaResol) {fParameters[14] = SigmaResol;}\r
- void SetNResol(Double_t NResol) {fParameters[15] = NResol;}\r
-\r
- void SetAllParameters(const Double_t* parameters);\r
- void SetIntegral(Double_t integral) { fIntegral = integral; }\r
- void SetIntegralFunB(Double_t integral) { fintxFunB = integral; }\r
- void SetIntegralBkgPos(Double_t integral) { fintxDecayTimeBkgPos = integral; }\r
- void SetIntegralBkgNeg(Double_t integral) { fintxDecayTimeBkgNeg = integral; }\r
- void SetIntegralBkgSym(Double_t integral) { fintxDecayTimeBkgSym = integral; }\r
- void SetIntegralMassSig(Double_t integral) { fintmMassSig = integral; }\r
- void SetIntegralRes(Double_t integral) { fintxRes = integral; }\r
- void SetIntegralMassBkg(Double_t integral) { fintmMassBkg = integral; }\r
-\r
- void SetCsiMC(const TH1F* MCtemplate) {fhCsiMC = (TH1F*)MCtemplate->Clone("fhCsiMC");}\r
-\r
- void SetResolutionConstants();\r
- void SetMassWndHigh(Double_t limit) { fMassWndHigh = TDatabasePDG::Instance()->GetParticle(443)->Mass() + limit ;}\r
- void SetMassWndLow(Double_t limit) { fMassWndLow = TDatabasePDG::Instance()->GetParticle(443)->Mass() - limit ;}\r
- void SetCrystalBallFunction(Bool_t okCB) {fCrystalBallParam = okCB;}\r
-\r
- void ConvertFromSpherical() { fFPlus = TMath::Power((fParameters[0]*TMath::Cos(fParameters[1])),2.);\r
- fFMinus = TMath::Power((fParameters[0]*TMath::Sin(fParameters[1])*TMath::Sin(fParameters[2])),2.);\r
- fFSym = TMath::Power((fParameters[0]*TMath::Sin(fParameters[1])*TMath::Cos(fParameters[2])),2.);} \r
-\r
- void ComputeIntegral(); \r
-\r
- void ReadMCtemplates(Int_t BinNum);\r
-\r
- void PrintStatus();\r
-\r
- private: \r
- //\r
- Double_t fParameters[16]; /* par[0] = fRadius; \r
- par[1] = fTheta;\r
- par[2] = fPhi;\r
- par[3] = fOneOvLamPlus;\r
- par[4] = fOneOvLamMinus;\r
- par[5] = fOneOvLamSym;\r
- par[6] = fMassBkgSlope;\r
- par[7] = fFractionJpsiFromBeauty;\r
- par[8] = fFsig;\r
- par[9] = fCrystalBallMmean;\r
- par[10] = fCrystalBallNexp;\r
- par[11] = fCrystalBallSigma;\r
- par[12] = fCrystalBallAlpha;\r
- par[13] = fCrystalBallNorm;\r
- par[14] = fSigmaResol;\r
- par[15] = fNResol; */\r
-\r
-\r
- Double_t fFPlus; // parameters of the log-likelihood function\r
- Double_t fFMinus; // Slopes of the x distributions of the background \r
- Double_t fFSym; // functions \r
-\r
- Double_t fIntegral; // integral values of log-likelihood function terms\r
- Double_t fintxFunB;\r
- Double_t fintxDecayTimeBkgPos;\r
- Double_t fintxDecayTimeBkgNeg;\r
- Double_t fintxDecayTimeBkgSym;\r
- Double_t fintmMassSig;\r
- Double_t fintxRes;\r
- Double_t fintmMassBkg;\r
- TH1F *fhCsiMC; // X distribution used as MC template for JPSI from B\r
- Double_t fResolutionConstants[5]; // constants for the parametrized resolution function R(X)\r
- Double_t fMassWndHigh; // JPSI Mass window higher limit\r
- Double_t fMassWndLow; // JPSI Mass window lower limit\r
- Bool_t fCrystalBallParam; // Boolean to switch to Crystall Ball parameterisation\r
-\r
- ////\r
-\r
- Double_t EvaluateCDFfunc(Double_t x, Double_t m) const ;\r
-\r
- Double_t EvaluateCDFfuncNorm(Double_t x, Double_t m) const ;\r
-\r
- ////\r
-\r
- Double_t EvaluateCDFfuncSignalPart(Double_t x, Double_t m) const ; // Signal part \r
-\r
- Double_t EvaluateCDFDecayTimeSigDistr(Double_t x) const ;\r
-\r
- Double_t EvaluateCDFInvMassSigDistr(Double_t m) const ;\r
-\r
- Double_t EvaluateCDFfuncBkgPart(Double_t x,Double_t m) const ; // Background part\r
-\r
- Double_t EvaluateCDFDecayTimeBkgDistr(Double_t x) const ;\r
- \r
- Double_t EvaluateCDFInvMassBkgDistr(Double_t m) const ;\r
-\r
- ////\r
-\r
- Double_t FunB(Double_t x) const ;\r
-\r
- Double_t FunP(Double_t x) const ;\r
-\r
- Double_t CsiMC(Double_t x) const ;\r
-\r
- Double_t FunBkgPos(Double_t x) const ;\r
-\r
- Double_t FunBkgNeg(Double_t x) const ;\r
-\r
- Double_t FunBkgSym(Double_t x) const ;\r
-\r
- Double_t ResolutionFunc(Double_t x) const ;\r
- //\r
- ClassDef (AliBtoJPSItoEleCDFfitFCN,1); // Unbinned log-likelihood fit \r
-\r
-};\r
-\r
-#endif\r
+++ /dev/null
-/**************************************************************************
- * 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 <TVirtualFitter.h>
-#include <TStopwatch.h>
-#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);
-}
-
+++ /dev/null
-#ifndef ALIBTOJPSITOELECDFFITHANDLER_H\r
-#define ALIBTOJPSITOELECDFFITHANDLER_H\r
-/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
- * See cxx source for full Copyright notice */\r
-\r
-//___________________________________________________________________________\r
-// Class AliBtoJPSItoEleCDFfitHandler\r
-// Class to perform unbinned log-likelihood fit\r
-// \r
-// Origin: C. Di Giglio\r
-// Contact: Carmelo.Digiglio@ba.infn.it; Giuseppe.Bruno@ba.infn.it\r
-//____________________________________________________________________________\r
-\r
-#include <TNamed.h>\r
-#include <TBits.h>\r
-#include "TMath.h"\r
-#include "AliBtoJPSItoEleCDFfitFCN.h"\r
-#include "AliLog.h"\r
-\r
-class TBits; \r
-\r
-class AliBtoJPSItoEleCDFfitHandler : public TNamed {\r
- public:\r
- //\r
- AliBtoJPSItoEleCDFfitHandler();\r
- AliBtoJPSItoEleCDFfitHandler& operator= (const AliBtoJPSItoEleCDFfitHandler& c);\r
- AliBtoJPSItoEleCDFfitHandler(const AliBtoJPSItoEleCDFfitHandler& c);\r
- AliBtoJPSItoEleCDFfitHandler(Double_t* decaytime, Double_t* invariantmass, Int_t ncand);\r
- ~AliBtoJPSItoEleCDFfitHandler(); \r
- //Double_t Up() const { return fUp*TMath::Sqrt(fLikely->GetIntegral()); }\r
- Double_t Up() const { return fUp; }\r
- void SetErrorDef(Double_t up) {fUp = up;}\r
- void SetPrintStatus(Bool_t okPrintStatus) { fPrintStatus = okPrintStatus; } \r
- Bool_t GetPrintStatus() { return fPrintStatus ; }\r
- void SetParamStartValues (Double_t*);\r
- Double_t* GetStartParamValues() { return fParamStartValues; }\r
- TBits GetFixedParamList() { return fIsParamFixed; }\r
- void FixParam(UInt_t param, Bool_t value) { fIsParamFixed.SetBitNumber(param,value); }\r
- void FixAllParam(Bool_t value) { for(UInt_t par=0;par<16;par++) fIsParamFixed.SetBitNumber(par,value); }\r
- Bool_t IsParamFixed(UInt_t param) { return fIsParamFixed.TestBitNumber(param); }\r
- void SetResolutionConstants();\r
- void SetCrystalBallFunction(Bool_t okCB);\r
- void SetMassWndHigh(Double_t limit);\r
- void SetMassWndLow(Double_t limit);\r
-\r
- Double_t operator()(const Double_t* par) const ;\r
- void CdfFCN(Int_t & /* npar */, Double_t * /* gin */, Double_t &f, Double_t *par, Int_t /* iflag */);\r
-\r
- Double_t* Decaytime() const { return fX; }\r
- Double_t* InvariantMass() const { return fM; }\r
- AliBtoJPSItoEleCDFfitFCN* LikelihoodPointer() const { return fLikely; }\r
- Int_t DoMinimization();\r
-\r
- private:\r
- //\r
- TBits fIsParamFixed; //array of bits: 0 = param free; 1 = param fixed;\r
- Bool_t fPrintStatus; //flag to enable the prit out of the algorithm at each step\r
- Double_t fParamStartValues[16]; //array of parameters input value\r
- Double_t fUp; //error definition \r
- Double_t* fX; //pseudo-proper decay time X\r
- Double_t* fM; //invariant mass M\r
- AliBtoJPSItoEleCDFfitFCN* fLikely; //Log likelihood function\r
- Int_t fNcand; //number of candidates\r
- //\r
- ClassDef(AliBtoJPSItoEleCDFfitHandler,1);\r
-\r
-}; \r
-#endif\r
+++ /dev/null
-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;
-}
+++ /dev/null
-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;
-}
+++ /dev/null
-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;
-}
}
+ // 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;
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()<<endl;
- hMass->Fill(d->InvMassD0(),0.5);
- hMass->Fill(d->InvMassD0bar(),0.5);
- hCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
- hSecVtxZ->Fill(d->GetSecVtxZ());
- //cout<<d->GetSecVtxX() <<endl;
-
// get daughter AOD tracks
AliAODTrack *trk0 = (AliAODTrack*)d->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"<<endl;
- }
- //cout<<"pt of positive track: "<<trk0->Pt()<<endl;
-
- // make a AliNeutralTrackParam from the D0
- // and calculate impact parameters to primary vertex
- AliNeutralTrackParam trackD0(d);
- //cout<<"pt of D0: "<<d->Pt()<<" (AliAODRecoDecay); "<<trackD0.Pt()<<" (track)"<<endl;
- //trackD0.Print();
- Double_t dz[2],covdz[3];
- trackD0.PropagateToDCA(vtx1,aod->GetMagneticField(),1000.,dz,covdz);
- //cout<<"D0 impact parameter rphi: "<<dz[0]<<" +- "<<TMath::Sqrt(covdz[0])<<endl;
- }
- if(unsetvtx) d->UnsetOwnPrimaryVtx();
+
+ if(trk0->GetStatus()) printf("ok %d\n",iD0toKpi);
}
nTotDstar += nDstar;
cout<<"Number of D*->D0pi: "<<nDstar<<endl;
+ AliRDHFCutsDStartoKpipi *cuts = new AliRDHFCutsDStartoKpipi("CutsDStartoKpipi");
+
+
+ AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
+ //esdTrackCuts->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: "<<trk->Pt()<<endl;
- //cout<<"pt of D0: "<<d->Pt()<<endl;
- }
- if(unsetvtx) {
- c->UnsetOwnPrimaryVtx();
- c->Get2Prong()->UnsetOwnPrimaryVtx();
- }
+ if(cuts->IsSelected(c,AliRDHFCuts::kTracks)) printf("passed\n");
+
}
Int_t nVtxsHF = arrayVerticesHF->GetEntriesFast();
nTotHF += nVtxsHF;
cout<<"Number of heavy-flavour vertices: "<<nVtxsHF<<endl;
- for (Int_t iVtx = 0; iVtx < nVtxsHF; iVtx++) {
+ for (Int_t iVtx = 0; iVtx <0; iVtx++) {
AliAODVertex *vtxHF = (AliAODVertex*)arrayVerticesHF->UncheckedAt(iVtx);
// print info
//cout << iVtx << ": vertex z position: " << vtxHF->GetZ() << endl;
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff