Changes requested in report #61429: PID: Separating response functions from ESD ...

[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDStar.cxx

/**************************************************************************\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
 * appear 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
/////////////////////////////////////////////////////////////\r
//\r
// AliAnalysisTaskSE for D* candidates invariant mass histogram\r
// and comparison with the MC truth.\r
//\r
// Authors: Y.Wang, yifei@physi.uni-heidelberg.de\r
/////////////////////////////////////////////////////////////\r
\r
#include <Riostream.h>\r
#include <TClonesArray.h>\r
#include <TList.h>\r
#include <TH1F.h>\r
\r
\r
#include "AliPID.h"\r
#include "AliTPCPIDResponse.h"\r
#include "AliAnalysisManager.h"\r
#include "AliAODHandler.h"\r
#include "AliESDtrack.h"\r
#include "AliAODEvent.h"\r
#include "AliAODVertex.h"\r
#include "AliAODTrack.h"\r
#include "AliAODMCHeader.h"\r
#include "AliAODMCParticle.h"\r
#include "AliAODRecoDecayHF2Prong.h"\r
#include "AliAODRecoCascadeHF.h"\r
#include "AliAnalysisVertexingHF.h"\r
#include "AliAnalysisTaskSE.h"\r
#include "AliAnalysisTaskSEDStar.h"\r
\r
ClassImp(AliAnalysisTaskSEDStar)\r
\r
\r
//________________________________________________________________________\r
AliAnalysisTaskSEDStar::AliAnalysisTaskSEDStar():\r
AliAnalysisTaskSE(),\r
fOutput(0), \r
fhistMass(0),\r
fhistSgn(0),\r
fhistBkg(0),\r
fReadMC(0),\r
fPID(1),\r
fNSigma(3),\r
fVHF(0)\r
{\r
  // Default constructor\r
}\r
\r
//________________________________________________________________________\r
AliAnalysisTaskSEDStar::AliAnalysisTaskSEDStar(const char *name):\r
AliAnalysisTaskSE(name),\r
fOutput(0), \r
fhistMass(0),\r
fhistSgn(0),\r
fhistBkg(0),\r
fReadMC(0),\r
fPID(1),\r
fNSigma(3),\r
fVHF(0)\r
{\r
  // Default constructor\r
\r
  // Output slot #1 writes into a TList container\r
  DefineOutput(1,TList::Class());  //My private output\r
}\r
\r
//________________________________________________________________________\r
AliAnalysisTaskSEDStar::~AliAnalysisTaskSEDStar()\r
{\r
  // Destructor\r
\r
  if (fhistMass){\r
    delete fhistMass;\r
    fhistMass=0;\r
  }\r
\r
  if (fhistSgn){\r
    delete fhistSgn;\r
    fhistSgn=0;\r
  }\r
\r
  if (fhistBkg){\r
    delete fhistBkg;\r
    fhistBkg=0;\r
  }\r
\r
  if (fOutput) {\r
    delete fOutput;\r
    fOutput = 0;\r
  }\r
  if (fVHF) {\r
    delete fVHF;\r
    fVHF = 0;\r
  }\r
\r
}  \r
\r
//________________________________________________________________________\r
void AliAnalysisTaskSEDStar::Init()\r
{\r
  // Initialization\r
\r
  if(fDebug > 1) printf("AnalysisTaskSEDStar::Init() \n");\r
\r
//  gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/ConfigVertexingHF.C");\r
  gROOT->LoadMacro("ConfigVertexingHF.C");\r
\r
  fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");  \r
  // set dedidcated cuts\r
  //fVHF->SetD0fromDstarCuts(0.7,0.03,0.8,0.06,0.06,0.05,0.05,-0.0002,0.6); //2.p-p vertex reconstructed\r
  fVHF->SetD0fromDstarCuts(0.3,999999.,1.1,0.,0.,999999.,999999.,999999.,0.);\r
  fVHF->SetDstarCuts(0.3, 0.1, 0.05, 100000000000.0, 0.5);\r
  fVHF->PrintStatus();\r
\r
  return;\r
}\r
\r
//________________________________________________________________________\r
void AliAnalysisTaskSEDStar::UserCreateOutputObjects()\r
{\r
\r
  // Create the output container\r
  //\r
  if(fDebug > 1) printf("AnalysisTaskSEDStar::UserCreateOutputObjects() \n");\r
\r
  // Several histograms are more conveniently managed in a TList\r
  fOutput = new TList();\r
  fOutput->SetOwner();\r
\r
    fhistMass = new TH1F("fhistMass","D^{*+}-D^{0} invariant mass; M [GeV]; Entries",200,0,0.3);\r
    fhistSgn = new TH1F("fhistSgn", "D^{*+}-D^{0} invariant mass - MC; M [GeV]; Entries",200,0,0.3);\r
    fhistBkg = new TH1F("fhistBkg", "D^{*+}-Background invariant mass - MC; M [GeV]; Entries",200,0,0.3);\r
 \r
  fOutput->Add(fhistMass);\r
  fOutput->Add(fhistSgn);\r
  fOutput->Add(fhistBkg);\r
\r
\r
  return;\r
}\r
\r
//________________________________________________________________________\r
void AliAnalysisTaskSEDStar::UserExec(Option_t */*option*/)\r
{\r
  // Execute analysis for current event:\r
  // heavy flavor candidates association to MC truth\r
  \r
  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());\r
 \r
  // load D*->D0 pi candidates                                                 \r
  TClonesArray *inputArray=0;\r
  \r
  if(!aod && 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
    aod = 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
      inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");\r
    }\r
  } else {\r
    inputArray=(TClonesArray*)aod->GetList()->FindObject("Dstar");\r
  }\r
