AliPhysics  9fe175b (9fe175b)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalLight.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2016, 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 #include <RVersion.h>
16 #include <TClonesArray.h>
17 #include <TList.h>
18 #include <TObject.h>
19 #include <TH1F.h>
20 #include <TProfile.h>
21 #include <TSystem.h>
22 #include <TFile.h>
23 #include <TChain.h>
24 #include <TKey.h>
25 
26 #include "AliStack.h"
27 #include "AliAODEvent.h"
28 #include "AliAnalysisManager.h"
29 #include "AliCentrality.h"
30 #include "AliEMCALGeometry.h"
31 #include "AliESDEvent.h"
32 #include "AliEmcalParticle.h"
33 #include "AliEventplane.h"
34 #include "AliInputEventHandler.h"
35 #include "AliLog.h"
36 #include "AliMCParticle.h"
37 #include "AliVCluster.h"
38 #include "AliVEventHandler.h"
39 #include "AliVParticle.h"
40 #include "AliAODTrack.h"
41 #include "AliVCaloTrigger.h"
42 #include "AliGenPythiaEventHeader.h"
43 #include "AliAODMCHeader.h"
44 #include "AliMCEvent.h"
45 #include "AliEMCALTriggerPatchInfo.h"
46 
47 #include "AliMultSelection.h"
48 
50 
52 
56 
61  AliAnalysisTaskSE(),
62  fForceBeamType(kNA),
63  fGeneralHistograms(kFALSE),
64  fCreateHisto(kTRUE),
65  fNeedEmcalGeom(kTRUE),
66  fNcentBins(4),
67  fUseNewCentralityEstimation(kFALSE),
68  fIsPythia(kFALSE),
69  fCaloCellsName(),
70  fCaloTriggersName(),
71  fCaloTriggerPatchInfoName(),
72  fCentEst("V0M"),
73  fParticleCollArray(),
74  fClusterCollArray(),
75  fTriggerSelectionBitMap(0),
76  fMinCent(-999),
77  fMaxCent(-999),
78  fMinVz(-999),
79  fMaxVz(999),
80  fZvertexDiff(0.5),
81  fMinPtTrack(0),
82  fMinNTrack(0),
83  fMinPtTrackInEmcal(0),
84  fSelectPtHardBin(-999),
85  fAcceptedTriggerClasses(),
86  fRejectedTriggerClasses(),
87  fInitialized(kFALSE),
88  fDataType(kAOD),
89  fGeom(0),
90  fCaloCells(0),
91  fCaloTriggers(0),
92  fTriggerPatchInfo(0),
93  fCent(0),
94  fCentBin(-1),
95  fEPV0(-1.0),
96  fEPV0A(-1.0),
97  fEPV0C(-1.0),
98  fNVertCont(0),
99  fNVertSPDCont(0),
100  fFiredTriggerBitMap(0),
101  fFiredTriggerClasses(),
102  fBeamType(kNA),
103  fPythiaHeader(0),
104  fPtHard(0),
105  fPtHardBin(0),
106  fNTrials(0),
107  fXsection(0),
108  fOutput(0),
109  fHistEventCount(0),
110  fHistTrialsAfterSel(0),
111  fHistEventsAfterSel(0),
112  fHistXsectionAfterSel(0),
113  fHistTrials(0),
114  fHistEvents(0),
115  fHistXsection(0),
116  fHistPtHard(0),
117  fHistCentrality(0),
118  fHistZVertex(0),
119  fHistEventPlane(0),
120  fHistEventRejection(0),
121  fHistTriggerClasses(0)
122 {
123  fVertex[0] = 0;
124  fVertex[1] = 0;
125  fVertex[2] = 0;
126  fVertexSPD[0] = 0;
127  fVertexSPD[1] = 0;
128  fVertexSPD[2] = 0;
129 
130  fParticleCollArray.SetOwner(kTRUE);
131  fClusterCollArray.SetOwner(kTRUE);
132 }
133 
145  AliAnalysisTaskSE(name),
146  fForceBeamType(kNA),
147  fGeneralHistograms(kFALSE),
148  fCreateHisto(kTRUE),
149  fNeedEmcalGeom(kTRUE),
150  fNcentBins(4),
151  fUseNewCentralityEstimation(kFALSE),
152  fIsPythia(kFALSE),
153  fCaloCellsName(),
154  fCaloTriggersName(),
155  fCaloTriggerPatchInfoName(),
156  fCentEst("V0M"),
157  fParticleCollArray(),
158  fClusterCollArray(),
159  fTriggerSelectionBitMap(0),
160  fMinCent(-999),
161  fMaxCent(-999),
162  fMinVz(-999),
163  fMaxVz(999),
164  fZvertexDiff(0.5),
165  fMinPtTrack(0),
166  fMinNTrack(0),
167  fMinPtTrackInEmcal(0),
168  fSelectPtHardBin(-999),
169  fAcceptedTriggerClasses(),
170  fRejectedTriggerClasses(),
171  fInitialized(kFALSE),
172  fDataType(kAOD),
173  fGeom(0),
174  fCaloCells(0),
175  fCaloTriggers(0),
176  fTriggerPatchInfo(0),
177  fCent(0),
178  fCentBin(-1),
179  fEPV0(-1.0),
180  fEPV0A(-1.0),
181  fEPV0C(-1.0),
182  fNVertCont(0),
183  fNVertSPDCont(0),
184  fFiredTriggerBitMap(0),
185  fFiredTriggerClasses(),
186  fBeamType(kNA),
187  fPythiaHeader(0),
188  fPtHard(0),
189  fPtHardBin(0),
190  fNTrials(0),
191  fXsection(0),
192  fOutput(0),
193  fHistEventCount(0),
194  fHistTrialsAfterSel(0),
195  fHistEventsAfterSel(0),
196  fHistXsectionAfterSel(0),
197  fHistTrials(0),
198  fHistEvents(0),
199  fHistXsection(0),
200  fHistPtHard(0),
201  fHistCentrality(0),
202  fHistZVertex(0),
203  fHistEventPlane(0),
204  fHistEventRejection(0),
205  fHistTriggerClasses(0)
206 {
207  fVertex[0] = 0;
208  fVertex[1] = 0;
209  fVertex[2] = 0;
210  fVertexSPD[0] = 0;
211  fVertexSPD[1] = 0;
212  fVertexSPD[2] = 0;
213  fParticleCollArray.SetOwner(kTRUE);
214  fClusterCollArray.SetOwner(kTRUE);
215 
216  if (fCreateHisto) {
217  DefineOutput(1, TList::Class());
218  }
219 }
220 
225 {
226 }
227 
248 {
249  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
250  if (mgr) {
251  AliVEventHandler *evhand = mgr->GetInputEventHandler();
252  if (evhand) {
253  if (evhand->InheritsFrom("AliESDInputHandler")) {
254  fDataType = kESD;
255  }
256  else {
257  fDataType = kAOD;
258  }
259  }
260  else {
261  AliError("Event handler not found!");
262  }
263  }
264  else {
265  AliError("Analysis manager not found!");
266  }
267 
268  if (!fCreateHisto)
269  return;
270 
271  OpenFile(1);
272  fOutput = new TList();
273  fOutput->SetOwner();
274 
275  if (fForceBeamType == kpp)
276  fNcentBins = 1;
277 
278  if (!fGeneralHistograms)
279  return;
280 
281  if (fIsPythia) {
282  fHistTrialsAfterSel = new TH1F("fHistTrialsAfterSel", "fHistTrialsAfterSel", 11, 0, 11);
283  fHistTrialsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
284  fHistTrialsAfterSel->GetYaxis()->SetTitle("trials");
286 
287  fHistEventsAfterSel = new TH1F("fHistEventsAfterSel", "fHistEventsAfterSel", 11, 0, 11);
288  fHistEventsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
289  fHistEventsAfterSel->GetYaxis()->SetTitle("total events");
291 
292  fHistXsectionAfterSel = new TProfile("fHistXsectionAfterSel", "fHistXsectionAfterSel", 11, 0, 11);
293  fHistXsectionAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
294  fHistXsectionAfterSel->GetYaxis()->SetTitle("xsection");
296 
297  fHistTrials = new TH1F("fHistTrials", "fHistTrials", 11, 0, 11);
298  fHistTrials->GetXaxis()->SetTitle("p_{T} hard bin");
299  fHistTrials->GetYaxis()->SetTitle("trials");
300  fOutput->Add(fHistTrials);
301 
302  fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
303  fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
304  fHistEvents->GetYaxis()->SetTitle("total events");
305  fOutput->Add(fHistEvents);
306 
307  fHistXsection = new TProfile("fHistXsection", "fHistXsection", 11, 0, 11);
308  fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
309  fHistXsection->GetYaxis()->SetTitle("xsection");
310  fOutput->Add(fHistXsection);
311 
312  const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
313  const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
314 
315  for (Int_t i = 1; i < 12; i++) {
316  fHistTrialsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
317  fHistEventsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
318 
319  fHistTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
320  fHistXsection->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
321  fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
322  }
323 
324  fHistPtHard = new TH1F("fHistPtHard", "fHistPtHard", 250, 0, 1000);
325  fHistPtHard->GetXaxis()->SetTitle("p_{T,hard} (GeV/c)");
326  fHistPtHard->GetYaxis()->SetTitle("counts");
327  fOutput->Add(fHistPtHard);
328  }
329 
330  fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
331  fHistZVertex->GetXaxis()->SetTitle("z");
332  fHistZVertex->GetYaxis()->SetTitle("counts");
333  fOutput->Add(fHistZVertex);
334 
335  if (fForceBeamType != kpp) {
336  fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 200, 0, 100);
337  fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
338  fHistCentrality->GetYaxis()->SetTitle("counts");
339  fOutput->Add(fHistCentrality);
340 
341  fHistEventPlane = new TH1F("fHistEventPlane","Event plane", 120, -TMath::Pi(), TMath::Pi());
342  fHistEventPlane->GetXaxis()->SetTitle("event plane");
343  fHistEventPlane->GetYaxis()->SetTitle("counts");
344  fOutput->Add(fHistEventPlane);
345  }
346 
347  fHistEventRejection = new TH1F("fHistEventRejection","Reasons to reject event",20,0,20);
348 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
349  fHistEventRejection->SetBit(TH1::kCanRebin);
350 #else
351  fHistEventRejection->SetCanExtend(TH1::kAllAxes);
352 #endif
353  fHistEventRejection->GetXaxis()->SetBinLabel(1,"PhysSel");
354  fHistEventRejection->GetXaxis()->SetBinLabel(2,"trigger");
355  fHistEventRejection->GetXaxis()->SetBinLabel(3,"trigTypeSel");
356  fHistEventRejection->GetXaxis()->SetBinLabel(4,"Cent");
357  fHistEventRejection->GetXaxis()->SetBinLabel(5,"vertex contr.");
358  fHistEventRejection->GetXaxis()->SetBinLabel(6,"Vz");
359  fHistEventRejection->GetXaxis()->SetBinLabel(7,"VzSPD");
360  fHistEventRejection->GetXaxis()->SetBinLabel(8,"trackInEmcal");
361  fHistEventRejection->GetXaxis()->SetBinLabel(9,"minNTrack");
362  fHistEventRejection->GetXaxis()->SetBinLabel(10,"VtxSel2013pA");
363  fHistEventRejection->GetXaxis()->SetBinLabel(11,"PileUp");
364  fHistEventRejection->GetXaxis()->SetBinLabel(12,"EvtPlane");
365  fHistEventRejection->GetXaxis()->SetBinLabel(13,"SelPtHardBin");
366  fHistEventRejection->GetXaxis()->SetBinLabel(14,"Bkg evt");
367  fHistEventRejection->GetYaxis()->SetTitle("counts");
369 
370  fHistTriggerClasses = new TH1F("fHistTriggerClasses","fHistTriggerClasses",3,0,3);
371 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
372  fHistTriggerClasses->SetBit(TH1::kCanRebin);
373 #else
374  fHistTriggerClasses->SetCanExtend(TH1::kAllAxes);
375 #endif
377 
378  fHistEventCount = new TH1F("fHistEventCount","fHistEventCount",2,0,2);
379  fHistEventCount->GetXaxis()->SetBinLabel(1,"Accepted");
380  fHistEventCount->GetXaxis()->SetBinLabel(2,"Rejected");
381  fHistEventCount->GetYaxis()->SetTitle("counts");
382  fOutput->Add(fHistEventCount);
383 
384  PostData(1, fOutput);
385 }
386 
401 {
402  if (fIsPythia) {
406  fHistPtHard->Fill(fPtHard);
407  }
408 
409  fHistZVertex->Fill(fVertex[2]);
410 
411  if (fForceBeamType != kpp) {
412  fHistCentrality->Fill(fCent);
413  fHistEventPlane->Fill(fEPV0);
414  }
415 
416  TObjArray* triggerClasses = InputEvent()->GetFiredTriggerClasses().Tokenize(" ");
417  TIter next(triggerClasses);
418  TObjString* triggerClass = 0;
419  while ((triggerClass = static_cast<TObjString*>(next()))) {
420  fHistTriggerClasses->Fill(triggerClass->GetString(), 1);
421  }
422  delete triggerClasses;
423  triggerClasses = 0;
424 
425  return kTRUE;
426 }
427 
448 {
449  if (!fInitialized)
450  ExecOnce();
451 
452  if (!fInitialized)
453  return;
454 
455  if (!RetrieveEventObjects())
456  return;
457 
458  if (IsEventSelected()) {
459  if (fGeneralHistograms) fHistEventCount->Fill("Accepted",1);
460  }
461  else {
462  if (fGeneralHistograms) fHistEventCount->Fill("Rejected",1);
463  return;
464  }
465 
467  if (!FillGeneralHistograms())
468  return;
469  }
470 
471  if (!Run())
472  return;
473 
474  if (fCreateHisto) {
475  if (!FillHistograms())
476  return;
477  }
478 
479  if (fCreateHisto && fOutput) {
480  // information for this iteration of the UserExec in the container
481  PostData(1, fOutput);
482  }
483 }
484 
496 Bool_t AliAnalysisTaskEmcalLight::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
497 {
498 
499  TString file(currFile);
500  fXsec = 0;
501  fTrials = 1;
502 
503  if (file.Contains(".zip#")) {
504  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
505  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
506  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
507  file.Replace(pos+1,pos2-pos1,"");
508  } else {
509  // not an archive take the basename....
510  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
511  }
512  AliDebug(1,Form("File name: %s",file.Data()));
513 
514  // Get the pt hard bin
515  TString strPthard(file);
516 
517  strPthard.Remove(strPthard.Last('/'));
518  strPthard.Remove(strPthard.Last('/'));
519  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
520  strPthard.Remove(0,strPthard.Last('/')+1);
521  if (strPthard.IsDec())
522  pthard = strPthard.Atoi();
523  else
524  AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
525 
526  // problem that we cannot really test the existance of a file in a archive so we have to live with open error message from root
527  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
528 
529  if (!fxsec) {
530  // next trial fetch the histgram file
531  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
532  if (!fxsec) {
533  // not a severe condition but inciate that we have no information
534  return kFALSE;
535  } else {
536  // find the tlist we want to be independtent of the name so use the Tkey
537  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
538  if (!key) {
539  fxsec->Close();
540  return kFALSE;
541  }
542  TList *list = dynamic_cast<TList*>(key->ReadObj());
543  if (!list) {
544  fxsec->Close();
545  return kFALSE;
546  }
547  fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
548  fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
549  fxsec->Close();
550  }
551  } else { // no tree pyxsec.root
552  TTree *xtree = (TTree*)fxsec->Get("Xsection");
553  if (!xtree) {
554  fxsec->Close();
555  return kFALSE;
556  }
557  UInt_t ntrials = 0;
558  Double_t xsection = 0;
559  xtree->SetBranchAddress("xsection",&xsection);
560  xtree->SetBranchAddress("ntrials",&ntrials);
561  xtree->GetEntry(0);
562  fTrials = ntrials;
563  fXsec = xsection;
564  fxsec->Close();
565  }
566  return kTRUE;
567 }
568 
583 {
585  return kTRUE;
586 
587  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
588  if (!tree) {
589  AliError(Form("%s - UserNotify: No current tree!",GetName()));
590  return kFALSE;
591  }
592 
593  Float_t xsection = 0;
594  Float_t trials = 0;
595  Int_t pthardbin = 0;
596 
597  TFile *curfile = tree->GetCurrentFile();
598  if (!curfile) {
599  AliError(Form("%s - UserNotify: No current file!",GetName()));
600  return kFALSE;
601  }
602 
603  TChain *chain = dynamic_cast<TChain*>(tree);
604  if (chain) tree = chain->GetTree();
605 
606  Int_t nevents = tree->GetEntriesFast();
607 
608  Bool_t res = PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
609 
610  if (!res) return kTRUE;
611 
612  fHistTrials->Fill(pthardbin, trials);
613  fHistXsection->Fill(pthardbin, xsection);
614  fHistEvents->Fill(pthardbin, nevents);
615 
616  return kTRUE;
617 }
618 
630 {
631  if (!InputEvent()) {
632  AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
633  return;
634  }
635 
636  if (fNeedEmcalGeom) {
637  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(InputEvent()->GetRunNumber());
638  if (!fGeom) {
639  AliFatal(Form("%s: Can not get EMCal geometry instance. If you do not need the EMCal geometry, disable it by setting task->SetNeedEmcalGeometry(kFALSE).", GetName()));
640  return;
641  }
642  }
643 
644  //Load all requested track branches - each container knows name already
645  for (Int_t i =0; i<fParticleCollArray.GetEntriesFast(); i++) {
646  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
647  cont->SetArray(InputEvent());
648  }
649 
650  //Load all requested cluster branches - each container knows name already
651  for (Int_t i =0; i<fClusterCollArray.GetEntriesFast(); i++) {
652  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
653  cont->SetArray(InputEvent());
654  }
655 
656  if (!fCaloCellsName.IsNull() && !fCaloCells) {
657  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
658  if (!fCaloCells) {
659  AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
660  return;
661  }
662  }
663 
664  if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
665  fCaloTriggers = dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
666  if (!fCaloTriggers) {
667  AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data()));
668  return;
669  }
670  }
671 
672  if (!fCaloTriggerPatchInfoName.IsNull() && !fTriggerPatchInfo) {
673  fTriggerPatchInfo = GetArrayFromEvent(fCaloTriggerPatchInfoName.Data(),"AliEMCALTriggerPatchInfo");
674  if (!fTriggerPatchInfo) {
675  AliError(Form("%s: Could not retrieve calo trigger patch info %s!", GetName(), fCaloTriggerPatchInfoName.Data()));
676  return;
677  }
678 
679  }
680 
681  fInitialized = kTRUE;
682 }
683 
690 {
691  if (fForceBeamType != kNA)
692  return fForceBeamType;
693 
694  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
695  if (esd) {
696  const AliESDRun *run = esd->GetESDRun();
697  TString beamType = run->GetBeamType();
698  if (beamType == "p-p")
699  return kpp;
700  else if (beamType == "A-A")
701  return kAA;
702  else if (beamType == "p-A")
703  return kpA;
704  else
705  return kNA;
706  } else {
707  Int_t runNumber = InputEvent()->GetRunNumber();
708  if ((runNumber >= 136851 && runNumber <= 139517) || // LHC10h
709  (runNumber >= 166529 && runNumber <= 170593)) { // LHC11h
710  return kAA;
711  } else if ((runNumber>=188365 && runNumber <= 188366) || // LHC12g
712  (runNumber >= 195344 && runNumber <= 196608)) { // LHC13b-f
713  return kpA;
714  } else {
715  return kpp;
716  }
717  }
718 }
719 
742 {
744  if (fGeneralHistograms) fHistEventRejection->Fill("PhysSel",1);
745  return kFALSE;
746  }
747 
748  Bool_t acceptedTrgClassFound = kFALSE;
749  if (fAcceptedTriggerClasses.GetEntriesFast() > 0) {
750  TIter acceptedTrigger(&fAcceptedTriggerClasses);
751  TObject* obj = 0;
752  while ((obj = acceptedTrigger())) {
753  if (fFiredTriggerClasses.Contains(obj->GetName())) {
754  acceptedTrgClassFound = kTRUE;
755  break;
756  }
757  }
758 
759  if (!acceptedTrgClassFound) {
760  if (fGeneralHistograms) fHistEventRejection->Fill("Trg class (acc)",1);
761  return kFALSE;
762  }
763  }
764 
765  if (fRejectedTriggerClasses.GetEntriesFast() > 0) {
766  TIter rejectedTrigger(&fRejectedTriggerClasses);
767  TObject* obj = 0;
768  while ((obj = rejectedTrigger())) {
769  if (fFiredTriggerClasses.Contains(obj->GetName())) {
770  if (fGeneralHistograms) fHistEventRejection->Fill("Trg class (rej)",1);
771  return kFALSE;
772  }
773  }
774  }
775 
776  if ((fMinCent != -999) && (fMaxCent != -999)) {
777  if (fCent < fMinCent || fCent > fMaxCent) {
778  if (fGeneralHistograms) fHistEventRejection->Fill("Cent",1);
779  return kFALSE;
780  }
781  }
782 
783  if ((fMinVz > -998.) && (fMaxVz < 998.)) {
784  if (fNVertCont == 0 ) {
785  if (fGeneralHistograms) fHistEventRejection->Fill("vertex contr.",1);
786  return kFALSE;
787  }
788  Double_t vz = fVertex[2];
789  if (vz < fMinVz || vz > fMaxVz) {
790  if (fGeneralHistograms) fHistEventRejection->Fill("Vz",1);
791  return kFALSE;
792  }
793 
794  if (fNVertSPDCont > 0 && fZvertexDiff < 999) {
795  Double_t vzSPD = fVertexSPD[2];
796  Double_t dvertex = TMath::Abs(vz-vzSPD);
797  //if difference larger than fZvertexDiff
798  if (dvertex > fZvertexDiff) {
799  if (fGeneralHistograms) fHistEventRejection->Fill("VzSPD",1);
800  return kFALSE;
801  }
802  }
803  }
804 
805  if (fMinPtTrackInEmcal > 0 && fGeom && GetParticleContainer(0)) {
806  Bool_t trackInEmcalOk = kFALSE;
807  Int_t ntracks = GetParticleContainer(0)->GetNParticles();
808  for (Int_t i = 0; i < ntracks; i++) {
809  AliVParticle *track = GetParticleContainer(0)->GetAcceptParticle(i);
810  if (!track)
811  continue;
812 
813  Double_t phiMin = fGeom->GetArm1PhiMin() * TMath::DegToRad();
814  Double_t phiMax = fGeom->GetArm1PhiMax() * TMath::DegToRad();
815  Int_t runNumber = InputEvent()->GetRunNumber();
816  if (runNumber>=177295 && runNumber<=197470) { //small SM masked in 2012 and 2013
817  phiMin = 1.4;
818  phiMax = TMath::Pi();
819  }
820 
821  if (track->Eta() < fGeom->GetArm1EtaMin() || track->Eta() > fGeom->GetArm1EtaMax() || track->Phi() < phiMin || track->Phi() > phiMax)
822  continue;
823  if (track->Pt() > fMinPtTrackInEmcal) {
824  trackInEmcalOk = kTRUE;
825  break;
826  }
827  }
828  if (!trackInEmcalOk) {
829  if (fGeneralHistograms) fHistEventRejection->Fill("trackInEmcal",1);
830  return kFALSE;
831  }
832  }
833 
834  if (fMinNTrack > 0) {
835  Int_t nTracksAcc = 0;
836  Int_t ntracks = GetParticleContainer(0)->GetNParticles();
837  for (Int_t i = 0; i < ntracks; i++) {
838  AliVParticle *track = GetParticleContainer(0)->GetAcceptParticle(i);
839  if (!track)
840  continue;
841  if (track->Pt() > fMinPtTrack) {
842  nTracksAcc++;
843  if (nTracksAcc >= fMinNTrack)
844  break;
845  }
846  }
847  if (nTracksAcc<fMinNTrack) {
848  if (fGeneralHistograms) fHistEventRejection->Fill("minNTrack",1);
849  return kFALSE;
850  }
851  }
852 
853  if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin) {
854  if (fGeneralHistograms) fHistEventRejection->Fill("SelPtHardBin",1);
855  return kFALSE;
856  }
857 
858  return kTRUE;
859 }
860 
869 TClonesArray *AliAnalysisTaskEmcalLight::GetArrayFromEvent(const char *name, const char *clname)
870 {
871  TClonesArray *arr = 0;
872  TString sname(name);
873  if (!sname.IsNull()) {
874  arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
875  if (!arr) {
876  AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
877  return 0;
878  }
879  } else {
880  return 0;
881  }
882 
883  if (!clname)
884  return arr;
885 
886  TString objname(arr->GetClass()->GetName());
887  TClass cls(objname);
888  if (!cls.InheritsFrom(clname)) {
889  AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
890  GetName(), cls.GetName(), name, clname));
891  return 0;
892  }
893  return arr;
894 }
895 
901 {
902  fVertex[0] = 0;
903  fVertex[1] = 0;
904  fVertex[2] = 0;
905  fNVertCont = 0;
906 
907  fVertexSPD[0] = 0;
908  fVertexSPD[1] = 0;
909  fVertexSPD[2] = 0;
910  fNVertSPDCont = 0;
911 
912  fFiredTriggerClasses = InputEvent()->GetFiredTriggerClasses();
913 
914  if (fDataType == kESD) {
915  fFiredTriggerBitMap = static_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())->IsEventSelected();
916  }
917  else {
918  fFiredTriggerBitMap = static_cast<AliVAODHeader*>(InputEvent()->GetHeader())->GetOfflineTrigger();
919  }
920 
921  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
922  if (vert) {
923  vert->GetXYZ(fVertex);
924  fNVertCont = vert->GetNContributors();
925  }
926 
927  const AliVVertex *vertSPD = InputEvent()->GetPrimaryVertexSPD();
928  if (vertSPD) {
929  vertSPD->GetXYZ(fVertexSPD);
930  fNVertSPDCont = vertSPD->GetNContributors();
931  }
932 
934 
935  if (fBeamType == kAA || fBeamType == kpA ) {
937  AliMultSelection *MultSelection = static_cast<AliMultSelection*>(InputEvent()->FindListObject("MultSelection"));
938  if (MultSelection) {
939  fCent = MultSelection->GetMultiplicityPercentile(fCentEst.Data());
940  }
941  else {
942  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
943  }
944  }
945  else { // old centrality estimation < 2015
946  AliCentrality *aliCent = InputEvent()->GetCentrality();
947  if (aliCent) {
948  fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
949  }
950  else {
951  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
952  }
953  }
954 
955  if (fNcentBins==4) {
956  if (fCent >= 0 && fCent < 10) fCentBin = 0;
957  else if (fCent >= 10 && fCent < 30) fCentBin = 1;
958  else if (fCent >= 30 && fCent < 50) fCentBin = 2;
959  else if (fCent >= 50 && fCent <= 100) fCentBin = 3;
960  else {
961  AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
962  fCentBin = fNcentBins-1;
963  }
964  }
965  else if (fNcentBins==5) { // for PbPb 2015
966  if (fCent >= 0 && fCent < 10) fCentBin = 0;
967  else if (fCent >= 10 && fCent < 30) fCentBin = 1;
968  else if (fCent >= 30 && fCent < 50) fCentBin = 2;
969  else if (fCent >= 50 && fCent <= 90) fCentBin = 3;
970  else if (fCent > 90) {
971  fCent = 99;
972  fCentBin = 4;
973  }
974  else {
975  AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
976  fCentBin = fNcentBins-1;
977  }
978  }
979  else {
980  Double_t centWidth = (fMaxCent-fMinCent)/(Double_t)fNcentBins;
981  if(centWidth>0.) {
982  fCentBin = TMath::FloorNint(fCent/centWidth);
983  }
984  else {
985  fCentBin = 0;
986  }
987  if (fCentBin>=fNcentBins) {
988  AliWarning(Form("%s: fCentBin too large: cent = %f fCentBin = %d. Assuming 99", GetName(),fCent,fCentBin));
989  fCentBin = fNcentBins-1;
990  }
991  }
992 
993  AliEventplane *aliEP = InputEvent()->GetEventplane();
994  if (aliEP) {
995  fEPV0 = aliEP->GetEventplane("V0" ,InputEvent());
996  fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
997  fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
998  } else {
999  AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
1000  }
1001  }
1002  else {
1003  fCent = 99;
1004  fCentBin = 0;
1005  }
1006 
1007  if (fIsPythia) {
1008 
1009  if (MCEvent()) {
1010  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
1011  if (!fPythiaHeader) {
1012  // Check if AOD
1013  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
1014 
1015  if (aodMCH) {
1016  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
1017  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
1018  if (fPythiaHeader) break;
1019  }
1020  }
1021  }
1022  }
1023 
1024  if (fPythiaHeader) {
1025  fPtHard = fPythiaHeader->GetPtHard();
1026 
1027  const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
1028  const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
1029  for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
1030  if (fPtHard >= ptHardLo[fPtHardBin] && fPtHard < ptHardHi[fPtHardBin])
1031  break;
1032  }
1033 
1034  fXsection = fPythiaHeader->GetXsection();
1035  fNTrials = fPythiaHeader->Trials();
1036  }
1037  }
1038 
1039  AliEmcalContainer* cont = 0;
1040 
1041  TIter nextPartColl(&fParticleCollArray);
1042  while ((cont = static_cast<AliEmcalContainer*>(nextPartColl()))) cont->NextEvent();
1043 
1044  TIter nextClusColl(&fClusterCollArray);
1045  while ((cont = static_cast<AliParticleContainer*>(nextClusColl()))) cont->NextEvent();
1046 
1047  return kTRUE;
1048 }
1049 
1058 {
1059  if (TString(n).IsNull()) return 0;
1060 
1062 
1063  fParticleCollArray.Add(cont);
1064 
1065  return cont;
1066 }
1067 
1076 {
1077  if (TString(n).IsNull()) return 0;
1078 
1079  AliTrackContainer* cont = new AliTrackContainer(n);
1080 
1081  fParticleCollArray.Add(cont);
1082 
1083  return cont;
1084 }
1085 
1094 {
1095  if (TString(n).IsNull()) return 0;
1096 
1098 
1099  fParticleCollArray.Add(cont);
1100 
1101  return cont;
1102 }
1103 
1112 {
1113  if (TString(n).IsNull()) return 0;
1114 
1116 
1117  fClusterCollArray.Add(cont);
1118 
1119  return cont;
1120 }
1121 
1128 {
1129  if (i<0 || i>fParticleCollArray.GetEntriesFast()) return 0;
1130  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
1131  return cont;
1132 }
1133 
1140 {
1141  if (i<0 || i>fClusterCollArray.GetEntriesFast()) return 0;
1142  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
1143  return cont;
1144 }
1145 
1152 {
1153  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.FindObject(name));
1154  return cont;
1155 }
1156 
1163 {
1164  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.FindObject(name));
1165  return cont;
1166 }
1167 
1173 void AliAnalysisTaskEmcalLight::AddObjectToEvent(TObject *obj, Bool_t attempt)
1174 {
1175  if (!(InputEvent()->FindListObject(obj->GetName()))) {
1176  InputEvent()->AddObject(obj);
1177  }
1178  else {
1179  if (!attempt) {
1180  AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1181  }
1182  }
1183 }
1184 
1192 Bool_t AliAnalysisTaskEmcalLight::IsTrackInEmcalAcceptance(AliVParticle* part, Double_t edges) const
1193 {
1194 
1195  if (!fGeom) {
1196  AliWarning(Form("%s - AliAnalysisTaskEmcalBase::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
1197  return kFALSE;
1198  }
1199 
1200  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
1201  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
1202 
1203  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
1204  return kTRUE;
1205  }
1206  else {
1207  return kFALSE;
1208  }
1209 }
1210 
1212 {
1213  axis->SetBinLabel(1, "NullObject");
1214  axis->SetBinLabel(2, "Pt");
1215  axis->SetBinLabel(3, "Acceptance");
1216  axis->SetBinLabel(4, "MCLabel");
1217  axis->SetBinLabel(5, "BitMap");
1218  axis->SetBinLabel(6, "HF cut");
1219  axis->SetBinLabel(7, "Bit6");
1220  axis->SetBinLabel(8, "NotHybridTrack");
1221  axis->SetBinLabel(9, "MCFlag");
1222  axis->SetBinLabel(10, "MCGenerator");
1223  axis->SetBinLabel(11, "ChargeCut");
1224  axis->SetBinLabel(12, "MinDistanceTPCSectorEdge");
1225  axis->SetBinLabel(13, "Bit12");
1226  axis->SetBinLabel(14, "IsEMCal");
1227  axis->SetBinLabel(15, "Time");
1228  axis->SetBinLabel(16, "Energy");
1229  axis->SetBinLabel(17, "ExoticCut");
1230  axis->SetBinLabel(18, "Bit17");
1231  axis->SetBinLabel(19, "Area");
1232  axis->SetBinLabel(20, "AreaEmc");
1233  axis->SetBinLabel(21, "ZLeadingCh");
1234  axis->SetBinLabel(22, "ZLeadingEmc");
1235  axis->SetBinLabel(23, "NEF");
1236  axis->SetBinLabel(24, "MinLeadPt");
1237  axis->SetBinLabel(25, "MaxTrackPt");
1238  axis->SetBinLabel(26, "MaxClusterPt");
1239  axis->SetBinLabel(27, "Flavour");
1240  axis->SetBinLabel(28, "TagStatus");
1241  axis->SetBinLabel(29, "MinNConstituents");
1242  axis->SetBinLabel(30, "Bit29");
1243  axis->SetBinLabel(31, "Bit30");
1244  axis->SetBinLabel(32, "Bit31");
1245 }
1246 
1253 Double_t AliAnalysisTaskEmcalLight::GetParallelFraction(AliVParticle* part1, AliVParticle* part2)
1254 {
1255  TVector3 vect1(part1->Px(), part1->Py(), part1->Pz());
1256  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1257  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1258  return z;
1259 }
1260 
1267 Double_t AliAnalysisTaskEmcalLight::GetParallelFraction(const TVector3& vect1, AliVParticle* part2)
1268 {
1269  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1270  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1271  return z;
1272 }
1273 
1282 void AliAnalysisTaskEmcalLight::GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
1283 {
1284  phidiff = 999;
1285  etadiff = 999;
1286 
1287  if (!t||!v) return;
1288 
1289  Double_t veta = t->GetTrackEtaOnEMCal();
1290  Double_t vphi = t->GetTrackPhiOnEMCal();
1291 
1292  Float_t pos[3] = {0};
1293  v->GetPosition(pos);
1294  TVector3 cpos(pos);
1295  Double_t ceta = cpos.Eta();
1296  Double_t cphi = cpos.Phi();
1297  etadiff=veta-ceta;
1298  phidiff=TVector2::Phi_mpi_pi(vphi-cphi);
1299 }
1300 
1307 {
1308  Byte_t ret = 0;
1309  if (t->TestBit(BIT(22)) && !t->TestBit(BIT(23)))
1310  ret = 1;
1311  else if (!t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1312  ret = 2;
1313  else if (t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1314  ret = 3;
1315  return ret;
1316 }
1317 
1327 Byte_t AliAnalysisTaskEmcalLight::GetTrackType(const AliAODTrack *aodTrack, UInt_t filterBit1, UInt_t filterBit2)
1328 {
1329 
1330  Int_t res = 0;
1331 
1332  if (aodTrack->TestFilterBit(filterBit1)) {
1333  res = 0;
1334  }
1335  else if (aodTrack->TestFilterBit(filterBit2)) {
1336  if ((aodTrack->GetStatus()&AliVTrack::kITSrefit)!=0) {
1337  res = 1;
1338  }
1339  else {
1340  res = 2;
1341  }
1342  }
1343  else {
1344  res = 3;
1345  }
1346 
1347  return res;
1348 }
Double_t fVertexSPD[3]
!event Svertex
TString fCaloTriggersName
name of calo triggers collection
TObjArray fAcceptedTriggerClasses
list of accepted trigger classes
Int_t fNcentBins
how many centrality bins
AliEMCALGeometry * fGeom
!emcal geometry
EBeamType_t fBeamType
!event beam type
TClonesArray * GetArrayFromEvent(const char *name, const char *clname=0)
TObjArray fClusterCollArray
cluster collection array
Float_t fXsection
!x-section from pythia header
Double_t fEPV0A
!event plane V0A
Double_t fEPV0
!event plane V0
TH1 * fHistEventCount
!incoming and selected events
TObjArray fRejectedTriggerClasses
list of accepted trigger classes
TObjArray fParticleCollArray
particle/track collection array
Int_t GetNParticles() const
Container with name, TClonesArray and cuts for particles.
Double_t fZvertexDiff
upper limit for distance between primary and SPD vertex
TSystem * gSystem
AliTrackContainer * AddTrackContainer(const char *n)
TString fCaloTriggerPatchInfoName
trigger patch info array name
Double_t fMinVz
min vertex for event selection
TList * list
TProfile * fHistXsectionAfterSel
!x section from pythia header
TString fCentEst
name of the centrality estimator
UInt_t fTriggerSelectionBitMap
trigger selection bit map
TH1 * fHistCentrality
!event centrality distribution
static Byte_t GetTrackType(const AliVTrack *t)
TH1 * fHistTriggerClasses
!number of events in each trigger class
TString fCaloCellsName
name of calo cell collection
TH1 * fHistEvents
!total number of events per pt hard bin
Bool_t IsTrackInEmcalAcceptance(AliVParticle *part, Double_t edges=0.9) const
Double_t fMinCent
min centrality for event selection
Double_t fMinPtTrack
cut on track pt in event selection
Container for particles within the EMCAL framework.
static void GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
void AddObjectToEvent(TObject *obj, Bool_t attempt=kFALSE)
static Double_t fgkEMCalDCalPhiDivide
phi value used to distinguish between DCal and EMCal
EBeamType_t fForceBeamType
forced beam type
AliClusterContainer * AddClusterContainer(const char *n)
Bool_t fCreateHisto
whether or not create histograms
TH1 * fHistEventRejection
!book keep reasons for rejecting event
TH1 * fHistTrials
!trials from pyxsec.root
TString fFiredTriggerClasses
!trigger classes fired by the current event
Base task in the EMCAL framework (lighter version of AliAnalysisTaskEmcal)
Double_t fVertex[3]
!event vertex
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Int_t fPtHardBin
!event pt hard bin
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
Bool_t fIsPythia
if it is a PYTHIA production
Bool_t fInitialized
!whether or not the task has been already initialized
TProfile * fHistXsection
!x section from pyxsec.root
Int_t fCentBin
!event centrality bin
Int_t fSelectPtHardBin
select one pt hard bin for analysis
TH1 * fHistEventsAfterSel
!total number of events per pt hard bin after selection
Double_t fMaxCent
max centrality for event selection
ULong_t fFiredTriggerBitMap
!bit map of fired triggers
Int_t fNVertSPDCont
!event SPD vertex number of contributors
Double_t fMinPtTrackInEmcal
min pt track in emcal
Bool_t PythiaInfoFromFile(const char *currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
Int_t fMinNTrack
minimum nr of tracks in event with pT>fTrackPtCut
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
AliParticleContainer * GetParticleContainer(Int_t i=0) const
TFile * file
Int_t nevents[nsamples]
Int_t fNVertCont
!event vertex number of contributors
AliMCParticleContainer * AddMCParticleContainer(const char *n)
Bool_t fGeneralHistograms
whether or not it should fill some general histograms
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
AliGenPythiaEventHeader * fPythiaHeader
!event Pythia header
TH1 * fHistPtHard
!pt hard distribution
TH1 * fHistTrialsAfterSel
!total number of trials per pt hard bin after selection
TH1 * fHistEventPlane
!event plane distribution
EBeamType_t
Switch for the beam type.
AliVCaloTrigger * fCaloTriggers
!calo triggers
Container structure for EMCAL clusters.
Container for MC-true particles within the EMCAL framework.
TH1 * fHistZVertex
!z vertex position
AliParticleContainer * AddParticleContainer(const char *n)
Double_t fMaxVz
max vertex for event selection
Double_t fCent
!event centrality
Double_t fEPV0C
!event plane V0C
Bool_t fUseNewCentralityEstimation
Use new centrality estimation (for 2015 data)
TClonesArray * fTriggerPatchInfo
!trigger patch info array
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
static Double_t GetParallelFraction(AliVParticle *part1, AliVParticle *part2)
EDataType_t fDataType
!data type (ESD or AOD)