PWG3/vertexingHF/AliAnalysisTaskSEHFQA.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 //
  18 // AliAnalysisTaskSE for HF quality assurance
  19 //
  20 // Author: Chiara Bianchin, chiara.bianchin@pd.infn.it
  21 /////////////////////////////////////////////////////////////
  22
  23 #include <Riostream.h>
  24 #include <TClonesArray.h>
  25 #include <TCanvas.h>
  26 #include <TList.h>
  27 #include <TH1F.h>
  28 #include <TH2F.h>
  29 #include <TDatabasePDG.h>
  30
  31 #include <AliAnalysisDataSlot.h>
  32 #include <AliAnalysisDataContainer.h>
  33 #include "AliAnalysisManager.h"
  34 #include "AliESDtrack.h"
  35 #include "AliVertexerTracks.h"
  36 #include "AliPID.h"
  37 #include "AliTPCPIDResponse.h"
  38 #include "AliAODHandler.h"
  39 #include "AliAODEvent.h"
  40 #include "AliAODVertex.h"
  41 #include "AliAODTrack.h"
  42 #include "AliAODMCParticle.h"
  43 #include "AliAODRecoDecayHF2Prong.h"
  44 #include "AliAODRecoCascadeHF.h"
  45 #include "AliAnalysisVertexingHF.h"
  46 #include "AliAnalysisTaskSE.h"
  47 #include "AliRDHFCuts.h"
  48 #include "AliRDHFCutsDplustoKpipi.h"
  49 #include "AliRDHFCutsD0toKpipipi.h"
  50 #include "AliRDHFCutsDstoKKpi.h"
  51 #include "AliRDHFCutsDStartoKpipi.h"
  52 #include "AliRDHFCutsD0toKpi.h"
  53 #include "AliRDHFCutsLctopKpi.h"
  54
  55 #include "AliAnalysisTaskSEHFQA.h"
  56
  57
  58 ClassImp(AliAnalysisTaskSEHFQA)
  59
  60 //____________________________________________________________________________
  61
  62 AliAnalysisTaskSEHFQA::AliAnalysisTaskSEHFQA():AliAnalysisTaskSE(),
  63   fNEntries(0x0),
  64   fOutputPID(0x0),
  65   fOutputTrack(0x0),
  66   fDecayChannel(AliAnalysisTaskSEHFQA::kD0toKpi),
  67   fCuts(0x0)
  68 {
  69   //default constructor
  70 }
  71
  72 //____________________________________________________________________________
  73 AliAnalysisTaskSEHFQA::AliAnalysisTaskSEHFQA(const char *name, AliAnalysisTaskSEHFQA::DecChannel ch,AliRDHFCuts* cuts):
  74 AliAnalysisTaskSE(name),
  75 fNEntries(0x0),
  76 fOutputPID(0x0),
  77 fOutputTrack(0x0),
  78 fDecayChannel(ch),
  79 fCuts(0x0)
  80 {
  81   //constructor
  82
  83   //SetCutObject(cuts);
  84   fCuts=cuts;
  85
  86   // Output slot #1 writes into a TH1F container (number of events)
  87   DefineOutput(1,TH1F::Class());  //My private output
  88   // Output slot #2 writes into a TList container (PID)
  89   DefineOutput(2,TList::Class());  //My private output
  90   // Output slot #3 writes into a TList container (Tracks)
  91   DefineOutput(3,TList::Class());  //My private output
  92   // Output slot #4 writes into a AliRDHFCuts container (cuts)
  93   switch(fDecayChannel){
  94   case 0:
  95     DefineOutput(4,AliRDHFCutsDplustoKpipi::Class());  //My private output
  96     break;
  97   case 1:
  98     DefineOutput(4,AliRDHFCutsD0toKpi::Class());  //My private output
  99     break;
 100   case 2:
 101     DefineOutput(4,AliRDHFCutsDStartoKpipi::Class());  //My private output
 102     break;
 103   case 3:
 104     DefineOutput(4,AliRDHFCutsDstoKKpi::Class());  //My private output
 105     break;
 106   case 4:
 107     DefineOutput(4,AliRDHFCutsD0toKpipipi::Class());  //My private output
 108     break;
 109    case 5:
 110     DefineOutput(4,AliRDHFCutsLctopKpi::Class());  //My private output
 111     break;
 112  }
 113 }
 114
 115 //___________________________________________________________________________
 116 AliAnalysisTaskSEHFQA::~AliAnalysisTaskSEHFQA()
 117 {
 118   //destructor
 119   if(fNEntries){
 120     delete fNEntries;
 121     fNEntries=0;
 122   }
 123   if(fOutputPID){
 124     delete fOutputPID;
 125     fOutputPID=0;
 126   }
 127   if(fOutputTrack){
 128     delete fOutputTrack;
 129     fOutputTrack=0;
 130   }
 131 }
 132
 133 //___________________________________________________________________________
 134 void AliAnalysisTaskSEHFQA::Init(){
 135
 136   //initialization
 137   if(fDebug > 1) printf("AnalysisTaskSEHFQA::Init() \n");
 138
 139   switch(fDecayChannel){
 140   case 0:
 141     {
 142       AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
 143       // Post the data
 144       PostData(4,copycut);
 145     }
 146     break;
 147   case 1:
 148     {
 149       AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
 150       // Post the data
 151       PostData(4,copycut);
 152     }
 153     break;
 154   case 2:
 155     {
 156       AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
 157       // Post the data
 158       PostData(4,copycut);
 159     }
 160     break;
 161   case 3:
 162     {
 163       AliRDHFCutsDstoKKpi* copycut=new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
 164       // Post the data
 165       PostData(4,copycut);
 166     }
 167     break;
 168   case 4:
 169     {
 170       AliRDHFCutsD0toKpipipi* copycut=new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
 171       // Post the data
 172       PostData(4,copycut);
 173     }
 174     break;
 175   case 5:
 176     {
 177       AliRDHFCutsLctopKpi* copycut=new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fCuts)));
 178       // Post the data
 179       PostData(4,copycut);
 180     }
 181     break;
 182
 183   default:
 184     return;
 185   }
 186
 187
 188
 189 }
 190
 191 //___________________________________________________________________________
 192 void AliAnalysisTaskSEHFQA::UserCreateOutputObjects()
 193 {
 194
 195   //create the output container
 196   if(fDebug > 1) printf("AnalysisTaskSEHFQA::UserCreateOutputObjects() \n");
 197
 198   //count events
 199
 200   fNEntries=new TH1F(GetOutputSlot(1)->GetContainer()->GetName(), "Counts the number of events", 6,-0.5,5.5);
 201   fNEntries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 202   fNEntries->GetXaxis()->SetBinLabel(2,"Pile-up Rej");
 203   fNEntries->GetXaxis()->SetBinLabel(3,"No VertexingHF");
 204   fNEntries->GetXaxis()->SetBinLabel(4,"nCandidates(AnCuts)");
 205   fNEntries->GetXaxis()->SetBinLabel(5,"EventsWithGoodVtx");
 206   fNEntries->GetXaxis()->SetBinLabel(6,"N. of 0SMH");
 207   fNEntries->GetXaxis()->SetNdivisions(1,kFALSE);
 208
 209   //PID
 210
 211   fOutputPID=new TList();
 212   fOutputPID->SetOwner();
 213   fOutputPID->SetName(GetOutputSlot(2)->GetContainer()->GetName());
 214
 215   //TOF pid
 216   TString hname="hTOFsig";
 217   TH1F* hTOFsig=new TH1F(hname.Data(),"Distribution of TOF signal;TOF time [ps];Entries", 100, -2.e3,40.e3);
 218
 219   hname="hTOFtime";
 220   TH1F* hTOFtime=new TH1F(hname.Data(),"Distribution of TOF time Kaon;TOF time(Kaon) [ps];Entries", 1000, 0.,50000.);
 221
 222   hname="hTOFtimeKaonHyptime";
 223   TH2F* hTOFtimeKaonHyptime=new TH2F(hname.Data(),"TOFtime - timeHypothesisForKaon;p[GeV/c];TOFtime - timeHypothesisForKaon [ps]",200,0.,4.,1000,-20000.,20000.);
 224
 225   hname="hTOFtimeKaonHyptimeAC";
 226   TH2F* hTOFtimeKaonHyptimeAC=new TH2F(hname.Data(),"TOFtime - timeHypothesisForKaon;p[GeV/c];TOFtime - timeHypothesisForKaon [ps]",200,0.,4.,1000,-20000.,20000.);
 227
 228   hname="hTOFsigmaK160";
 229   TH2F* hTOFsigmaK160=new TH2F(hname.Data(),"(TOFsignal-timeK)/160 ps;p[GeV/c];(TOFsignal-timeK)/160 ps",200,0.,4.,100,-5,5);
 230
 231   hname="hTOFsigmaPion160";
 232   TH2F* hTOFsigmaPion160=new TH2F(hname.Data(),"(TOFsignal-time#pi)/160 ps;p[GeV/c];(TOFsignal-time#pi)/160 ps",200,0.,4.,100,-5,5);
 233
 234   hname="hTOFsigmaProton160";
 235   TH2F* hTOFsigmaProton160=new TH2F(hname.Data(),"(TOFsignal-timep)/160 ps;p[GeV/c];(TOFsignal-time p)/160 ps",200,0.,4.,100,-5,5);
 236
 237   hname="hTOFsigmaK120";
 238   TH2F* hTOFsigmaK120=new TH2F(hname.Data(),"(TOFsignal-timeK)/120 ps;p[GeV/c];(TOFsignal-timeK)/120 ps",200,0.,4.,100,-5,5);
 239
 240   hname="hTOFsigmaPion120";
 241   TH2F* hTOFsigmaPion120=new TH2F(hname.Data(),"(TOFsignal-time#pi)/120 ps;p[GeV/c];(TOFsignal-time#pi)/120 ps",200,0.,4.,100,-5,5);
 242
 243   hname="hTOFsigmaProton120";
 244   TH2F* hTOFsigmaProton120=new TH2F(hname.Data(),"(TOFsignal-timep)/120 ps;p[GeV/c];(TOFsignal-time p)/120 ps",200,0.,4.,100,-5,5);
 245
 246   //TPC pid
 247   hname="hTPCsig";
 248   TH1F* hTPCsig=new TH1F(hname.Data(),"Distribution of TPC signal;TPC sig;Entries", 100, 35.,100.);
 249
 250   hname="hTPCsigvsp";
 251   TH2F* hTPCsigvsp=new TH2F(hname.Data(),"TPCsig vs p;TPC p[GeV/c];TPCsig",200,0.,4.,1000,35.,100.);
 252
 253   hname="hTPCsigvspAC";
 254   TH2F* hTPCsigvspAC=new TH2F(hname.Data(),"TPCsig vs p;TPCp[GeV/c];TPCsig",200,0.,4.,1000,35.,100.);
 255
 256   hname="hTPCsigmaK";
 257   TH2F* hTPCsigmaK=new TH2F(hname.Data(),"TPC Sigma for K as a function of momentum;p[GeV/c];Sigma Kaon",200,0.,4.,200,-5,5);
 258
 259   hname="hTPCsigmaPion";
 260   TH2F* hTPCsigmaPion=new TH2F(hname.Data(),"TPC Sigma for #pi as a function of momentum;p[GeV/c];Sigma #pi",200,0.,4.,200,-5,5);
 261
 262   hname="hTPCsigmaProton";
 263   TH2F* hTPCsigmaProton=new TH2F(hname.Data(),"TPC Sigma for proton as a function of momentum;p[GeV/c];Sigma Proton",200,0.,4.,200,-5,5);
 264
 265   fOutputPID->Add(hTOFsig);
 266   fOutputPID->Add(hTPCsig);
 267   fOutputPID->Add(hTOFtime);
 268   fOutputPID->Add(hTOFtimeKaonHyptime);
 269   fOutputPID->Add(hTOFtimeKaonHyptimeAC);
 270   fOutputPID->Add(hTOFsigmaK160);
 271   fOutputPID->Add(hTOFsigmaPion160);
 272   fOutputPID->Add(hTOFsigmaProton160);
 273   fOutputPID->Add(hTOFsigmaK120);
 274   fOutputPID->Add(hTOFsigmaPion120);
 275   fOutputPID->Add(hTOFsigmaProton120);
 276   fOutputPID->Add(hTPCsigvsp);
 277   fOutputPID->Add(hTPCsigvspAC);
 278   fOutputPID->Add(hTPCsigmaK);
 279   fOutputPID->Add(hTPCsigmaPion);
 280   fOutputPID->Add(hTPCsigmaProton);
 281
 282   //quality of the tracks
 283
 284   fOutputTrack=new TList();
 285   fOutputTrack->SetOwner();
 286   fOutputTrack->SetName(GetOutputSlot(3)->GetContainer()->GetName());
 287
 288   hname="hnClsITS";
 289   TH1F* hnClsITS=new TH1F(hname.Data(),"Distribution of number of ITS clusters;nITScls;Entries",7,-0.5,6.5);
 290
 291   hname="hnClsITS-SA";
 292   TH1F* hnClsITSSA=new TH1F(hname.Data(),"Distribution of number of ITS clusters(ITS-SA);nITScls;Entries",7,-0.5,6.5);
 293
 294   hname="hnClsSPD";
 295   TH1F* hnClsSPD=new TH1F(hname.Data(),"Distribution of number of SPD clusters;nSPDcls;Entries",3,-0.5,2.5);
 296
 297   hname="hptGoodTr";
 298   TH1F* hptGoodTr=new TH1F(hname.Data(),"Pt distribution of 'good' tracks;p_{t}[GeV];Entries/0.05 GeV/c",400,0.,20.);
 299   hptGoodTr->SetTitleOffset(1.3,"Y");
 300
 301   hname="hdistrGoodTr";
 302   TH1F* hdistrGoodTr=new TH1F(hname.Data(),"Distribution of number of 'good' tracks per event;no.good-tracks/ev;Entries",4000,-0.5,3999.5);
 303   hdistrGoodTr->SetTitleOffset(1.3,"Y");
 304
 305   hname="hNtracklets";
 306   TH1F* hNtracklets=new TH1F(hname.Data(),"Number of tracklets;ntracklets;Entries",5000,-0.5,4999.5);
 307
 308   hname="hMult";
 309   TH1F* hMult=new TH1F(hname.Data(),"Multiplicity;multiplicity;Entries",10000,-0.5,9999.5);
 310
 311   hname="hd0";
 312   TH1F* hd0=new TH1F(hname.Data(),"Impact parameter distribution of 'good' tracks;d_{0}[cm];Entries/10^{3} cm",200,-0.1,0.1);
 313
 314   fOutputTrack->Add(hnClsITS);
 315   fOutputTrack->Add(hnClsITSSA);
 316   fOutputTrack->Add(hnClsSPD);
 317   fOutputTrack->Add(hptGoodTr);
 318   fOutputTrack->Add(hdistrGoodTr);
 319   fOutputTrack->Add(hNtracklets);
 320   fOutputTrack->Add(hMult);
 321   fOutputTrack->Add(hd0);
 322
 323   // Post the data
 324   PostData(1,fNEntries);
 325   PostData(2,fOutputPID);
 326   PostData(3,fOutputTrack);
 327   PostData(4,fCuts);
 328 }
 329
 330 //___________________________________________________________________________
 331 void AliAnalysisTaskSEHFQA::UserExec(Option_t */*option*/)
 332 {
 333   // Execute analysis for current event
 334
 335   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 336   if(fDebug>2) printf("Analysing decay %d\n",fDecayChannel);
 337   // Post the data already here
 338   PostData(1,fNEntries);
 339   PostData(2,fOutputPID);
 340   PostData(3,fOutputTrack);
 341   PostData(4,fCuts);
 342
 343   TClonesArray *arrayProng =0;
 344   Int_t pdg=0;
 345   if(!aod && AODEvent() && IsStandardAOD()) {
 346     // In case there is an AOD handler writing a standard AOD, use the AOD
 347     // event in memory rather than the input (ESD) event.
 348     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 349     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 350     // have to taken from the AOD event hold by the AliAODExtension
 351     AliAODHandler* aodHandler = (AliAODHandler*)
 352       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 353     if(aodHandler->GetExtensions()) {
 354
 355       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 356       AliAODEvent *aodFromExt = ext->GetAOD();
 357
 358
 359
 360       switch(fDecayChannel){
 361       case 0:
 362         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 363         pdg=411;
 364         break;
 365       case 1:
 366         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 367         pdg=421;
 368         break;
 369       case 2:
 370         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 371         pdg=413;
 372         break;
 373       case 3:
 374         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 375         pdg=431;
 376         break;
 377       case 4:
 378         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
 379         pdg=421;
 380         break;
 381       case 5:
 382         arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 383         pdg=4122;
 384         break;
 385       }
 386     }
 387   } else {
 388     switch(fDecayChannel){
 389     case 0:
 390       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 391       pdg=411;
 392       break;
 393     case 1:
 394       arrayProng=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
 395       pdg=421;
 396       break;
 397     case 2:
 398       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Dstar");
 399       pdg=413;
 400       break;
 401     case 3:
 402       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 403       pdg=431;
 404       break;
 405     case 4:
 406       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm4Prong");
 407       pdg=421;
 408       break;
 409     case 5:
 410       arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 411       pdg=4122;
 412       break;
 413     }
 414   }
 415   Bool_t isSimpleMode=kFALSE;
 416   if(!arrayProng) {
 417     AliInfo("Branch not found! The output will contain only trak related histograms\n");
 418     isSimpleMode=kTRUE;
 419     fNEntries->Fill(2);
 420   }
 421
 422   // fix for temporary bug in ESDfilter
 423   // the AODs with null vertex pointer didn't pass the PhysSel
 424   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 425
 426   // count event
 427   fNEntries->Fill(0);
 428
 429   //count events with good vertex
 430   // AOD primary vertex
 431   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 432   TString primTitle = vtx1->GetTitle();
 433   if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) fNEntries->Fill(4);
 434
 435   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
 436   TString trigclass=aod->GetFiredTriggerClasses();
 437   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNEntries->Fill(5);
 438
 439
 440   //select event
 441   if(!fCuts->IsEventSelected(aod)) {
 442     // rejected for pileup
 443     if(fCuts->GetWhyRejection()==1) fNEntries->Fill(1);
 444     return;
 445   }
 446
 447   Int_t ntracks=0;
 448   Int_t isGoodTrack=0;
 449
 450   if(aod) ntracks=aod->GetNTracks();
 451
 452
 453   ((TH1F*)fOutputTrack->FindObject("hNtracklets"))->Fill(aod->GetTracklets()->GetNumberOfTracklets());
 454   ((TH1F*)fOutputTrack->FindObject("hMult"))->Fill(aod->GetHeader()->GetRefMultiplicity());
 455
 456
 457   //loop on tracks in the event
 458   for (Int_t k=0;k<ntracks;k++){
 459     AliAODTrack* track=aod->GetTrack(k);
 460     AliAODPidHF* pidHF=fCuts->GetPidHF();
 461     AliAODPid *pid = track->GetDetPid();
 462     if(!pid)  {if (fDebug>1)cout<<"No AliAODPid found"<<endl; continue;}
 463
 464     Double_t times[AliPID::kSPECIES];
 465     pid->GetIntegratedTimes(times);
 466
 467     //check TOF
 468     if(pidHF && pidHF->CheckStatus(track,"TOF")){
 469       ((TH1F*)fOutputPID->FindObject("hTOFtime"))->Fill(times[AliPID::kProton]);
 470       ((TH2F*)fOutputPID->FindObject("hTOFtimeKaonHyptime"))->Fill(track->P(),pid->GetTOFsignal()-times[3]); //3 is kaon
 471       ((TH1F*)fOutputPID->FindObject("hTOFsig"))->Fill(pid->GetTOFsignal());
 472       if (pid->GetTOFsignal()< 0) ((TH1F*)fOutputPID->FindObject("hTOFsig"))->Fill(-1);
 473
 474       ((TH2F*)fOutputPID->FindObject("hTOFsigmaK160"))->Fill(track->P(),(pid->GetTOFsignal()-times[AliPID::kKaon])/160);
 475       ((TH2F*)fOutputPID->FindObject("hTOFsigmaPion160"))->Fill(track->P(),(pid->GetTOFsignal()-times[AliPID::kPion])/160);
 476       ((TH2F*)fOutputPID->FindObject("hTOFsigmaProton160"))->Fill(track->P(),(pid->GetTOFsignal()-times[AliPID::kProton])/160);
 477       ((TH2F*)fOutputPID->FindObject("hTOFsigmaK120"))->Fill(track->P(),(pid->GetTOFsignal()-times[AliPID::kKaon])/120);
 478       ((TH2F*)fOutputPID->FindObject("hTOFsigmaPion120"))->Fill(track->P(),(pid->GetTOFsignal()-times[AliPID::kPion])/120);
 479       ((TH2F*)fOutputPID->FindObject("hTOFsigmaProton120"))->Fill(track->P(),(pid->GetTOFsignal()-times[AliPID::kProton])/120);
 480
 481     }//if TOF status
 482
 483     if(pidHF && pidHF->CheckStatus(track,"TPC")){
 484       Double_t alephParameters[5];
 485       //this is recommended for LHC10d
 486       alephParameters[0] = 1.34490e+00/50.;
 487       alephParameters[1] =  2.69455e+01;
 488       alephParameters[2] =  TMath::Exp(-2.97552e+01);
 489       alephParameters[3] = 2.35339e+00;
 490       alephParameters[4] = 5.98079e+00;
 491       /*
 492       //this is recommended for LHC10bc
 493       alephParameters[0] = 0.0283086/0.97;
 494       alephParameters[1] = 2.63394e+01;
 495       alephParameters[2] = 5.04114e-11;
 496       alephParameters[3] = 2.12543e+00;
 497       alephParameters[4] = 4.88663e+00;
 498       */
 499       AliTPCPIDResponse* tpcres=new AliTPCPIDResponse();
 500       tpcres->SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2],alephParameters[3],alephParameters[4]);
 501       Double_t TPCp=pid->GetTPCmomentum();
 502       Double_t TPCsignal=pid->GetTPCsignal();
 503       ((TH1F*)fOutputPID->FindObject("hTPCsig"))->Fill(TPCsignal);
 504       ((TH1F*)fOutputPID->FindObject("hTPCsigvsp"))->Fill(TPCp,TPCsignal);
 505       ((TH2F*)fOutputPID->FindObject("hTPCsigmaK"))->Fill(TPCp,tpcres->GetNumberOfSigmas(TPCp,TPCsignal,track->GetTPCNcls(),AliPID::kKaon));
 506       ((TH2F*)fOutputPID->FindObject("hTPCsigmaPion"))->Fill(TPCp,tpcres->GetNumberOfSigmas(TPCp,TPCsignal,track->GetTPCNcls(),AliPID::kPion));
 507       ((TH2F*)fOutputPID->FindObject("hTPCsigmaProton"))->Fill(TPCp,tpcres->GetNumberOfSigmas(TPCp,TPCsignal,track->GetTPCNcls(),AliPID::kProton));
 508       delete tpcres;
 509     }//if TPC status
 510
 511     //check clusters of the tracks
 512     Int_t nclsTot=0,nclsSPD=0;
 513
 514     for(Int_t l=0;l<6;l++) {
 515       if(TESTBIT(track->GetITSClusterMap(),l)) {
 516         nclsTot++; if(l<2) nclsSPD++;
 517       }
 518     }
 519     ((TH1F*)fOutputTrack->FindObject("hnClsITS"))->Fill(nclsTot);
 520     ((TH1F*)fOutputTrack->FindObject("hnClsSPD"))->Fill(nclsSPD);
 521
 522     if(!(track->GetStatus()&AliESDtrack::kTPCin) && track->GetStatus()&AliESDtrack::kITSrefit && !(track->GetStatus()&AliESDtrack::kITSpureSA)){//tracks retrieved in the ITS and not reconstructed in the TPC
 523       ((TH1F*)fOutputTrack->FindObject("hnClsITS-SA"))->Fill(nclsTot);
 524     }
 525
 526     if(isSimpleMode){
 527
 528       if (track->Pt()>0.3 &&
 529           track->GetStatus()&AliESDtrack::kTPCrefit &&
 530           track->GetStatus()&AliESDtrack::kITSrefit &&
 531           /*nclsTot>3 &&*/
 532           nclsSPD>0) {//fill hist good tracks
 533
 534         ((TH1F*)fOutputTrack->FindObject("hptGoodTr"))->Fill(track->Pt());
 535
 536         isGoodTrack++;
 537
 538         ((TH1F*)fOutputTrack->FindObject("hdistrGoodTr"))->Fill(isGoodTrack);
 539       }
 540     }//simple mode: no IsSelected on tracks: use "manual" cuts
 541
 542   } //end loop on tracks
 543
 544   if(!isSimpleMode){
 545   // loop over candidates
 546   Int_t nCand = arrayProng->GetEntriesFast();
 547   Int_t ndaugh=3;
 548   if(fDecayChannel==AliAnalysisTaskSEHFQA::kD0toKpi) ndaugh=2;
 549   if(fDecayChannel==AliAnalysisTaskSEHFQA::kD0toKpipipi) ndaugh=4;
 550
 551   for (Int_t iCand = 0; iCand < nCand; iCand++) {
 552     AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayProng->UncheckedAt(iCand);
 553
 554     for(Int_t id=0;id<ndaugh;id++){
 555
 556       //other histograms to be filled when the cut object is given
 557       AliAODTrack* track=(AliAODTrack*)d->GetDaughter(id);
 558
 559       //track quality
 560
 561       if (fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(pdg)) && fCuts->IsSelected(d,AliRDHFCuts::kTracks,aod)) {
 562
 563         ((TH1F*)fOutputTrack->FindObject("hptGoodTr"))->Fill(track->Pt());
 564         isGoodTrack++;
 565
 566         ((TH1F*)fOutputTrack->FindObject("hdistrGoodTr"))->Fill(isGoodTrack);
 567
 568         ((TH1F*)fOutputTrack->FindObject("hd0"))->Fill(d->Getd0Prong(id));
 569
 570         if (fCuts->IsSelected(d,AliRDHFCuts::kAll,aod)){
 571
 572           AliAODPid *pid = track->GetDetPid();
 573           AliAODPidHF* pidHF=fCuts->GetPidHF();
 574           Double_t times[5];
 575           pid->GetIntegratedTimes(times);
 576           if(pidHF && pidHF->CheckStatus(track,"TOF")) ((TH2F*)fOutputPID->FindObject("hTOFtimeKaonHyptimeAC"))->Fill(track->P(),pid->GetTOFsignal()-times[AliPID::kKaon]);
 577           if(pidHF && pidHF->CheckStatus(track,"TPC")) ((TH2F*)fOutputPID->FindObject("hTPCsigvspAC"))->Fill(pid->GetTPCmomentum(),pid->GetTPCsignal());
 578
 579           fNEntries->Fill(3);
 580         } //end analysis cuts
 581       } //end acceptance and track cuts
 582     } //end loop on tracks in the candidate
 583   } //end loop on candidates
 584   }
 585 }
 586
 587   //____________________________________________________________________________
 588 void AliAnalysisTaskSEHFQA::Terminate(Option_t */*option*/){
 589   //terminate analysis
 590
 591   fNEntries = dynamic_cast<TH1F*>(GetOutputData(1));
 592   if(!fNEntries){
 593     printf("ERROR: %s not available\n",GetOutputSlot(1)->GetContainer()->GetName());
 594     return;
 595   }
 596
 597   fOutputPID = dynamic_cast<TList*> (GetOutputData(2));
 598   if (!fOutputPID) {
 599     printf("ERROR: %s not available\n",GetOutputSlot(2)->GetContainer()->GetName());
 600     return;
 601   }
 602
 603   fOutputTrack = dynamic_cast<TList*> (GetOutputData(3));
 604   if (!fOutputTrack) {
 605     printf("ERROR: %s not available\n",GetOutputSlot(3)->GetContainer()->GetName());
 606     return;
 607   }
 608
 609 }
 610