\r
\r
  if(!inputArray) {\r
    printf("AliAnalysisTaskSEDStar::UserExec: input branch not found!\n");\r
    return;\r
  }\r
  \r
  // AOD primary vertex\r
  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();\r
  //vtx1->Print();\r
  \r
  TClonesArray *mcArray = 0;\r
  AliAODMCHeader *mcHeader = 0;\r
\r
  if(fReadMC) {\r
    // load MC particles\r
    mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
    if(!mcArray) {\r
      printf("AliAnalysisTaskSEDStar::UserExec: MC particles branch not found!\n");\r
      return;\r
    }\r
  \r
    // load MC header\r
    mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());\r
    if(!mcHeader) {\r
      printf("AliAnalysisTaskSEDStar::UserExec: MC header branch not found!\n");\r
      return;\r
    }\r
  }\r
  \r
  // loop over D*->D0 pi candidates\r
  Int_t nInDstar = inputArray->GetEntriesFast();\r
  if(fDebug > 1) printf("Number of D*->D0 pi: %d\n",nInDstar);\r
///////\r
  \r
  for (Int_t iDstar = 0; iDstar < nInDstar; iDstar++) {\r
    AliAODRecoCascadeHF *d = (AliAODRecoCascadeHF*)inputArray->UncheckedAt(iDstar);\r
\r
    Bool_t unsetvtx=kFALSE;\r
    if(!d->GetOwnPrimaryVtx()) {\r
      d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables\r
      d->Get2Prong()->SetOwnPrimaryVtx(vtx1);\r
      unsetvtx=kTRUE;\r
    }\r
    \r
    if(d->SelectDstar(fVHF->GetDstarCuts(),fVHF->GetD0fromDstarCuts(),kTRUE)) {//selected\r
      Double_t invmassDelta = d->DeltaInvMass();\r
      //TVector3 p3Trk0(d->PxProng(0),d->PyProng(0),d->PzProng(0)); // pi_s\r
      //TVector3 p3Trk1(d->PxProng(1),d->PyProng(1),d->PzProng(1)); // D0\r
      //Double_t CosOpenAngle = p3Trk0.Dot(p3Trk1)/(p3Trk0.Mag()*p3Trk1.Mag());\r
\r
      //PID of D0 daughters\r
      AliAODTrack *pos = (AliAODTrack*)d->Get2Prong()->GetDaughter(0);\r
      AliAODTrack *neg = (AliAODTrack*)d->Get2Prong()->GetDaughter(1);\r
\r
      if (fPID) {\r
        if(fDebug > 1) printf("AnalysisTaskSEDStar::TPCPIDon \n");\r
        if(fDebug > 1) printf("AnalysisTaskSEDStar::NSigmaTPC: %f\n", fNSigma);\r
\r
        if (d->Charge()>0){\r
          if(!SelectTPCPID(pos, 2, fNSigma)) continue;//pion+\r
          if(!SelectTPCPID(neg, 3, fNSigma)) continue;//kaon-\r
        }else{\r
          if(!SelectTPCPID(pos, 3, fNSigma)) continue;//kaon+\r
          if(!SelectTPCPID(neg, 2, fNSigma)) continue;//pion-\r
        }\r
      }\r
\r
      if(fReadMC) {\r
      Int_t labDstar = d->MatchToMC(413,421,mcArray); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)\r
      if(labDstar>=0) {\r
        //AliAODMCParticle *partDstar = (AliAODMCParticle*)mcArray->At(labDstar);\r
        //AliAODMCParticle *partPis = (AliAODMCParticle*)mcArray->At(partDstar->GetDaughter(1));\r
        //AliAODMCParticle *partD0 = (AliAODMCParticle*)mcArray->At(partDstar->GetDaughter(0));\r
        //AliAODMCParticle *partD0daughter = (AliAODMCParticle*)mcArray->At(partD0->GetDaughter(0));\r
        fhistSgn->Fill(invmassDelta);\r
\r
      }\r
      else {//background\r
        fhistBkg->Fill(invmassDelta);\r
      }\r
      }\r
      //no MC info, just cut selection\r
      fhistMass->Fill(invmassDelta);\r
\r
    } //else cout<<"NOT SELECTED"<<endl;\r
