Update of the Xiaoming code for pp900 first muon analysis. Fixing wanirngs and violti...
[u/mrichter/AliRoot.git] / PWG3 / muon / AliMuonsHFHeader.cxx
CommitLineData
fd1d0cb9 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
4292b3b6 24#include <TMath.h>
fd1d0cb9 25#include <TH1.h>
26#include <TH1F.h>
27#include <TH2F.h>
28#include <TList.h>
29#include <TObjArray.h>
30#include <TObjString.h>
4292b3b6 31
fd1d0cb9 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"
4292b3b6 42#include "AliMuonsHFHeader.h"
43
fd1d0cb9 44class TNamed;
45class AliESDVertex;
46
4292b3b6 47ClassImp(AliMuonsHFHeader)
48
fd1d0cb9 49const TString AliMuonsHFHeader::fgkStdBranchName("MuEvsH");
50Int_t AliMuonsHFHeader::fgAnaMode = 0;
51Bool_t AliMuonsHFHeader::fgIsMC = kFALSE;
52Bool_t AliMuonsHFHeader::fgIsEventSelected = kFALSE;
53Double_t AliMuonsHFHeader::fgCuts[3] = { -999999., 999999., 999999.};
54
4292b3b6 55//_____________________________________________________________________________
56AliMuonsHFHeader::AliMuonsHFHeader() :
57TNamed(),
58fTriggerMask(0),
fd1d0cb9 59fFiredTrigger(0),
60fNFiredTrigger(0),
61fIsPhysicsTriggered(kFALSE),
62fIsPhysicsAccepted(kFALSE),
63fEventType(0),
64fUnrecoVertex(kFALSE),
4292b3b6 65fNContributors(0),
fd1d0cb9 66fUnrecoVtxSPD(kFALSE),
67fNContributorsSPD(0),
68fNTrackletsSPD(0),
69fCentrality(0.)
4292b3b6 70{
71 //
72 // default constructor
73 //
fd1d0cb9 74 for (Int_t i=3; i--;) fVtx[i] = 0.;
75 for (Int_t i=3; i--;) fVtxSPD[i] = 0.;
76}
77
78//_____________________________________________________________________________
79AliMuonsHFHeader::AliMuonsHFHeader(const AliMuonsHFHeader &src) :
80TNamed(),
81fTriggerMask(src.fTriggerMask),
82fFiredTrigger(src.fFiredTrigger),
83fNFiredTrigger(src.fNFiredTrigger),
84fIsPhysicsTriggered(src.fIsPhysicsTriggered),
85fIsPhysicsAccepted(src.fIsPhysicsAccepted),
86fEventType(src.fEventType),
87fUnrecoVertex(src.fUnrecoVertex),
88fNContributors(src.fNContributors),
89fUnrecoVtxSPD(src.fUnrecoVtxSPD),
90fNContributorsSPD(src.fNContributorsSPD),
91fNTrackletsSPD(src.fNTrackletsSPD),
92fCentrality(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//_____________________________________________________________________________
102AliMuonsHFHeader& 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;
4292b3b6 124}
125
126//_____________________________________________________________________________
127AliMuonsHFHeader::~AliMuonsHFHeader()
128{
129 //
130 // default destructor
131 //
132}
133
134//_____________________________________________________________________________
fd1d0cb9 135void AliMuonsHFHeader::SetEvent(AliAODEvent *event)
4292b3b6 136{
fd1d0cb9 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();
4292b3b6 150 return;
151}
152
153//_____________________________________________________________________________
fd1d0cb9 154void AliMuonsHFHeader::SetEvent(AliESDEvent *event)
4292b3b6 155{
fd1d0cb9 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);
4292b3b6 163 fNContributors = vertex->GetNContributors();
fd1d0cb9 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//_____________________________________________________________________________
182void 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//_____________________________________________________________________________
221void 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//_____________________________________________________________________________
232Bool_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//_____________________________________________________________________________
247void 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//_____________________________________________________________________________
259void 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//_____________________________________________________________________________
272void 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//_____________________________________________________________________________
288void 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//_____________________________________________________________________________
317void 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//_____________________________________________________________________________
343void 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//_____________________________________________________________________________
371void 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//_____________________________________________________________________________
381void 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//_____________________________________________________________________________
391void 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//_____________________________________________________________________________
409void 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//_____________________________________________________________________________
431void 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//_____________________________________________________________________________
452void 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//_____________________________________________________________________________
464void 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);
4292b3b6 471 return;
472}