-/**************************************************************************
- * 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 *
- * appeuear 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. *
- **************************************************************************/
-
-/* $Id$ */
-
-//
-//
-// Base class for DStar Analysis
-//
-//
-// The D* spectra study is done in pt bins:
-// [0,0.5] [0.5,1] [1,2] [2,3] [3,4] [4,5] [5,6] [6,7] [7,8],
-// [8,10],[10,12], [12,16], [16,20] and [20,24]
-//
-// Cuts arew centralized in AliRDHFCutsDStartoKpipi
-// Side Band and like sign background are implemented in the macro
-//
-//-----------------------------------------------------------------------
-//
-// Author A.Grelli
-// ERC-QGP Utrecht University - a.grelli@uu.nl,
-// Author Y.Wang
-// University of Heidelberg - yifei@physi.uni-heidelberg.de
-// Author C.Ivan
-// ERC-QGP Utrecht University - c.ivan@uu.nl,
-//
-//-----------------------------------------------------------------------
-
-#include <TSystem.h>
-#include <TParticle.h>
-#include <TH1I.h>
-#include "TROOT.h"
-#include <TDatabasePDG.h>
-#include <AliAnalysisDataSlot.h>
-#include <AliAnalysisDataContainer.h>
-#include "AliRDHFCutsDStartoKpipi.h"
-#include "AliStack.h"
-#include "AliMCEvent.h"
-#include "AliAnalysisManager.h"
-#include "AliAODMCHeader.h"
-#include "AliAODHandler.h"
-#include "AliLog.h"
-#include "AliAODVertex.h"
-#include "AliAODRecoDecay.h"
-#include "AliAODRecoDecayHF.h"
-#include "AliAODRecoCascadeHF.h"
-#include "AliAODRecoDecayHF2Prong.h"
-#include "AliAnalysisVertexingHF.h"
-#include "AliESDtrack.h"
-#include "AliAODMCParticle.h"
-#include "AliAnalysisTaskSE.h"
-#include "AliAnalysisTaskSEDStarSpectra.h"
-#include "AliNormalizationCounter.h"
-
-ClassImp(AliAnalysisTaskSEDStarSpectra)
-
-//__________________________________________________________________________
-AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra():
- AliAnalysisTaskSE(),
- fEvents(0),
- fAnalysis(0),
- fD0Window(0),
- fPeakWindow(0),
- fUseMCInfo(kFALSE),
- fDoSearch(kFALSE),
- fOutput(0),
- fOutputAll(0),
- fOutputPID(0),
- fNSigma(3),
- fCuts(0),
- fCEvents(0),
- fTrueDiff2(0),
- fDeltaMassD1(0),
- fCounter(0)
-{
- //
- // Default ctor
- //
-}
-//___________________________________________________________________________
-AliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra(const Char_t* name, AliRDHFCutsDStartoKpipi* cuts) :
- AliAnalysisTaskSE(name),
- fEvents(0),
- fAnalysis(0),
- fD0Window(0),
- fPeakWindow(0),
- fUseMCInfo(kFALSE),
- fDoSearch(kFALSE),
- fOutput(0),
- fOutputAll(0),
- fOutputPID(0),
- fNSigma(3),
- fCuts(0),
- fCEvents(0),
- fTrueDiff2(0),
- fDeltaMassD1(0),
- fCounter(0)
-{
- //
- // Constructor. Initialization of Inputs and Outputs
- //
- Info("AliAnalysisTaskSEDStarSpectra","Calling Constructor");
-
- fCuts=cuts;
-
- DefineOutput(1,TList::Class()); //conters
- DefineOutput(2,TList::Class()); //All Entries output
- DefineOutput(3,TList::Class()); //3sigma PID output
- DefineOutput(4,AliRDHFCutsDStartoKpipi::Class()); //My private output
- DefineOutput(5,AliNormalizationCounter::Class()); // normalization
-}
-
-//___________________________________________________________________________
-AliAnalysisTaskSEDStarSpectra::~AliAnalysisTaskSEDStarSpectra() {
- //
- // destructor
- //
- Info("~AliAnalysisTaskSEDStarSpectra","Calling Destructor");
-
- if (fOutput) {
- delete fOutput;
- fOutput = 0;
- }
- if (fOutputAll) {
- delete fOutputAll;
- fOutputAll = 0;
- }
- if (fOutputPID) {
- delete fOutputPID;
- fOutputPID = 0;
- }
- if (fCuts) {
- delete fCuts;
- fCuts = 0;
- }
- if(fCEvents){
- delete fCEvents;
- fCEvents =0;
- }
- if(fDeltaMassD1){
- delete fDeltaMassD1;
- fDeltaMassD1 =0;
- }
-}
-//_________________________________________________
-void AliAnalysisTaskSEDStarSpectra::Init(){
- //
- // Initialization
- //
-
- if(fDebug > 1) printf("AnalysisTaskSEDStarSpectra::Init() \n");
- AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
- // Post the data
- PostData(4,copyfCuts);
-
- return;
-}
-
-//_________________________________________________
-void AliAnalysisTaskSEDStarSpectra::UserExec(Option_t *)
-{
- // user exec
- if (!fInputEvent) {
- Error("UserExec","NO EVENT FOUND!");
- return;
- }
-
- fEvents++;
-
- AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
- TClonesArray *arrayDStartoD0pi=0;
-
- fCEvents->Fill(1);
-
- if(!aodEvent && 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.
- aodEvent = 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();
- arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
- }
- } else {
- arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
- }
-
- // fix for temporary bug in ESDfilter
- // the AODs with null vertex pointer didn't pass the PhysSel
- if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
-
- fCEvents->Fill(2);
-
- fCounter->StoreEvent(aodEvent,fUseMCInfo);
-
- if(!fCuts->IsEventSelected(aodEvent)) {
- if(fCuts->GetWhyRejection()==6) // rejected for Z vertex
- fCEvents->Fill(6);
- return;
- }
-
- // Load the event
- AliInfo(Form("Event %d",fEvents));
- if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
-
- // counters for efficiencies
- Int_t icountReco = 0;
-
- //D* and D0 prongs needed to MatchToMC method
- Int_t pdgDgDStartoD0pi[2]={421,211};
- Int_t pdgDgD0toKpi[2]={321,211};
-
- // AOD primary vertex
- AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
- if(!vtx1) return;
- if(vtx1->GetNContributors()<1) return;
-
- if (!arrayDStartoD0pi){
- AliInfo("Could not find array of HF vertices, skipping the event");
- return;
- }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
-
- Int_t nSelectedAna =0;
- Int_t nSelectedProd =0;
-
- // loop over the tracks to search for candidates soft pion
- for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
-
- // D* candidates and D0 from D*
- AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
- if(!dstarD0pi->GetSecondaryVtx()) continue;
- AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
- if (!theD0particle) continue;
-
- Int_t isDStar = 0;
- TClonesArray *mcArray = 0; // fix coverity
-
- // mc analysis
- if(fUseMCInfo){
- //MC array need for maching
- mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
- if (!mcArray) {
- AliError("Could not find Monte-Carlo in AOD");
- return;
- }
- // find associated MC particle for D* ->D0toKpi
- Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
- if(mcLabel>=0) isDStar = 1;
- }
-
- Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt());
-
- // set the D0 search window bin by bin - useful to calculate side band bkg
- if (ptbin==0){
- if(fAnalysis==1){
- fD0Window=0.035;
- fPeakWindow=0.03;
- }else{
- fD0Window=0.020;
- fPeakWindow=0.0018;
- }
- }
- if (ptbin==1){
- if(fAnalysis==1){
- fD0Window=0.035;
- fPeakWindow=0.03;
- }else{
- fD0Window=0.020;
- fPeakWindow=0.0018;
- }
- }
- if (ptbin==2){
- if(fAnalysis==1){
- fD0Window=0.035;
- fPeakWindow=0.03;
- }else{
- fD0Window=0.020;
- fPeakWindow=0.0018;
- }
- }
- if (ptbin==3){
- if(fAnalysis==1){
- fD0Window=0.035;
- fPeakWindow=0.03;
- }else{
- fD0Window=0.022;
- fPeakWindow=0.0016;
- }
- }
- if (ptbin==4){
- if(fAnalysis==1){
- fD0Window=0.035;
- fPeakWindow=0.03;
- }else{
- fD0Window=0.026;
- fPeakWindow=0.0014;
- }
- }
- if (ptbin==5){
- if(fAnalysis==1){
- fD0Window=0.045;
- fPeakWindow=0.03;
- }else{
- fD0Window=0.026;
- fPeakWindow=0.0014;
- }
- }
- if (ptbin==6){
- if(fAnalysis==1){
- fD0Window=0.045;
- fPeakWindow=0.03;
- }else{
- fD0Window=0.026;
- fPeakWindow=0.006;
- }
- }
- if (ptbin==7){
- if(fAnalysis==1){
- fD0Window=0.055;
- fPeakWindow=0.03;
- }else{
- fD0Window=0.026;
- fPeakWindow=0.006;
- }
- }
- if (ptbin>7){
- if(fAnalysis==1){
- fD0Window=0.074;
- fPeakWindow=0.03;
- }else{
- fD0Window=0.026;
- fPeakWindow=0.006;
- }
- }
- fCEvents->Fill(9);
-
- nSelectedProd++;
- nSelectedAna++;
-
- Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
-
- // here for lead-lead
- Double_t invmassD0 = dstarD0pi->InvMassD0();
- if (TMath::Abs(invmassD0-mPDGD0)>fD0Window) continue;
-
- Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
- Double_t invmassDelta = dstarD0pi->DeltaInvMass();
-
- if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>fPeakWindow) continue;
-
- Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
- if(!isTkSelected) continue;
-
- if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
-
- Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
- if (!isSelected) continue;
-
- // fill PID
- FillSpectrum(dstarD0pi,isDStar,fCuts,fOutputPID);
- SideBandBackground(dstarD0pi,fCuts,fOutputPID);
- WrongSignForDStar(dstarD0pi,fCuts,fOutputPID);
-
- // fill no pid
- fCuts->SetUsePID(kFALSE);
- FillSpectrum(dstarD0pi,isDStar,fCuts,fOutputAll);
- SideBandBackground(dstarD0pi,fCuts,fOutputAll);
- WrongSignForDStar(dstarD0pi,fCuts,fOutputAll);
- fCuts->SetUsePID(kTRUE);
-
- // rare D search ------
- if(fDoSearch){
- TLorentzVector LorentzTrack1(0,0,0,0); // lorentz 4 vector
- TLorentzVector LorentzTrack2(0,0,0,0); // lorentz 4 vector
-
- for (Int_t i=0; i<aodEvent->GetNTracks(); i++){
-
- AliAODTrack* aodTrack = aodEvent->GetTrack(i);
-
- if(dstarD0pi->Charge() == aodTrack->Charge()) continue;
- if((!(aodTrack->GetStatus()&AliESDtrack::kITSrefit)|| (!(aodTrack->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
- if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>0.02) continue;
-
- //build the D1 mass
- Double_t mass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
-
- LorentzTrack1.SetPxPyPzE( dstarD0pi->Px(),dstarD0pi->Py(), dstarD0pi->Pz(), dstarD0pi->E(413) );
- LorentzTrack2.SetPxPyPzE( aodTrack->Px(),aodTrack->Py(), aodTrack->Pz(),aodTrack->E(mass) );
-
- //D1 mass
- Double_t d1mass = ((LorentzTrack1+LorentzTrack2).M());
- //mass difference - at 0.4117 and 0.4566
- fDeltaMassD1->Fill(d1mass-dstarD0pi->InvMassDstarKpipi());
- }
- }
-
- if(isDStar == 1) {
- fTrueDiff2->Fill(dstarD0pi->Pt(),dstarD0pi->DeltaInvMass());
- }
-
- }
-
- fCounter->StoreCandidates(aodEvent,nSelectedProd,kTRUE);
- fCounter->StoreCandidates(aodEvent,nSelectedAna,kFALSE);
-
- AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
-
- PostData(1,fOutput);
- PostData(2,fOutputAll);
- PostData(3,fOutputPID);
- PostData(5,fCounter);
-
-}
-//________________________________________ terminate ___________________________
-void AliAnalysisTaskSEDStarSpectra::Terminate(Option_t*)
-{
- // The Terminate() function is the last function to be called during
- // a query. It always runs on the client, it can be used to present
- // the results graphically or save the results to file.
-
- //Info("Terminate","");
- AliAnalysisTaskSE::Terminate();
-
- fOutput = dynamic_cast<TList*> (GetOutputData(1));
- if (!fOutput) {
- printf("ERROR: fOutput not available\n");
- return;
- }
-
- fCEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fCEvents"));
- fDeltaMassD1 = dynamic_cast<TH1F*>(fOutput->FindObject("fDeltaMassD1"));
- fTrueDiff2 = dynamic_cast<TH2F*>(fOutput->FindObject("fTrueDiff2"));
-
- fOutputAll = dynamic_cast<TList*> (GetOutputData(1));
- if (!fOutputAll) {
- printf("ERROR: fOutputAll not available\n");
- return;
- }
- fOutputPID = dynamic_cast<TList*> (GetOutputData(2));
- if (!fOutputPID) {
- printf("ERROR: fOutputPID not available\n");
- return;
- }
-
-
- return;
-}
-//___________________________________________________________________________
-void AliAnalysisTaskSEDStarSpectra::UserCreateOutputObjects() {
- // output
- Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
-
- //slot #1
- //OpenFile(1);
- fOutput = new TList();
- fOutput->SetOwner();
- fOutput->SetName("chist0");
-
- fOutputAll = new TList();
- fOutputAll->SetOwner();
- fOutputAll->SetName("listAll");
-
- fOutputPID = new TList();
- fOutputPID->SetOwner();
- fOutputPID->SetName("listPID");
-
- // define histograms
- DefineHistograms();
-
- //Counter for Normalization
- TString normName="NormalizationCounter";
- AliAnalysisDataContainer *cont = GetOutputSlot(4)->GetContainer();
- if(cont)normName=(TString)cont->GetName();
- fCounter = new AliNormalizationCounter(normName.Data());
- fCounter->SetRejectPileUp(fCuts->GetOptPileUp());
-
- PostData(1,fOutput);
- PostData(2,fOutputAll);
- PostData(3,fOutputPID);
-
- return;
-}
-//___________________________________ hiostograms _______________________________________
-void AliAnalysisTaskSEDStarSpectra::DefineHistograms(){
-
- fCEvents = new TH1F("fCEvents","conter",11,0,11);
- fCEvents->SetStats(kTRUE);
- fCEvents->GetXaxis()->SetTitle("1");
- fCEvents->GetYaxis()->SetTitle("counts");
- fOutput->Add(fCEvents);
-
- fTrueDiff2 = new TH2F("DiffDstar_pt","True Reco diff vs pt",200,0,15,900,0,0.3);
- fOutput->Add(fTrueDiff2);
-
- fDeltaMassD1 = new TH1F("DeltaMassD1","delta mass d1",600,0,0.8);
- fOutput->Add(fDeltaMassD1);
-
- const Int_t nhist=14;
- TString nameMass=" ", nameSgn=" ", nameBkg=" ";
-
- for(Int_t i=-2;i<nhist;i++){
- nameMass="histDeltaMass_";
- nameMass+=i+1;
- nameSgn="histDeltaSgn_";
- nameSgn+=i+1;
- nameBkg="histDeltaBkg_";
- nameBkg+=i+1;
-
- if (i==-2) {
- nameMass="histDeltaMass";
- nameSgn="histDeltaSgn";
- nameBkg="histDeltaBkg";
- }
-
- TH1F* spectrumMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} invariant mass; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
- TH1F* spectrumSgn = new TH1F(nameSgn.Data(), "D^{*}-D^{0} Signal invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
- TH1F* spectrumBkg = new TH1F(nameBkg.Data(), "D^{*}-D^{0} Background invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
-
- nameMass="histD0Mass_";
- nameMass+=i+1;
- nameSgn="histD0Sgn_";
- nameSgn+=i+1;
- nameBkg="histD0Bkg_";
- nameBkg+=i+1;
-
- if (i==-2) {
- nameMass="histD0Mass";
- nameSgn="histD0Sgn";
- nameBkg="histD0Bkg";
- }
-
- TH1F* spectrumD0Mass = new TH1F(nameMass.Data(),"D^{0} invariant mass; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
- TH1F* spectrumD0Sgn = new TH1F(nameSgn.Data(), "D^{0} Signal invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
- TH1F* spectrumD0Bkg = new TH1F(nameBkg.Data(), "D^{0} Background invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
-
- nameMass="histDstarMass_";
- nameMass+=i+1;
- nameSgn="histDstarSgn_";
- nameSgn+=i+1;
- nameBkg="histDstarBkg_";
- nameBkg+=i+1;
-
- if (i==-2) {
- nameMass="histDstarMass";
- nameSgn="histDstarSgn";
- nameBkg="histDstarBkg";
- }
-
- TH1F* spectrumDstarMass = new TH1F(nameMass.Data(),"D^{*} invariant mass; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
- TH1F* spectrumDstarSgn = new TH1F(nameSgn.Data(), "D^{*} Signal invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
- TH1F* spectrumDstarBkg = new TH1F(nameBkg.Data(), "D^{*} Background invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
-
- nameMass="histSideBandMass_";
- nameMass+=i+1;
- if (i==-2) {
- nameMass="histSideBandMass";
- }
-
- TH1F* spectrumSideBandMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} sideband mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
-
- nameMass="histWrongSignMass_";
- nameMass+=i+1;
- if (i==-2) {
- nameMass="histWrongSignMass";
- }
-
- TH1F* spectrumWrongSignMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} wrongsign mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
-
-
- spectrumMass->Sumw2();
- spectrumSgn->Sumw2();
- spectrumBkg->Sumw2();
-
- spectrumMass->SetLineColor(6);
- spectrumSgn->SetLineColor(2);
- spectrumBkg->SetLineColor(4);
-
- spectrumMass->SetMarkerStyle(20);
- spectrumSgn->SetMarkerStyle(20);
- spectrumBkg->SetMarkerStyle(20);
- spectrumMass->SetMarkerSize(0.6);
- spectrumSgn->SetMarkerSize(0.6);
- spectrumBkg->SetMarkerSize(0.6);
- spectrumMass->SetMarkerColor(6);
- spectrumSgn->SetMarkerColor(2);
- spectrumBkg->SetMarkerColor(4);
-
- spectrumD0Mass->Sumw2();
- spectrumD0Sgn->Sumw2();
- spectrumD0Bkg->Sumw2();
-
- spectrumD0Mass->SetLineColor(6);
- spectrumD0Sgn->SetLineColor(2);
- spectrumD0Bkg->SetLineColor(4);
-
- spectrumD0Mass->SetMarkerStyle(20);
- spectrumD0Sgn->SetMarkerStyle(20);
- spectrumD0Bkg->SetMarkerStyle(20);
- spectrumD0Mass->SetMarkerSize(0.6);
- spectrumD0Sgn->SetMarkerSize(0.6);
- spectrumD0Bkg->SetMarkerSize(0.6);
- spectrumD0Mass->SetMarkerColor(6);
- spectrumD0Sgn->SetMarkerColor(2);
- spectrumD0Bkg->SetMarkerColor(4);
-
- spectrumDstarMass->Sumw2();
- spectrumDstarSgn->Sumw2();
- spectrumDstarBkg->Sumw2();
-
- spectrumDstarMass->SetLineColor(6);
- spectrumDstarSgn->SetLineColor(2);
- spectrumDstarBkg->SetLineColor(4);
-
- spectrumDstarMass->SetMarkerStyle(20);
- spectrumDstarSgn->SetMarkerStyle(20);
- spectrumDstarBkg->SetMarkerStyle(20);
- spectrumDstarMass->SetMarkerSize(0.6);
- spectrumDstarSgn->SetMarkerSize(0.6);
- spectrumDstarBkg->SetMarkerSize(0.6);
- spectrumDstarMass->SetMarkerColor(6);
- spectrumDstarSgn->SetMarkerColor(2);
- spectrumDstarBkg->SetMarkerColor(4);
-
- spectrumSideBandMass->Sumw2();
- spectrumSideBandMass->SetLineColor(4);
- spectrumSideBandMass->SetMarkerStyle(20);
- spectrumSideBandMass->SetMarkerSize(0.6);
- spectrumSideBandMass->SetMarkerColor(4);
-
- spectrumWrongSignMass->Sumw2();
- spectrumWrongSignMass->SetLineColor(4);
- spectrumWrongSignMass->SetMarkerStyle(20);
- spectrumWrongSignMass->SetMarkerSize(0.6);
- spectrumWrongSignMass->SetMarkerColor(4);
-
- TH1F* allMass = (TH1F*)spectrumMass->Clone();
- TH1F* allSgn = (TH1F*)spectrumSgn->Clone();
- TH1F* allBkg = (TH1F*)spectrumBkg->Clone();
-
- TH1F* pidMass = (TH1F*)spectrumMass->Clone();
- TH1F* pidSgn = (TH1F*)spectrumSgn->Clone();
- TH1F* pidBkg = (TH1F*)spectrumBkg->Clone();
-
- fOutputAll->Add(allMass);
- fOutputAll->Add(allSgn);
- fOutputAll->Add(allBkg);
-
- fOutputPID->Add(pidMass);
- fOutputPID->Add(pidSgn);
- fOutputPID->Add(pidBkg);
-
- TH1F* allD0Mass = (TH1F*)spectrumD0Mass->Clone();
- TH1F* allD0Sgn = (TH1F*)spectrumD0Sgn->Clone();
- TH1F* allD0Bkg = (TH1F*)spectrumD0Bkg->Clone();
-
- TH1F* pidD0Mass = (TH1F*)spectrumD0Mass->Clone();
- TH1F* pidD0Sgn = (TH1F*)spectrumD0Sgn->Clone();
- TH1F* pidD0Bkg = (TH1F*)spectrumD0Bkg->Clone();
-
- fOutputAll->Add(allD0Mass);
- fOutputAll->Add(allD0Sgn);
- fOutputAll->Add(allD0Bkg);
-
- fOutputPID->Add(pidD0Mass);
- fOutputPID->Add(pidD0Sgn);
- fOutputPID->Add(pidD0Bkg);
-
- TH1F* allDstarMass = (TH1F*)spectrumDstarMass->Clone();
- TH1F* allDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
- TH1F* allDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
-
- TH1F* pidDstarMass = (TH1F*)spectrumDstarMass->Clone();
- TH1F* pidDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
- TH1F* pidDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
-
- fOutputAll->Add(allDstarMass);
- fOutputAll->Add(allDstarSgn);
- fOutputAll->Add(allDstarBkg);
-
- fOutputPID->Add(pidDstarMass);
- fOutputPID->Add(pidDstarSgn);
- fOutputPID->Add(pidDstarBkg);
-
- TH1F* allSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
- TH1F* pidSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
-
- fOutputAll->Add(allSideBandMass);
- fOutputPID->Add(pidSideBandMass);
-
- TH1F* allWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
- TH1F* pidWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
-
- fOutputAll->Add(allWrongSignMass);
- fOutputPID->Add(pidWrongSignMass);
-
- }
-
- // pt spectra
- nameMass="ptMass";
- nameSgn="ptSgn";
- nameBkg="ptBkg";
-
- TH1F* ptspectrumMass = new TH1F(nameMass.Data(),"D^{*} p_{T}; p_{T} [GeV]; Entries",200,0,10);
- TH1F* ptspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);
- TH1F* ptspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);
-
- ptspectrumMass->Sumw2();
- ptspectrumSgn->Sumw2();
- ptspectrumBkg->Sumw2();
-
- ptspectrumMass->SetLineColor(6);
- ptspectrumSgn->SetLineColor(2);
- ptspectrumBkg->SetLineColor(4);
-
- ptspectrumMass->SetMarkerStyle(20);
- ptspectrumSgn->SetMarkerStyle(20);
- ptspectrumBkg->SetMarkerStyle(20);
- ptspectrumMass->SetMarkerSize(0.6);
- ptspectrumSgn->SetMarkerSize(0.6);
- ptspectrumBkg->SetMarkerSize(0.6);
- ptspectrumMass->SetMarkerColor(6);
- ptspectrumSgn->SetMarkerColor(2);
- ptspectrumBkg->SetMarkerColor(4);
-
- TH1F* ptallMass = (TH1F*)ptspectrumMass->Clone();
- TH1F* ptallSgn = (TH1F*)ptspectrumSgn->Clone();
- TH1F* ptallBkg = (TH1F*)ptspectrumBkg->Clone();
-
- TH1F* ptpidMass = (TH1F*)ptspectrumMass->Clone();
- TH1F* ptpidSgn = (TH1F*)ptspectrumSgn->Clone();
- TH1F* ptpidBkg = (TH1F*)ptspectrumBkg->Clone();
-
- fOutputAll->Add(ptallMass);
- fOutputAll->Add(ptallSgn);
- fOutputAll->Add(ptallBkg);
-
- fOutputPID->Add(ptpidMass);
- fOutputPID->Add(ptpidSgn);
- fOutputPID->Add(ptpidBkg);
-
- // eta spectra
- nameMass="etaMass";
- nameSgn="etaSgn";
- nameBkg="etaBkg";
-
- TH1F* etaspectrumMass = new TH1F(nameMass.Data(),"D^{*} #eta; #eta; Entries",200,-1,1);
- TH1F* etaspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal #eta - MC; #eta; Entries",200,-1,1);
- TH1F* etaspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background #eta - MC; #eta; Entries",200,-1,1);
-
- etaspectrumMass->Sumw2();
- etaspectrumSgn->Sumw2();
- etaspectrumBkg->Sumw2();
-
- etaspectrumMass->SetLineColor(6);
- etaspectrumSgn->SetLineColor(2);
- etaspectrumBkg->SetLineColor(4);
-
- etaspectrumMass->SetMarkerStyle(20);
- etaspectrumSgn->SetMarkerStyle(20);
- etaspectrumBkg->SetMarkerStyle(20);
- etaspectrumMass->SetMarkerSize(0.6);
- etaspectrumSgn->SetMarkerSize(0.6);
- etaspectrumBkg->SetMarkerSize(0.6);
- etaspectrumMass->SetMarkerColor(6);
- etaspectrumSgn->SetMarkerColor(2);
- etaspectrumBkg->SetMarkerColor(4);
-
- TH1F* etaallMass = (TH1F*)etaspectrumMass->Clone();
- TH1F* etaallSgn = (TH1F*)etaspectrumSgn->Clone();
- TH1F* etaallBkg = (TH1F*)etaspectrumBkg->Clone();
-
- TH1F* etapidMass = (TH1F*)etaspectrumMass->Clone();
- TH1F* etapidSgn = (TH1F*)etaspectrumSgn->Clone();
- TH1F* etapidBkg = (TH1F*)etaspectrumBkg->Clone();
-
- fOutputAll->Add(etaallMass);
- fOutputAll->Add(etaallSgn);
- fOutputAll->Add(etaallBkg);
-
- fOutputPID->Add(etapidMass);
- fOutputPID->Add(etapidSgn);
- fOutputPID->Add(etapidBkg);
-
- return;
-}
-//________________________________________________________________________
-void AliAnalysisTaskSEDStarSpectra::FillSpectrum(AliAODRecoCascadeHF *part, Int_t isDStar, AliRDHFCutsDStartoKpipi *cuts, TList *listout){
- //
- // Fill histos for D* spectrum
- //
-
- Int_t ptbin=cuts->PtBin(part->Pt());
- Double_t invmassD0 = part->InvMassD0();
-
- Double_t pt = part->Pt();
- Double_t eta = part->Eta();
-
- Double_t invmassDelta = part->DeltaInvMass();
- Double_t invmassDstar = part->InvMassDstarKpipi();
-
- TString fillthis="";
- Bool_t massInRange=kFALSE;
-
- Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
- Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
-
- // delta M(Kpipi)-M(Kpi)
- if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<fPeakWindow) massInRange=kTRUE;
-
- if(fUseMCInfo) {
- if(isDStar==1) {
- fillthis="histD0Sgn_";
- fillthis+=ptbin;
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
- fillthis="histD0Sgn";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
- fillthis="histDstarSgn_";
- fillthis+=ptbin;
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
- fillthis="histDstarSgn";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
- fillthis="histDeltaSgn_";
- fillthis+=ptbin;
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
- fillthis="histDeltaSgn";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
- if (massInRange) {
- fillthis="ptSgn";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
- fillthis="etaSgn";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
- }
- }
- else {//background
- fillthis="histD0Bkg_";
- fillthis+=ptbin;
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
- fillthis="histD0Bkg";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
- fillthis="histDstarBkg_";
- fillthis+=ptbin;
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
- fillthis="histDstarBkg";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
- fillthis="histDeltaBkg_";
- fillthis+=ptbin;
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
- fillthis="histDeltaBkg";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
- if (massInRange) {
- fillthis="ptBkg";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
- fillthis="etaBkg";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
- }
- }
- }
- //no MC info, just cut selection
- fillthis="histD0Mass_";
- fillthis+=ptbin;
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
- fillthis="histD0Mass";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
- fillthis="histDstarMass_";
- fillthis+=ptbin;
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
- fillthis="histDstarMass";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
- fillthis="histDeltaMass_";
- fillthis+=ptbin;
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
- fillthis="histDeltaMass";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
-
- if (massInRange) {
- fillthis="ptMass";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
- fillthis="etaMass";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
- }
-
- return;
-}
-//______________________________ side band background for D*___________________________________
-void AliAnalysisTaskSEDStarSpectra::SideBandBackground(AliAODRecoCascadeHF *part, AliRDHFCutsDStartoKpipi *cuts, TList *listout){
-
- // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
- // (expected detector resolution) on the left and right frm the D0 mass. Each band
- // has a width of ~5 sigmas. Two band needed for opening angle considerations
-
- Int_t ptbin=cuts->PtBin(part->Pt());
-
- Bool_t massInRange=kFALSE;
-
- // select the side bands intervall
- Double_t invmassD0 = part->InvMassD0();
- if(TMath::Abs(invmassD0-1.865)>4*fD0Window && TMath::Abs(invmassD0-1.865)<8*fD0Window){
-
- // for pt and eta
- Double_t invmassDelta = part->DeltaInvMass();
- if (TMath::Abs(invmassDelta-0.14557)<fPeakWindow) massInRange=kTRUE;
-
- TString fillthis="";
- fillthis="histSideBandMass_";
- fillthis+=ptbin;
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
- fillthis="histSideBandMass";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
-
- }
-}
-//________________________________________________________________________________________________________________
-void AliAnalysisTaskSEDStarSpectra::WrongSignForDStar(AliAODRecoCascadeHF *part, AliRDHFCutsDStartoKpipi *cuts, TList *listout){
- //
- // assign the wrong charge to the soft pion to create background
- //
- Int_t ptbin=cuts->PtBin(part->Pt());
-
- AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)part->Get2Prong();
-
- Int_t okD0WrongSign,okD0barWrongSign;
- Double_t wrongMassD0=0.;
-
- Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
- if (!isSelected){
- return;
- }
-
- okD0WrongSign = 1;
- okD0barWrongSign = 1;
-
- //if is D*+ than assume D0bar
- if(part->Charge()>0 && (isSelected ==1)) {
- okD0WrongSign = 0;
- }
- if(part->Charge()<0 && (isSelected ==2)){
- okD0barWrongSign = 0;
- }
-
- // assign the wrong mass in case the cuts return both D0 and D0bar
- if(part->Charge()>0 && (isSelected ==3)) {
- okD0WrongSign = 0;
- } else if(part->Charge()<0 && (isSelected ==3)){
- okD0barWrongSign = 0;
- }
-
- //wrong D0 inv mass
- if(okD0WrongSign!=0){
- wrongMassD0 = theD0particle->InvMassD0();
- }else if(okD0WrongSign==0){
- wrongMassD0 = theD0particle->InvMassD0bar();
- }
-
- if(TMath::Abs(wrongMassD0-1.865)<fD0Window){
-
- // wrong D* inv mass
- Double_t e[3];
- if (part->Charge()>0){
- e[0]=theD0particle->EProng(0,321);
- e[1]=theD0particle->EProng(1,211);
- }else{
- e[0]=theD0particle->EProng(0,211);
- e[1]=theD0particle->EProng(1,321);
- }
- e[2]=part->EProng(0,211);
-
- Double_t esum = e[0]+e[1]+e[2];
- Double_t pds = part->P();
-
- Double_t wrongMassDstar = TMath::Sqrt(esum*esum-pds*pds);
-
- TString fillthis="";
- fillthis="histWrongSignMass_";
- fillthis+=ptbin;
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
- fillthis="histWrongSignMass";
- ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
-
- }
-}
+/**************************************************************************
\r * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
\r * *
\r * Author: The ALICE Off-line Project. *
\r * Contributors are mentioned in the code where appropriate. *
\r * *
\r * Permission to use, copy, modify and distribute this software and its *
\r * documentation strictly for non-commercial purposes is hereby granted *
\r * without fee, provided that the above copyright notice appears in all *
\r * copies and that both the copyright notice and this permission notice *
\r * appeuear in the supporting documentation. The authors make no claims *
\r * about the suitability of this software for any purpose. It is *
\r * provided "as is" without express or implied warranty. *
\r **************************************************************************/
\r\r/* $Id$ */
\r\r//
\r//
\r// Base class for DStar Analysis
\r//
\r//
\r// The D* spectra study is done in pt bins:
\r// [0,0.5] [0.5,1] [1,2] [2,3] [3,4] [4,5] [5,6] [6,7] [7,8],
\r// [8,10],[10,12], [12,16], [16,20] and [20,24]
\r//
\r// Cuts arew centralized in AliRDHFCutsDStartoKpipi
\r// Side Band and like sign background are implemented in the macro
\r//
\r//-----------------------------------------------------------------------
\r//
\r// Author A.Grelli
\r// ERC-QGP Utrecht University - a.grelli@uu.nl,
\r// Author Y.Wang
\r// University of Heidelberg - yifei@physi.uni-heidelberg.de
\r// Author C.Ivan
\r// ERC-QGP Utrecht University - c.ivan@uu.nl,
\r//
\r//-----------------------------------------------------------------------
\r\r#include <TSystem.h>
\r#include <TParticle.h>
\r#include <TH1I.h>
\r#include "TROOT.h"
\r#include <TDatabasePDG.h>
\r#include <AliAnalysisDataSlot.h>
\r#include <AliAnalysisDataContainer.h>
\r#include "AliRDHFCutsDStartoKpipi.h"
\r#include "AliStack.h"
\r#include "AliMCEvent.h"
\r#include "AliAnalysisManager.h"
\r#include "AliAODMCHeader.h"
\r#include "AliAODHandler.h"
\r#include "AliLog.h"
\r#include "AliAODVertex.h"
\r#include "AliAODRecoDecay.h"
\r#include "AliAODRecoDecayHF.h"
\r#include "AliAODRecoCascadeHF.h"
\r#include "AliAODRecoDecayHF2Prong.h"
\r#include "AliAnalysisVertexingHF.h"
\r#include "AliESDtrack.h"
\r#include "AliAODMCParticle.h"
\r#include "AliAnalysisTaskSE.h"
\r#include "AliAnalysisTaskSEDStarSpectra.h"
\r#include "AliNormalizationCounter.h"
\r\rClassImp(AliAnalysisTaskSEDStarSpectra)
\r\r//__________________________________________________________________________
\rAliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra():
\r AliAnalysisTaskSE(),
\r fEvents(0),
\r fAnalysis(0),
\r fD0Window(0),
\r fPeakWindow(0),
\r fUseMCInfo(kFALSE),
\r fDoSearch(kFALSE),
\r fOutput(0),
\r fOutputAll(0),
\r fOutputPID(0),
\r fNSigma(3),
\r fCuts(0),
\r fCEvents(0),
\r fTrueDiff2(0),
\r fDeltaMassD1(0),
\r fCounter(0)
\r{
\r //
\r // Default ctor
\r //
\r}
\r//___________________________________________________________________________
\rAliAnalysisTaskSEDStarSpectra::AliAnalysisTaskSEDStarSpectra(const Char_t* name, AliRDHFCutsDStartoKpipi* cuts) :
\r AliAnalysisTaskSE(name),
\r fEvents(0),
\r fAnalysis(0),
\r fD0Window(0),
\r fPeakWindow(0),
\r fUseMCInfo(kFALSE),
\r fDoSearch(kFALSE),
\r fOutput(0),
\r fOutputAll(0),
\r fOutputPID(0),
\r fNSigma(3),
\r fCuts(0),
\r fCEvents(0),
\r fTrueDiff2(0),
\r fDeltaMassD1(0),
\r fCounter(0)
\r{
\r //
\r // Constructor. Initialization of Inputs and Outputs
\r //
\r Info("AliAnalysisTaskSEDStarSpectra","Calling Constructor");
\r\r fCuts=cuts;
\r\r DefineOutput(1,TList::Class()); //conters
\r DefineOutput(2,TList::Class()); //All Entries output
\r DefineOutput(3,TList::Class()); //3sigma PID output
\r DefineOutput(4,AliRDHFCutsDStartoKpipi::Class()); //My private output
\r DefineOutput(5,AliNormalizationCounter::Class()); // normalization
\r}
\r\r//___________________________________________________________________________
\rAliAnalysisTaskSEDStarSpectra::~AliAnalysisTaskSEDStarSpectra() {
\r //
\r // destructor
\r //
\r Info("~AliAnalysisTaskSEDStarSpectra","Calling Destructor");
\r \r if (fOutput) {
\r delete fOutput;
\r fOutput = 0;
\r }
\r if (fOutputAll) {
\r delete fOutputAll;
\r fOutputAll = 0;
\r }
\r if (fOutputPID) {
\r delete fOutputPID;
\r fOutputPID = 0;
\r }
\r if (fCuts) {
\r delete fCuts;
\r fCuts = 0;
\r }
\r if(fCEvents){
\r delete fCEvents;
\r fCEvents =0;
\r }
\r if(fDeltaMassD1){
\r delete fDeltaMassD1;
\r fDeltaMassD1 =0;
\r }
\r}
\r//_________________________________________________
\rvoid AliAnalysisTaskSEDStarSpectra::Init(){
\r //
\r // Initialization
\r //
\r\r if(fDebug > 1) printf("AnalysisTaskSEDStarSpectra::Init() \n");
\r AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
\r // Post the data
\r PostData(4,copyfCuts);
\r\r return;
\r}
\r\r//_________________________________________________
\rvoid AliAnalysisTaskSEDStarSpectra::UserExec(Option_t *)
\r{
\r // user exec
\r if (!fInputEvent) {
\r Error("UserExec","NO EVENT FOUND!");
\r return;
\r }
\r\r fEvents++;
\r\r AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
\r TClonesArray *arrayDStartoD0pi=0;
\r\r fCEvents->Fill(1);
\r\r if(!aodEvent && AODEvent() && IsStandardAOD()) {
\r // In case there is an AOD handler writing a standard AOD, use the AOD
\r // event in memory rather than the input (ESD) event.
\r aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
\r // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
\r // have to taken from the AOD event hold by the AliAODExtension
\r AliAODHandler* aodHandler = (AliAODHandler*)
\r ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
\r if(aodHandler->GetExtensions()) {
\r AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
\r AliAODEvent *aodFromExt = ext->GetAOD();
\r arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
\r }
\r } else {
\r arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
\r }
\r\r // fix for temporary bug in ESDfilter
\r // the AODs with null vertex pointer didn't pass the PhysSel
\r if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
\r \r fCEvents->Fill(2);
\r\r fCounter->StoreEvent(aodEvent,fCuts,fUseMCInfo);
\r\r if(!fCuts->IsEventSelected(aodEvent)) {
\r if(fCuts->GetWhyRejection()==6) // rejected for Z vertex
\r fCEvents->Fill(6);
\r return;
\r }
\r\r // Load the event
\r AliInfo(Form("Event %d",fEvents));
\r if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
\r\r // counters for efficiencies
\r Int_t icountReco = 0;
\r \r //D* and D0 prongs needed to MatchToMC method
\r Int_t pdgDgDStartoD0pi[2]={421,211};
\r Int_t pdgDgD0toKpi[2]={321,211};
\r\r // AOD primary vertex
\r AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
\r if(!vtx1) return;
\r if(vtx1->GetNContributors()<1) return;
\r\r if (!arrayDStartoD0pi){
\r AliInfo("Could not find array of HF vertices, skipping the event");
\r return;
\r }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
\r\r Int_t nSelectedAna =0;
\r Int_t nSelectedProd =0;
\r\r // loop over the tracks to search for candidates soft pion
\r for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
\r\r // D* candidates and D0 from D*
\r AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
\r if(!dstarD0pi->GetSecondaryVtx()) continue;
\r AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
\r if (!theD0particle) continue;
\r \r Int_t isDStar = 0;
\r TClonesArray *mcArray = 0; // fix coverity
\r\r // mc analysis
\r if(fUseMCInfo){
\r //MC array need for maching
\r mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
\r if (!mcArray) {
\r AliError("Could not find Monte-Carlo in AOD");
\r return;
\r }
\r // find associated MC particle for D* ->D0toKpi
\r Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
\r if(mcLabel>=0) isDStar = 1;
\r }
\r \r Int_t ptbin=fCuts->PtBin(dstarD0pi->Pt());
\r \r // set the D0 search window bin by bin - useful to calculate side band bkg
\r if (ptbin==0){
\r if(fAnalysis==1){
\r fD0Window=0.035;
\r fPeakWindow=0.03;
\r }else{
\r fD0Window=0.020;
\r fPeakWindow=0.0018;
\r }
\r }
\r if (ptbin==1){
\r if(fAnalysis==1){
\r fD0Window=0.035;
\r fPeakWindow=0.03;
\r }else{
\r fD0Window=0.020;
\r fPeakWindow=0.0018;
\r }
\r }
\r if (ptbin==2){
\r if(fAnalysis==1){
\r fD0Window=0.035;
\r fPeakWindow=0.03;
\r }else{
\r fD0Window=0.020;
\r fPeakWindow=0.0018;
\r }
\r }
\r if (ptbin==3){
\r if(fAnalysis==1){
\r fD0Window=0.035;
\r fPeakWindow=0.03;
\r }else{
\r fD0Window=0.022;
\r fPeakWindow=0.0016;
\r }
\r }
\r if (ptbin==4){
\r if(fAnalysis==1){
\r fD0Window=0.035;
\r fPeakWindow=0.03;
\r }else{
\r fD0Window=0.026;
\r fPeakWindow=0.0014;
\r }
\r }
\r if (ptbin==5){
\r if(fAnalysis==1){
\r fD0Window=0.045;
\r fPeakWindow=0.03;
\r }else{
\r fD0Window=0.026;
\r fPeakWindow=0.0014;
\r }
\r }
\r if (ptbin==6){
\r if(fAnalysis==1){
\r fD0Window=0.045;
\r fPeakWindow=0.03;
\r }else{
\r fD0Window=0.026;
\r fPeakWindow=0.006;
\r }
\r }
\r if (ptbin==7){
\r if(fAnalysis==1){
\r fD0Window=0.055;
\r fPeakWindow=0.03;
\r }else{
\r fD0Window=0.026;
\r fPeakWindow=0.006;
\r }
\r }
\r if (ptbin>7){
\r if(fAnalysis==1){
\r fD0Window=0.074;
\r fPeakWindow=0.03;
\r }else{
\r fD0Window=0.026;
\r fPeakWindow=0.006;
\r }
\r }
\r fCEvents->Fill(9);
\r \r nSelectedProd++;
\r nSelectedAna++;
\r \r Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
\r \r // here for lead-lead
\r Double_t invmassD0 = dstarD0pi->InvMassD0();
\r if (TMath::Abs(invmassD0-mPDGD0)>fD0Window) continue;
\r \r Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
\r Double_t invmassDelta = dstarD0pi->DeltaInvMass();
\r \r if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>fPeakWindow) continue;
\r \r Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
\r if(!isTkSelected) continue;
\r \r if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
\r \r Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
\r if (!isSelected) continue;
\r\r // fill PID
\r FillSpectrum(dstarD0pi,isDStar,fCuts,fOutputPID);
\r SideBandBackground(dstarD0pi,fCuts,fOutputPID);
\r WrongSignForDStar(dstarD0pi,fCuts,fOutputPID);
\r\r // fill no pid
\r fCuts->SetUsePID(kFALSE);
\r FillSpectrum(dstarD0pi,isDStar,fCuts,fOutputAll);
\r SideBandBackground(dstarD0pi,fCuts,fOutputAll);
\r WrongSignForDStar(dstarD0pi,fCuts,fOutputAll);
\r fCuts->SetUsePID(kTRUE);
\r\r // rare D search ------
\r if(fDoSearch){
\r TLorentzVector LorentzTrack1(0,0,0,0); // lorentz 4 vector
\r TLorentzVector LorentzTrack2(0,0,0,0); // lorentz 4 vector
\r \r for (Int_t i=0; i<aodEvent->GetNTracks(); i++){
\r \r AliAODTrack* aodTrack = aodEvent->GetTrack(i);
\r \r if(dstarD0pi->Charge() == aodTrack->Charge()) continue;
\r if((!(aodTrack->GetStatus()&AliESDtrack::kITSrefit)|| (!(aodTrack->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
\r if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>0.02) continue;
\r \r //build the D1 mass
\r Double_t mass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
\r\r LorentzTrack1.SetPxPyPzE( dstarD0pi->Px(),dstarD0pi->Py(), dstarD0pi->Pz(), dstarD0pi->E(413) );
\r LorentzTrack2.SetPxPyPzE( aodTrack->Px(),aodTrack->Py(), aodTrack->Pz(),aodTrack->E(mass) );
\r \r //D1 mass
\r Double_t d1mass = ((LorentzTrack1+LorentzTrack2).M());
\r //mass difference - at 0.4117 and 0.4566
\r fDeltaMassD1->Fill(d1mass-dstarD0pi->InvMassDstarKpipi());
\r }
\r }
\r\r if(isDStar == 1) {
\r fTrueDiff2->Fill(dstarD0pi->Pt(),dstarD0pi->DeltaInvMass());
\r }
\r \r }
\r \r fCounter->StoreCandidates(aodEvent,nSelectedProd,kTRUE);
\r fCounter->StoreCandidates(aodEvent,nSelectedAna,kFALSE);
\r\r AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
\r \r PostData(1,fOutput);
\r PostData(2,fOutputAll);
\r PostData(3,fOutputPID);
\r PostData(5,fCounter);
\r\r}
\r//________________________________________ terminate ___________________________
\rvoid AliAnalysisTaskSEDStarSpectra::Terminate(Option_t*)
\r{
\r // The Terminate() function is the last function to be called during
\r // a query. It always runs on the client, it can be used to present
\r // the results graphically or save the results to file.
\r \r //Info("Terminate","");
\r AliAnalysisTaskSE::Terminate();
\r \r fOutput = dynamic_cast<TList*> (GetOutputData(1));
\r if (!fOutput) {
\r printf("ERROR: fOutput not available\n");
\r return;
\r }
\r \r fCEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fCEvents"));
\r fDeltaMassD1 = dynamic_cast<TH1F*>(fOutput->FindObject("fDeltaMassD1"));
\r fTrueDiff2 = dynamic_cast<TH2F*>(fOutput->FindObject("fTrueDiff2"));
\r\r fOutputAll = dynamic_cast<TList*> (GetOutputData(1));
\r if (!fOutputAll) {
\r printf("ERROR: fOutputAll not available\n");
\r return;
\r }
\r fOutputPID = dynamic_cast<TList*> (GetOutputData(2));
\r if (!fOutputPID) {
\r printf("ERROR: fOutputPID not available\n");
\r return;
\r }
\r \r\r return;
\r}
\r//___________________________________________________________________________
\rvoid AliAnalysisTaskSEDStarSpectra::UserCreateOutputObjects() {
\r // output
\r Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
\r \r //slot #1
\r //OpenFile(1);
\r fOutput = new TList();
\r fOutput->SetOwner();
\r fOutput->SetName("chist0");
\r\r fOutputAll = new TList();
\r fOutputAll->SetOwner();
\r fOutputAll->SetName("listAll");
\r\r fOutputPID = new TList();
\r fOutputPID->SetOwner();
\r fOutputPID->SetName("listPID");
\r \r // define histograms
\r DefineHistograms();
\r\r //Counter for Normalization
\r fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
\r\r PostData(1,fOutput);
\r PostData(2,fOutputAll);
\r PostData(3,fOutputPID);
\r\r return;
\r}
\r//___________________________________ hiostograms _______________________________________
\rvoid AliAnalysisTaskSEDStarSpectra::DefineHistograms(){
\r\r fCEvents = new TH1F("fCEvents","conter",11,0,11);
\r fCEvents->SetStats(kTRUE);
\r fCEvents->GetXaxis()->SetTitle("1");
\r fCEvents->GetYaxis()->SetTitle("counts");
\r fOutput->Add(fCEvents);
\r\r fTrueDiff2 = new TH2F("DiffDstar_pt","True Reco diff vs pt",200,0,15,900,0,0.3);
\r fOutput->Add(fTrueDiff2);
\r\r fDeltaMassD1 = new TH1F("DeltaMassD1","delta mass d1",600,0,0.8);
\r fOutput->Add(fDeltaMassD1);
\r\r const Int_t nhist=14;
\r TString nameMass=" ", nameSgn=" ", nameBkg=" ";
\r\r for(Int_t i=-2;i<nhist;i++){
\r nameMass="histDeltaMass_";
\r nameMass+=i+1;
\r nameSgn="histDeltaSgn_";
\r nameSgn+=i+1;
\r nameBkg="histDeltaBkg_";
\r nameBkg+=i+1;
\r \r if (i==-2) {
\r nameMass="histDeltaMass";
\r nameSgn="histDeltaSgn";
\r nameBkg="histDeltaBkg";
\r }
\r \r TH1F* spectrumMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} invariant mass; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
\r TH1F* spectrumSgn = new TH1F(nameSgn.Data(), "D^{*}-D^{0} Signal invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
\r TH1F* spectrumBkg = new TH1F(nameBkg.Data(), "D^{*}-D^{0} Background invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",200,0.1,0.2);
\r \r nameMass="histD0Mass_";
\r nameMass+=i+1;
\r nameSgn="histD0Sgn_";
\r nameSgn+=i+1;
\r nameBkg="histD0Bkg_";
\r nameBkg+=i+1;
\r \r if (i==-2) {
\r nameMass="histD0Mass";
\r nameSgn="histD0Sgn";
\r nameBkg="histD0Bkg";
\r }
\r\r TH1F* spectrumD0Mass = new TH1F(nameMass.Data(),"D^{0} invariant mass; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
\r TH1F* spectrumD0Sgn = new TH1F(nameSgn.Data(), "D^{0} Signal invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
\r TH1F* spectrumD0Bkg = new TH1F(nameBkg.Data(), "D^{0} Background invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
\r\r nameMass="histDstarMass_";
\r nameMass+=i+1;
\r nameSgn="histDstarSgn_";
\r nameSgn+=i+1;
\r nameBkg="histDstarBkg_";
\r nameBkg+=i+1;
\r\r if (i==-2) {
\r nameMass="histDstarMass";
\r nameSgn="histDstarSgn";
\r nameBkg="histDstarBkg";
\r }
\r\r TH1F* spectrumDstarMass = new TH1F(nameMass.Data(),"D^{*} invariant mass; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
\r TH1F* spectrumDstarSgn = new TH1F(nameSgn.Data(), "D^{*} Signal invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
\r TH1F* spectrumDstarBkg = new TH1F(nameBkg.Data(), "D^{*} Background invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
\r\r nameMass="histSideBandMass_";
\r nameMass+=i+1;
\r if (i==-2) {
\r nameMass="histSideBandMass";
\r }
\r \r TH1F* spectrumSideBandMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} sideband mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
\r\r nameMass="histWrongSignMass_";
\r nameMass+=i+1;
\r if (i==-2) {
\r nameMass="histWrongSignMass";
\r }
\r \r TH1F* spectrumWrongSignMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} wrongsign mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
\r\r\r spectrumMass->Sumw2();
\r spectrumSgn->Sumw2();
\r spectrumBkg->Sumw2();
\r \r spectrumMass->SetLineColor(6);
\r spectrumSgn->SetLineColor(2);
\r spectrumBkg->SetLineColor(4);
\r \r spectrumMass->SetMarkerStyle(20);
\r spectrumSgn->SetMarkerStyle(20);
\r spectrumBkg->SetMarkerStyle(20);
\r spectrumMass->SetMarkerSize(0.6);
\r spectrumSgn->SetMarkerSize(0.6);
\r spectrumBkg->SetMarkerSize(0.6);
\r spectrumMass->SetMarkerColor(6);
\r spectrumSgn->SetMarkerColor(2);
\r spectrumBkg->SetMarkerColor(4);
\r\r spectrumD0Mass->Sumw2();
\r spectrumD0Sgn->Sumw2();
\r spectrumD0Bkg->Sumw2();
\r\r spectrumD0Mass->SetLineColor(6);
\r spectrumD0Sgn->SetLineColor(2);
\r spectrumD0Bkg->SetLineColor(4);
\r\r spectrumD0Mass->SetMarkerStyle(20);
\r spectrumD0Sgn->SetMarkerStyle(20);
\r spectrumD0Bkg->SetMarkerStyle(20);
\r spectrumD0Mass->SetMarkerSize(0.6);
\r spectrumD0Sgn->SetMarkerSize(0.6);
\r spectrumD0Bkg->SetMarkerSize(0.6);
\r spectrumD0Mass->SetMarkerColor(6);
\r spectrumD0Sgn->SetMarkerColor(2);
\r spectrumD0Bkg->SetMarkerColor(4);
\r\r spectrumDstarMass->Sumw2();
\r spectrumDstarSgn->Sumw2();
\r spectrumDstarBkg->Sumw2();
\r\r spectrumDstarMass->SetLineColor(6);
\r spectrumDstarSgn->SetLineColor(2);
\r spectrumDstarBkg->SetLineColor(4);
\r\r spectrumDstarMass->SetMarkerStyle(20);
\r spectrumDstarSgn->SetMarkerStyle(20);
\r spectrumDstarBkg->SetMarkerStyle(20);
\r spectrumDstarMass->SetMarkerSize(0.6);
\r spectrumDstarSgn->SetMarkerSize(0.6);
\r spectrumDstarBkg->SetMarkerSize(0.6);
\r spectrumDstarMass->SetMarkerColor(6);
\r spectrumDstarSgn->SetMarkerColor(2);
\r spectrumDstarBkg->SetMarkerColor(4);
\r\r spectrumSideBandMass->Sumw2();
\r spectrumSideBandMass->SetLineColor(4);
\r spectrumSideBandMass->SetMarkerStyle(20);
\r spectrumSideBandMass->SetMarkerSize(0.6);
\r spectrumSideBandMass->SetMarkerColor(4);
\r\r spectrumWrongSignMass->Sumw2();
\r spectrumWrongSignMass->SetLineColor(4);
\r spectrumWrongSignMass->SetMarkerStyle(20);
\r spectrumWrongSignMass->SetMarkerSize(0.6);
\r spectrumWrongSignMass->SetMarkerColor(4);
\r\r TH1F* allMass = (TH1F*)spectrumMass->Clone();
\r TH1F* allSgn = (TH1F*)spectrumSgn->Clone();
\r TH1F* allBkg = (TH1F*)spectrumBkg->Clone();
\r\r TH1F* pidMass = (TH1F*)spectrumMass->Clone();
\r TH1F* pidSgn = (TH1F*)spectrumSgn->Clone();
\r TH1F* pidBkg = (TH1F*)spectrumBkg->Clone();
\r\r fOutputAll->Add(allMass);
\r fOutputAll->Add(allSgn);
\r fOutputAll->Add(allBkg);
\r\r fOutputPID->Add(pidMass);
\r fOutputPID->Add(pidSgn);
\r fOutputPID->Add(pidBkg);
\r\r TH1F* allD0Mass = (TH1F*)spectrumD0Mass->Clone();
\r TH1F* allD0Sgn = (TH1F*)spectrumD0Sgn->Clone();
\r TH1F* allD0Bkg = (TH1F*)spectrumD0Bkg->Clone();
\r\r TH1F* pidD0Mass = (TH1F*)spectrumD0Mass->Clone();
\r TH1F* pidD0Sgn = (TH1F*)spectrumD0Sgn->Clone();
\r TH1F* pidD0Bkg = (TH1F*)spectrumD0Bkg->Clone();
\r\r fOutputAll->Add(allD0Mass);
\r fOutputAll->Add(allD0Sgn);
\r fOutputAll->Add(allD0Bkg);
\r\r fOutputPID->Add(pidD0Mass);
\r fOutputPID->Add(pidD0Sgn);
\r fOutputPID->Add(pidD0Bkg);
\r \r TH1F* allDstarMass = (TH1F*)spectrumDstarMass->Clone();
\r TH1F* allDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
\r TH1F* allDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
\r \r TH1F* pidDstarMass = (TH1F*)spectrumDstarMass->Clone();
\r TH1F* pidDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
\r TH1F* pidDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
\r \r fOutputAll->Add(allDstarMass);
\r fOutputAll->Add(allDstarSgn);
\r fOutputAll->Add(allDstarBkg);
\r\r fOutputPID->Add(pidDstarMass);
\r fOutputPID->Add(pidDstarSgn);
\r fOutputPID->Add(pidDstarBkg);
\r\r TH1F* allSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
\r TH1F* pidSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
\r\r fOutputAll->Add(allSideBandMass);
\r fOutputPID->Add(pidSideBandMass);
\r \r TH1F* allWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
\r TH1F* pidWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
\r\r fOutputAll->Add(allWrongSignMass);
\r fOutputPID->Add(pidWrongSignMass);
\r \r }
\r \r // pt spectra
\r nameMass="ptMass";
\r nameSgn="ptSgn";
\r nameBkg="ptBkg";
\r \r TH1F* ptspectrumMass = new TH1F(nameMass.Data(),"D^{*} p_{T}; p_{T} [GeV]; Entries",200,0,10);
\r TH1F* ptspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);
\r TH1F* ptspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background p_{T} - MC; p_{T} [GeV]; Entries",200,0,10);
\r \r ptspectrumMass->Sumw2();
\r ptspectrumSgn->Sumw2();
\r ptspectrumBkg->Sumw2();
\r \r ptspectrumMass->SetLineColor(6);
\r ptspectrumSgn->SetLineColor(2);
\r ptspectrumBkg->SetLineColor(4);
\r \r ptspectrumMass->SetMarkerStyle(20);
\r ptspectrumSgn->SetMarkerStyle(20);
\r ptspectrumBkg->SetMarkerStyle(20);
\r ptspectrumMass->SetMarkerSize(0.6);
\r ptspectrumSgn->SetMarkerSize(0.6);
\r ptspectrumBkg->SetMarkerSize(0.6);
\r ptspectrumMass->SetMarkerColor(6);
\r ptspectrumSgn->SetMarkerColor(2);
\r ptspectrumBkg->SetMarkerColor(4);
\r \r TH1F* ptallMass = (TH1F*)ptspectrumMass->Clone();
\r TH1F* ptallSgn = (TH1F*)ptspectrumSgn->Clone();
\r TH1F* ptallBkg = (TH1F*)ptspectrumBkg->Clone();
\r \r TH1F* ptpidMass = (TH1F*)ptspectrumMass->Clone();
\r TH1F* ptpidSgn = (TH1F*)ptspectrumSgn->Clone();
\r TH1F* ptpidBkg = (TH1F*)ptspectrumBkg->Clone();
\r \r fOutputAll->Add(ptallMass);
\r fOutputAll->Add(ptallSgn);
\r fOutputAll->Add(ptallBkg);
\r \r fOutputPID->Add(ptpidMass);
\r fOutputPID->Add(ptpidSgn);
\r fOutputPID->Add(ptpidBkg);
\r \r // eta spectra
\r nameMass="etaMass";
\r nameSgn="etaSgn";
\r nameBkg="etaBkg";
\r \r TH1F* etaspectrumMass = new TH1F(nameMass.Data(),"D^{*} #eta; #eta; Entries",200,-1,1);
\r TH1F* etaspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal #eta - MC; #eta; Entries",200,-1,1);
\r TH1F* etaspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background #eta - MC; #eta; Entries",200,-1,1);
\r \r etaspectrumMass->Sumw2();
\r etaspectrumSgn->Sumw2();
\r etaspectrumBkg->Sumw2();
\r \r etaspectrumMass->SetLineColor(6);
\r etaspectrumSgn->SetLineColor(2);
\r etaspectrumBkg->SetLineColor(4);
\r \r etaspectrumMass->SetMarkerStyle(20);
\r etaspectrumSgn->SetMarkerStyle(20);
\r etaspectrumBkg->SetMarkerStyle(20);
\r etaspectrumMass->SetMarkerSize(0.6);
\r etaspectrumSgn->SetMarkerSize(0.6);
\r etaspectrumBkg->SetMarkerSize(0.6);
\r etaspectrumMass->SetMarkerColor(6);
\r etaspectrumSgn->SetMarkerColor(2);
\r etaspectrumBkg->SetMarkerColor(4);
\r \r TH1F* etaallMass = (TH1F*)etaspectrumMass->Clone();
\r TH1F* etaallSgn = (TH1F*)etaspectrumSgn->Clone();
\r TH1F* etaallBkg = (TH1F*)etaspectrumBkg->Clone();
\r \r TH1F* etapidMass = (TH1F*)etaspectrumMass->Clone();
\r TH1F* etapidSgn = (TH1F*)etaspectrumSgn->Clone();
\r TH1F* etapidBkg = (TH1F*)etaspectrumBkg->Clone();
\r \r fOutputAll->Add(etaallMass);
\r fOutputAll->Add(etaallSgn);
\r fOutputAll->Add(etaallBkg);
\r \r fOutputPID->Add(etapidMass);
\r fOutputPID->Add(etapidSgn);
\r fOutputPID->Add(etapidBkg);
\r \r return;
\r}
\r//________________________________________________________________________
\rvoid AliAnalysisTaskSEDStarSpectra::FillSpectrum(AliAODRecoCascadeHF *part, Int_t isDStar, AliRDHFCutsDStartoKpipi *cuts, TList *listout){
\r //
\r // Fill histos for D* spectrum
\r //
\r\r Int_t ptbin=cuts->PtBin(part->Pt());
\r Double_t invmassD0 = part->InvMassD0();
\r \r Double_t pt = part->Pt();
\r Double_t eta = part->Eta();
\r \r Double_t invmassDelta = part->DeltaInvMass();
\r Double_t invmassDstar = part->InvMassDstarKpipi();
\r \r TString fillthis="";
\r Bool_t massInRange=kFALSE;
\r \r Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
\r Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
\r \r // delta M(Kpipi)-M(Kpi)
\r if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<fPeakWindow) massInRange=kTRUE;
\r \r if(fUseMCInfo) {
\r if(isDStar==1) {
\r fillthis="histD0Sgn_";
\r fillthis+=ptbin;
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
\r fillthis="histD0Sgn";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
\r fillthis="histDstarSgn_";
\r fillthis+=ptbin;
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
\r fillthis="histDstarSgn";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
\r fillthis="histDeltaSgn_";
\r fillthis+=ptbin;
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
\r fillthis="histDeltaSgn";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
\r if (massInRange) {
\r fillthis="ptSgn";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
\r fillthis="etaSgn";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
\r }
\r }
\r else {//background
\r fillthis="histD0Bkg_";
\r fillthis+=ptbin;
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
\r fillthis="histD0Bkg";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
\r fillthis="histDstarBkg_";
\r fillthis+=ptbin;
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
\r fillthis="histDstarBkg";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
\r fillthis="histDeltaBkg_";
\r fillthis+=ptbin;
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
\r fillthis="histDeltaBkg";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
\r if (massInRange) {
\r fillthis="ptBkg";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
\r fillthis="etaBkg";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
\r }
\r }
\r }
\r //no MC info, just cut selection
\r fillthis="histD0Mass_";
\r fillthis+=ptbin;
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
\r fillthis="histD0Mass";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
\r fillthis="histDstarMass_";
\r fillthis+=ptbin;
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
\r fillthis="histDstarMass";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDstar);
\r fillthis="histDeltaMass_";
\r fillthis+=ptbin;
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
\r fillthis="histDeltaMass";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
\r \r if (massInRange) {
\r fillthis="ptMass";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(pt);
\r fillthis="etaMass";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(eta);
\r }
\r \r return;
\r}
\r//______________________________ side band background for D*___________________________________
\rvoid AliAnalysisTaskSEDStarSpectra::SideBandBackground(AliAODRecoCascadeHF *part, AliRDHFCutsDStartoKpipi *cuts, TList *listout){
\r\r // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
\r // (expected detector resolution) on the left and right frm the D0 mass. Each band
\r // has a width of ~5 sigmas. Two band needed for opening angle considerations
\r\r Int_t ptbin=cuts->PtBin(part->Pt());
\r \r Bool_t massInRange=kFALSE;
\r \r // select the side bands intervall
\r Double_t invmassD0 = part->InvMassD0();
\r if(TMath::Abs(invmassD0-1.865)>4*fD0Window && TMath::Abs(invmassD0-1.865)<8*fD0Window){
\r \r // for pt and eta
\r Double_t invmassDelta = part->DeltaInvMass();
\r if (TMath::Abs(invmassDelta-0.14557)<fPeakWindow) massInRange=kTRUE;
\r \r TString fillthis="";
\r fillthis="histSideBandMass_";
\r fillthis+=ptbin;
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
\r fillthis="histSideBandMass";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassDelta);
\r \r }
\r}
\r//________________________________________________________________________________________________________________
\rvoid AliAnalysisTaskSEDStarSpectra::WrongSignForDStar(AliAODRecoCascadeHF *part, AliRDHFCutsDStartoKpipi *cuts, TList *listout){
\r //
\r // assign the wrong charge to the soft pion to create background
\r //
\r Int_t ptbin=cuts->PtBin(part->Pt());
\r \r AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)part->Get2Prong();
\r\r Int_t okD0WrongSign,okD0barWrongSign;
\r Double_t wrongMassD0=0.;
\r \r Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
\r if (!isSelected){
\r return;
\r }
\r\r okD0WrongSign = 1;
\r okD0barWrongSign = 1;
\r \r //if is D*+ than assume D0bar
\r if(part->Charge()>0 && (isSelected ==1)) {
\r okD0WrongSign = 0;
\r }
\r if(part->Charge()<0 && (isSelected ==2)){
\r okD0barWrongSign = 0;
\r }
\r \r // assign the wrong mass in case the cuts return both D0 and D0bar
\r if(part->Charge()>0 && (isSelected ==3)) {
\r okD0WrongSign = 0;
\r } else if(part->Charge()<0 && (isSelected ==3)){
\r okD0barWrongSign = 0;
\r }
\r \r //wrong D0 inv mass
\r if(okD0WrongSign!=0){
\r wrongMassD0 = theD0particle->InvMassD0();
\r }else if(okD0WrongSign==0){
\r wrongMassD0 = theD0particle->InvMassD0bar();
\r }
\r \r if(TMath::Abs(wrongMassD0-1.865)<fD0Window){
\r \r // wrong D* inv mass
\r Double_t e[3];
\r if (part->Charge()>0){
\r e[0]=theD0particle->EProng(0,321);
\r e[1]=theD0particle->EProng(1,211);
\r }else{
\r e[0]=theD0particle->EProng(0,211);
\r e[1]=theD0particle->EProng(1,321);
\r }
\r e[2]=part->EProng(0,211);
\r \r Double_t esum = e[0]+e[1]+e[2];
\r Double_t pds = part->P();
\r\r Double_t wrongMassDstar = TMath::Sqrt(esum*esum-pds*pds);
\r\r TString fillthis="";
\r fillthis="histWrongSignMass_";
\r fillthis+=ptbin;
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
\r fillthis="histWrongSignMass";
\r ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
\r \r }
\r}
\r\ No newline at end of file
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff