Update of the Xiaoming code for pp900 first muon analysis. Fixing wanirngs and violti...
[u/mrichter/AliRoot.git] / 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 }