\r
    if(unsetvtx) {\r
      d->UnsetOwnPrimaryVtx();\r
      d->Get2Prong()->UnsetOwnPrimaryVtx();\r
    }\r
  }  \r
\r
  // Post the data\r
  PostData(1,fOutput);\r
\r
\r
  return;\r
}\r
//________________________________________________________________________\r
Bool_t AliAnalysisTaskSEDStar::SelectTPCPID(AliAODTrack *trk, Int_t pid, Double_t nsig){//pid(0-4): {e,mu,pi,K,p}\r
  Bool_t flag=kTRUE;\r
  if ((trk->GetStatus()&AliESDtrack::kTPCpid )==0) return flag;\r
  AliAODPid *detpid = trk->GetDetPid();\r
  static AliTPCPIDResponse TPCpid;\r
  Double_t nsigma = TPCpid.GetNumberOfSigmas(trk->P(),detpid->GetTPCsignal(),trk->GetTPCClusterMap().CountBits(),(AliPID::EParticleType)pid);\r
  if (TMath::Abs(nsigma)>nsig) flag=kFALSE;\r
  return flag;\r
}\r
\r
//________________________________________________________________________\r
void AliAnalysisTaskSEDStar::Terminate(Option_t */*option*/)\r
{\r
  // Terminate analysis\r
  //\r
  if(fDebug > 1) printf("AnalysisTaskSEDStar: Terminate() \n");\r
\r
  fOutput = dynamic_cast<TList*> (GetOutputData(1));\r
  if (!fOutput) {     \r
    printf("ERROR: fOutput not available\n");\r
    return;\r
  }\r
\r
  \r
  fhistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fhistMass"));\r
  fhistSgn = dynamic_cast<TH1F*>(fOutput->FindObject("fhistSgn"));\r
  fhistBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fhistBkg"));\r
\r
  return;\r
}\r
\r