PWG3/vertexingHF/AliCFHeavyFlavourTask.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 // Example of task
  18 // which provides standard way of calculating acceptance and efficiency
  19 // between different steps of the procedure.
  20 // The ouptut of the task is a AliCFContainer from which the efficiencies
  21 // can be calculated
  22 // Applied on AliAODRecoDecayHF2Prong particles (D0->Kpi)
  23 //-----------------------------------------------------------------------
  24 // Author : C. Zampolli starting from an example from R. Vernet
  25 //-----------------------------------------------------------------------
  26
  27 #include <TH1I.h>
  28 #include <TCanvas.h>
  29 #include <TStyle.h>
  30
  31 #include "AliCFHeavyFlavourTask.h"
  32 #include "AliCFContainer.h"
  33 #include "AliLog.h"
  34 #include "AliAODEvent.h"
  35 #include "AliAODRecoDecayHF2Prong.h"
  36 #include "AliAODMCParticle.h"
  37 #include "AliCFManager.h"
  38
  39 //__________________________________________________________________________
  40 AliCFHeavyFlavourTask::AliCFHeavyFlavourTask() :
  41         AliAnalysisTaskSE(),
  42         fPDG(0),
  43         fCFManager(0x0),
  44         fHistEventsProcessed(0x0),
  45         fCountMC(0),
  46         fEvents(0)
  47 {
  48         //
  49         //Default ctor
  50         //
  51 }
  52 //___________________________________________________________________________
  53 AliCFHeavyFlavourTask::AliCFHeavyFlavourTask(const Char_t* name) :
  54         AliAnalysisTaskSE(name),
  55         fPDG(0),
  56         fCFManager(0x0),
  57         fHistEventsProcessed(0x0),
  58         fCountMC(0),
  59         fEvents(0)
  60 {
  61         //
  62         // Constructor. Initialization of Inputs and Outputs
  63         //
  64         /*
  65           DefineInput(0) and DefineOutput(0)
  66           are taken care of by AliAnalysisTaskSE constructor
  67         */
  68         DefineOutput(1,TH1I::Class());
  69         DefineOutput(2,AliCFContainer::Class());
  70 }
  71
  72 //___________________________________________________________________________
  73 AliCFHeavyFlavourTask& AliCFHeavyFlavourTask::operator=(const AliCFHeavyFlavourTask& c)
  74 {
  75         //
  76         // Assignment operator
  77         //
  78         if (this!=&c) {
  79                 AliAnalysisTaskSE::operator=(c) ;
  80                 fPDG      = c.fPDG;
  81                 fCFManager  = c.fCFManager;
  82                 fHistEventsProcessed = c.fHistEventsProcessed;
  83         }
  84         return *this;
  85 }
  86
  87 //___________________________________________________________________________
  88 AliCFHeavyFlavourTask::AliCFHeavyFlavourTask(const AliCFHeavyFlavourTask& c) :
  89         AliAnalysisTaskSE(c),
  90         fPDG(c.fPDG),
  91         fCFManager(c.fCFManager),
  92         fHistEventsProcessed(c.fHistEventsProcessed),
  93         fCountMC(c.fCountMC),
  94         fEvents(c.fEvents)
  95 {
  96         //
  97         // Copy Constructor
  98         //
  99 }
 100
 101 //___________________________________________________________________________
 102 AliCFHeavyFlavourTask::~AliCFHeavyFlavourTask() {
 103         //
 104         //destructor
 105         //
 106         if (fCFManager)           delete fCFManager ;
 107         if (fHistEventsProcessed) delete fHistEventsProcessed ;
 108 }
 109
 110 //_________________________________________________
 111 void AliCFHeavyFlavourTask::UserExec(Option_t *)
 112 {
 113         //
 114         // Main loop function
 115         //
 116
 117         if (!fInputEvent) {
 118                 Error("UserExec","NO EVENT FOUND!");
 119                 return;
 120         }
 121
 122         fEvents++;
 123         if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
 124         AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 125         fCFManager->SetEventInfo(aodEvent);
 126
 127         // MC-event selection
 128         Double_t containerInput[2] ;
 129
 130         // loop on the MC event
 131         TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 132         if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
 133         Int_t icountMC = 0;
 134
 135         for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
 136                 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
 137                 if (!mcPart) AliWarning("Particle not found in tree");
 138
 139                 // check the MC-level cuts
 140                 if (!fCFManager->CheckParticleCuts(0,mcPart)) continue; // 0 stands for MC level
 141
 142                 // check whether the D0 decays in pi+K
 143                 // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 144                 // could use a cut in the CF, but not clear how to define a TDecayChannel
 145                 // implemented for the time being as a cut on the number of daughters - see later when
 146                 // getting the daughters
 147
 148                 // getting the daughters:
 149                 // AliMCParticle::GetDaughter(0) returns the index of the 1st daughter
 150                 // AliMCParticle::GetDaughter(1) returns the index of the last daughter
 151                 // so when a particle has 2 daughters, the difference must be 1
 152                 Int_t daughter0 = mcPart->GetDaughter(0);
 153                 Int_t daughter1 = mcPart->GetDaughter(1);
 154                 if (daughter0 == 0 || daughter1 == 0) {
 155                         AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
 156                         continue;
 157                 }
 158                 if (TMath::Abs(daughter1 - daughter0 != 1)) {
 159                         AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
 160                         continue;
 161                 }
 162                 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
 163                 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
 164                 if (!mcPartDaughter0 || !mcPartDaughter1) {
 165                         AliWarning("At least one Daughter Particle not found in tree, skipping");
 166                         continue;
 167                 }
 168
 169                 // fill the container for Gen-level selection
 170                 containerInput[0] = mcPart->Pt();
 171                 containerInput[1] = mcPart->Y();
 172                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
 173                 icountMC++;
 174         }
 175
 176         AliDebug(2,Form("Found %i MC particles that are D0!!",icountMC));
 177         // Now go to rec level
 178         fCountMC += icountMC;
 179         // load D0->Kpi candidates
 180         TClonesArray *arrayD0toKpi = (TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
 181
 182         if (!arrayD0toKpi) AliError("Could not find array of D0->Kpi candidates");
 183
 184         AliDebug(2, Form("Found %d vertices",arrayD0toKpi->GetEntriesFast()));
 185
 186         for (Int_t iD0 = 0; iD0<arrayD0toKpi->GetEntriesFast(); iD0++) {
 187
 188                 AliAODRecoDecayHF2Prong* rd = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->At(iD0);
 189
 190                 // cuts can't be applied to RecoDecays particles
 191                 // if (!fCFManager->CheckParticleCuts(1, rd)) continue;  // 1 stands for AOD level
 192
 193                 // check if associated MC v0 passes the cuts
 194                 // first get the label
 195                 Int_t mcLabel = rd->MatchToMC(421,mcArray);
 196                 if (mcLabel == -1)
 197                         {
 198                                 AliDebug(2,"No MC particle found");
 199                         }
 200                 else {
 201                         AliAODMCParticle* mcVtxHF = (AliAODMCParticle*)mcArray->At(mcLabel);
 202                         if (!mcVtxHF) {
 203                                 AliWarning("Could not find associated MC in AOD MC tree");
 204                                 continue;
 205                         }
 206
 207                         AliDebug(2, Form("pdg code from MC = %d",TMath::Abs(mcVtxHF->GetPdgCode())));
 208                         if (!fCFManager->CheckParticleCuts(0, mcVtxHF)) {  // 0 stands for MC level
 209                                 continue;
 210                         }
 211
 212                         //fill the container
 213                         Double_t pt = rd->Pt();
 214                         Double_t rapidity = rd->YD0();
 215
 216                         AliDebug(2, Form("Filling the container with pt = %f and rapidity = %f", pt, rapidity));
 217                         containerInput[0] = pt ;
 218                         containerInput[1] = rapidity ;
 219                         fCFManager->GetParticleContainer()->Fill(containerInput,1) ;
 220                 }
 221         }
 222
 223         fHistEventsProcessed->Fill(0);
 224         // PostData(0) is taken care of by AliAnalysisTaskSE
 225         PostData(1,fHistEventsProcessed) ;
 226         PostData(2,fCFManager->GetParticleContainer()) ;
 227 }
 228
 229
 230 //___________________________________________________________________________
 231 void AliCFHeavyFlavourTask::Terminate(Option_t*)
 232 {
 233         // The Terminate() function is the last function to be called during
 234         // a query. It always runs on the client, it can be used to present
 235         // the results graphically or save the results to file.
 236
 237         AliAnalysisTaskSE::Terminate();
 238
 239         AliInfo(Form("Found %i MC particles that are D0 in MC in %d events",fCountMC,fEvents));
 240
 241         //draw some example plots....
 242
 243         AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
 244
 245         // projecting the containers to obtain histograms
 246         // first argument = variable, second argument = step
 247         TH1D* h00 =   cont->ShowProjection(0,0) ;  // pt, MC
 248         TH1D* h01 =   cont->ShowProjection(0,1) ;  // pt, Reco
 249
 250         TH1D* h10 =   cont->ShowProjection(1,0) ;  // rapidity, MC
 251         TH1D* h11 =   cont->ShowProjection(1,1) ;  // rapidity, Reco
 252
 253         Double_t max1 = h00->GetMaximum();
 254         Double_t max2 = h10->GetMaximum();
 255
 256         // MC histos
 257         h00->SetTitle("p_{T} (GeV/c)");
 258         h10->SetTitle("rapidity");
 259         h00->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 260         h10->GetXaxis()->SetTitle("rapidity");
 261         h00->GetYaxis()->SetRangeUser(0,max1*1.2);
 262         h10->GetYaxis()->SetRangeUser(0,max2*1.2);
 263         h00->SetMarkerStyle(20) ;
 264         h10->SetMarkerStyle(24) ;
 265         h00->SetMarkerColor(2);
 266         h10->SetMarkerColor(2);
 267
 268         // Reco histos
 269         h01->SetTitle("p_{T} (GeV/c)");
 270         h11->SetTitle("rapidity");
 271         h01->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 272         h11->GetXaxis()->SetTitle("rapidity");
 273         h01->GetYaxis()->SetRangeUser(0,max1*1.2);
 274         h11->GetYaxis()->SetRangeUser(0,max2*1.2);
 275         h01->SetMarkerStyle(20) ;
 276         h11->SetMarkerStyle(24) ;
 277         h01->SetMarkerColor(4);
 278         h11->SetMarkerColor(4);
 279
 280         gStyle->SetCanvasColor(0);
 281         gStyle->SetFrameFillColor(0);
 282         gStyle->SetTitleFillColor(0);
 283         gStyle->SetStatColor(0);
 284
 285         TCanvas * c =new TCanvas("c","",1000,600);
 286         c->Divide(2,2);
 287
 288         c->cd(1);
 289         h00->Draw("p");
 290         c->cd(2);
 291         h01->Draw("p");
 292         c->cd(3);
 293         h10->Draw("p");
 294         c->cd(4);
 295         h11->Draw("p");
 296         c->cd();
 297
 298         c->SaveAs("plots.eps");
 299 }
 300
 301 //___________________________________________________________________________
 302 void AliCFHeavyFlavourTask::UserCreateOutputObjects() {
 303         //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 304         //TO BE SET BEFORE THE EXECUTION OF THE TASK
 305         //
 306         Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 307
 308         //slot #1
 309         OpenFile(1);
 310         fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
 311 }
 312