PWG3/hfe/AliHFEelecbackground.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 // First implementation of a class
  18 // to reject tagged electron coming from conversion, pi0 and eta
  19 // by calculating the e+e- invariant mass
  20 // of the tagged electron with other tracks
  21 // after looser cuts for the partner.
  22 // PostProcess should extract the background yield
  23 // If running with MC, it can be compared to the expected background yield
  24 //
  25 // Authors:
  26 //   Raphaelle Bailhache <rbailhache@ikf.uni-frankfurt.de > <R.Bailhache@gsi.de >
  27 //
  28 //
  29
  30 #include <THnSparse.h>
  31 #include <TParticle.h>
  32 #include <TFile.h>
  33 #include <TCanvas.h>
  34 #include <TAxis.h>
  35 #include <TH1D.h>
  36 #include <TString.h>
  37 #include <TLegend.h>
  38
  39 #include <AliESDEvent.h>
  40 #include <AliAODEvent.h>
  41 #include <AliESDtrack.h>
  42 #include <AliAODTrack.h>
  43 #include "AliHFEelecbackground.h"
  44 #include <AliMCEvent.h>
  45 #include <AliStack.h>
  46
  47
  48 ClassImp(AliHFEelecbackground)
  49
  50   const Double_t    AliHFEelecbackground::fgkMe = 0.00051099892;
  51
  52 //___________________________________________________________________________________________
  53 AliHFEelecbackground::AliHFEelecbackground():
  54   fESD1(0x0)
  55   ,fAOD1(0x0)
  56   ,fMCEvent(0x0)
  57   ,fBz(0)
  58   ,fkVertex(0x0)
  59   ,fUseMCPID(kFALSE)
  60   ,fPtESD(0.0)
  61   ,fIndexTrack(0)
  62   ,fIndexTrackPart(0)
  63   ,fPdg(0)
  64   ,fLabMother(-1)
  65   ,fIsFrom(-1)
  66   ,fMotherGamma(-1)
  67   ,fMotherPi0(-1)
  68   ,fMotherC(-1)
  69   ,fMotherB(-1)
  70   ,fMotherEta(-1)
  71   ,fIsPartner(kFALSE)
  72   ,fList(0x0)
  73   ,fListPostProcess(0x0)
  74 {
  75   //
  76   // Default constructor
  77   //
  78
  79 }
  80
  81 //_______________________________________________________________________________________________
  82 AliHFEelecbackground::AliHFEelecbackground(const AliHFEelecbackground &p):
  83   TObject(p)
  84   ,fESD1(0x0)
  85   ,fAOD1(0x0)
  86   ,fMCEvent(0x0)
  87   ,fBz(p.fBz)
  88   ,fkVertex(p.fkVertex)
  89   ,fUseMCPID(p.fUseMCPID)
  90   ,fPtESD(p.fPtESD)
  91   ,fIndexTrack(0)
  92   ,fIndexTrackPart(0)
  93   ,fPdg(0)
  94   ,fLabMother(-1)
  95   ,fIsFrom(-1)
  96   ,fMotherGamma(-1)
  97   ,fMotherPi0(-1)
  98   ,fMotherC(-1)
  99   ,fMotherB(-1)
 100   ,fMotherEta(-1)
 101   ,fIsPartner(kFALSE)
 102   ,fList(0x0)
 103   ,fListPostProcess(0x0)
 104 {
 105   //
 106   // Copy constructor
 107   //
 108 }
 109
 110 //_______________________________________________________________________________________________
 111 AliHFEelecbackground&
 112 AliHFEelecbackground::operator=(const AliHFEelecbackground &)
 113 {
 114   //
 115   // Assignment operator
 116   //
 117
 118   AliInfo("Not yet implemented.");
 119   return *this;
 120 }
 121
 122 //_______________________________________________________________________________________________
 123 AliHFEelecbackground::~AliHFEelecbackground()
 124 {
 125   //
 126   // Destructor
 127   //
 128
 129   if(fList){
 130     fList->Clear();
 131     delete fList;
 132   }
 133
 134   if(fListPostProcess){
 135     fListPostProcess->Clear();
 136     delete fListPostProcess;
 137   }
 138 }
 139 //___________________________________________________________________________________________
 140 Bool_t AliHFEelecbackground::Load(const Char_t *filename)
 141 {
 142   //
 143   // Generic container loader
 144   //
 145
 146   if(!TFile::Open(filename)){
 147     return kFALSE;
 148   }
 149   TList *o = 0x0;
 150   if(!(o = (TList*)gFile->Get("Results"))){
 151     return kFALSE;
 152   }
 153   TList *oe = 0x0;
 154   if(!(oe = (TList*)dynamic_cast<TList *>(o->FindObject("HFEelecbackground")))){
 155     return kFALSE;
 156   }
 157   fList = (TList*)oe->Clone("HFEelecbackground");
 158   gFile->Close();
 159   return kTRUE;
 160 }
 161 //_______________________________________________________________________________________________
 162 void AliHFEelecbackground::CreateHistograms(TList* const qaList)
 163 {
 164   //
 165   // create histograms
 166   //
 167   if(!qaList) return;
 168
 169   fList = qaList;
 170   fList->SetName("HFEelecbackground");
 171
 172   // bins
 173
 174   Int_t nBinsPt = 25;
 175   Double_t minPt = 0.001;
 176   Double_t maxPt = 10.0;
 177
 178   Int_t nBinsInv = 50;
 179   Double_t minInv = 0.0;
 180   Double_t maxInv = 0.2;
 181
 182   Int_t nBinsOp = 50;
 183   Double_t minOp = 0.0;
 184   Double_t maxOp = 2;
 185
 186   Double_t *binLimLogPt = new Double_t[nBinsPt+1];
 187   Double_t *binLimPt    = new Double_t[nBinsPt+1];
 188   for(Int_t i=0; i<=nBinsPt; i++) binLimLogPt[i]=(Double_t)TMath::Log10(minPt) + (TMath::Log10(maxPt)-TMath::Log10(minPt))/nBinsPt*(Double_t)i ;
 189   for(Int_t i=0; i<=nBinsPt; i++) binLimPt[i]=(Double_t)TMath::Power(10,binLimLogPt[i]);
 190
 191   Double_t *binLimInv    = new Double_t[nBinsInv+1];
 192   for(Int_t i=0; i<=nBinsInv; i++) binLimInv[i]=(Double_t)minInv  + (maxInv-minInv)  /nBinsInv*(Double_t)i ;
 193
 194   Double_t *binLimOp    = new Double_t[nBinsOp+1];
 195   for(Int_t i=0; i<=nBinsOp; i++) binLimOp[i]=(Double_t)minOp  + (maxOp-minOp) /nBinsOp*(Double_t)i ;
 196
 197
 198   const Int_t nvarO = 3; // ptrectaggede, ptrecmother, openingangle or invmass
 199
 200   Int_t iBinOInv[nvarO];
 201   iBinOInv[0]=nBinsPt;
 202   iBinOInv[1]=nBinsPt;
 203   iBinOInv[2]=nBinsInv;
 204
 205   Int_t iBinOOp[nvarO];
 206   iBinOOp[0]=nBinsPt;
 207   iBinOOp[1]=nBinsPt;
 208   iBinOOp[2]=nBinsOp;
 209
 210   //
 211   //
 212   //
 213
 214   THnSparseF *openinganglepp = new THnSparseF("openinganglepp","",nvarO,iBinOOp);
 215   openinganglepp->SetBinEdges(0,binLimPt);
 216   openinganglepp->SetBinEdges(1,binLimPt);
 217   openinganglepp->SetBinEdges(2,binLimOp);
 218   openinganglepp->Sumw2();
 219
 220   THnSparseF *openinganglenn = new THnSparseF("openinganglenn","",nvarO,iBinOOp);
 221   openinganglenn->SetBinEdges(0,binLimPt);
 222   openinganglenn->SetBinEdges(1,binLimPt);
 223   openinganglenn->SetBinEdges(2,binLimOp);
 224   openinganglenn->Sumw2();
 225
 226   THnSparseF *openingangless = new THnSparseF("openingangless","",nvarO,iBinOOp);
 227   openingangless->SetBinEdges(0,binLimPt);
 228   openingangless->SetBinEdges(1,binLimPt);
 229   openingangless->SetBinEdges(2,binLimOp);
 230   openingangless->Sumw2();
 231
 232   THnSparseF *openingangler = new THnSparseF("openingangler","",nvarO,iBinOOp);
 233   openingangler->SetBinEdges(0,binLimPt);
 234   openingangler->SetBinEdges(1,binLimPt);
 235   openingangler->SetBinEdges(2,binLimOp);
 236   openingangler->Sumw2();
 237
 238   THnSparseF *openingangleos = new THnSparseF("openingangleos","",nvarO,iBinOOp);
 239   openingangleos->SetBinEdges(0,binLimPt);
 240   openingangleos->SetBinEdges(1,binLimPt);
 241   openingangleos->SetBinEdges(2,binLimOp);
 242   openingangleos->Sumw2();
 243
 244   THnSparseF *openinganglepi0=0x0;
 245   THnSparseF *openingangleeta=0x0;
 246   THnSparseF *openinganglegamma=0x0;
 247   THnSparseF *openingangleC=0x0;
 248   THnSparseF *openingangleB=0x0;
 249
 250   if(HasMCData()) {
 251
 252     openinganglepi0 = new THnSparseF("openinganglepi0","",nvarO,iBinOOp);
 253     openinganglepi0->SetBinEdges(0,binLimPt);
 254     openinganglepi0->SetBinEdges(1,binLimPt);
 255     openinganglepi0->SetBinEdges(2,binLimOp);
 256     openinganglepi0->Sumw2();
 257
 258     openingangleeta = new THnSparseF("openingangleeta","",nvarO,iBinOOp);
 259     openingangleeta->SetBinEdges(0,binLimPt);
 260     openingangleeta->SetBinEdges(1,binLimPt);
 261     openingangleeta->SetBinEdges(2,binLimOp);
 262     openingangleeta->Sumw2();
 263
 264     openinganglegamma = new THnSparseF("openinganglegamma","",nvarO,iBinOOp);
 265     openinganglegamma->SetBinEdges(0,binLimPt);
 266     openinganglegamma->SetBinEdges(1,binLimPt);
 267     openinganglegamma->SetBinEdges(2,binLimOp);
 268     openinganglegamma->Sumw2();
 269
 270     openingangleC = new THnSparseF("openingangleC","",nvarO,iBinOOp);
 271     openingangleC->SetBinEdges(0,binLimPt);
 272     openingangleC->SetBinEdges(1,binLimPt);
 273     openingangleC->SetBinEdges(2,binLimOp);
 274     openingangleC->Sumw2();
 275
 276     openingangleB = new THnSparseF("openingangleB","",nvarO,iBinOOp);
 277     openingangleB->SetBinEdges(0,binLimPt);
 278     openingangleB->SetBinEdges(1,binLimPt);
 279     openingangleB->SetBinEdges(2,binLimOp);
 280     openingangleB->Sumw2();
 281
 282   }
 283
 284   //
 285
 286   THnSparseF *invmasspp = new THnSparseF("invmasspp","",nvarO,iBinOInv);
 287   invmasspp->SetBinEdges(0,binLimPt);
 288   invmasspp->SetBinEdges(1,binLimPt);
 289   invmasspp->SetBinEdges(2,binLimInv);
 290   invmasspp->Sumw2();
 291
 292   THnSparseF *invmassnn = new THnSparseF("invmassnn","",nvarO,iBinOInv);
 293   invmassnn->SetBinEdges(0,binLimPt);
 294   invmassnn->SetBinEdges(1,binLimPt);
 295   invmassnn->SetBinEdges(2,binLimInv);
 296   invmassnn->Sumw2();
 297
 298   THnSparseF *invmassss = new THnSparseF("invmassss","",nvarO,iBinOInv);
 299   invmassss->SetBinEdges(0,binLimPt);
 300   invmassss->SetBinEdges(1,binLimPt);
 301   invmassss->SetBinEdges(2,binLimInv);
 302   invmassss->Sumw2();
 303
 304   THnSparseF *invmassr = new THnSparseF("invmassr","",nvarO,iBinOInv);
 305   invmassr->SetBinEdges(0,binLimPt);
 306   invmassr->SetBinEdges(1,binLimPt);
 307   invmassr->SetBinEdges(2,binLimInv);
 308   invmassr->Sumw2();
 309
 310   THnSparseF *invmassos = new THnSparseF("invmassos","",nvarO,iBinOInv);
 311   invmassos->SetBinEdges(0,binLimPt);
 312   invmassos->SetBinEdges(1,binLimPt);
 313   invmassos->SetBinEdges(2,binLimInv);
 314   invmassos->Sumw2();
 315
 316   THnSparseF *invmasspi0=0x0;
 317   THnSparseF *invmasseta=0x0;
 318   THnSparseF *invmassgamma=0x0;
 319   THnSparseF *invmassC=0x0;
 320   THnSparseF *invmassB=0x0;
 321
 322   if(HasMCData()) {
 323
 324     invmasspi0 = new THnSparseF("invmasspi0","",nvarO,iBinOInv);
 325     invmasspi0->SetBinEdges(0,binLimPt);
 326     invmasspi0->SetBinEdges(1,binLimPt);
 327     invmasspi0->SetBinEdges(2,binLimInv);
 328     invmasspi0->Sumw2();
 329
 330     invmasseta = new THnSparseF("invmasseta","",nvarO,iBinOInv);
 331     invmasseta->SetBinEdges(0,binLimPt);
 332     invmasseta->SetBinEdges(1,binLimPt);
 333     invmasseta->SetBinEdges(2,binLimInv);
 334     invmasseta->Sumw2();
 335
 336     invmassgamma = new THnSparseF("invmassgamma","",nvarO,iBinOInv);
 337     invmassgamma->SetBinEdges(0,binLimPt);
 338     invmassgamma->SetBinEdges(1,binLimPt);
 339     invmassgamma->SetBinEdges(2,binLimInv);
 340     invmassgamma->Sumw2();
 341
 342     invmassC = new THnSparseF("invmassC","",nvarO,iBinOInv);
 343     invmassC->SetBinEdges(0,binLimPt);
 344     invmassC->SetBinEdges(1,binLimPt);
 345     invmassC->SetBinEdges(2,binLimInv);
 346     invmassC->Sumw2();
 347
 348     invmassB = new THnSparseF("invmassB","",nvarO,iBinOInv);
 349     invmassB->SetBinEdges(0,binLimPt);
 350     invmassB->SetBinEdges(1,binLimPt);
 351     invmassB->SetBinEdges(2,binLimInv);
 352     invmassB->Sumw2();
 353
 354   }
 355
 356   //
 357   //
 358   //
 359
 360   fList->AddAt(openinganglepp,kPp);
 361   fList->AddAt(openinganglenn,kNn);
 362   fList->AddAt(openingangless,kSs);
 363   fList->AddAt(openingangler,kR);
 364   fList->AddAt(openingangleos,kOs);
 365   if(HasMCData()) {
 366     fList->AddAt(openinganglepi0,2*kNSignComb+kElectronFromPi0);
 367     fList->AddAt(openingangleeta,2*kNSignComb+kElectronFromEta);
 368     fList->AddAt(openinganglegamma,2*kNSignComb+kElectronFromGamma);
 369     fList->AddAt(openingangleC,2*kNSignComb+kElectronFromC);
 370     fList->AddAt(openingangleB,2*kNSignComb+kElectronFromB);
 371   }
 372
 373   fList->AddAt(invmasspp,kNSignComb+kPp);
 374   fList->AddAt(invmassnn,kNSignComb+kNn);
 375   fList->AddAt(invmassss,kNSignComb+kSs);
 376   fList->AddAt(invmassr,kNSignComb+kR);
 377   fList->AddAt(invmassos,kNSignComb+kOs);
 378   if(HasMCData()) {
 379     fList->AddAt(invmasspi0,2*kNSignComb+kNMCInfo+kElectronFromPi0);
 380     fList->AddAt(invmasseta,2*kNSignComb+kNMCInfo+kElectronFromEta);
 381     fList->AddAt(invmassgamma,2*kNSignComb+kNMCInfo+kElectronFromGamma);
 382     fList->AddAt(invmassC,2*kNSignComb+kNMCInfo+kElectronFromC);
 383     fList->AddAt(invmassB,2*kNSignComb+kNMCInfo+kElectronFromB);
 384   }
 385
 386 }
 387 //_______________________________________________________________________________________________
 388 void AliHFEelecbackground::PairAnalysis(AliESDtrack* const track, AliESDtrack* const trackPart)
 389 {
 390   //
 391   // calculate (tagged e-partner) dca, opening angle, invariant mass
 392   //
 393
 394   ////////////////////////
 395   // Partner track cut
 396   ////////////////////////
 397   if(!SingleTrackCut(trackPart)) return;
 398
 399   ////////////////////////
 400   // Take label
 401   ////////////////////////
 402   Int_t indexTrack = TMath::Abs(track->GetLabel());
 403   Int_t indexTrackPart = TMath::Abs(trackPart->GetLabel());
 404
 405   /////////////////////////
 406   // If MC data
 407   ////////////////////////
 408
 409   if(HasMCData()) {
 410
 411     // Take info track if not already done
 412     if(indexTrack!= fIndexTrack) {
 413       fIsFrom = -1;
 414
 415       fPdg = GetPdg(indexTrack);
 416       fLabMother = GetLabMother(indexTrack);
 417
 418       fMotherGamma = IsMotherGamma(indexTrack);
 419       if((fMotherGamma != -1) && ((TMath::Abs(fPdg)) == 11)) fIsFrom = kElectronFromGamma;
 420       fMotherPi0 = IsMotherPi0(indexTrack);
 421       if((fMotherPi0 != -1) && ((TMath::Abs(fPdg)) == 11)) fIsFrom = kElectronFromPi0;
 422       fMotherC = IsMotherC(indexTrack);
 423       if((fMotherC != -1) && ((TMath::Abs(fPdg)) == 11)) fIsFrom = kElectronFromC;
 424       fMotherB = IsMotherB(indexTrack);
 425       if((fMotherB != -1) && ((TMath::Abs(fPdg)) == 11)) fIsFrom = kElectronFromB;
 426       fMotherEta = IsMotherEta(indexTrack);
 427       if((fMotherEta != -1) && ((TMath::Abs(fPdg)) == 11)) fIsFrom = kElectronFromEta;
 428
 429     }
 430
 431     // Look at trackPart
 432     Int_t pdgPart = GetPdg(indexTrackPart);
 433     if(TMath::Abs(pdgPart) == 11) {
 434       Int_t labMotherPart = GetLabMother(indexTrackPart);
 435       if((labMotherPart == fLabMother) && (indexTrack != indexTrackPart) && (TMath::Abs(pdgPart)==11) && (fPdg*pdgPart < 0) && (fLabMother >=0) && (fLabMother < (((AliStack *)fMCEvent->Stack())->GetNtrack()))) fIsPartner = kTRUE;
 436     }
 437
 438   }
 439
 440   //////////////////////
 441   // Sign
 442   /////////////////////
 443   Int_t sign = -1;
 444   if((track->Charge() > 0.0) && (trackPart->Charge() > 0.0)) sign = kPp;
 445   if((track->Charge() < 0.0) && (trackPart->Charge() < 0.0)) sign = kNn;
 446   if(((track->Charge() > 0.0) && (trackPart->Charge() < 0.0)) || ((track->Charge() < 0.0) && (trackPart->Charge() > 0.0))) sign = kOs;
 447
 448   //////////////////////
 449   // DCA
 450   /////////////////////
 451
 452   Double_t xthis,xp;
 453   Double_t dca = track->GetDCA(trackPart,fBz,xthis,xp);
 454   if(TMath::Abs(dca) > 3.0) return;
 455
 456   /////////////////////////////
 457   // Propagate
 458   ////////////////////////////
 459
 460   Double_t norradius = TMath::Sqrt(fkVertex->GetX()*fkVertex->GetX()+fkVertex->GetY()*fkVertex->GetY());
 461
 462   AliESDtrack *trackCopy = new AliESDtrack(*track);
 463   AliESDtrack *trackPartCopy = new AliESDtrack(*trackPart);
 464   Bool_t propagateok = kTRUE;
 465   if((!(trackPartCopy->PropagateTo(norradius,fBz))) || (!(trackCopy->PropagateTo(norradius,fBz)))) propagateok = kFALSE;
 466   if(!propagateok) {
 467     if(trackCopy) delete trackCopy;
 468     if(trackPartCopy) delete trackPartCopy;
 469     return;
 470   }
 471
 472   ///////////////////////////////////
 473   // Calcul mother variables
 474   ///////////////////////////////////
 475   Double_t results[5];
 476   Double_t resultsr[5];
 477
 478   CalculateMotherVariable(trackCopy,trackPartCopy,&results[0]);
 479   CalculateMotherVariableR(trackCopy,trackPartCopy,&resultsr[0]);
 480
 481   /////////////////////////////////////
 482   // Fill
 483   /////////////////////////////////////
 484
 485   FillOutput(results, resultsr, sign);
 486
 487   if(trackCopy) delete trackCopy;
 488   if(trackPartCopy) delete trackPartCopy;
 489
 490 }
 491 //_____________________________________________________________________________________
 492 void AliHFEelecbackground::CalculateMotherVariable(AliESDtrack* const track, AliESDtrack* const trackpart, Double_t *results)
 493 {
 494   //
 495   // variables mother and take the pt of the track
 496   //
 497   // results contain: ptmother, invmass, etamother, phimother, openingangle
 498   //
 499
 500   TVector3 v3Dtagged;
 501   TVector3 v3Dpart;
 502
 503   Double_t *pxyz = new Double_t[3];
 504   track->PxPyPz(pxyz);
 505   v3Dtagged.SetXYZ(pxyz[0],pxyz[1],pxyz[2]);
 506   fPtESD = TMath::Sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]);
 507
 508   Double_t *pxyzpart = new Double_t[3];
 509   trackpart->PxPyPz(pxyzpart);
 510   v3Dpart.SetXYZ(pxyzpart[0],pxyzpart[1],pxyzpart[2]);
 511
 512
 513
 514   TVector3 motherrec = v3Dtagged + v3Dpart;
 515
 516   Double_t etaESDmother = motherrec.Eta();
 517   Double_t ptESDmother  = motherrec.Pt();
 518   Double_t phiESDmother = motherrec.Phi();
 519   if(phiESDmother > TMath::Pi()) phiESDmother = phiESDmother - (2*TMath::Pi());
 520
 521
 522   // openinganglepropagated
 523   Double_t openingangle = v3Dtagged.Angle(v3Dpart);
 524
 525   // invmass
 526   Double_t pESD      = TMath::Sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]+pxyz[2]*pxyz[2]);
 527   Double_t pESDpart  = TMath::Sqrt(pxyzpart[0]*pxyzpart[0]+pxyzpart[1]*pxyzpart[1]+pxyzpart[2]*pxyzpart[2]);
 528
 529   // e propagate
 530   Double_t eESD     = TMath::Sqrt(pESD*pESD+fgkMe*fgkMe);
 531   Double_t eESDpart = TMath::Sqrt(pESDpart*pESDpart+fgkMe*fgkMe);
 532
 533   Double_t invmass = TMath::Sqrt((eESD+eESDpart)*(eESD+eESDpart)-(motherrec.Px()*motherrec.Px()+motherrec.Py()*motherrec.Py()+motherrec.Pz()*motherrec.Pz()));
 534
 535   if(!results) {
 536     results = new Double_t[5];
 537   }
 538   results[0] = ptESDmother;
 539   results[1] = etaESDmother;
 540   results[2] = phiESDmother;
 541   results[3] = invmass;
 542   results[4] = openingangle;
 543
 544 }
 545 //_____________________________________________________________________________________
 546 void AliHFEelecbackground::CalculateMotherVariableR(AliESDtrack* const track, AliESDtrack* const trackpart, Double_t *results)
 547 {
 548   //
 549   // variables mother
 550   //
 551   // results contain: ptmother, invmass, etamother, phimother, openingangle
 552   //
 553
 554   TVector3 v3Dtagged;
 555   TVector3 v3Dpart;
 556
 557   Double_t *pxyz = new Double_t[3];
 558   track->PxPyPz(pxyz);
 559   v3Dtagged.SetXYZ(pxyz[0],pxyz[1],pxyz[2]);
 560   Double_t *pxyzpart = new Double_t[3];
 561   trackpart->PxPyPz(pxyzpart);
 562   v3Dpart.SetXYZ(pxyzpart[0],pxyzpart[1],pxyzpart[2]);
 563
 564   // rotate the partner
 565   v3Dpart.RotateZ(TMath::Pi());
 566   v3Dpart.GetXYZ(pxyzpart);
 567
 568
 569   TVector3 motherrec = v3Dtagged + v3Dpart;
 570
 571   Double_t etaESDmother = motherrec.Eta();
 572   Double_t ptESDmother  = motherrec.Pt();
 573   Double_t phiESDmother = motherrec.Phi();
 574   if(phiESDmother > TMath::Pi()) phiESDmother = phiESDmother - (2*TMath::Pi());
 575
 576
 577   // openinganglepropagated
 578   Double_t openingangle = v3Dtagged.Angle(v3Dpart);
 579   //printf("Openingangle %f\n",openingangle);
 580
 581   // invmass
 582   Double_t pESD      = TMath::Sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]+pxyz[2]*pxyz[2]);
 583   Double_t pESDpart  = TMath::Sqrt(pxyzpart[0]*pxyzpart[0]+pxyzpart[1]*pxyzpart[1]+pxyzpart[2]*pxyzpart[2]);
 584   // e propagate
 585   Double_t eESD     = TMath::Sqrt(pESD*pESD+fgkMe*fgkMe);
 586   Double_t eESDpart = TMath::Sqrt(pESDpart*pESDpart+fgkMe*fgkMe);
 587
 588   Double_t invmass = TMath::Sqrt((eESD+eESDpart)*(eESD+eESDpart)-(motherrec.Px()*motherrec.Px()+motherrec.Py()*motherrec.Py()+motherrec.Pz()*motherrec.Pz()));
 589
 590   if(!results) {
 591     results = new Double_t[5];
 592   }
 593   results[0] = ptESDmother;
 594   results[1] = etaESDmother;
 595   results[2] = phiESDmother;
 596   results[3] = invmass;
 597   results[4] = openingangle;
 598
 599
 600 }
 601 //_________________________________________________________________________________
 602 void AliHFEelecbackground::FillOutput(Double_t *results, Double_t *resultsr, Int_t sign)
 603 {
 604   //
 605   // Fill the invariant mass and opening angle distributions
 606   //
 607
 608   Double_t co[3];
 609   co[0] = fPtESD;
 610
 611   switch(sign){
 612
 613       case kPp:
 614         co[1] = results[0];
 615         co[2] = TMath::Abs(results[4]);
 616         (dynamic_cast<THnSparseF *>(fList->At(kPp)))->Fill(co);
 617         (dynamic_cast<THnSparseF *>(fList->At(kSs)))->Fill(co);
 618
 619
 620         co[2] = results[3];
 621         if(TMath::Abs(results[4]) < 0.8){
 622           (dynamic_cast<THnSparseF *>(fList->At(kPp+kNSignComb)))->Fill(co);
 623           (dynamic_cast<THnSparseF *>(fList->At(kSs+kNSignComb)))->Fill(co);
 624         }
 625
 626         break;
 627
 628       case kNn:
 629         co[1] = results[0];
 630         co[2] = TMath::Abs(results[4]);
 631         (dynamic_cast<THnSparseF *>(fList->At(kNn)))->Fill(co);
 632         (dynamic_cast<THnSparseF *>(fList->At(kSs)))->Fill(co);
 633
 634         co[2] = results[3];
 635         if(TMath::Abs(results[4]) < 0.8){
 636           (dynamic_cast<THnSparseF *>(fList->At(kNn+kNSignComb)))->Fill(co);
 637           (dynamic_cast<THnSparseF *>(fList->At(kSs+kNSignComb)))->Fill(co);
 638         }
 639         break;
 640
 641       case kOs:
 642         co[1] = results[0];
 643         co[2] = TMath::Abs(results[4]);
 644         (dynamic_cast<THnSparseF *>(fList->At(kOs)))->Fill(co);
 645         if(HasMCData()) {
 646           if((fIsFrom == kElectronFromPi0) && (fIsPartner)) (dynamic_cast<THnSparseF *>(fList->At(2*kNSignComb+kElectronFromPi0)))->Fill(co);
 647           if((fIsFrom == kElectronFromEta) && (fIsPartner)) (dynamic_cast<THnSparseF *>(fList->At(2*kNSignComb+kElectronFromEta)))->Fill(co);
 648           if((fIsFrom == kElectronFromGamma) && (fIsPartner)) (dynamic_cast<THnSparseF *>(fList->At(2*kNSignComb+kElectronFromGamma)))->Fill(co);
 649           if((fIsFrom == kElectronFromC) && (fIsPartner)) (dynamic_cast<THnSparseF *>(fList->At(2*kNSignComb+kElectronFromC)))->Fill(co);
 650           if((fIsFrom == kElectronFromB) && (fIsPartner)) (dynamic_cast<THnSparseF *>(fList->At(2*kNSignComb+kElectronFromB)))->Fill(co);
 651         }
 652
 653         co[2] = results[3];
 654         if(TMath::Abs(results[4]) < 0.8){
 655           (dynamic_cast<THnSparseF *>(fList->At(kOs+kNSignComb)))->Fill(co);
 656           if(HasMCData()) {
 657             if((fIsFrom == kElectronFromPi0) && (fIsPartner)) (dynamic_cast<THnSparseF *>(fList->At(2*kNSignComb+kElectronFromPi0+kNMCInfo)))->Fill(co);
 658             if((fIsFrom == kElectronFromEta) && (fIsPartner)) (dynamic_cast<THnSparseF *>(fList->At(2*kNSignComb+kElectronFromEta+kNMCInfo)))->Fill(co);
 659             if((fIsFrom == kElectronFromGamma) && (fIsPartner)) (dynamic_cast<THnSparseF *>(fList->At(2*kNSignComb+kElectronFromGamma+kNMCInfo)))->Fill(co);
 660             if((fIsFrom == kElectronFromC) && (fIsPartner)) (dynamic_cast<THnSparseF *>(fList->At(2*kNSignComb+kElectronFromC+kNMCInfo)))->Fill(co);
 661             if((fIsFrom == kElectronFromB) && (fIsPartner)) (dynamic_cast<THnSparseF *>(fList->At(2*kNSignComb+kElectronFromB+kNMCInfo)))->Fill(co);
 662           }
 663         }
 664
 665         // rotated
 666         co[1] = resultsr[0];
 667         co[2] = TMath::Abs(resultsr[4]);
 668
 669         (dynamic_cast<THnSparseF *>(fList->At(kR)))->Fill(co);
 670
 671         co[2] = resultsr[3];
 672         if(TMath::Abs(resultsr[4]) < 0.8){
 673           (dynamic_cast<THnSparseF *>(fList->At(kR+kNSignComb)))->Fill(co);
 674         }
 675         break;
 676
 677   default:
 678
 679     break;
 680
 681   }
 682 }
 683 //_______________________________________________________________________________________________
 684 Bool_t AliHFEelecbackground::SingleTrackCut(AliESDtrack* const track) const
 685 {
 686   //
 687   // Return minimum quality for the partner
 688   //
 689
 690   UInt_t status = track->GetStatus();
 691
 692   if(((status & AliESDtrack::kTPCin)==0) && (status & AliESDtrack::kITSin)) {
 693
 694     Int_t nbcl = track->GetITSclusters(0);
 695     if(nbcl > 1) return kTRUE;
 696     else return kFALSE;
 697
 698   }
 699
 700   if(status & AliESDtrack::kTPCin) {
 701
 702     if(status & AliESDtrack::kTPCrefit)  return kTRUE;
 703     else return kFALSE;
 704
 705   }
 706
 707   return kFALSE;
 708
 709 }
 710 //______________________________________________________________________________________________
 711 void AliHFEelecbackground::SetEvent(AliESDEvent* const ESD)
 712 {
 713   //
 714   // Set the AliESD Event, the magnetic field and the primary vertex
 715   //
 716
 717   fESD1=ESD;
 718   fBz=fESD1->GetMagneticField();
 719   fkVertex = fESD1->GetPrimaryVertex();
 720
 721 }
 722 //____________________________________________________________________________________________________________
 723 Int_t AliHFEelecbackground::IsMotherGamma(Int_t tr) {
 724
 725   //
 726   // Return the lab of gamma mother or -1 if not gamma
 727   //
 728
 729   AliStack* stack = fMCEvent->Stack();
 730   if((tr < 0) || (tr >= stack->GetNtrack())) return -1;
 731
 732   // Take mother
 733   TParticle * particle = stack->Particle(tr);
 734   if(!particle) return -1;
 735   Int_t imother   = particle->GetFirstMother();
 736   if((imother < 0) || (imother >= stack->GetNtrack())) return -1;
 737   TParticle * mother = stack->Particle(imother);
 738   if(!mother) return -1;
 739
 740   // Check gamma
 741   Int_t pdg = mother->GetPdgCode();
 742   if(TMath::Abs(pdg) == 22) return imother;
 743   if(TMath::Abs(pdg) == 11) {
 744     return IsMotherGamma(imother);
 745   }
 746   return -1;
 747
 748 }
 749 //
 750 //____________________________________________________________________________________________________________
 751 Int_t AliHFEelecbackground::IsMotherPi0(Int_t tr) {
 752
 753   //
 754   // Return the lab of pi0 mother or -1 if not pi0
 755   //
 756
 757   AliStack* stack = fMCEvent->Stack();
 758   if((tr < 0) || (tr >= stack->GetNtrack())) return -1;
 759
 760   // Take mother
 761   TParticle * particle = stack->Particle(tr);
 762   if(!particle) return -1;
 763   Int_t imother   = particle->GetFirstMother();
 764   if((imother < 0) || (imother >= stack->GetNtrack())) return -1;
 765   TParticle * mother = stack->Particle(imother);
 766   if(!mother) return -1;
 767
 768   // Check gamma
 769   Int_t pdg = mother->GetPdgCode();
 770   if(TMath::Abs(pdg) == 111) return imother;
 771   if(TMath::Abs(pdg) == 11) {
 772     return IsMotherPi0(imother);
 773   }
 774   return -1;
 775
 776 }
 777 //____________________________________________________________________________________________________________
 778 Int_t AliHFEelecbackground::IsMotherEta(Int_t tr) {
 779
 780   //
 781   // Return the lab of pi0 mother or -1 if not pi0
 782   //
 783
 784   AliStack* stack = fMCEvent->Stack();
 785   if((tr < 0) || (tr >= stack->GetNtrack())) return -1;
 786
 787   // Take mother
 788   TParticle * particle = stack->Particle(tr);
 789   if(!particle) return -1;
 790   Int_t imother   = particle->GetFirstMother();
 791   if((imother < 0) || (imother >= stack->GetNtrack())) return -1;
 792   TParticle * mother = stack->Particle(imother);
 793   if(!mother) return -1;
 794
 795   // Check gamma
 796   Int_t pdg = mother->GetPdgCode();
 797   if(TMath::Abs(pdg) == 221) return imother;
 798   if(TMath::Abs(pdg) == 11) {
 799     return IsMotherEta(imother);
 800   }
 801   return -1;
 802
 803 }
 804 //____________________________________________________________________________________________________________
 805 Int_t AliHFEelecbackground::IsMotherC(Int_t tr) {
 806
 807   //
 808   // Return the lab of signal mother or -1 if not signal
 809   //
 810
 811   AliStack* stack = fMCEvent->Stack();
 812   if((tr < 0) || (tr >= stack->GetNtrack())) return -1;
 813
 814   // Take mother
 815   TParticle * particle = stack->Particle(tr);
 816   if(!particle) return -1;
 817   Int_t imother   = particle->GetFirstMother();
 818   if((imother < 0) || (imother >= stack->GetNtrack())) return -1;
 819   TParticle * mother = stack->Particle(imother);
 820   if(!mother) return -1;
 821
 822   // Check gamma
 823   Int_t pdg = mother->GetPdgCode();
 824   if((TMath::Abs(pdg)==411) || (TMath::Abs(pdg)==421) || (TMath::Abs(pdg)==431) || (TMath::Abs(pdg)==4122) || (TMath::Abs(pdg)==4132) || (TMath::Abs(pdg)==4232) || (TMath::Abs(pdg)==43320)) return imother;
 825   if(TMath::Abs(pdg) == 11) {
 826     return IsMotherC(imother);
 827   }
 828   return -1;
 829
 830 }
 831 //____________________________________________________________________________________________________________
 832 Int_t AliHFEelecbackground::IsMotherB(Int_t tr) {
 833
 834   //
 835   // Return the lab of signal mother or -1 if not signal
 836   //
 837
 838   AliStack* stack = fMCEvent->Stack();
 839   if((tr < 0) || (tr >= stack->GetNtrack())) return -1;
 840
 841   // Take mother
 842   TParticle * particle = stack->Particle(tr);
 843   if(!particle) return -1;
 844   Int_t imother   = particle->GetFirstMother();
 845   if((imother < 0) || (imother >= stack->GetNtrack())) return -1;
 846   TParticle * mother = stack->Particle(imother);
 847   if(!mother) return -1;
 848
 849   // Check gamma
 850   Int_t pdg = mother->GetPdgCode();
 851   if((TMath::Abs(pdg)==511) || (TMath::Abs(pdg)==521) || (TMath::Abs(pdg)==531) || (TMath::Abs(pdg)==5122) || (TMath::Abs(pdg)==5132) || (TMath::Abs(pdg)==5232) || (TMath::Abs(pdg)==53320)) return imother;
 852   if(TMath::Abs(pdg) == 11) {
 853     return IsMotherB(imother);
 854   }
 855   return -1;
 856
 857 }
 858 //____________________________________________________________________________________________________________
 859 Int_t AliHFEelecbackground::GetPdg(Int_t tr) {
 860
 861   //
 862   // Simply pdg code
 863   //
 864
 865   AliStack* stack = fMCEvent->Stack();
 866   if((tr < 0) || (tr >= stack->GetNtrack())) return -1;
 867
 868   // MC Information
 869   TParticle * particle = stack->Particle(tr);
 870   if(!particle) return -1;
 871   Int_t pdg = particle->GetPdgCode();
 872
 873   return pdg;
 874
 875 }
 876 //____________________________________________________________________________________________________________
 877 Int_t AliHFEelecbackground::GetLabMother(Int_t tr) {
 878
 879   //
 880   // Simply lab mother
 881   //
 882
 883   AliStack* stack = fMCEvent->Stack();
 884   if((tr < 0) || (tr >= stack->GetNtrack())) return -1;
 885
 886   // MC Information
 887   TParticle * particle = stack->Particle(tr);
 888   if(!particle) return -1;
 889   Int_t imother = particle->GetFirstMother();
 890
 891   return imother;
 892
 893 }
 894 //_______________________________________________________________________________________________
 895 void AliHFEelecbackground::PostProcess()
 896 {
 897   //
 898   // Post process the histos and extract the background pt spectra
 899   //
 900
 901   if(!fList) return;
 902
 903   // invariant mass input spectra
 904   THnSparseF *invmassss = dynamic_cast<THnSparseF *>(fList->FindObject("invmassss"));
 905   THnSparseF *invmassr  = dynamic_cast<THnSparseF *>(fList->FindObject("invmassr"));
 906   THnSparseF *invmassos = dynamic_cast<THnSparseF *>(fList->FindObject("invmassos"));
 907   THnSparseF *invmassgamma = dynamic_cast<THnSparseF *>(fList->FindObject("invmassgamma"));
 908   THnSparseF *invmasspi0 = dynamic_cast<THnSparseF *>(fList->FindObject("invmasspi0"));
 909   THnSparseF *invmasseta = dynamic_cast<THnSparseF *>(fList->FindObject("invmasseta"));
 910   THnSparseF *invmassC = dynamic_cast<THnSparseF *>(fList->FindObject("invmassC"));
 911   THnSparseF *invmassB = dynamic_cast<THnSparseF *>(fList->FindObject("invmassB"));
 912
 913   TAxis *ptaxisinvmass = invmassss->GetAxis(0);
 914   Int_t  nbinsptinvmass = ptaxisinvmass->GetNbins();
 915
 916   // outputs
 917   TH1D **invmassosptproj = new TH1D*[nbinsptinvmass];
 918   TH1D **invmassssptproj = new TH1D*[nbinsptinvmass];
 919   TH1D **invmassrptproj = new TH1D*[nbinsptinvmass];
 920   TH1D **invmassdiffptproj = new TH1D*[nbinsptinvmass];
 921   TH1D **invmassgammaptproj = new TH1D*[nbinsptinvmass];
 922   TH1D **invmasspi0ptproj = new TH1D*[nbinsptinvmass];
 923   TH1D **invmassetaptproj = new TH1D*[nbinsptinvmass];
 924   TH1D **invmassCptproj = new TH1D*[nbinsptinvmass];
 925   TH1D **invmassBptproj = new TH1D*[nbinsptinvmass];
 926
 927   TH1D *yieldPtFound = (TH1D *) invmassss->Projection(0);
 928   yieldPtFound->SetName("Found yield");
 929   yieldPtFound->Reset();
 930
 931   TH1D *yieldPtSourcesMC = 0x0;
 932   if(invmasspi0 && invmasseta && invmassgamma) {
 933     yieldPtSourcesMC = (TH1D *) invmassss->Projection(0);
 934     yieldPtSourcesMC->SetName("Found yield");
 935     yieldPtSourcesMC->Reset();
 936   }
 937
 938   TH1D *yieldPtSignalCutMC = 0x0;
 939   if(invmassC && invmassB) {
 940     yieldPtSignalCutMC = (TH1D *) invmassss->Projection(0);
 941     yieldPtSignalCutMC->SetName("Found yield");
 942     yieldPtSignalCutMC->Reset();
 943   }
 944
 945   // canvas
 946   Int_t nbrow = (Int_t) (nbinsptinvmass/5);
 947   TString namecanvas("Invmassnamecanvas");
 948   TCanvas * canvas =new TCanvas(namecanvas,namecanvas,800,800);
 949   canvas->Divide(5,nbrow+1);
 950
 951
 952   // loop on pt bins
 953   for(Int_t k=1; k <= nbinsptinvmass; k++){
 954
 955     Double_t lowedge = ptaxisinvmass->GetBinLowEdge(k);
 956     Double_t upedge  = ptaxisinvmass->GetBinUpEdge(k);
 957
 958     ((TAxis *)invmassss->GetAxis(0))->SetRange(k,k);
 959     ((TAxis *)invmassr->GetAxis(0))->SetRange(k,k);
 960     ((TAxis *)invmassos->GetAxis(0))->SetRange(k,k);
 961
 962     invmassosptproj[k-1] = invmassos->Projection(2);
 963     invmassssptproj[k-1] = invmassss->Projection(2);
 964     invmassrptproj[k-1]  = invmassr->Projection(2);
 965     invmassgammaptproj[k-1] = 0x0;
 966     invmasspi0ptproj[k-1] = 0x0;
 967     invmassetaptproj[k-1] = 0x0;
 968     invmassCptproj[k-1] = 0x0;
 969     invmassBptproj[k-1] = 0x0;
 970     if(invmassgamma) invmassgammaptproj[k-1] = invmassgamma->Projection(2);
 971     if(invmasspi0) invmasspi0ptproj[k-1] = invmasspi0->Projection(2);
 972     if(invmasseta) invmassetaptproj[k-1] = invmasseta->Projection(2);
 973     if(invmassC) invmassCptproj[k-1] = invmassC->Projection(2);
 974     if(invmassB) invmassBptproj[k-1] = invmassB->Projection(2);
 975
 976     invmassdiffptproj[k-1] = (TH1D *) invmassosptproj[k-1]->Clone();
 977     TString name("Invmassdiffptbin");
 978     name += k;
 979     invmassdiffptproj[k-1]->SetName(name);
 980     invmassdiffptproj[k-1]->Add(invmassssptproj[k-1],-1.0);
 981
 982     TString namee("p_{T} tagged from ");
 983     namee += lowedge;
 984     namee += " GeV/c to ";
 985     namee += upedge;
 986     namee += " GeV/c";
 987
 988     invmassosptproj[k-1]->SetTitle((const char*)namee);
 989     invmassssptproj[k-1]->SetTitle((const char*)namee);
 990     invmassrptproj[k-1]->SetTitle((const char*)namee);
 991     invmassdiffptproj[k-1]->SetTitle((const char*)namee);
 992     if(invmassgammaptproj[k-1]) invmassgammaptproj[k-1]->SetTitle((const char*)namee);
 993     if(invmasspi0ptproj[k-1]) invmasspi0ptproj[k-1]->SetTitle((const char*)namee);
 994     if(invmassetaptproj[k-1]) invmassetaptproj[k-1]->SetTitle((const char*)namee);
 995     if(invmassCptproj[k-1]) invmassCptproj[k-1]->SetTitle((const char*)namee);
 996     if(invmassBptproj[k-1]) invmassBptproj[k-1]->SetTitle((const char*)namee);
 997
 998
 999
1000     invmassosptproj[k-1]->SetStats(0);
1001     invmassssptproj[k-1]->SetStats(0);
1002     invmassrptproj[k-1]->SetStats(0);
1003     invmassdiffptproj[k-1]->SetStats(0);
1004     if(invmassgammaptproj[k-1]) invmassgammaptproj[k-1]->SetStats(0);
1005     if(invmasspi0ptproj[k-1]) invmasspi0ptproj[k-1]->SetStats(0);
1006     if(invmassetaptproj[k-1]) invmassetaptproj[k-1]->SetStats(0);
1007     if(invmassCptproj[k-1]) invmassCptproj[k-1]->SetStats(0);
1008     if(invmassBptproj[k-1]) invmassBptproj[k-1]->SetStats(0);
1009
1010     Double_t yieldf = invmassdiffptproj[k-1]->Integral();
1011     if(invmassetaptproj[k-1] && invmasspi0ptproj[k-1] && invmassgammaptproj[k-1] && invmassCptproj[k-1] && invmassBptproj[k-1]) {
1012       Double_t yieldg = invmassetaptproj[k-1]->Integral() + invmasspi0ptproj[k-1]->Integral() + invmassgammaptproj[k-1]->Integral();
1013       yieldPtSourcesMC->SetBinContent(k,yieldg);
1014
1015       Double_t yieldsignal = invmassCptproj[k-1]->Integral() + invmassBptproj[k-1]->Integral();
1016       yieldPtSignalCutMC->SetBinContent(k,yieldsignal);
1017     }
1018
1019     yieldPtFound->SetBinContent(k,yieldf);
1020
1021     canvas->cd(k);
1022     invmassosptproj[k-1]->Draw();
1023     invmassssptproj[k-1]->Draw("same");
1024     invmassdiffptproj[k-1]->Draw("same");
1025     invmassrptproj[k-1]->Draw("same");
1026     TLegend *legiv = new TLegend(0.4,0.6,0.89,0.89);
1027     legiv->AddEntry(invmassosptproj[k-1],"Opposite signs","p");
1028     legiv->AddEntry(invmassssptproj[k-1],"Same signs","p");
1029     legiv->AddEntry(invmassdiffptproj[k-1],"(Opposite - Same) signs","p");
1030     legiv->AddEntry(invmassrptproj[k-1],"rotated","p");
1031     if(invmassgammaptproj[k-1]) legiv->AddEntry(invmassgammaptproj[k-1],"e^{+}e^{-} from #gamma","p");
1032     if(invmasspi0ptproj[k-1]) legiv->AddEntry(invmasspi0ptproj[k-1],"e^{+}e^{-} from #pi^{0}","p");
1033     if(invmassetaptproj[k-1]) legiv->AddEntry(invmassetaptproj[k-1],"e^{+}e^{-} from #eta","p");
1034     legiv->Draw("same");
1035
1036   }
1037
1038   yieldPtFound->SetStats(0);
1039   if(yieldPtSourcesMC) yieldPtSourcesMC->SetStats(0);
1040   if(yieldPtSignalCutMC) yieldPtSignalCutMC->SetStats(0);
1041
1042   TCanvas * canvasfin =new TCanvas("results","results",800,800);
1043   canvasfin->cd(1);
1044   yieldPtFound->Draw();
1045   if(yieldPtSourcesMC && yieldPtSignalCutMC) {
1046     yieldPtSourcesMC->Draw("same");
1047     yieldPtSignalCutMC->Draw("same");
1048     TLegend *lega = new TLegend(0.4,0.6,0.89,0.89);
1049     lega->AddEntry(yieldPtFound,"Contributions found","p");
1050     lega->AddEntry(yieldPtSourcesMC,"Contributions of e^{+}e^{-} from #gamma, #pi^{0} and #eta","p");
1051     lega->AddEntry(yieldPtSignalCutMC,"Contributions of e^{+}e^{-} from C and B","p");
1052     lega->Draw("same");
1053   }
1054
1055
1056
1057   if(!fListPostProcess) fListPostProcess = new TList();
1058   fListPostProcess->SetName("ListPostProcess");
1059
1060   for(Int_t k=0; k < nbinsptinvmass; k++){
1061     fListPostProcess->AddAt(invmassosptproj[k],kOos+kNOutput*k);
1062     fListPostProcess->AddAt(invmassssptproj[k],kOss+kNOutput*k);
1063     fListPostProcess->AddAt(invmassrptproj[k],kOr+kNOutput*k);
1064     fListPostProcess->AddAt(invmassdiffptproj[k],kOdiff+kNOutput*k);
1065     if(invmassgammaptproj[k]) fListPostProcess->AddAt(invmassgammaptproj[k],kNOutput*nbinsptinvmass+kNMCInfo*k+kElectronFromGamma);
1066     if(invmasspi0ptproj[k]) fListPostProcess->AddAt(invmasspi0ptproj[k],kNOutput*nbinsptinvmass+kNMCInfo*k+kElectronFromPi0);
1067     if(invmassetaptproj[k]) fListPostProcess->AddAt(invmassetaptproj[k],kNOutput*nbinsptinvmass+kNMCInfo*k+kElectronFromEta);
1068     if(invmassCptproj[k]) fListPostProcess->AddAt(invmassCptproj[k],kNOutput*nbinsptinvmass+kNMCInfo*k+kElectronFromC);
1069     if(invmassBptproj[k]) fListPostProcess->AddAt(invmassBptproj[k],kNOutput*nbinsptinvmass+kNMCInfo*k+kElectronFromB);
1070   }
1071
1072   fListPostProcess->AddAt(yieldPtFound,kNOutput*nbinsptinvmass+kNMCInfo*nbinsptinvmass);
1073   if(yieldPtSourcesMC) fListPostProcess->AddAt(yieldPtSourcesMC,kNOutput*nbinsptinvmass+kNMCInfo*nbinsptinvmass+1);
1074   if(yieldPtSignalCutMC) fListPostProcess->AddAt(yieldPtSignalCutMC,kNOutput*nbinsptinvmass+kNMCInfo*nbinsptinvmass+2);
1075
1076   // delete dynamic array
1077   delete[] invmassosptproj;
1078   delete[] invmassssptproj;
1079   delete[] invmassrptproj;
1080   delete[] invmassdiffptproj;
1081   delete[] invmassgammaptproj;
1082   delete[] invmasspi0ptproj;
1083   delete[] invmassetaptproj;
1084   delete[] invmassCptproj;
1085   delete[] invmassBptproj;
1086
1087 }
1088 //_______________________________________________________________________________________________
1089 void AliHFEelecbackground::Plot() const
1090 {
1091   //
1092   // Plot the output
1093   //
1094
1095   if(!fList) return;
1096
1097   // opening angle
1098   THnSparseF *openinganglepp = dynamic_cast<THnSparseF *>(fList->FindObject("openinganglepp"));
1099   THnSparseF *openinganglenn = dynamic_cast<THnSparseF *>(fList->FindObject("openinganglenn"));
1100   THnSparseF *openingangless = dynamic_cast<THnSparseF *>(fList->FindObject("openingangless"));
1101   THnSparseF *openingangler  = dynamic_cast<THnSparseF *>(fList->FindObject("openingangler"));
1102   THnSparseF *openingangleos = dynamic_cast<THnSparseF *>(fList->FindObject("openingangleos"));
1103
1104   THnSparseF *openinganglegamma = dynamic_cast<THnSparseF *>(fList->FindObject("openinganglegamma"));
1105   THnSparseF *openinganglepi0 = dynamic_cast<THnSparseF *>(fList->FindObject("openinganglepi0"));
1106   THnSparseF *openingangleC = dynamic_cast<THnSparseF *>(fList->FindObject("openingangleC"));
1107   THnSparseF *openingangleB = dynamic_cast<THnSparseF *>(fList->FindObject("openingangleB"));
1108   THnSparseF *openingangleeta = dynamic_cast<THnSparseF *>(fList->FindObject("openingangleeta"));
1109
1110
1111   // invariant mass
1112   THnSparseF *invmasspp = dynamic_cast<THnSparseF *>(fList->FindObject("invmasspp"));
1113   THnSparseF *invmassnn = dynamic_cast<THnSparseF *>(fList->FindObject("invmassnn"));
1114   THnSparseF *invmassss = dynamic_cast<THnSparseF *>(fList->FindObject("invmassss"));
1115   THnSparseF *invmassr  = dynamic_cast<THnSparseF *>(fList->FindObject("invmassr"));
1116   THnSparseF *invmassos = dynamic_cast<THnSparseF *>(fList->FindObject("invmassos"));
1117
1118   THnSparseF *invmassgamma = dynamic_cast<THnSparseF *>(fList->FindObject("invmassgamma"));
1119   THnSparseF *invmasspi0 = dynamic_cast<THnSparseF *>(fList->FindObject("invmasspi0"));
1120   THnSparseF *invmassC = dynamic_cast<THnSparseF *>(fList->FindObject("invmassC"));
1121   THnSparseF *invmassB = dynamic_cast<THnSparseF *>(fList->FindObject("invmassB"));
1122   THnSparseF *invmasseta = dynamic_cast<THnSparseF *>(fList->FindObject("invmasseta"));
1123
1124   // Projection over all pt
1125   TH1D *openingangleppproj = openinganglepp->Projection(2);
1126   TH1D *openinganglennproj = openinganglenn->Projection(2);
1127   TH1D *openinganglessproj = openingangless->Projection(2);
1128   TH1D *openinganglerproj  = openingangler->Projection(2);
1129   TH1D *openingangleosproj = openingangleos->Projection(2);
1130
1131   TH1D *openinganglegammaproj = 0x0;
1132   TH1D *openinganglepi0proj = 0x0;
1133   TH1D *openingangleCproj = 0x0;
1134   TH1D *openingangleBproj = 0x0;
1135   TH1D *openingangleetaproj = 0x0;
1136   if(openinganglegamma) openinganglegammaproj = openinganglegamma->Projection(2);
1137   if(openinganglepi0) openinganglepi0proj = openinganglepi0->Projection(2);
1138   if(openingangleC) openingangleCproj = openingangleC->Projection(2);
1139   if(openingangleB) openingangleBproj = openingangleB->Projection(2);
1140   if(openingangleeta) openingangleetaproj = openingangleeta->Projection(2);
1141
1142
1143   TH1D *invmassppproj = invmasspp->Projection(2);
1144   TH1D *invmassnnproj = invmassnn->Projection(2);
1145   TH1D *invmassssproj = invmassss->Projection(2);
1146   TH1D *invmassrproj  = invmassr->Projection(2);
1147   TH1D *invmassosproj = invmassos->Projection(2);
1148
1149   TH1D *invmassgammaproj = 0x0;
1150   TH1D *invmasspi0proj = 0x0;
1151   TH1D *invmassCproj = 0x0;
1152   TH1D *invmassBproj = 0x0;
1153   TH1D *invmassetaproj = 0x0;
1154   if(invmassgamma) invmassgammaproj = invmassgamma->Projection(2);
1155   if(invmasspi0) invmasspi0proj = invmasspi0->Projection(2);
1156   if(invmassC) invmassCproj = invmassC->Projection(2);
1157   if(invmassB) invmassBproj = invmassB->Projection(2);
1158   if(invmasseta) invmassetaproj = invmasseta->Projection(2);
1159
1160   openingangleppproj->SetStats(0);
1161   openinganglennproj->SetStats(0);
1162   openinganglessproj->SetStats(0);
1163   openinganglerproj->SetStats(0);
1164   openingangleosproj->SetStats(0);
1165   if(openinganglegammaproj) openinganglegammaproj->SetStats(0);
1166   if(openinganglepi0proj) openinganglepi0proj->SetStats(0);
1167   if(openingangleCproj) openingangleCproj->SetStats(0);
1168   if(openingangleBproj) openingangleBproj->SetStats(0);
1169   if(openingangleetaproj) openingangleetaproj->SetStats(0);
1170
1171   invmassppproj->SetStats(0);
1172   invmassnnproj->SetStats(0);
1173   invmassssproj->SetStats(0);
1174   invmassrproj->SetStats(0);
1175   invmassosproj->SetStats(0);
1176   if(invmassgammaproj) invmassgammaproj->SetStats(0);
1177   if(invmasspi0proj) invmasspi0proj->SetStats(0);
1178   if(invmassCproj) invmassCproj->SetStats(0);
1179   if(invmassBproj) invmassBproj->SetStats(0);
1180   if(invmassetaproj) invmassetaproj->SetStats(0);
1181
1182   openingangleppproj->SetTitle("");
1183   openinganglennproj->SetTitle("");
1184   openinganglessproj->SetTitle("");
1185   openinganglerproj->SetTitle("");
1186   openingangleosproj->SetTitle("");
1187   if(openinganglegammaproj) openinganglegammaproj->SetTitle("");
1188   if(openinganglepi0proj) openinganglepi0proj->SetTitle("");
1189   if(openingangleCproj) openingangleCproj->SetTitle("");
1190   if(openingangleBproj) openingangleBproj->SetTitle("");
1191   if(openingangleetaproj) openingangleetaproj->SetTitle("");
1192
1193   invmassppproj->SetTitle("");
1194   invmassnnproj->SetTitle("");
1195   invmassssproj->SetTitle("");
1196   invmassrproj->SetTitle("");
1197   invmassosproj->SetTitle("");
1198   if(invmassgammaproj) invmassgammaproj->SetTitle("");
1199   if(invmasspi0proj) invmasspi0proj->SetTitle("");
1200   if(invmassCproj) invmassCproj->SetTitle("");
1201   if(invmassBproj) invmassBproj->SetTitle("");
1202   if(invmassetaproj) invmassetaproj->SetTitle("");
1203
1204   // Draw histograms for opening angle
1205   TCanvas * copeningangle =new TCanvas("openingangle","Openingangle",800,800);
1206   copeningangle->cd();
1207   openingangleppproj->Draw();
1208   openinganglennproj->Draw("same");
1209   openinganglessproj->Draw("same");
1210   openinganglerproj->Draw("same");
1211   openingangleosproj->Draw("same");
1212   if(openinganglegammaproj) openinganglegammaproj->Draw("same");
1213   if(openinganglepi0proj) openinganglepi0proj->Draw("same");
1214   if(openingangleCproj) openingangleCproj->Draw("same");
1215   if(openingangleBproj) openingangleBproj->Draw("same");
1216   if(openingangleetaproj) openingangleetaproj->Draw("same");
1217   TLegend *lego = new TLegend(0.4,0.6,0.89,0.89);
1218   lego->AddEntry(openingangleppproj,"positive-positive","p");
1219   lego->AddEntry(openinganglennproj,"negative-negative","p");
1220   lego->AddEntry(openinganglessproj,"same-sign","p");
1221   lego->AddEntry(openinganglerproj,"rotated","p");
1222   lego->AddEntry(openingangleosproj,"positive-negative","p");
1223   if(openinganglegammaproj) lego->AddEntry(openinganglegammaproj,"e^{+}e^{-} from #gamma","p");
1224   if(openinganglepi0proj) lego->AddEntry(openinganglepi0proj,"e^{+}e^{-} from #pi^{0}","p");
1225   if(openingangleCproj) lego->AddEntry(openingangleCproj,"e^{+}e^{-} from c","p");
1226   if(openingangleBproj) lego->AddEntry(openingangleBproj,"e^{+}e^{-} from b","p");
1227   if(openingangleetaproj) lego->AddEntry(openingangleetaproj,"e^{+}e^{-} from #eta","p");
1228   lego->Draw("same");
1229
1230   // Draw histograms for opening angle
1231   TCanvas * cinvmass =new TCanvas("invmass","Invmass",800,800);
1232   cinvmass->cd();
1233   invmassppproj->Draw();
1234   invmassnnproj->Draw("same");
1235   invmassssproj->Draw("same");
1236   invmassrproj->Draw("same");
1237   invmassosproj->Draw("same");
1238   if(invmassgammaproj) invmassgammaproj->Draw("same");
1239   if(invmasspi0proj) invmasspi0proj->Draw("same");
1240   if(invmassCproj) invmassCproj->Draw("same");
1241   if(invmassBproj) invmassBproj->Draw("same");
1242   if(invmassetaproj) invmassetaproj->Draw("same");
1243   TLegend *legi = new TLegend(0.4,0.6,0.89,0.89);
1244   legi->AddEntry(invmassppproj,"positive-positive","p");
1245   legi->AddEntry(invmassnnproj,"negative-negative","p");
1246   legi->AddEntry(invmassssproj,"same-sign","p");
1247   legi->AddEntry(invmassrproj,"rotated","p");
1248   legi->AddEntry(invmassosproj,"positive-negative","p");
1249   if(invmassgammaproj) legi->AddEntry(invmassgammaproj,"e^{+}e^{-} from #gamma","p");
1250   if(invmasspi0proj) legi->AddEntry(invmasspi0proj,"e^{+}e^{-} from #pi^{0}","p");
1251   if(invmassCproj) legi->AddEntry(invmassCproj,"e^{+}e^{-} from c","p");
1252   if(invmassBproj) legi->AddEntry(invmassBproj,"e^{+}e^{-} from b","p");
1253   if(invmassetaproj) legi->AddEntry(invmassetaproj,"e^{+}e^{-} from #eta","p");
1254   legi->Draw("same");
1255
1256 }