PWGHF/correlationHF/AliDxHFECorrelation.cxx

   1 // $Id$
   2
   3 //**************************************************************************
   4 //* This file is property of and copyright by the ALICE Project            *
   5 //* ALICE Experiment at CERN, All rights reserved.                         *
   6 //*                                                                        *
   7 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
   8 //*                  Sedat Altinpinar <Sedat.Altinpinar@cern.ch>           *
   9 //*                  Hege Erdal       <hege.erdal@gmail.com>               *
  10 //*                                                                        *
  11 //* Permission to use, copy, modify and distribute this software and its   *
  12 //* documentation strictly for non-commercial purposes is hereby granted   *
  13 //* without fee, provided that the above copyright notice appears in all   *
  14 //* copies and that both the copyright notice and this permission notice   *
  15 //* appear in the supporting documentation. The authors make no claims     *
  16 //* about the suitability of this software for any purpose. It is          *
  17 //* provided "as is" without express or implied warranty.                  *
  18 //**************************************************************************
  19
  20 /// @file   AliDxHFECorrelation.cxx
  21 /// @author Sedat Altinpinar, Hege Erdal, Matthias Richter
  22 /// @date   2012-04-25
  23 /// @brief  Worker class for D0-HF electron correlation
  24 ///
  25
  26 #include "AliDxHFECorrelation.h"
  27 #include "AliVParticle.h"
  28 #include "AliLog.h"
  29 //#include "AliAnalysisCuts.h"         // required dependency libANALYSISalice.so
  30 //#include "AliFlowTrackSimple.h"      // required dependency libPWGflowBase.so
  31 //#include "AliFlowCandidateTrack.h"   // required dependency libPWGflowTasks.so
  32 //#include "AliCFContainer.h"          // required dependency libCORRFW.so
  33 #include "AliAODRecoDecayHF2Prong.h"   // libPWGHFvertexingHF
  34 #include "AliRDHFCutsD0toKpi.h"
  35 #include "TObjArray.h"
  36 #include "TH1D.h"
  37 #include "TH2D.h"
  38 #include "THnSparse.h"
  39 #include "TMath.h"
  40 #include "TFile.h"
  41 #include "TCanvas.h"
  42 #include "TDatabasePDG.h"
  43 #include "TLorentzVector.h"
  44 #include <iostream>
  45 #include <cerrno>
  46 #include <memory>
  47
  48 using namespace std;
  49
  50 ClassImp(AliDxHFECorrelation)
  51
  52 AliDxHFECorrelation::AliDxHFECorrelation(const char* name)
  53   : TNamed(name?name:"AliDxHFECorrelation", "")
  54   , fHistograms(NULL)
  55   , fControlObjects(NULL)
  56   , fCorrProperties(NULL)
  57   , fhEventControlCorr(NULL)
  58   , fCuts(NULL)
  59   , fUseMC(kFALSE)
  60 {
  61   // default constructor
  62   //
  63   //
  64
  65 }
  66
  67 const char* AliDxHFECorrelation::fgkEventControlBinNames[]={
  68   "nEventsAll",
  69   "nEventsSelected",
  70   "nEventsD0"
  71   "nEventsD0e"
  72 };
  73
  74 // TODO: maybe delete PtD0 in the future, note: there are more places!
  75 const char* AliDxHFECorrelation::fgkCorrControlBinNames[]={
  76   "D0InvMass",
  77   "PtD0",
  78   "PhiD0",
  79   "PtBinD0",
  80   "dPhi",
  81   "Pt electron",
  82 };
  83
  84 AliDxHFECorrelation::~AliDxHFECorrelation()
  85 {
  86   // destructor
  87   //
  88   //
  89   if (fHistograms) delete fHistograms;
  90   fHistograms=NULL;
  91
  92   // NOTE: fControlObjects owns the object, and they are deleted in the
  93   // destructor of TList
  94   if (fControlObjects) delete fControlObjects;
  95   fControlObjects=NULL;
  96   fCorrProperties=NULL;
  97   fhEventControlCorr=NULL;
  98
  99   // NOTE: the external object is deleted elsewhere
 100   fCuts=NULL;
 101 }
 102
 103 int AliDxHFECorrelation::Init()
 104 {
 105   AliInfo("Setting up THnSparse");
 106   // Are using THnSparse instead of histograms to store information
 107   // At the moment use old setup for THnSparse, change to new later
 108
 109   /*
 110     TString name;
 111     const int thnSize = 7;
 112     const double pi=TMath::Pi();
 113     //                           0      1     2      3      4      5    6
 114     //                        mass      Pt   Phi    ePt    ePhi  eEta  DeltaPhi
 115     int    thnBins[thnSize] = {   200,   1000,  100,   100,   100,  100,  180  };
 116     double thnMin [thnSize] = {  1.5648,   0,    0,     0,    0.0, -1.0,  0.0  };
 117     double thnMax [thnSize] = {  2.1648, 100,  (2*pi), 100, (2*pi), 1.0, (2*pi)};
 118
 119     name.Form("%s info", GetName());
 120     std::auto_ptr<THnSparseF> CorrProperties(new THnSparseF(name, name, thnSize, thnBins, thnMin, thnMax));
 121
 122     if (CorrProperties.get()==NULL) {
 123     return -ENOMEM;
 124     }
 125     int axis=0;
 126     CorrProperties->GetAxis(axis++)->SetTitle("D0 Inv Mass");
 127     CorrProperties->GetAxis(axis++)->SetTitle("D0 Pt");
 128     CorrProperties->GetAxis(axis++)->SetTitle("D0 Phi");
 129     CorrProperties->GetAxis(axis++)->SetTitle("electron Pt");
 130     CorrProperties->GetAxis(axis++)->SetTitle("electron Phi");
 131     CorrProperties->GetAxis(axis++)->SetTitle("electron Eta");
 132     CorrProperties->GetAxis(axis++)->SetTitle("#Delta#Phi D0 -HFE"); */
 133
 134
 135   // TODO: think about removing PtD0, but remember to change this in all places!
 136   static const int sizeEventdphi = 6;
 137   Double_t minPhi= -TMath::Pi()/2;
 138   Double_t maxPhi= 3*TMath::Pi()/2;
 139   //                                        0         1      2      3       4       5
 140   //                                      D0invmass PtD0   PhiD0  PtbinD0  dphi   Pte
 141   int    binsEventdphi[sizeEventdphi] = {   200,    1000,    100,     21,  100,    1000};
 142   double minEventdphi [sizeEventdphi] = { 1.5648,    0,      0,      0,    minPhi,   0 };
 143   double maxEventdphi [sizeEventdphi] = {  2.1648,  100, 2*(TMath::Pi()), 20,  maxPhi, 100 };
 144
 145   TString name;
 146   name.Form("%s info", GetName());
 147   std::auto_ptr<THnSparseF> CorrProperties(new THnSparseF(name, name, sizeEventdphi, binsEventdphi,minEventdphi , maxEventdphi));
 148   if (CorrProperties.get()==NULL) {
 149     return -ENOMEM;
 150   }
 151   int iLabel=0;
 152
 153   for (iLabel=0; iLabel<sizeEventdphi; iLabel++)
 154     CorrProperties->GetAxis(iLabel)->SetTitle(fgkCorrControlBinNames[iLabel]);
 155
 156   //----------------------------------------------
 157   // Histogram for storing event information
 158
 159   std::auto_ptr<TH1D> hEventControl(new TH1D("hEventControlCorr", "hEventControlCorr", 10, 0, 10));
 160   if (!hEventControl.get()) {
 161     return -ENOMEM;
 162   }
 163
 164   for (iLabel=0; iLabel<kNEventControlLabels; iLabel++)
 165     hEventControl->GetXaxis()->SetBinLabel(iLabel, fgkEventControlBinNames[iLabel]);
 166
 167   fCorrProperties=CorrProperties.release();
 168   AddControlObject(fCorrProperties);
 169   fhEventControlCorr=hEventControl.release();
 170   AddControlObject(fhEventControlCorr);
 171
 172   return 0;
 173 }
 174
 175 int AliDxHFECorrelation::AddControlObject(TObject* pObj)
 176 {
 177   AliInfo("Adding object");
 178   /// add control object to list, the base class becomes owner of the object
 179   if (!pObj) return -EINVAL;
 180   if (!fControlObjects) {
 181     fControlObjects=new TList;
 182     if (!fControlObjects) return -ENOMEM;
 183     fControlObjects->SetOwner();
 184   }
 185   if (fControlObjects->FindObject(pObj->GetName())) {
 186     AliError(Form("ignoring duplicate object '%s' of type %s", pObj->GetName(), pObj->ClassName()));
 187     return -EEXIST;
 188   }
 189   fControlObjects->Add(pObj);
 190   return 0;
 191 }
 192
 193 int AliDxHFECorrelation::HistogramEventProperties(int bin)
 194 {
 195   /// histogram event properties
 196   if (!fhEventControlCorr) return 0;
 197   fhEventControlCorr->Fill(bin);
 198
 199   return 0;
 200 }
 201
 202 int AliDxHFECorrelation::Fill(const TObjArray* triggerCandidates, const TObjArray* associatedTracks)
 203 {
 204   /// fill ThnSparse from array of AliVParticle objects
 205   if (!triggerCandidates || !associatedTracks) return -EINVAL;
 206   if (!fControlObjects) {
 207     Init();
 208   }
 209   if (!fControlObjects) {
 210     AliError("Initialisation failed, can not fill THnSparse");
 211   }
 212
 213   const double Pii=TMath::Pi();
 214
 215   TIter itrigger(triggerCandidates);
 216   TObject* otrigger=NULL;
 217   int ctrigger=-1;
 218
 219   // For the moment this is very specific to D0-electron correlation. Should be
 220   // changed to be more specific.
 221   while ((otrigger=itrigger())!=NULL) {
 222     // loop over trigger D0 particle
 223     ctrigger++;
 224     AliVParticle* ptrigger=dynamic_cast<AliVParticle*>(otrigger);
 225     AliAODRecoDecayHF2Prong *d0 = dynamic_cast<AliAODRecoDecayHF2Prong*>(ptrigger);
 226
 227     if (!ptrigger) continue;
 228
 229     TIter iElectron(associatedTracks);
 230     TObject* oElectron=NULL;
 231     int cElectron=-1;
 232     while ((oElectron=iElectron())!=NULL) {
 233       // loop over electrons
 234
 235       cElectron++;
 236       AliVParticle* pElectron=dynamic_cast<AliVParticle*>(oElectron);
 237       if (!pElectron) continue;
 238
 239       //Calculating dPhi using TLorentzVectors
 240       Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 241       TLorentzVector D0vector(0.,0.,0.,0.);
 242       TLorentzVector evector(0.,0.,0.,0.);
 243       D0vector.SetXYZM(d0->Px(),d0->Py(),d0->Pz(),mPDG);
 244       evector.SetXYZM(pElectron->Px(),pElectron->Py(),pElectron->Pz(),0.000511);
 245       Double_t DeltaPhi=D0vector.DeltaPhi(evector);
 246       if(DeltaPhi<-TMath::PiOver2()) DeltaPhi=DeltaPhi+(2*Pii);
 247
 248       /*Double_t CorrelationArray[]={d0->InvMassD0(),d0->Pt(),d0->Phi(),
 249         pElectron->Pt(),pElectron->Phi(), pElectron->Eta(),
 250         DeltaPhi};*/
 251
 252       // TODO: think about a method to retrieve the pt bin from the
 253       // selection object
 254       Int_t ptbin=fCuts->PtBin(d0->Pt());
 255       Double_t CorrelationArray[]={d0->InvMassD0(),d0->Pt(),d0->Phi(),ptbin,
 256                                    DeltaPhi, pElectron->Pt() };
 257       fCorrProperties->Fill(CorrelationArray);
 258
 259     } // loop over associated tracks
 260   } // loop over trigger particle
 261
 262   return 0;
 263 }
 264
 265 void AliDxHFECorrelation::Clear(Option_t * /*option*/)
 266 {
 267   /// overloaded from TObject: cleanup
 268
 269   // nothing to be done so far
 270   return TObject::Clear();
 271 }
 272
 273 void AliDxHFECorrelation::Print(Option_t */*option*/) const
 274 {
 275   /// overloaded from TObject: print info
 276   cout << "====================================================================" << endl;
 277   TObject::Print();
 278   if (fHistograms) {
 279     fHistograms->Print();
 280   }
 281 }
 282
 283 void AliDxHFECorrelation::Draw(Option_t */*option*/)
 284 {
 285   /// overloaded from TObject: draw histograms
 286 }
 287
 288 TObject* AliDxHFECorrelation::FindObject(const char *name) const
 289 {
 290   /// overloaded from TObject: find object by name
 291   if (fControlObjects) {
 292     return fControlObjects->FindObject(name);
 293   }
 294   return NULL;
 295 }
 296
 297 TObject* AliDxHFECorrelation::FindObject(const TObject *obj) const
 298 {
 299   /// overloaded from TObject: find object by pointer
 300   if (fControlObjects) {
 301     return fControlObjects->FindObject(obj);
 302   }
 303   return NULL;
 304 }
 305
 306 void AliDxHFECorrelation::SaveAs(const char *filename, Option_t */*option*/) const
 307 {
 308   /// overloaded from TObject: save to file
 309   std::auto_ptr<TFile> output(TFile::Open(filename, "RECREATE"));
 310   if (!output.get() || output->IsZombie()) {
 311     AliError(Form("can not open file %s from writing", filename));
 312     return;
 313   }
 314   output->cd();
 315   if (fControlObjects) fControlObjects->Write();
 316   output->Close();
 317 }
 318
 319 AliDxHFECorrelation& AliDxHFECorrelation::operator+=(const AliDxHFECorrelation& other)
 320 {
 321   /// add histograms from another instance
 322   // TODO - need to change this to ThnSparse?
 323   if (!fHistograms || !other.fHistograms) return *this;
 324
 325   for (int i=0; i<kNofHistograms; i++) {
 326     if (fHistograms->At(i)==NULL || other.fHistograms->At(i)==NULL) continue;
 327     TH1* target=reinterpret_cast<TH1*>(fHistograms->At(i));
 328     TH1* source=reinterpret_cast<TH1*>(other.fHistograms->At(i));
 329     if (!target || !source) continue;
 330     TString name(fHistograms->At(i)->GetName());
 331     if (name.CompareTo(target->GetName())!=0) {
 332       AliWarning(Form("skipping incompatible objects at position %d: %s vs %s", i, source->GetName(), target->GetName()));
 333       continue;
 334     }
 335     if (source->IsA()!=target->IsA()) {
 336       AliWarning(Form("skipping incompatible classes at position %d: %s vs %s", i, source->ClassName(), target->ClassName()));
 337       continue;
 338     }
 339     target->Add(source);
 340   }
 341   return *this;
 342 }