PWGJE/FlavourJetTasks/AliAnalysisTaskSEDmesonsFilterCJ.cxx

   1 // $Id$
   2 /**************************************************************************
   3  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16 //
  17 //
  18 // Task for selecting D mesons to be used as an input for D-Jet correlations
  19 //
  20 //-----------------------------------------------------------------------
  21 // Authors:
  22 // C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
  23 // A.Grelli (Utrecht University) a.grelli@uu.nl
  24 // Xiaoming Zhang (LBNL)  XMZhang@lbl.gov
  25 //-----------------------------------------------------------------------
  26
  27 #include <TDatabasePDG.h>
  28 #include <TParticle.h>
  29 #include <TVector3.h>
  30 #include "TROOT.h"
  31 #include <TH3F.h>
  32
  33 #include "AliRDHFCutsDStartoKpipi.h"
  34 #include "AliRDHFCutsD0toKpi.h"
  35 #include "AliAODMCHeader.h"
  36 #include "AliAODHandler.h"
  37 #include "AliAnalysisManager.h"
  38 #include "AliLog.h"
  39 #include "AliAODVertex.h"
  40 #include "AliAODRecoDecay.h"
  41 #include "AliAODRecoCascadeHF.h"
  42 #include "AliAODRecoDecayHF2Prong.h"
  43 #include "AliESDtrack.h"
  44 #include "AliAODMCParticle.h"
  45 #include "AliAnalysisTaskSEDmesonsFilterCJ.h"
  46
  47 ClassImp(AliAnalysisTaskSEDmesonsFilterCJ)
  48
  49 //_____________________________________________________________________________
  50 AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ() :
  51   AliAnalysisTaskSE(),
  52   fUseMCInfo(kFALSE),
  53   fCandidateType(0),
  54   fCandidateName(""),
  55   fPDGmother(0),
  56   fNProngs(0),
  57   fBranchName(""),
  58   fOutput(0),
  59 //fOutputCandidates(0),
  60   fCuts(0),
  61   fMinMass(0.),
  62   fMaxMass(0.),
  63   fCandidateArray(0)
  64 //fIsSelectedArray(0)
  65 {
  66 //
  67 // Default constructor
  68 //
  69
  70   for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
  71   for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
  72 }
  73
  74 //_____________________________________________________________________________
  75 AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ(const char *name, AliRDHFCuts *cuts, ECandidateType candtype) :
  76   AliAnalysisTaskSE(name),
  77   fUseMCInfo(kFALSE),
  78   fCandidateType(candtype),
  79   fCandidateName(""),
  80   fPDGmother(0),
  81   fNProngs(0),
  82 //fPDGdaughters(),
  83   fBranchName(""),
  84   fOutput(0),
  85 //fOutputCandidates(0),
  86   fCuts(cuts),
  87   fMinMass(0.),
  88   fMaxMass(0.),
  89   fCandidateArray(0)
  90 //fIsSelectedArray(0)
  91 {
  92 //
  93 // Constructor. Initialization of Inputs and Outputs
  94 //
  95
  96   Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");
  97
  98   for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
  99   for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
 100
 101   const Int_t nptbins = fCuts->GetNPtBins();
 102   Float_t defaultSigmaD013[13] = { 0.012, 0.012, 0.012, 0.015, 0.015, 0.018, 0.018, 0.020, 0.020, 0.030, 0.030, 0.037, 0.040 };
 103
 104   switch (fCandidateType) {
 105   case 0 :
 106     fCandidateName = "D0";
 107     fPDGmother = 421;
 108     fNProngs = 2;
 109     fPDGdaughters[0] = 211;  // pi
 110     fPDGdaughters[1] = 321;  // K
 111     fPDGdaughters[2] = 0;    // empty
 112     fPDGdaughters[3] = 0;    // empty
 113     fBranchName = "D0toKpi";
 114     break;
 115   case 1 :
 116     fCandidateName = "Dstar";
 117     fPDGmother = 413;
 118     fNProngs = 3;
 119     fPDGdaughters[1] = 211; // pi soft
 120     fPDGdaughters[0] = 421; // D0
 121     fPDGdaughters[2] = 211; // pi fromD0
 122     fPDGdaughters[3] = 321; // K from D0
 123     fBranchName = "Dstar";
 124
 125     if (nptbins<=13) {
 126       for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];
 127     } else {
 128       AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
 129     }
 130     break;
 131   default :
 132     printf("%d not accepted!!\n",fCandidateType);
 133     break;
 134   }
 135
 136   if (fCandidateType==kD0toKpi) SetMassLimits(0.15, fPDGmother);
 137   if (fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
 138
 139   DefineOutput(1, TList::Class());       // histos
 140   DefineOutput(2, AliRDHFCuts::Class()); // my cuts
 141 }
 142
 143 //_____________________________________________________________________________
 144 AliAnalysisTaskSEDmesonsFilterCJ::~AliAnalysisTaskSEDmesonsFilterCJ()
 145 {
 146 //
 147 // destructor
 148 //
 149
 150   Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");
 151
 152   if (fOutput) { delete fOutput; fOutput = 0; }
 153   if (fCuts)   { delete fCuts;   fCuts   = 0; }
 154   if (fCandidateArray)  { delete fCandidateArray;  fCandidateArray  = 0; }
 155 //if (fIsSelectedArray) { delete fIsSelectedArray; fIsSelectedArray = 0; }
 156 }
 157
 158 //_____________________________________________________________________________
 159 void AliAnalysisTaskSEDmesonsFilterCJ::Init()
 160 {
 161 //
 162 // Initialization
 163 //
 164
 165   if(fDebug>1) printf("AnalysisTaskSEDmesonsForJetCorrelations::Init() \n");
 166
 167   switch (fCandidateType) {
 168     case 0: {
 169               AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
 170               copyfCutsDzero->SetName("AnalysisCutsDzero");
 171               PostData(2, copyfCutsDzero);  // Post the data
 172             } break;
 173     case 1: {
 174               AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
 175               copyfCutsDstar->SetName("AnalysisCutsDStar");
 176               PostData(2, copyfCutsDstar); // Post the cuts
 177             } break;
 178     default: return;
 179   }
 180
 181   return;
 182 }
 183
 184 //_____________________________________________________________________________
 185 void AliAnalysisTaskSEDmesonsFilterCJ::UserCreateOutputObjects()
 186 {
 187 //
 188 // output
 189 //
 190
 191   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 192
 193 //OpenFile(1);
 194   fOutput = new TList(); fOutput->SetOwner();
 195   DefineHistoForAnalysis(); // define histograms
 196
 197 //fOutputCandidates= new TList();
 198 //fOutputCandidates->SetOwner();
 199
 200   fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);
 201   fCandidateArray->SetName(Form("fCandidateArray%s",fCandidateName.Data()));
 202
 203 //fOutputCandidates->Add(fCandidateArray);
 204 //fIsSelectedArray = new TClonesArray("Int_t",0);
 205 //fIsSelectedArray->SetName(Form("fIsSelectedArray%s",suffix.Data()));
 206 //fOutputCandidates->Add(fIsSelectedArray);
 207
 208   PostData(1, fOutput);
 209 //PostData(3, fOutputCandidates);
 210   return;
 211 }
 212
 213 //_____________________________________________________________________________
 214 void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *)
 215 {
 216 //
 217 // user exec
 218 //
 219
 220   // add cadidate branch
 221   fCandidateArray->Delete();
 222   if (!(InputEvent()->FindListObject(Form("fCandidateArray%s",fCandidateName.Data())))) InputEvent()->AddObject(fCandidateArray);
 223
 224   // Load the event
 225   AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 226
 227   TClonesArray *arrayDStartoD0pi = 0;
 228   if (!aodEvent && AODEvent() && IsStandardAOD()) {
 229
 230     // In case there is an AOD handler writing a standard AOD, use the AOD
 231     // event in memory rather than the input (ESD) event.
 232     aodEvent = dynamic_cast<AliAODEvent*>(AODEvent());
 233
 234     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 235     // have to taken from the AOD event hold by the AliAODExtension
 236     AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 237     if(aodHandler->GetExtensions()) {
 238       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 239       AliAODEvent *aodFromExt = ext->GetAOD();
 240       arrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
 241     }
 242   } else {
 243     arrayDStartoD0pi = (TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
 244   }
 245
 246   if (!arrayDStartoD0pi) {
 247     AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
 248     return;
 249   } else {
 250     AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
 251   }
 252
 253   TClonesArray* mcArray = 0x0;
 254   if (fUseMCInfo) {
 255     mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 256     if (!mcArray) {
 257       printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
 258       return;
 259     }
 260   }
 261
 262   //Histograms
 263   TH1I* hstat = (TH1I*)fOutput->FindObject("hstat");
 264   TH2F* hInvMassptD = (TH2F*)fOutput->FindObject("hInvMassptD");
 265
 266   TH1F* hPtPion=0x0;
 267   if (fCandidateType==kDstartoKpipi) hPtPion = (TH1F*)fOutput->FindObject("hPtPion");
 268   hstat->Fill(0);
 269
 270   // fix for temporary bug in ESDfilter
 271   // the AODs with null vertex pointer didn't pass the PhysSel
 272   if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
 273
 274   //Event selection
 275   Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
 276 //TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
 277   if (!iseventselected) return;
 278   hstat->Fill(1);
 279
 280   const Int_t nD = arrayDStartoD0pi->GetEntriesFast();
 281   hstat->Fill(2, nD);
 282
 283 //Int_t icountReco = 0;  // counters for efficiencies
 284
 285   //D* and D0 prongs needed to MatchToMC method
 286   Int_t pdgDgDStartoD0pi[2] = { 421, 211 };  // D0,pi
 287   Int_t pdgDgD0toKpi[2] = { 321, 211 };      // K, pi
 288
 289   //D0 from D0 bar
 290   Int_t pdgdaughtersD0[2] = { 211, 321 };     // pi,K
 291   Int_t pdgdaughtersD0bar[2] = { 321, 211 };  // K,pi
 292
 293   Int_t iCand =0;
 294   Int_t isSelected = 0;
 295   AliAODRecoDecayHF *charmCand = 0;
 296
 297   Int_t mcLabel = -9999;
 298   Int_t pdgMeson = 413;
 299   if (fCandidateType==kD0toKpi) pdgMeson = 421;
 300
 301   for (Int_t icharm=0; icharm<nD; icharm++) {   //loop over D candidates
 302     charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates
 303     if (!charmCand) return;
 304
 305
 306     if (fUseMCInfo) { // Look in MC, try to simulate the z
 307       if (fCandidateType==kDstartoKpipi) {
 308         AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
 309         mcLabel = temp->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
 310       }
 311
 312       if (fCandidateType==kD0toKpi) mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);
 313
 314       if (mcLabel<=0) continue;
 315     }
 316
 317     // region of interest + cuts
 318     if (!fCuts->IsInFiducialAcceptance(charmCand->Pt(),charmCand->Y(pdgMeson))) continue;
 319
 320     isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected
 321     if (!isSelected) continue;
 322
 323     new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand);
 324 //  new ((*fIsSelectedArray)[iCand]) Int_t(isSelected);
 325     iCand++;
 326
 327     Double_t masses[2];
 328     Double_t ptD = charmCand->Pt();
 329     if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi
 330
 331       //softpion from D* decay
 332       AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
 333       AliAODTrack *track2 = (AliAODTrack*)temp->GetBachelor();
 334       hstat->Fill(3);
 335
 336       // select D* in the D0 window.
 337       // In the cut object window is loose to allow side bands
 338       Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 339
 340       // retrieve the corresponding pt bin for the candidate
 341       // and set the expected D0 width (x3)
 342 //    static const Int_t n = fCuts->GetNPtBins();
 343       Int_t bin = fCuts->PtBin(ptD);
 344       if(bin<0 || bin>=fCuts->GetNPtBins()) {
 345         AliError(Form("Pt %.3f out of bounds",ptD));
 346         continue;
 347       }
 348
 349       AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin]));
 350       if ((temp->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) {
 351         masses[0] = temp->DeltaInvMass(); //D*
 352         masses[1] = 0.; //dummy for D*
 353
 354         //D*  delta mass
 355         hInvMassptD->Fill(masses[0], ptD); // 2 D slice for pt bins
 356
 357         // D* pt and soft pion pt for good candidates
 358         Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
 359         Double_t invmassDelta = temp->DeltaInvMass();
 360         if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt());
 361       }
 362     } //Dstar specific
 363
 364     if (fCandidateType==kD0toKpi) { //D0->Kpi
 365
 366       //needed quantities
 367       masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
 368       masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
 369       hstat->Fill(3);
 370
 371       // mass vs pt
 372       if (isSelected==1 || isSelected==3) hInvMassptD->Fill(masses[0],ptD);
 373       if (isSelected>=2) hInvMassptD->Fill(masses[1],ptD);
 374     } //D0 specific
 375
 376     charmCand = 0;
 377   } // end of D cand loop
 378
 379 //AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
 380
 381 //PostData(1,fOutput);
 382 //PostData(3,fOutputCandidates);
 383
 384   return;
 385 }
 386
 387 //_____________________________________________________________________________
 388 void AliAnalysisTaskSEDmesonsFilterCJ::Terminate(Option_t*)
 389 {
 390 // The Terminate() function is the last function to be called during
 391 // a query. It always runs on the client, it can be used to present
 392 // the results graphically or save the results to file.
 393
 394   Info("Terminate"," terminate");
 395   AliAnalysisTaskSE::Terminate();
 396
 397   fOutput = dynamic_cast<TList*>(GetOutputData(1));
 398   if (!fOutput) {
 399     printf("ERROR: fOutput not available\n");
 400     return;
 401   }
 402
 403   return;
 404 }
 405
 406 //_____________________________________________________________________________
 407 void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t range, Int_t pdg)
 408 {
 409 //
 410 // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
 411 //
 412
 413   Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
 414
 415   // compute the Delta mass for the D*
 416   if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
 417 //cout << "mass ---------------" << endl;
 418
 419   fMinMass = mass - range;
 420   fMaxMass = mass + range;
 421
 422   AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
 423   if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) AliFatal("Wrong mass limits!\n");
 424
 425   return;
 426 }
 427
 428 //_____________________________________________________________________________
 429 void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t lowlimit, Double_t uplimit)
 430 {
 431 //
 432 // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
 433 //
 434
 435   if (uplimit>lowlimit) {
 436     fMinMass = lowlimit;
 437     fMaxMass = uplimit;
 438   } else {
 439     printf("Error! Lower limit larger than upper limit!\n");
 440     fMinMass = uplimit - uplimit*0.2;
 441     fMaxMass = uplimit;
 442   }
 443
 444   return;
 445 }
 446
 447 //_____________________________________________________________________________
 448 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar(Int_t nptbins, Float_t *width)
 449 {
 450 //
 451 // AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar
 452 //
 453
 454   if (nptbins>30) {
 455     AliInfo("Maximum number of bins allowed is 30!");
 456     return kFALSE;
 457   }
 458
 459   if (!width) return kFALSE;
 460   for (Int_t ipt=0; ipt<nptbins; ipt++) fSigmaD0[ipt]=width[ipt];
 461
 462   return kTRUE;
 463 }
 464
 465
 466 //_____________________________________________________________________________
 467 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()
 468 {
 469 //
 470 // AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis
 471 //
 472
 473   // Statistics
 474   TH1I* hstat = new TH1I("hstat","Statistics",6,-0.5,5.5);
 475   hstat->GetXaxis()->SetBinLabel(1, "N ev anal");
 476   hstat->GetXaxis()->SetBinLabel(2, "N ev sel");
 477   hstat->GetXaxis()->SetBinLabel(3, "N cand");
 478   hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
 479   hstat->GetXaxis()->SetBinLabel(5, "N jets");
 480   hstat->GetXaxis()->SetBinLabel(6, "N cand in jet");
 481 /*if(fUseMCInfo) {
 482     hstat->GetXaxis()->SetBinLabel(7,"N D");
 483     hstat->GetXaxis()->SetBinLabel(8,"N D in jet");
 484   }*/
 485   hstat->SetNdivisions(1);
 486   fOutput->Add(hstat);
 487
 488   // Invariant mass related histograms
 489   const Int_t nbinsmass = 200;
 490   TH2F *hInvMass = new TH2F("hInvMassptD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, 100, 0., 50.);
 491   hInvMass->SetStats(kTRUE);
 492   hInvMass->GetXaxis()->SetTitle("mass (GeV/c)");
 493   hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 494   fOutput->Add(hInvMass);
 495
 496   if (fCandidateType==kDstartoKpipi) {
 497     TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.);
 498     hPtPion->SetStats(kTRUE);
 499     hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 500     hPtPion->GetYaxis()->SetTitle("entries");
 501     fOutput->Add(hPtPion);
 502   }
 503
 504   return kTRUE;
 505 }