PWG3/vertexingHF/AliAnalysisTaskSEDStar.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 \r
  16 /////////////////////////////////////////////////////////////\r
  17 //\r
  18 // AliAnalysisTaskSE for D* candidates invariant mass histogram\r
  19 // and comparison with the MC truth.\r
  20 //\r
  21 // Authors: Y.Wang, yifei@physi.uni-heidelberg.de\r
  22 /////////////////////////////////////////////////////////////\r
  23 \r
  24 #include <Riostream.h>\r
  25 #include <TClonesArray.h>\r
  26 #include <TList.h>\r
  27 #include <TH1F.h>\r
  28 \r
  29 \r
  30 #include "AliPID.h"\r
  31 #include "AliTPCPIDResponse.h"\r
  32 #include "AliAnalysisManager.h"\r
  33 #include "AliAODHandler.h"\r
  34 #include "AliESDtrack.h"\r
  35 #include "AliAODEvent.h"\r
  36 #include "AliAODVertex.h"\r
  37 #include "AliAODTrack.h"\r
  38 #include "AliAODMCHeader.h"\r
  39 #include "AliAODMCParticle.h"\r
  40 #include "AliAODRecoDecayHF2Prong.h"\r
  41 #include "AliAODRecoCascadeHF.h"\r
  42 #include "AliAnalysisVertexingHF.h"\r
  43 #include "AliAnalysisTaskSE.h"\r
  44 #include "AliAnalysisTaskSEDStar.h"\r
  45 \r
  46 ClassImp(AliAnalysisTaskSEDStar)\r
  47 \r
  48 \r
  49 //________________________________________________________________________\r
  50 AliAnalysisTaskSEDStar::AliAnalysisTaskSEDStar():\r
  51 AliAnalysisTaskSE(),\r
  52 fOutput(0), \r
  53 fhistMass(0),\r
  54 fhistSgn(0),\r
  55 fhistBkg(0),\r
  56 fReadMC(0),\r
  57 fPID(1),\r
  58 fNSigma(3),\r
  59 fVHF(0)\r
  60 {\r
  61   // Default constructor\r
  62 }\r
  63 \r
  64 //________________________________________________________________________\r
  65 AliAnalysisTaskSEDStar::AliAnalysisTaskSEDStar(const char *name):\r
  66 AliAnalysisTaskSE(name),\r
  67 fOutput(0), \r
  68 fhistMass(0),\r
  69 fhistSgn(0),\r
  70 fhistBkg(0),\r
  71 fReadMC(0),\r
  72 fPID(1),\r
  73 fNSigma(3),\r
  74 fVHF(0)\r
  75 {\r
  76   // Default constructor\r
  77 \r
  78   // Output slot #1 writes into a TList container\r
  79   DefineOutput(1,TList::Class());  //My private output\r
  80 }\r
  81 \r
  82 //________________________________________________________________________\r
  83 AliAnalysisTaskSEDStar::~AliAnalysisTaskSEDStar()\r
  84 {\r
  85   // Destructor\r
  86 \r
  87   if (fhistMass){\r
  88     delete fhistMass;\r
  89     fhistMass=0;\r
  90   }\r
  91 \r
  92   if (fhistSgn){\r
  93     delete fhistSgn;\r
  94     fhistSgn=0;\r
  95   }\r
  96 \r
  97   if (fhistBkg){\r
  98     delete fhistBkg;\r
  99     fhistBkg=0;\r
 100   }\r
 101 \r
 102   if (fOutput) {\r
 103     delete fOutput;\r
 104     fOutput = 0;\r
 105   }\r
 106   if (fVHF) {\r
 107     delete fVHF;\r
 108     fVHF = 0;\r
 109   }\r
 110 \r
 111 }  \r
 112 \r
 113 //________________________________________________________________________\r
 114 void AliAnalysisTaskSEDStar::Init()\r
 115 {\r
 116   // Initialization\r
 117 \r
 118   if(fDebug > 1) printf("AnalysisTaskSEDStar::Init() \n");\r
 119 \r
 120 //  gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/ConfigVertexingHF.C");\r
 121   gROOT->LoadMacro("ConfigVertexingHF.C");\r
 122 \r
 123   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");  \r
 124   // set dedidcated cuts\r
 125   //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
 126   fVHF->SetD0fromDstarCuts(0.3,999999.,1.1,0.,0.,999999.,999999.,999999.,0.);\r
 127   fVHF->SetDstarCuts(0.3, 0.1, 0.05, 100000000000.0, 0.5);\r
 128   fVHF->PrintStatus();\r
 129 \r
 130   return;\r
 131 }\r
 132 \r
 133 //________________________________________________________________________\r
 134 void AliAnalysisTaskSEDStar::UserCreateOutputObjects()\r
 135 {\r
 136 \r
 137   // Create the output container\r
 138   //\r
 139   if(fDebug > 1) printf("AnalysisTaskSEDStar::UserCreateOutputObjects() \n");\r
 140 \r
 141   // Several histograms are more conveniently managed in a TList\r
 142   fOutput = new TList();\r
 143   fOutput->SetOwner();\r
 144 \r
 145     fhistMass = new TH1F("fhistMass","D^{*+}-D^{0} invariant mass; M [GeV]; Entries",200,0,0.3);\r
 146     fhistSgn = new TH1F("fhistSgn", "D^{*+}-D^{0} invariant mass - MC; M [GeV]; Entries",200,0,0.3);\r
 147     fhistBkg = new TH1F("fhistBkg", "D^{*+}-Background invariant mass - MC; M [GeV]; Entries",200,0,0.3);\r
 148  \r
 149   fOutput->Add(fhistMass);\r
 150   fOutput->Add(fhistSgn);\r
 151   fOutput->Add(fhistBkg);\r
 152 \r
 153 \r
 154   return;\r
 155 }\r
 156 \r
 157 //________________________________________________________________________\r
 158 void AliAnalysisTaskSEDStar::UserExec(Option_t */*option*/)\r
 159 {\r
 160   // Execute analysis for current event:\r
 161   // heavy flavor candidates association to MC truth\r
 162   \r
 163   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());\r
 164  \r
 165   // load D*->D0 pi candidates                                                 \r
 166   TClonesArray *inputArray=0;\r
 167   \r
 168   if(!aod && AODEvent() && IsStandardAOD()) {\r
 169     // In case there is an AOD handler writing a standard AOD, use the AOD \r
 170     // event in memory rather than the input (ESD) event.    \r
 171     aod = dynamic_cast<AliAODEvent*> (AODEvent());\r
 172     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)\r
 173     // have to taken from the AOD event hold by the AliAODExtension\r
 174     AliAODHandler* aodHandler = (AliAODHandler*) \r
 175       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());\r
 176     if(aodHandler->GetExtensions()) {\r
 177       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");\r
 178       AliAODEvent *aodFromExt = ext->GetAOD();\r
 179       inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");\r
 180     }\r
 181   } else {\r
 182     inputArray=(TClonesArray*)aod->GetList()->FindObject("Dstar");\r
 183   }\r
 184 \r
 185 \r
 186   if(!inputArray) {\r
 187     printf("AliAnalysisTaskSEDStar::UserExec: input branch not found!\n");\r
 188     return;\r
 189   }\r
 190   \r
 191   // AOD primary vertex\r
 192   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();\r
 193   //vtx1->Print();\r
 194   \r
 195   TClonesArray *mcArray = 0;\r
 196   AliAODMCHeader *mcHeader = 0;\r
 197 \r
 198   if(fReadMC) {\r
 199     // load MC particles\r
 200     mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
 201     if(!mcArray) {\r
 202       printf("AliAnalysisTaskSEDStar::UserExec: MC particles branch not found!\n");\r
 203       return;\r
 204     }\r
 205   \r
 206     // load MC header\r
 207     mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());\r
 208     if(!mcHeader) {\r
 209       printf("AliAnalysisTaskSEDStar::UserExec: MC header branch not found!\n");\r
 210       return;\r
 211     }\r
 212   }\r
 213 \r
 214   Int_t pdgDg[2]={421,211},pdgD0Dg[2]={321,211};\r
 215   \r
 216   // loop over D*->D0 pi candidates\r
 217   Int_t nInDstar = inputArray->GetEntriesFast();\r
 218   if(fDebug > 1) printf("Number of D*->D0 pi: %d\n",nInDstar);\r
 219 ///////\r
 220   \r
 221   for (Int_t iDstar = 0; iDstar < nInDstar; iDstar++) {\r
 222     AliAODRecoCascadeHF *d = (AliAODRecoCascadeHF*)inputArray->UncheckedAt(iDstar);\r
 223 \r
 224     Bool_t unsetvtx=kFALSE;\r
 225     if(!d->GetOwnPrimaryVtx()) {\r
 226       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables\r
 227       d->Get2Prong()->SetOwnPrimaryVtx(vtx1);\r
 228       unsetvtx=kTRUE;\r
 229     }\r
 230     \r
 231     if(d->SelectDstar(fVHF->GetDstarCuts(),fVHF->GetD0fromDstarCuts(),kTRUE)) {//selected\r
 232       Double_t invmassDelta = d->DeltaInvMass();\r
 233       //TVector3 p3Trk0(d->PxProng(0),d->PyProng(0),d->PzProng(0)); // pi_s\r
 234       //TVector3 p3Trk1(d->PxProng(1),d->PyProng(1),d->PzProng(1)); // D0\r
 235       //Double_t CosOpenAngle = p3Trk0.Dot(p3Trk1)/(p3Trk0.Mag()*p3Trk1.Mag());\r
 236 \r
 237       //PID of D0 daughters\r
 238       AliAODTrack *pos = (AliAODTrack*)d->Get2Prong()->GetDaughter(0);\r
 239       AliAODTrack *neg = (AliAODTrack*)d->Get2Prong()->GetDaughter(1);\r
 240 \r
 241       if (fPID) {\r
 242         if(fDebug > 1) printf("AnalysisTaskSEDStar::TPCPIDon \n");\r
 243         if(fDebug > 1) printf("AnalysisTaskSEDStar::NSigmaTPC: %f\n", fNSigma);\r
 244 \r
 245         if (d->Charge()>0){\r
 246           if(!SelectTPCPID(pos, 2, fNSigma)) continue;//pion+\r
 247           if(!SelectTPCPID(neg, 3, fNSigma)) continue;//kaon-\r
 248         }else{\r
 249           if(!SelectTPCPID(pos, 3, fNSigma)) continue;//kaon+\r
 250           if(!SelectTPCPID(neg, 2, fNSigma)) continue;//pion-\r
 251         }\r
 252       }\r
 253 \r
 254       if(fReadMC) {\r
 255       Int_t labDstar = d->MatchToMC(413,421,pdgDg,pdgD0Dg,mcArray); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)\r
 256       if(labDstar>=0) {\r
 257         //AliAODMCParticle *partDstar = (AliAODMCParticle*)mcArray->At(labDstar);\r
 258         //AliAODMCParticle *partPis = (AliAODMCParticle*)mcArray->At(partDstar->GetDaughter(1));\r
 259         //AliAODMCParticle *partD0 = (AliAODMCParticle*)mcArray->At(partDstar->GetDaughter(0));\r
 260         //AliAODMCParticle *partD0daughter = (AliAODMCParticle*)mcArray->At(partD0->GetDaughter(0));\r
 261         fhistSgn->Fill(invmassDelta);\r
 262 \r
 263       }\r
 264       else {//background\r
 265         fhistBkg->Fill(invmassDelta);\r
 266       }\r
 267       }\r
 268       //no MC info, just cut selection\r
 269       fhistMass->Fill(invmassDelta);\r
 270 \r
 271     } //else cout<<"NOT SELECTED"<<endl;\r
 272 \r
 273     if(unsetvtx) {\r
 274       d->UnsetOwnPrimaryVtx();\r
 275       d->Get2Prong()->UnsetOwnPrimaryVtx();\r
 276     }\r
 277   }  \r
 278 \r
 279   // Post the data\r
 280   PostData(1,fOutput);\r
 281 \r
 282 \r
 283   return;\r
 284 }\r
 285 //________________________________________________________________________\r
 286 Bool_t AliAnalysisTaskSEDStar::SelectTPCPID(AliAODTrack *trk, Int_t pid, Double_t nsig){//pid(0-4): {e,mu,pi,K,p}\r
 287   Bool_t flag=kTRUE;\r
 288   if ((trk->GetStatus()&AliESDtrack::kTPCpid )==0) return flag;\r
 289   AliAODPid *detpid = trk->GetDetPid();\r
 290   static AliTPCPIDResponse TPCpid;\r
 291   Double_t nsigma = TPCpid.GetNumberOfSigmas(trk->P(),detpid->GetTPCsignal(),trk->GetTPCClusterMap().CountBits(),(AliPID::EParticleType)pid);\r
 292   if (TMath::Abs(nsigma)>nsig) flag=kFALSE;\r
 293   return flag;\r
 294 }\r
 295 \r
 296 //________________________________________________________________________\r
 297 void AliAnalysisTaskSEDStar::Terminate(Option_t */*option*/)\r
 298 {\r
 299   // Terminate analysis\r
 300   //\r
 301   if(fDebug > 1) printf("AnalysisTaskSEDStar: Terminate() \n");\r
 302 \r
 303   fOutput = dynamic_cast<TList*> (GetOutputData(1));\r
 304   if (!fOutput) {     \r
 305     printf("ERROR: fOutput not available\n");\r
 306     return;\r
 307   }\r
 308 \r
 309   \r
 310   fhistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fhistMass"));\r
 311   fhistSgn = dynamic_cast<TH1F*>(fOutput->FindObject("fhistSgn"));\r
 312   fhistBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fhistBkg"));\r
 313 \r
 314   return;\r
 315 }\r
 316 \r