PWGJE/FlavourJetTasks/AliAnalysisTaskSEDmesonsFilterCJ.cxx

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