PWG3/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 /////////////////////////////////////////////////////////////
  17 //
  18 // class used to extract and store info at event level
  19 //
  20 // Author: X-M. Zhang, zhang@clermont.in2p3.fr
  21 //                     zhangxm@iopp.ccnu.edu.cn
  22 /////////////////////////////////////////////////////////////
  23
  24 #include <TMath.h>
  25 #include <TH1.h>
  26 #include <TH1F.h>
  27 #include <TH2F.h>
  28 #include <TList.h>
  29 #include <TObjArray.h>
  30 #include <TObjString.h>
  31
  32 #include "AliTriggerAnalysis.h"
  33 #include "AliBackgroundSelection.h"
  34 #include "AliAODVertex.h"
  35 #include "AliAODEvent.h"
  36 #include "AliESDEvent.h"
  37 #include "AliMultiplicity.h"
  38 #include "AliMuonInfoStoreRD.h"
  39 #include "AliMuonInfoStoreMC.h"
  40 #include "AliDimuInfoStoreRD.h"
  41 #include "AliDimuInfoStoreMC.h"
  42 #include "AliMuonsHFHeader.h"
  43
  44 class TNamed;
  45 class AliESDVertex;
  46
  47 ClassImp(AliMuonsHFHeader)
  48
  49 const TString AliMuonsHFHeader::fgkStdBranchName("MuEvsH");
  50 Int_t         AliMuonsHFHeader::fgAnaMode         = 0;
  51 Bool_t        AliMuonsHFHeader::fgIsMC            = kFALSE;
  52 Bool_t        AliMuonsHFHeader::fgIsEventSelected = kFALSE;
  53 Double_t      AliMuonsHFHeader::fgCuts[3] = { -999999., 999999., 999999.};
  54
  55 //_____________________________________________________________________________
  56 AliMuonsHFHeader::AliMuonsHFHeader() :
  57 TNamed(),
  58 fTriggerMask(0),
  59 fFiredTrigger(0),
  60 fNFiredTrigger(0),
  61 fIsPhysicsTriggered(kFALSE),
  62 fIsPhysicsAccepted(kFALSE),
  63 fEventType(0),
  64 fUnrecoVertex(kFALSE),
  65 fNContributors(0),
  66 fUnrecoVtxSPD(kFALSE),
  67 fNContributorsSPD(0),
  68 fNTrackletsSPD(0),
  69 fCentrality(0.)
  70 {
  71   //
  72   // default constructor
  73   //
  74   for (Int_t i=3; i--;) fVtx[i] = 0.;
  75   for (Int_t i=3; i--;) fVtxSPD[i] = 0.;
  76 }
  77
  78 //_____________________________________________________________________________
  79 AliMuonsHFHeader::AliMuonsHFHeader(const AliMuonsHFHeader &src) :
  80 TNamed(),
  81 fTriggerMask(src.fTriggerMask),
  82 fFiredTrigger(src.fFiredTrigger),
  83 fNFiredTrigger(src.fNFiredTrigger),
  84 fIsPhysicsTriggered(src.fIsPhysicsTriggered),
  85 fIsPhysicsAccepted(src.fIsPhysicsAccepted),
  86 fEventType(src.fEventType),
  87 fUnrecoVertex(src.fUnrecoVertex),
  88 fNContributors(src.fNContributors),
  89 fUnrecoVtxSPD(src.fUnrecoVtxSPD),
  90 fNContributorsSPD(src.fNContributorsSPD),
  91 fNTrackletsSPD(src.fNTrackletsSPD),
  92 fCentrality(src.fCentrality)
  93 {
  94   //
  95   // copy constructor
  96   //
  97   for (Int_t i=3; i--;) fVtx[i]    = src.fVtx[i];
  98   for (Int_t i=3; i--;) fVtxSPD[i] = src.fVtxSPD[i];
  99 }
 100
 101 //_____________________________________________________________________________
 102 AliMuonsHFHeader& AliMuonsHFHeader::operator=(const AliMuonsHFHeader &src)
 103 {
 104   //
 105   // assignment constructor
 106   //
 107   fTriggerMask        = src.fTriggerMask;
 108   fFiredTrigger       = src.fFiredTrigger;
 109   fNFiredTrigger      = src.fNFiredTrigger;
 110   fIsPhysicsTriggered = src.fIsPhysicsTriggered;
 111   fIsPhysicsAccepted  = src.fIsPhysicsAccepted;
 112   fEventType          = src.fEventType;
 113   fUnrecoVertex       = src.fUnrecoVertex;
 114   fNContributors      = src.fNContributors;
 115   fUnrecoVtxSPD       = src.fUnrecoVtxSPD;
 116   fNContributorsSPD   = src.fNContributorsSPD;
 117   fNTrackletsSPD      = src.fNTrackletsSPD;
 118   fCentrality         = src.fCentrality;
 119
 120   for (Int_t i=3; i--;) fVtx[i]      = src.fVtx[i];
 121   for (Int_t i=3; i--;) fVtxSPD[i]   = src.fVtxSPD[i];
 122
 123   return *this;
 124 }
 125
 126 //_____________________________________________________________________________
 127 AliMuonsHFHeader::~AliMuonsHFHeader()
 128 {
 129   //
 130   // default destructor
 131   //
 132 }
 133
 134 //_____________________________________________________________________________
 135 void AliMuonsHFHeader::SetEvent(AliAODEvent *event)
 136 {
 137   // extract event info from AOD event
 138
 139   fTriggerMask = event->GetTriggerMask();
 140
 141   AliAODVertex *vertex = event->GetPrimaryVertex();
 142   vertex->GetXYZ(fVtx);
 143   fNContributors = vertex->GetNContributors();
 144   fUnrecoVertex = (TMath::Abs(fVtx[0])<1e-6 && TMath::Abs(fVtx[1])<1e-6 &&
 145                    TMath::Abs(fVtx[2])<1e-6);
 146   this->SetTitle(vertex->GetTitle());
 147
 148   this->SetFiredTrigger(event->GetFiredTriggerClasses());
 149   this->EventSelection();
 150   return;
 151 }
 152
 153 //_____________________________________________________________________________
 154 void AliMuonsHFHeader::SetEvent(AliESDEvent *event)
 155 {
 156   // extract event info from ESD event
 157   // SPD vertex and Physics event selection are implimented
 158
 159   fTriggerMask = event->GetTriggerMask();
 160
 161   const AliESDVertex *vertex = event->GetPrimaryVertex();
 162   vertex->GetXYZ(fVtx);
 163   fNContributors = vertex->GetNContributors();
 164   fUnrecoVertex = (TMath::Abs(fVtx[0])<1e-6 && TMath::Abs(fVtx[1])<1e-6 &&
 165                    TMath::Abs(fVtx[2])<1e-6);
 166   this->SetTitle(vertex->GetTitle());
 167
 168   const AliESDVertex *vtxSPD = event->GetPrimaryVertexSPD();
 169   vtxSPD->GetXYZ(fVtxSPD);
 170   fNContributorsSPD = vtxSPD->GetNContributors();
 171   fUnrecoVtxSPD = (TMath::Abs(fVtxSPD[0])<1e-6 && TMath::Abs(fVtxSPD[1])<1e-6 &&
 172                    TMath::Abs(fVtxSPD[2])<1e-6);
 173   fNTrackletsSPD = event->GetMultiplicity()->GetNumberOfTracklets();
 174
 175   this->SetFiredTrigger(event->GetFiredTriggerClasses());
 176   this->PhysicsTriggerAna(event);
 177   this->EventSelection();
 178   return;
 179 }
 180
 181 //_____________________________________________________________________________
 182 void AliMuonsHFHeader::PhysicsTriggerAna(const AliESDEvent *esd)
 183 {
 184   // ESD event trigger condition analysis
 185   // according to the method in $ALICE_ROOT/ANALYSIS/AliPhysicsSelection.cxx
 186
 187   fEventType = esd->GetHeader()->GetEventType();
 188   fIsPhysicsTriggered = kFALSE;
 189   fIsPhysicsAccepted  = kFALSE;
 190
 191   AliTriggerAnalysis *triggerAna = new AliTriggerAnalysis();
 192   triggerAna->SetAnalyzeMC(fgIsMC);
 193   triggerAna->SetSPDGFOThreshhold(1);
 194
 195   Int_t  triggerHW  = triggerAna->SPDFiredChips(esd, 1);
 196   Bool_t isFiredv0A = triggerAna->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0A);
 197   Bool_t isFiredv0C = triggerAna->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0C);
 198   if (triggerHW==0 && !isFiredv0A && !isFiredv0C) {
 199     delete triggerAna;
 200     triggerAna = 0;
 201     return;
 202   }
 203
 204   Bool_t triggerBG = (triggerAna->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0ABG) ||
 205                       triggerAna->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0CBG));
 206   Int_t  isFiredSPD = triggerAna->SPDFiredChips(esd, 0);
 207   Bool_t triggerFD = ((isFiredSPD>1) || (isFiredSPD>0 && (isFiredv0A || isFiredv0C)) || (isFiredv0A && isFiredv0C));
 208   if ((!triggerBG) && triggerFD) fIsPhysicsTriggered = kTRUE;
 209   delete triggerAna;
 210   triggerAna = 0;
 211
 212   if (fIsPhysicsTriggered) {
 213     AliBackgroundSelection *bkgId = new AliBackgroundSelection();
 214     fIsPhysicsAccepted = bkgId->IsSelected(const_cast<AliESDEvent*>(esd));
 215     delete bkgId; bkgId=0;
 216   }
 217   return;
 218 }
 219
 220 //_____________________________________________________________________________
 221 void AliMuonsHFHeader::SetFiredTrigger(TString trigger)
 222 {
 223   // get info of fired trigger classes of event
 224
 225   fFiredTrigger = trigger;
 226   TObjArray *tokens = trigger.Tokenize(" ");
 227   fNFiredTrigger = tokens->GetEntries();
 228   return;
 229 }
 230
 231 //_____________________________________________________________________________
 232 Bool_t AliMuonsHFHeader::IsTriggerFired(TString trigger)
 233 {
 234   // check whether the trigger class "trigger" is fired in this event
 235
 236   if (fNFiredTrigger<=0) return kFALSE;
 237
 238   TObjArray *tokens = fFiredTrigger.Tokenize(" ");
 239   for (Int_t i=fNFiredTrigger; i--;) {
 240     TString fired = ((TObjString*)tokens->At(i))->String();
 241     if (fired.CompareTo(trigger)==0) return kTRUE;
 242   }
 243   return kFALSE;
 244 }
 245
 246 //_____________________________________________________________________________
 247 void AliMuonsHFHeader::EventSelection(TString triggerName)
 248 {
 249   // select event according to the "triggerName" & event selection cuts
 250
 251   fgIsEventSelected = kFALSE;
 252   if (!this->IsPhysicsAccepted()) return;
 253   if (!fgIsMC && !this->IsTriggerFired(triggerName)) return;
 254   this->EventSelection();
 255   return;
 256 }
 257
 258 //_____________________________________________________________________________
 259 void AliMuonsHFHeader::EventSelection()
 260 {
 261   // select event according to the event selection cuts
 262
 263   fgIsEventSelected = kFALSE;
 264   if (this->NVtxContributorsSPD()<fgCuts[0]) return;
 265   if (TMath::Abs(this->VzSPD())>fgCuts[1]) return;
 266   if (this->VtSPD()>fgCuts[2]) return;
 267   fgIsEventSelected = kTRUE;
 268   return;
 269 }
 270
 271 //_____________________________________________________________________________
 272 void AliMuonsHFHeader::CreateHistograms(TList *listEvent, TList *listMuon, TList *listDimu)
 273 {
 274   // create output histos of muon analysis according to the analysis mode & MC flag
 275
 276   this->CreateHistosEventH(listEvent);
 277   if (fgIsMC) {
 278     if (fgAnaMode!=2) this->CreateHistosMuonMC(listMuon);
 279     if (fgAnaMode!=1) this->CreateHistosDimuMC(listDimu);
 280   } else {
 281     if (fgAnaMode!=2) this->CreateHistosMuonRD(listMuon);
 282     if (fgAnaMode!=1) this->CreateHistosDimuRD(listDimu);
 283   }
 284   return;
 285 }
 286
 287 //_____________________________________________________________________________
 288 void AliMuonsHFHeader::CreateHistosEventH(TList *list)
 289 {
 290   // create histograms at event level
 291
 292   if (!list) list = new TList();
 293   list->SetOwner();
 294   Bool_t oldStatus = TH1::AddDirectoryStatus();
 295   TH1::AddDirectory(kFALSE);
 296
 297   const Int_t nHistos = 5;
 298   char    *name[nHistos] = { "hVztSPD", "hVzxSPD", "hVzySPD", "hVzNcontr", "hVtNcontr"  };
 299   Int_t  nbinsX[nHistos] = {      800 ,      800 ,      800 ,       800  ,        400   };
 300   Double_t xlow[nHistos] = {      -40.,      -40.,      -40.,       -40. ,          0.  };
 301   Double_t  xup[nHistos] = {       40.,       40.,       40.,        40. ,          4.  };
 302   Int_t  nbinsY[nHistos] = {      400 ,      600 ,      600 ,       202  ,        202   };
 303   Double_t ylow[nHistos] = {        0.,       -3.,       -3.,        -2.5,         -2.5 };
 304   Double_t  yup[nHistos] = {        4.,        3.,        3.,       199.5,        199.5 };
 305
 306   TH2F *histo = 0;
 307   for (Int_t i=0; i<nHistos; i++) {
 308     histo = new TH2F(name[i], name[i], nbinsX[i], xlow[i], xup[i], nbinsY[i], ylow[i], yup[i]);
 309     histo->Sumw2(); list->AddAt(histo, i); histo = 0;
 310   }
 311
 312   TH1::AddDirectory(oldStatus);
 313   return;
 314 }
 315
 316 //_____________________________________________________________________________
 317 void AliMuonsHFHeader::CreateHistosMuonRD(TList *list)
 318 {
 319   // create histograms for single muon
 320
 321   if (!list) list = new TList();
 322   list->SetOwner();
 323   Bool_t oldStatus = TH1::AddDirectoryStatus();
 324   TH1::AddDirectory(kFALSE);
 325
 326   const Int_t nHistos = 8;
 327   char    *name[nHistos] = {  "hP", "hPt", "hEta", "hDCA", "hTrg", "hCharge", "hUnfVtx" , "hEtaTimesDCA"};
 328   Int_t   nbins[nHistos] = { 1500 ,  300 ,   100 ,   500 ,    4  ,       3  ,       80  ,         10000 };
 329   Double_t xlow[nHistos] = {    0.,    0.,   -10.,     0.,   -0.5,      -1.5,      -40. ,             0.};
 330   Double_t  xup[nHistos] = {  150.,   30.,     0.,   500.,    3.5,       1.5,       40. ,          1000.};
 331
 332   TH1F *histo = 0;
 333   for (Int_t i=0; i<nHistos; i++) {
 334     histo = new TH1F(name[i], name[i], nbins[i], xlow[i], xup[i]);
 335     histo->Sumw2(); list->AddAt(histo, i); histo = 0;
 336   }
 337
 338   TH1::AddDirectory(oldStatus);
 339   return;
 340 }
 341
 342 //_____________________________________________________________________________
 343 void AliMuonsHFHeader::CreateHistosDimuRD(TList *list)
 344 {
 345   // create histograms for dimuon
 346
 347   if (!list) list = new TList();
 348   list->SetOwner();
 349   Bool_t oldStatus = TH1::AddDirectoryStatus();
 350   TH1::AddDirectory(kFALSE);
 351
 352   TH1F *histo = 0;
 353   const Int_t nHistos = 3;
 354   char    *name[nHistos] = {  "hP", "hPt", "hInvM" };
 355   Int_t   nbins[nHistos] = { 1500 ,  300 ,    300  };
 356   Double_t xlow[nHistos] = {    0.,    0.,      0. };
 357   Double_t  xup[nHistos] = {  150.,   30.,     30. };
 358   char *dimuName[3] = {"DimuNN", "DimuNP", "DimuPP"};
 359   for (Int_t i=0; i<3; i++) {
 360     for (Int_t j=0; j<nHistos; j++) {
 361       histo = new TH1F(Form("%s_%s",name[j],dimuName[i]), name[j], nbins[j], xlow[j], xup[j]);
 362       histo->Sumw2(); list->AddAt(histo,i*nHistos+j); histo = 0;
 363     }
 364   }
 365
 366   TH1::AddDirectory(oldStatus);
 367   return;
 368 }
 369
 370 //_____________________________________________________________________________
 371 void AliMuonsHFHeader::CreateHistosMuonMC(TList *list)
 372 {
 373   // create histograms for single muon with MC info
 374
 375   if (!list) list = new TList[AliMuonInfoStoreMC::NSources()];
 376   for (Int_t i=AliMuonInfoStoreMC::NSources(); i--;) this->CreateHistosMuonRD(&list[i]);
 377   return;
 378 }
 379
 380 //_____________________________________________________________________________
 381 void AliMuonsHFHeader::CreateHistosDimuMC(TList *list)
 382 {
 383   // create histograms for dimuon with MC info
 384
 385   if (!list) list = new TList[AliDimuInfoStoreMC::NSources()];
 386   for (Int_t i=AliDimuInfoStoreMC::NSources(); i--;) this->CreateHistosDimuRD(&list[i]);
 387   return;
 388 }
 389
 390 //_____________________________________________________________________________
 391 void AliMuonsHFHeader::FillHistosEventH(TList *list)
 392 {
 393   // fill histograms at event level according to event selection cuts
 394
 395   if (!list) return;
 396   if (!AliMuonsHFHeader::IsEventSelected()) return;
 397
 398   const Int_t nHistos = 5;
 399   Double_t vz = this->VzSPD();
 400   Double_t vt = this->VtSPD();
 401   Int_t    nc = this->NVtxContributorsSPD();
 402   Double_t distX[nHistos] = { vz,            vz,            vz, vz, vt };
 403   Double_t distY[nHistos] = { vt, this->VxSPD(), this->VySPD(), nc, nc };
 404   for (Int_t i=nHistos; i--;) ((TH2F*)list->At(i))->Fill(distX[i], distY[i]);
 405   return;
 406 }
 407
 408 //_____________________________________________________________________________
 409 void AliMuonsHFHeader::FillHistosMuonRD(TList *list, AliMuonInfoStoreRD* const muonStoreRD)
 410 {
 411   // fill histograms for single muon according to event & muon track selection cuts
 412
 413   if (!list) return;
 414   if (!AliMuonsHFHeader::IsEventSelected()) return;
 415   if (!muonStoreRD->MuonSelection()) return;
 416
 417   const Int_t nHistos = 8;
 418   Double_t dist[nHistos] = {muonStoreRD->Momentum().Mag(),
 419                             muonStoreRD->Momentum().Pt(),
 420                             muonStoreRD->Momentum().Eta(),
 421                             muonStoreRD->DCA(),
 422                             muonStoreRD->MatchTrigger(),
 423                             muonStoreRD->Charge(),
 424                             this->VzSPD(),
 425                             muonStoreRD->DCA()*TMath::Exp(-2.*muonStoreRD->Momentum().Eta())/1000.};
 426   for (Int_t i=nHistos; i--;) ((TH1F*)list->At(i))->Fill(dist[i]);
 427   return;
 428 }
 429
 430 //_____________________________________________________________________________
 431 void AliMuonsHFHeader::FillHistosDimuRD(TList *list, AliDimuInfoStoreRD* const dimuStoreRD)
 432 {
 433   // fill histograms for dimuon according to evnet & dimuon candidates selection cuts
 434
 435   if (!list) return;
 436   if (!AliMuonsHFHeader::IsEventSelected()) return;
 437   if (!dimuStoreRD->DimuSelection()) return;
 438
 439   Int_t theDimu = 0;
 440   if (dimuStoreRD->Charge()==0)     theDimu = 1;
 441   else if (dimuStoreRD->Charge()>0) theDimu = 2;
 442
 443   const Int_t nHistos = 3;
 444   Double_t dist[nHistos] = {dimuStoreRD->Momentum().Mag(),
 445                             dimuStoreRD->Momentum().Pt(),
 446                             dimuStoreRD->InvM()};
 447   for (Int_t i=nHistos; i--;) ((TH1F*)list->At(theDimu*nHistos+i))->Fill(dist[i]);
 448   return;
 449 }
 450
 451 //_____________________________________________________________________________
 452 void AliMuonsHFHeader::FillHistosMuonMC(TList *list, AliMuonInfoStoreMC* const muonStoreMC)
 453 {
 454   // fill histograms for single muon with MC info
 455
 456   Int_t srcMuon = muonStoreMC->MuonSource();
 457   this->FillHistosMuonRD(&list[6], (AliMuonInfoStoreRD*)muonStoreMC);
 458   if (srcMuon<0) this->FillHistosMuonRD(&list[5], (AliMuonInfoStoreRD*)muonStoreMC);
 459   else this->FillHistosMuonRD(&list[srcMuon], (AliMuonInfoStoreRD*)muonStoreMC);
 460   return;
 461 }
 462
 463 //_____________________________________________________________________________
 464 void AliMuonsHFHeader::FillHistosDimuMC(TList *list, AliDimuInfoStoreMC* const dimuStoreMC)
 465 {
 466   // fill histograms for dimuon with MC info
 467
 468   Int_t srcDimu = dimuStoreMC->DimuSource();
 469   this->FillHistosDimuRD(&list[6], (AliDimuInfoStoreRD*)dimuStoreMC);
 470   this->FillHistosDimuRD(&list[srcDimu], (AliDimuInfoStoreRD*)dimuStoreMC);
 471   return;
 472 }