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