PWG3/vertexingHF/AliAnalysisTaskSEDs.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-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 /* $Id: AliITSCorrMapSDD.cxx 32906 2009-06-12 16:56:53Z prino $ */
  17
  18 ///////////////////////////////////////////////////////////////////
  19 //                                                               //
  20 //  Analysis task to produce Ds candidates mass spectra          //
  21 // Origin: F.Prino, Torino, prino@to.infn.it                     //
  22 //                                                               //
  23 ///////////////////////////////////////////////////////////////////
  24
  25 #include <TClonesArray.h>
  26 #include <TNtuple.h>
  27 #include <TList.h>
  28 #include <TString.h>
  29 #include <TH1F.h>
  30 #include <TMath.h>
  31 #include <TDatabasePDG.h>
  32
  33 #include "AliAnalysisManager.h"
  34 #include "AliAODHandler.h"
  35 #include "AliAODEvent.h"
  36 #include "AliAODVertex.h"
  37 #include "AliAODTrack.h"
  38 #include "AliAODMCHeader.h"
  39 #include "AliAODMCParticle.h"
  40 #include "AliAODRecoDecayHF3Prong.h"
  41 #include "AliAnalysisVertexingHF.h"
  42 #include "AliAnalysisTaskSE.h"
  43 #include "AliAnalysisTaskSEDs.h"
  44
  45 ClassImp(AliAnalysisTaskSEDs)
  46
  47
  48 //________________________________________________________________________
  49 AliAnalysisTaskSEDs::AliAnalysisTaskSEDs():
  50   AliAnalysisTaskSE(),
  51   fOutput(0),
  52   fHistNEvents(0),
  53   fReadMC(kFALSE),
  54   fNPtBins(0),
  55   fMassRange(0.2),
  56   fVHF(0)
  57 {
  58   // Default constructor
  59 }
  60
  61 //________________________________________________________________________
  62 AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char *name):
  63   AliAnalysisTaskSE(name),
  64   fOutput(0),
  65   fHistNEvents(0),
  66   fReadMC(kFALSE),
  67   fNPtBins(0),
  68   fMassRange(0.2),
  69   fVHF(0)
  70 {
  71   // Default constructor
  72   // Output slot #1 writes into a TList container
  73   Double_t ptlim[6]={0.,1.,3.,5.,10.,9999999.};
  74   SetPtBins(5,ptlim);
  75
  76   DefineOutput(1,TList::Class());  //My private output
  77
  78 }
  79
  80 //________________________________________________________________________
  81 void AliAnalysisTaskSEDs::SetPtBins(Int_t n, Double_t* lim){
  82   // define pt bins for analysis
  83   if(n>kMaxPtBins){
  84     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
  85     fNPtBins=kMaxPtBins;
  86     fPtLimits[0]=0.;
  87     fPtLimits[1]=1.;
  88     fPtLimits[2]=3.;
  89     fPtLimits[3]=5.;
  90     fPtLimits[4]=10.;
  91     for(Int_t i=5; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
  92   }else{
  93     fNPtBins=n;
  94     for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
  95     for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
  96   }
  97   if(fDebug > 1){
  98     printf("Number of Pt bins = %d\n",fNPtBins);
  99     for(Int_t i=0; i<fNPtBins+1; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fPtLimits[i],fPtLimits[i+1]);
 100   }
 101 }
 102 //________________________________________________________________________
 103 AliAnalysisTaskSEDs::~AliAnalysisTaskSEDs()
 104 {
 105   // Destructor
 106   if (fOutput) {
 107     delete fOutput;
 108     fOutput = 0;
 109   }
 110   if (fVHF) {
 111     delete fVHF;
 112     fVHF = 0;
 113   }
 114 }
 115
 116 //________________________________________________________________________
 117 void AliAnalysisTaskSEDs::Init()
 118 {
 119   // Initialization
 120
 121   if(fDebug > 1) printf("AnalysisTaskSEDs::Init() \n");
 122
 123   gROOT->LoadMacro("ConfigVertexingHF.C");
 124
 125   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
 126   fVHF->PrintStatus();
 127
 128   return;
 129 }
 130
 131 //________________________________________________________________________
 132 void AliAnalysisTaskSEDs::UserCreateOutputObjects()
 133 {
 134   // Create the output container
 135   //
 136   if(fDebug > 1) printf("AnalysisTaskSEDs::UserCreateOutputObjects() \n");
 137
 138   // Several histograms are more conveniently managed in a TList
 139   fOutput = new TList();
 140   fOutput->SetOwner();
 141   fOutput->SetName("OutputHistos");
 142
 143   Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
 144   Double_t minMass=massDs-0.5*fMassRange;
 145   Double_t maxMass=massDs+0.5*fMassRange;
 146   TString hisname;
 147   Int_t index;
 148   for(Int_t i=0;i<fNPtBins;i++){
 149     index=GetHistoIndex(i);
 150     hisname.Form("hMassAllPt%d",i);
 151     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
 152     fMassHist[index]->Sumw2();
 153     hisname.Form("hMassAllPt%dCuts",i);
 154     fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
 155     fMassHistCuts[index]->Sumw2();
 156     hisname.Form("hCosPAllPt%d",i);
 157     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 158     fCosPHist[index]->Sumw2();
 159     hisname.Form("hDLenAllPt%d",i);
 160     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 161     fDLenHist[index]->Sumw2();
 162     hisname.Form("hDalitzAllPt%d",i);
 163     fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
 164     fDalitz[index]->Sumw2();
 165
 166     index=GetSignalHistoIndex(i);
 167     hisname.Form("hMassSigPt%d",i);
 168     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
 169     fMassHist[index]->Sumw2();
 170     hisname.Form("hMassSigPt%dCuts",i);
 171     fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
 172     fMassHistCuts[index]->Sumw2();
 173     hisname.Form("hCosPSigPt%d",i);
 174     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 175     fCosPHist[index]->Sumw2();
 176     hisname.Form("hDLenSigPt%d",i);
 177     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 178     fDLenHist[index]->Sumw2();
 179     hisname.Form("hDalitzSigPt%d",i);
 180     fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
 181     fDalitz[index]->Sumw2();
 182
 183     hisname.Form("hMassBkgPt%d",i);
 184     index=GetBackgroundHistoIndex(i);
 185     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
 186     fMassHist[index]->Sumw2();
 187     hisname.Form("hMassBkgPt%dCuts",i);
 188     fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
 189     fMassHistCuts[index]->Sumw2();
 190     hisname.Form("hCosPBkgPt%d",i);
 191     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 192     fCosPHist[index]->Sumw2();
 193     hisname.Form("hDLenBkgPt%d",i);
 194     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 195     fDLenHist[index]->Sumw2();
 196     hisname.Form("hDalitzBkgPt%d",i);
 197     fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
 198     fDalitz[index]->Sumw2();
 199   }
 200
 201   for(Int_t i=0; i<3*fNPtBins; i++){
 202     fOutput->Add(fMassHist[i]);
 203     fOutput->Add(fMassHistCuts[i]);
 204     fOutput->Add(fCosPHist[i]);
 205     fOutput->Add(fDLenHist[i]);
 206     fOutput->Add(fDalitz[i]);
 207   }
 208
 209   fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",4,-0.5,3.5);
 210   fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",4,-0.5,3.5);
 211   fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",4,-0.5,3.5);
 212   fChanHistCuts[0] = new TH1F("hChanAllCuts", "KKpi and piKK candidates",4,-0.5,3.5);
 213   fChanHistCuts[1] = new TH1F("hChanSigCuts", "KKpi and piKK candidates",4,-0.5,3.5);
 214   fChanHistCuts[2] = new TH1F("hChanBkgCuts", "KKpi and piKK candidates",4,-0.5,3.5);
 215   for(Int_t i=0;i<3;i++){
 216     fChanHist[i]->Sumw2();
 217     fChanHist[i]->SetMinimum(0);
 218     fChanHistCuts[i]->Sumw2();
 219     fChanHistCuts[i]->SetMinimum(0);
 220     fOutput->Add(fChanHist[i]);
 221     fOutput->Add(fChanHistCuts[i]);
 222   }
 223
 224   fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-1.5,1.5);
 225   fHistNEvents->Sumw2();
 226   fHistNEvents->SetMinimum(0);
 227   fOutput->Add(fHistNEvents);
 228
 229
 230   return;
 231 }
 232
 233 //________________________________________________________________________
 234 void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)
 235 {
 236   // Ds selection for current event, fill mass histos and selecetion variable histo
 237   // separate signal and backgound if fReadMC is activated
 238
 239   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 240   fHistNEvents->Fill(0); // count event
 241   // Post the data already here
 242   PostData(1,fOutput);
 243
 244
 245   TClonesArray *array3Prong = 0;
 246   if(!aod && AODEvent() && IsStandardAOD()) {
 247     // In case there is an AOD handler writing a standard AOD, use the AOD
 248     // event in memory rather than the input (ESD) event.
 249     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 250     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 251     // have to taken from the AOD event hold by the AliAODExtension
 252     AliAODHandler* aodHandler = (AliAODHandler*)
 253       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 254     if(aodHandler->GetExtensions()) {
 255       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 256       AliAODEvent *aodFromExt = ext->GetAOD();
 257       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 258     }
 259   } else {
 260     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 261   }
 262
 263   if(!array3Prong) {
 264     printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
 265     return;
 266   }
 267
 268
 269   TClonesArray *arrayMC=0;
 270   AliAODMCHeader *mcHeader=0;
 271
 272   // AOD primary vertex
 273   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 274   //    vtx1->Print();
 275
 276   // load MC particles
 277   if(fReadMC){
 278
 279     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 280     if(!arrayMC) {
 281       printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
 282       return;
 283     }
 284
 285   // load MC header
 286     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 287     if(!mcHeader) {
 288       printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
 289       return;
 290     }
 291   }
 292
 293   Int_t n3Prong = array3Prong->GetEntriesFast();
 294   if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
 295
 296
 297   Int_t pdgDstoKKpi[3]={321,321,211};
 298   Double_t cutsDs[14]={0.2,0.4,0.4,0.,0.,0.005,0.038,0.,0.,0.95,0.,0.1,0.004,0.035};
 299   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 300     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 301
 302
 303     Bool_t unsetvtx=kFALSE;
 304     if(!d->GetOwnPrimaryVtx()){
 305       d->SetOwnPrimaryVtx(vtx1);
 306       unsetvtx=kTRUE;
 307     }
 308     Int_t isKKpi,ispiKK;
 309     Int_t isPhi,isK0star;
 310     Int_t isKKpiTC,ispiKKTC;
 311     Int_t isPhiTC,isK0starTC;
 312     if(d->SelectDs(fVHF->GetDsCuts(),isKKpi,ispiKK,isPhi,isK0star)) {
 313       Double_t ptCand = d->Pt();
 314       Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,ptCand);
 315       Bool_t passTightCuts=d->SelectDs(cutsDs,isKKpiTC,ispiKKTC,isPhiTC,isK0starTC);
 316       Int_t labDs=-1;
 317       if(fReadMC){
 318         labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
 319       }
 320
 321       Double_t dlen=d->DecayLength();
 322       Double_t cosp=d->CosPointingAngle();
 323       Int_t index=GetHistoIndex(iPtBin);
 324       Int_t type=0;
 325       if(isKKpi) type+=1;
 326       if(ispiKK) type+=2;
 327       Int_t typeTC=0;
 328       if(isKKpiTC) typeTC+=1;
 329       if(ispiKKTC) typeTC+=2;
 330       fCosPHist[index]->Fill(cosp);
 331       fDLenHist[index]->Fill(dlen);
 332       fChanHist[0]->Fill(type);
 333       if(passTightCuts) fChanHistCuts[0]->Fill(typeTC);
 334       if(fReadMC){
 335         if(labDs>=0) {
 336           index=GetSignalHistoIndex(iPtBin);
 337           fCosPHist[index]->Fill(cosp);
 338           fDLenHist[index]->Fill(dlen);
 339           fChanHist[1]->Fill(type);
 340           if(passTightCuts) fChanHistCuts[1]->Fill(typeTC);
 341         }else{
 342           index=GetBackgroundHistoIndex(iPtBin);
 343           fCosPHist[index]->Fill(cosp);
 344           fDLenHist[index]->Fill(dlen);
 345           fChanHist[2]->Fill(type);
 346           if(passTightCuts) fChanHistCuts[2]->Fill(typeTC);
 347         }
 348       }
 349       if(isKKpi){
 350         index=GetHistoIndex(iPtBin);
 351         Double_t invMass=d->InvMassDsKKpi();
 352         fMassHist[index]->Fill(invMass);
 353         Double_t mass01=d->InvMass2Prongs(0,1,321,321);
 354         Double_t mass12=d->InvMass2Prongs(1,2,321,211);
 355         fDalitz[index]->Fill(mass01,mass12);
 356         if(passTightCuts && isKKpiTC) fMassHistCuts[index]->Fill(invMass);
 357         if(fReadMC){
 358           if(labDs>=0) {
 359             index=GetSignalHistoIndex(iPtBin);
 360             fMassHist[index]->Fill(invMass);
 361             fDalitz[index]->Fill(mass01,mass12);
 362             if(passTightCuts&& isKKpiTC) fMassHistCuts[index]->Fill(invMass);
 363           }else{
 364             index=GetBackgroundHistoIndex(iPtBin);
 365             fMassHist[index]->Fill(invMass);
 366             fDalitz[index]->Fill(mass01,mass12);
 367             if(passTightCuts&& isKKpiTC) fMassHistCuts[index]->Fill(invMass);
 368           }
 369         }
 370       }
 371       if(ispiKK){
 372         index=GetHistoIndex(iPtBin);
 373         Double_t invMass=d->InvMassDspiKK();
 374         fMassHist[index]->Fill(invMass);
 375         if(passTightCuts && ispiKKTC) fMassHistCuts[index]->Fill(invMass);
 376         if(fReadMC){
 377           if(labDs>=0) {
 378             index=GetSignalHistoIndex(iPtBin);
 379             fMassHist[index]->Fill(invMass);
 380             if(passTightCuts && ispiKKTC) fMassHistCuts[index]->Fill(invMass);
 381           }else{
 382             index=GetBackgroundHistoIndex(iPtBin);
 383             fMassHist[index]->Fill(invMass);
 384             if(passTightCuts && ispiKKTC) fMassHistCuts[index]->Fill(invMass);
 385           }
 386         }
 387       }
 388     }
 389     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 390   }
 391
 392
 393   PostData(1,fOutput);
 394   return;
 395 }
 396
 397 //_________________________________________________________________
 398
 399 void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
 400 {
 401   // Terminate analysis
 402   //
 403   if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
 404   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 405   if (!fOutput) {
 406     printf("ERROR: fOutput not available\n");
 407     return;
 408   }
 409   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
 410   fChanHist[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAll"));
 411   fChanHist[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSig"));
 412   fChanHist[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkg"));
 413   fChanHistCuts[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAllCuts"));
 414   fChanHistCuts[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSigCuts"));
 415   fChanHistCuts[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkgCuts"));
 416
 417
 418   TString hisname;
 419   Int_t index;
 420   for(Int_t i=0;i<fNPtBins;i++){
 421
 422     index=GetHistoIndex(i);
 423     hisname.Form("hMassAllPt%d",i);
 424     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 425     hisname.Form("hMassAllPt%dCuts",i);
 426     fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 427     hisname.Form("hCosPAllPt%d",i);
 428     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 429     hisname.Form("hDLenAllPt%d",i);
 430     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 431     hisname.Form("hDalitzAllPt%d",i);
 432     fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
 433
 434     index=GetSignalHistoIndex(i);
 435     hisname.Form("hMassSigPt%d",i);
 436     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 437     hisname.Form("hMassSigPt%dCuts",i);
 438     fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 439     hisname.Form("hCosPSigPt%d",i);
 440     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 441     hisname.Form("hDLenSigPt%d",i);
 442     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 443     hisname.Form("hDalitzSigPt%d",i);
 444     fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
 445
 446     index=GetBackgroundHistoIndex(i);
 447     hisname.Form("hMassBkgPt%d",i);
 448     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 449     hisname.Form("hMassBkgPt%dCuts",i);
 450     fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 451     hisname.Form("hCosPBkgPt%d",i);
 452     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 453     hisname.Form("hDLenBkgPt%d",i);
 454     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
 455     hisname.Form("hDalitzBkgPt%d",i);
 456     fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
 457
 458   }
 459   return;
 460 }