PWG/muon/AliMuonsHFHeader.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2006, 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$ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // class used to extract and store info at event level
  21 //
  22 // Author: X-M. Zhang, zhang@clermont.in2p3.fr
  23 //                     zhangxm@iopp.ccnu.edu.cn
  24 /////////////////////////////////////////////////////////////
  25
  26 #include <TMath.h>
  27 #include <TH1.h>
  28 #include <TH1D.h>
  29 #include <TList.h>
  30
  31 #include "AliInputEventHandler.h"
  32 #include "AliAODMCHeader.h"
  33 #include "AliVEvent.h"
  34 #include "AliAODEvent.h"
  35 #include "AliESDEvent.h"
  36 #include "AliVVertex.h"
  37 #include "AliCentrality.h"
  38 #include "AliEventplane.h"
  39
  40 #include "AliMuonTrackCuts.h"
  41 #include "AliMuonInfoStoreRD.h"
  42 #include "AliMuonInfoStoreMC.h"
  43 #include "AliDimuInfoStoreRD.h"
  44 #include "AliDimuInfoStoreMC.h"
  45 #include "AliMuonsHFHeader.h"
  46
  47 ClassImp(AliMuonsHFHeader)
  48
  49 const TString AliMuonsHFHeader::fgkStdBranchName("MuEvsH");
  50 Int_t         AliMuonsHFHeader::fgAnaMode = 0;
  51 Bool_t        AliMuonsHFHeader::fgIsMC    = kFALSE;
  52 Double_t      AliMuonsHFHeader::fgCuts[5] = { -999999., 999999., 999999., -999999., 999999. };
  53
  54 //_____________________________________________________________________________
  55 AliMuonsHFHeader::AliMuonsHFHeader() :
  56 TNamed(),
  57 fSelMask(AliVEvent::kAny),
  58 fIsMB(kFALSE),
  59 fIsMU(kFALSE),
  60 fIsPileupSPD(kFALSE),
  61 fVtxContrsN(0),
  62 fFiredTriggerClass(),
  63 fCentrality(-1.),
  64 fCentQA(-1),
  65 fEventPlane(0.)
  66 {
  67   //
  68   // default constructor
  69   //
  70   for (Int_t i=3; i--;) fVtx[i] = 0.;
  71   for (Int_t i=3; i--;) fVMC[i] = 0.;
  72 }
  73
  74 //_____________________________________________________________________________
  75 AliMuonsHFHeader::AliMuonsHFHeader(const AliMuonsHFHeader &src) :
  76 TNamed(),
  77 fSelMask(src.fSelMask),
  78 fIsMB(src.fIsMB),
  79 fIsMU(src.fIsMU),
  80 fIsPileupSPD(src.fIsPileupSPD),
  81 fVtxContrsN(src.fVtxContrsN),
  82 fFiredTriggerClass(src.fFiredTriggerClass),
  83 fCentrality(src.fCentrality),
  84 fCentQA(src.fCentQA),
  85 fEventPlane(src.fEventPlane)
  86 {
  87   //
  88   // copy constructor
  89   //
  90   for (Int_t i=3; i--;) fVtx[i] = src.fVtx[i];
  91   for (Int_t i=3; i--;) fVMC[i] = src.fVMC[i];
  92 }
  93
  94 //_____________________________________________________________________________
  95 AliMuonsHFHeader& AliMuonsHFHeader::operator=(const AliMuonsHFHeader &src)
  96 {
  97   //
  98   // assignment constructor
  99   //
 100
 101   if(&src==this) return *this;
 102
 103   fSelMask           = src.fSelMask;
 104   fIsMB              = src.fIsMB;
 105   fIsMU              = src.fIsMU;
 106   fIsPileupSPD       = src.fIsPileupSPD;
 107   fVtxContrsN        = src.fVtxContrsN;
 108   fFiredTriggerClass = src.fFiredTriggerClass;
 109   fCentrality        = src.fCentrality;
 110   fCentQA            = src.fCentQA;
 111   fEventPlane        = src.fEventPlane;
 112   for (Int_t i=3; i--;) fVtx[i] = src.fVtx[i];
 113   for (Int_t i=3; i--;) fVMC[i] = src.fVMC[i];
 114
 115   return *this;
 116 }
 117
 118 //_____________________________________________________________________________
 119 AliMuonsHFHeader::~AliMuonsHFHeader()
 120 {
 121   //
 122   // default destructor
 123   //
 124 }
 125
 126 //_____________________________________________________________________________
 127 void AliMuonsHFHeader::SetEventInfo(AliInputEventHandler* const handler, AliMCEvent* const eventMC)
 128 {
 129   // fill info at event level
 130
 131   AliVEvent *event = handler->GetEvent();
 132   AliAODEvent *aod = dynamic_cast<AliAODEvent*>(event);
 133   AliESDEvent *esd = dynamic_cast<AliESDEvent*>(event);
 134
 135   fSelMask = handler->IsEventSelected();
 136   if (aod) fFiredTriggerClass = aod->GetFiredTriggerClasses();
 137   if (esd) fFiredTriggerClass = esd->GetFiredTriggerClasses();
 138   fIsMB = fSelMask & AliVEvent::kMB;
 139   fIsMU = fSelMask & AliVEvent::kMUON;
 140
 141   const AliVVertex *vertex = event->GetPrimaryVertex();
 142   vertex->GetXYZ(fVtx);
 143   fVtxContrsN = vertex->GetNContributors();
 144   if (fgIsMC) {
 145     if (esd)   eventMC->GetPrimaryVertex()->GetXYZ(fVMC);
 146     if (aod) ((AliAODMCHeader*)aod->FindListObject(AliAODMCHeader::StdBranchName()))->GetVertex(fVMC);
 147   } this->SetTitle(vertex->GetTitle());
 148   fIsPileupSPD = (aod && !aod->GetTracklets()) ? event->IsPileupFromSPD(3,0.8,3.,2.,5.) : event->IsPileupFromSPDInMultBins();
 149
 150   AliCentrality *cent = event->GetCentrality();
 151   if (cent) {
 152     fCentrality = cent->GetCentralityPercentileUnchecked("V0M");
 153     fCentQA     = cent->GetQuality();
 154   }
 155
 156   AliEventplane *evnP = event->GetEventplane();
 157   if (evnP) fEventPlane = evnP->GetEventplane("Q");
 158 //if (evnP) fEventPlane = evnP->GetEventplane("V0A");
 159
 160   return;
 161 }
 162
 163 //_____________________________________________________________________________
 164 Bool_t AliMuonsHFHeader::IsSelected()
 165 {
 166   // select event according to the event selection cuts
 167   if (this->VtxContrsN()<fgCuts[0])       return kFALSE;
 168   if (TMath::Abs(this->Vz())>fgCuts[1])   return kFALSE;
 169   if (this->Vt()>fgCuts[2])               return kFALSE;
 170
 171   // centrality selection
 172   Float_t centr = this->Centrality();
 173   if (centr<fgCuts[3] || centr>fgCuts[4]) return kFALSE;
 174   return kTRUE;
 175 }
 176
 177 //_____________________________________________________________________________
 178 void AliMuonsHFHeader::CreateHistograms(TList *list)
 179 {
 180   // create output histos of muon analysis according to the analysis mode & MC flag
 181
 182   if (fgIsMC) {
 183     this->CreateHistosEvnH(list);
 184     if (fgAnaMode!=2) {
 185       TString sName[7] = { "Unidentified", "Hadron", "SecondaryMu", "PrimaryMu", "CharmMu", "BottomMu", "" };
 186       for (Int_t i=7; i--;) this->CreateHistosMuon(list, sName[i]);
 187     }
 188     if (fgAnaMode!=1) {
 189       TString sName[7] = { "Uncorr", "Resonance", "DDsame", "DDdiff", "BBsame", "BBdiff", "" };
 190       for (Int_t i=7; i--;) this->CreateHistosDimu(list, sName[i]);
 191     }
 192     return;
 193   }
 194
 195   this->CreateHistosEvnH(list,"MB"); this->CreateHistosEvnH(list,"MU");
 196   if (fgAnaMode!=2) { this->CreateHistosMuon(list,"MB"); this->CreateHistosMuon(list,"MU"); }
 197   if (fgAnaMode!=1) { this->CreateHistosDimu(list,"MB"); this->CreateHistosDimu(list,"MU"); }
 198   return;
 199 }
 200
 201 //_____________________________________________________________________________
 202 void AliMuonsHFHeader::CreateHistosEvnH(TList *list, TString sName)
 203 {
 204   // create histograms at event level
 205
 206   if (!list) list = new TList();
 207   list->SetOwner();
 208   Bool_t oldStatus = TH1::AddDirectoryStatus();
 209   TH1::AddDirectory(kFALSE);
 210
 211   const Int_t nhs    = 3;
 212   TString tName[nhs] = {  "Vz",  "Vt",  "VtxNcontr" };
 213   Int_t   nbins[nhs] = {  800 ,   40 ,        202   };
 214   Double_t xlow[nhs] = {  -40.,    0.,         -2.5 };
 215   Double_t  xup[nhs] = {   40.,    4.,        199.5 };
 216
 217   TH1D *histo = 0;
 218   for (Int_t i=0; i<nhs; i++) {
 219     char *hName = Form("h%s_%s", sName.Data(), tName[i].Data());
 220     histo = new TH1D(hName, hName, nbins[i], xlow[i], xup[i]);
 221     histo->Sumw2(); list->Add(histo); histo = 0;
 222   }
 223
 224   TH1::AddDirectory(oldStatus);
 225   return;
 226 }
 227
 228 //_____________________________________________________________________________
 229 void AliMuonsHFHeader::CreateHistosMuon(TList *list, TString sName)
 230 {
 231   // create histograms for single muon
 232
 233   if (!list) list = new TList();
 234   list->SetOwner();
 235   Bool_t oldStatus = TH1::AddDirectoryStatus();
 236   TH1::AddDirectory(kFALSE);
 237
 238   const Int_t nhs    = 7;
 239   TString tName[nhs] = {   "P",  "Pt",  "Eta",  "DCA",  "TrM",  "Charge", "Rabs" };
 240   Int_t   nbins[nhs] = { 1500 ,  300 ,    15 ,  1000 ,    4  ,       3  ,    48  };
 241   Double_t xlow[nhs] = {    0.,    0.,   -4.0,     0.,   -0.5,      -1.5,   17.6 };
 242   Double_t  xup[nhs] = {  150.,   30.,   -2.5,   500.,    3.5,       1.5,   80.0 };
 243
 244   TH1D *histo = 0;
 245   for (Int_t i=0; i<nhs; i++) {
 246     char *hName = Form("h%s_%s", sName.Data(), tName[i].Data());
 247     histo = new TH1D(hName, hName, nbins[i], xlow[i], xup[i]);
 248     histo->Sumw2(); list->Add(histo); histo = 0;
 249   }
 250
 251   TH1::AddDirectory(oldStatus);
 252   return;
 253 }
 254
 255 //_____________________________________________________________________________
 256 void AliMuonsHFHeader::CreateHistosDimu(TList *list, TString sName)
 257 {
 258   // create histograms for dimuon
 259
 260   if (!list) list = new TList();
 261   list->SetOwner();
 262   Bool_t oldStatus = TH1::AddDirectoryStatus();
 263   TH1::AddDirectory(kFALSE);
 264
 265   TH1D *histo = 0;
 266   const Int_t nhs    = 3;
 267   TString tName[nhs] = {   "P",  "Pt",  "InvM"   };
 268   Int_t   nbins[nhs] = { 1500 ,  300 ,    300    };
 269   Double_t xlow[nhs] = {    0.,    0.,      0.   };
 270   Double_t  xup[nhs] = {  150.,   30.,     30.   };
 271   TString dimuName[3] = { "DimuNN", "DimuNP", "DimuPP" };
 272   for (Int_t i=0; i<3; i++) {
 273     for (Int_t j=0; j<nhs; j++) {
 274       char *hName = Form("h%s_%s_%s", sName.Data(), dimuName[i].Data(), tName[j].Data());
 275       histo = new TH1D(hName, hName, nbins[j], xlow[j], xup[j]);
 276       histo->Sumw2(); list->Add(histo); histo = 0;
 277     }
 278   }
 279
 280   TH1::AddDirectory(oldStatus);
 281   return;
 282 }
 283
 284 //_____________________________________________________________________________
 285 void AliMuonsHFHeader::FillHistosEvnH(TList *list)
 286 {
 287   // fill histograms at event level according to event selection cuts
 288
 289   if (!list)               return;
 290   if (!this->IsSelected()) return;
 291
 292   const Int_t nhs    = 3;
 293   TString tName[nhs] = {       "Vz",       "Vt",        "VtxNcontr" };
 294   Double_t dist[nhs] = { this->Vz(), this->Vt(), this->VtxContrsN() };
 295   if (fgIsMC && (fSelMask & AliVEvent::kAny)) {
 296     for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h_%s",tName[i].Data())))->Fill(dist[i]);
 297   } else {
 298     if (fIsMB && (fSelMask & AliVEvent::kMB))   { for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s","MB",tName[i].Data())))->Fill(dist[i]); }
 299     if (fIsMU && (fSelMask & AliVEvent::kMUON)) { for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s","MU",tName[i].Data())))->Fill(dist[i]); }
 300   }
 301   return;
 302 }
 303
 304 //_____________________________________________________________________________
 305 void AliMuonsHFHeader::FillHistosMuon(TList *list, AliMuonInfoStoreRD* const infoStore, Int_t s)
 306 {
 307   // fill histograms for single muon according to event & muon track selection cuts
 308
 309   if (!list)                     return;
 310   if (!this->IsSelected())       return;
 311   if (!infoStore->IsSelected(0)) return;
 312
 313   const Int_t nhs    = 7;
 314   TString tName[nhs] = { "P", "Pt", "Eta", "DCA", "TrM", "Charge", "Rabs" };
 315   Double_t dist[nhs] = { infoStore->MomentumAtVtx().Mag(),
 316                          infoStore->MomentumAtVtx().Pt(),
 317                          infoStore->MomentumAtVtx().Eta(),
 318                          infoStore->DCA(),
 319                          infoStore->MatchTrigger(),
 320                          infoStore->Charge(),
 321                          infoStore->RabsEnd() };
 322
 323   if (fgIsMC && (fSelMask & AliVEvent::kAny)) {
 324     TString sName[7] = { "BottomMu", "CharmMu", "PrimaryMu", "SecondaryMu", "Hadron", "Unidentified", "" };
 325     for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s",sName[6].Data(),tName[i].Data())))->Fill(dist[i]);
 326     for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s",sName[s].Data(),tName[i].Data())))->Fill(dist[i]);
 327   } else {
 328     if (fIsMB && (fSelMask & AliVEvent::kMB))   { for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s","MB",tName[i].Data())))->Fill(dist[i]); }
 329     if (fIsMU && (fSelMask & AliVEvent::kMUON)) { for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s","MU",tName[i].Data())))->Fill(dist[i]); }
 330   }
 331
 332   return;
 333 }
 334
 335 //_____________________________________________________________________________
 336 void AliMuonsHFHeader::FillHistosDimu(TList *list, AliDimuInfoStoreRD* const infoStore, Int_t s)
 337 {
 338   // fill histograms for dimuon according to evnet & dimuon candidates selection cuts
 339
 340   if (!list)                     return;
 341   if (!this->IsSelected())       return;
 342   if (!infoStore->IsSelected(0)) return;
 343
 344   TString dimuName = "DimuNN";
 345   if (infoStore->Charge()==0)     dimuName = "DimuNP";
 346   else if (infoStore->Charge()>0) dimuName = "DimuPP";
 347
 348   const Int_t nhs    = 3;
 349   TString tName[nhs] = { "P", "Pt", "InvM" };
 350   Double_t dist[nhs] = { infoStore->Momentum().Mag(),
 351                          infoStore->Momentum().Pt(),
 352                          infoStore->InvM() };
 353
 354   if (fgIsMC && (fSelMask & AliVEvent::kAny)) {
 355     TString sName[7] = { "BBdiff", "BBsame", "DDdiff", "DDsame", "Resonance", "Uncorr", "" };
 356     for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s_%s",sName[6].Data(),dimuName.Data(),tName[i].Data())))->Fill(dist[i]);
 357     for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s_%s",sName[s].Data(),dimuName.Data(),tName[i].Data())))->Fill(dist[i]);
 358   } else {
 359     if (fIsMB && (fSelMask & AliVEvent::kMB)) {
 360       for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s_%s","MB",dimuName.Data(),tName[i].Data())))->Fill(dist[i]);
 361     }
 362     if (fIsMU && (fSelMask & AliVEvent::kMUON)) {
 363       for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("h%s_%s_%s","MU",dimuName.Data(),tName[i].Data())))->Fill(dist[i]);
 364     }
 365   }
 366
 367   return;
 368 }