AliPhysics  d9e9949 (d9e9949)
AliAnalysisTaskEmcalJetTriggerMatcher.cxx
Go to the documentation of this file.
1 // $Id$
2 /**************************************************************************
3  * Copyright(c) 1998-2015, ALICE Experiment at CERN, All rights reserved. *
4  * *
5  * Author: The ALICE Off-line Project. *
6  * Contributors are mentioned in the code where appropriate. *
7  * *
8  * Permission to use, copy, modify and distribute this software and its *
9  * documentation strictly for non-commercial purposes is hereby granted *
10  * without fee, provided that the above copyright notice appears in all *
11  * copies and that both the copyright notice and this permission notice *
12  * appear in the supporting documentation. The authors make no claims *
13  * about the suitability of this software for any purpose. It is *
14  * provided "as is" without express or implied warranty. *
15  **************************************************************************/
16 
18 //
19 // Task to match cluster from a jet to GA patch.....
20 //
21 //
22 // andrew.john.castro@cern.ch
24 
25 
26 
28 
29 // general ROOT includes
30 #include <TCanvas.h>
31 #include <TChain.h>
32 #include <TClonesArray.h>
33 #include <TH1F.h>
34 #include <TH2F.h>
35 #include <TH3F.h>
36 #include <THnSparse.h>
37 #include <TList.h>
38 #include <TLorentzVector.h>
39 #include <TParameter.h>
40 #include <TParticle.h>
41 #include <TTree.h>
42 #include <TVector3.h>
43 #include <TObjArray.h>
44 
45 // AliROOT includes
46 #include "AliAODEvent.h"
47 #include "AliESDEvent.h"
48 #include "AliAnalysisManager.h"
49 #include "AliAnalysisTask.h"
50 #include "AliCentrality.h"
51 #include "AliEmcalJet.h"
52 #include "AliAODJet.h"
53 #include "AliVCluster.h"
54 #include "AliVTrack.h"
55 #include <AliVEvent.h>
56 #include <AliVParticle.h>
57 #include "AliRhoParameter.h"
58 #include "AliLog.h"
59 #include "AliParticleContainer.h"
60 #include "AliClusterContainer.h"
61 #include "AliEmcalParticle.h"
62 #include "AliESDCaloCluster.h"
63 #include "AliEMCALTriggerPatchInfo.h"
64 
65 // event handler (and pico's) includes
66 #include <AliInputEventHandler.h>
67 #include <AliVEventHandler.h>
68 #include "AliESDInputHandler.h"
69 #include "AliPicoTrack.h"
70 
71 using std::cout;
72 using std::endl;
73 
75 
76 //________________________________________________________________________
78  AliAnalysisTaskEmcalJet("heavyF",kFALSE),
79  event(0),
80  fEventTrigEMCALL1Gamma1(0),
81  fEventTrigEMCALL1Gamma2(0),
82  fInputEvent(0x0),
83  fPhimin(-10), fPhimax(10),
84  fEtamin(-0.9), fEtamax(0.9),
85  fAreacut(0.0),
86  fJetHIpt(20.0),
87  fTrackEta(0.9),
88  fFillHists(0),
89  fMatchDist(0),
90  fTracksCont(0), fCaloClustersCont(0),
91  fESD(0), fAOD(0),
92  fMaxPatch(0),
93  fhnTriggerInfo(0), //trigger QA
94  fMainPatchType(kManual),
95  //fMainPatchType(kEmcalJet),
96  //fMainTrigCat(kTriggerLevel1Jet),
97  fMainTrigCat(kTriggerLevel1Gamma),
98  //fMainTrigCat(kTriggerRecalcJet), // Recalculated max trigger patch; does not need to be above trigger threshold
99  //fMainTrigCat(kTriggerRecalcGamma),
100  fMainTrigSimple(kFALSE),
101  fHistTriggerBitInfo(0),
102  fHistMaxTriggerBitInfo(0),
103  fHistEventSelection(0),
104  fHistRecalcGASize(0),
105  fHistRecalcGAEnergy(0),
106  fHistClusEvPatchE(0),
107  fHistdEtaPatchvdPhiPatch(0),
108  fHistRawJetPtvPatchE(0),
109  fHistMatchedClusJet(0),
110  fHistMatchdetadphi(0),
111  fHistdEtaPatchvdPhiPatchCMtoGeo(0),
112  fHistEventClusSpect(0),
113  fhnJetTrigger(0x0)
114 {
115  // Default constructor.
116  SetMakeGeneralHistograms(kTRUE);
117 
118 }
119 
120 //________________________________________________________________________
122  AliAnalysisTaskEmcalJet(name,kTRUE),
123  event(0),
124  fEventTrigEMCALL1Gamma1(0),
125  fEventTrigEMCALL1Gamma2(0),
126  fInputEvent(0x0),
127  fPhimin(-10), fPhimax(10),
128  fEtamin(-0.9), fEtamax(0.9),
129  fAreacut(0.0),
130  fJetHIpt(20.0),
131  fTrackEta(0.9),
132  fFillHists(0),
133  fMatchDist(0),
134  fTracksCont(0), fCaloClustersCont(0),
135  fESD(0), fAOD(0),
136  fMaxPatch(0),
137  fhnTriggerInfo(0), //trigger QA
138  fMainPatchType(kManual),
139  //fMainPatchType(kEmcalJet),
140  //fMainTrigCat(kTriggerLevel1Jet),
141  fMainTrigCat(kTriggerLevel1Gamma),
142  //fMainTrigCat(kTriggerRecalcJet), // Recalculated max trigger patch; does not need to be above trigger threshold
143  //fMainTrigCat(kTriggerRecalcGamma),
144  fMainTrigSimple(kFALSE),
145  fHistTriggerBitInfo(0),
146  fHistMaxTriggerBitInfo(0),
147  fHistEventSelection(0),
148  fHistRecalcGASize(0),
149  fHistRecalcGAEnergy(0),
150  fHistClusEvPatchE(0),
151  fHistdEtaPatchvdPhiPatch(0),
152  fHistRawJetPtvPatchE(0),
153  fHistMatchedClusJet(0),
154  fHistMatchdetadphi(0),
155  fHistdEtaPatchvdPhiPatchCMtoGeo(0),
156  fHistEventClusSpect(0),
157  fhnJetTrigger(0x0)
158 {
159 
161 
162  DefineInput(0,TChain::Class());
163  DefineOutput(1, TList::Class());
164 }
165 
166 //_______________________________________________________________________
168 {
169  // destructor
170  //
171  if (fOutput) {
172  delete fOutput;
173  fOutput = 0;
174  }
175 }
176 
177 //________________________________________________________________________
179 {
180  if (! fCreateHisto)
181  return;
183 
184  //no jets, just analysis tracks and clusters
187  fTracksCont->SetClassName("AliVTrack");
188  fCaloClustersCont->SetClassName("AliVCluster");
189  TString histname;
190  TString name;
191  TString title;
192 
193  if(fFillHists>0){
194  fHistClusEvPatchE = new TH2F("ClusterEvPatchE","Cluster Energy v Patch E; Cluster Energy (GeV); Trigger Patch Energy (GeV)",100,0,100,100,0,100);
195  fHistdEtaPatchvdPhiPatch = new TH2F("dEtaPatchdPhiPatch; #eta; #phi","dEtaPatchdPhiPatch; #Delta#eta; #Delta#phi",100,-0.1,0.1,100,-0.1,0.1);
196  fHistRawJetPtvPatchE = new TH2F("RawJetPtvTriggerE","RawJetPtvTriggerE; Jet Pt (GeV); Trigger Patch Energy (GeV)",100,0,100,100,0,100);
197  fHistTriggerBitInfo = new TH1F("TriggerBitInfo","TriggerBitInfo",15,0.5,15.5);
198  fHistMaxTriggerBitInfo = new TH1F("MaxPatchTriggerInfo","MaxPatchTriggerInfo",15,0.5,15.5);
199  fHistEventSelection = new TH1F("EventSelectionQA","EventSelectionQA",18,0.5,18.5);
200  fHistRecalcGASize = new TH2F("RecalcGAPatchSize","Patch_Size; #eta(Towers); #phi(Towers)",40,0,40,40,0,40);
201  fHistRecalcGAEnergy = new TH1F("ReclacGAEnergy","RecalcGAEnergy; Energy(GeV)",150,0,150);
202  fHistdEtaPatchvdPhiPatchCMtoGeo = new TH2F("dEtadPhi_GApatchCMtoGeo","dEtadPhi_GApatchCMtoGeo; #eta; #phi",100,-0.1,0.1,100,-0.1,0.1);
203  fHistEventClusSpect = new TH1F ("EventClusterSpectrum","ClusterEnergySpectrum",105,-5,100);
204 
205  fHistTriggerBitInfo->GetXaxis()->SetBinLabel(1,"IsLevel0");
206  fHistTriggerBitInfo->GetXaxis()->SetBinLabel(2,"IsJetlow");
207  fHistTriggerBitInfo->GetXaxis()->SetBinLabel(3,"IsJetHigh");
208  fHistTriggerBitInfo->GetXaxis()->SetBinLabel(4,"IsGammalow");
209  fHistTriggerBitInfo->GetXaxis()->SetBinLabel(5,"IsGammaHigh");
210  fHistTriggerBitInfo->GetXaxis()->SetBinLabel(6,"IsMainTrigger");
211  fHistTriggerBitInfo->GetXaxis()->SetBinLabel(7,"IsJetLowSimple");
212  fHistTriggerBitInfo->GetXaxis()->SetBinLabel(8,"IsJetHighSimple");
213  fHistTriggerBitInfo->GetXaxis()->SetBinLabel(9,"IsGammaLowSimple");
214  fHistTriggerBitInfo->GetXaxis()->SetBinLabel(10,"IsGammaHighSimple");
215  fHistTriggerBitInfo->GetXaxis()->SetBinLabel(11,"IsMainTriggerSimple");
216  fHistTriggerBitInfo->GetXaxis()->SetBinLabel(12,"IsOfflineSimple");
217  fHistTriggerBitInfo->GetXaxis()->SetBinLabel(13,"IsRecalcJet");
218  fHistTriggerBitInfo->GetXaxis()->SetBinLabel(14,"IsRecalcGamma");
219 
220  fHistMaxTriggerBitInfo->GetXaxis()->SetBinLabel(1,"IsLevel0");
221  fHistMaxTriggerBitInfo->GetXaxis()->SetBinLabel(2,"IsJetlow");
222  fHistMaxTriggerBitInfo->GetXaxis()->SetBinLabel(3,"IsJetHigh");
223  fHistMaxTriggerBitInfo->GetXaxis()->SetBinLabel(4,"IsGammalow");
224  fHistMaxTriggerBitInfo->GetXaxis()->SetBinLabel(5,"IsGammaHigh");
225  fHistMaxTriggerBitInfo->GetXaxis()->SetBinLabel(6,"IsMainTrigger");
226  fHistMaxTriggerBitInfo->GetXaxis()->SetBinLabel(7,"IsJetLowSimple");
227  fHistMaxTriggerBitInfo->GetXaxis()->SetBinLabel(8,"IsJetHighSimple");
228  fHistMaxTriggerBitInfo->GetXaxis()->SetBinLabel(9,"IsGammaLowSimple");
229  fHistMaxTriggerBitInfo->GetXaxis()->SetBinLabel(10,"IsGammaHighSimple");
230  fHistMaxTriggerBitInfo->GetXaxis()->SetBinLabel(11,"IsMainTriggerSimple");
231  fHistMaxTriggerBitInfo->GetXaxis()->SetBinLabel(12,"IsOfflineSimple");
232  fHistMaxTriggerBitInfo->GetXaxis()->SetBinLabel(13,"IsRecalcJet");
233  fHistMaxTriggerBitInfo->GetXaxis()->SetBinLabel(14,"IsRecalcGamma");
234 
235  fHistEventSelection->GetXaxis()->SetBinLabel(1,"AliVEvent::kMB");
236  fHistEventSelection->GetXaxis()->SetBinLabel(2,"AliVEvent::kINT7");
237  fHistEventSelection->GetXaxis()->SetBinLabel(3,"AliVEvent::kHighMult");
238  fHistEventSelection->GetXaxis()->SetBinLabel(4,"AliVEvent::kEMC1");
239  fHistEventSelection->GetXaxis()->SetBinLabel(5,"AliVEvent::kCINT5");
240  fHistEventSelection->GetXaxis()->SetBinLabel(6,"AliVEvent::kEMC7");
241  fHistEventSelection->GetXaxis()->SetBinLabel(7,"AliVEvent::kEMC8");
242  fHistEventSelection->GetXaxis()->SetBinLabel(8,"AliVEvent::kEMCEJE");
243  fHistEventSelection->GetXaxis()->SetBinLabel(9,"AliVEvent::kEMCEGA");
244  fHistEventSelection->GetXaxis()->SetBinLabel(10,"AliVEvent::kCentral");
245  fHistEventSelection->GetXaxis()->SetBinLabel(11,"AliVEvent::kSemiCentral");
246  fHistEventSelection->GetXaxis()->SetBinLabel(12,"AliVEvent::kZED");
247  fHistEventSelection->GetXaxis()->SetBinLabel(13,"AliVEvent::kINT8");
248  fHistEventSelection->GetXaxis()->SetBinLabel(14,"AliVEvent::kFastOnly");
249  fHistEventSelection->GetXaxis()->SetBinLabel(15,"AliVEvent::kAnyINT");
250  fHistEventSelection->GetXaxis()->SetBinLabel(16,"AliVEvent::kAny");
251  fHistEventSelection->GetXaxis()->SetBinLabel(17,"Event_Counter");
252 
253 
254  // PT bins used to be (2000, -100, 300)
255 
256  Int_t fgkNCentBins = 21;
257  Float_t kMinCent = 0.;
258  Float_t kMaxCent = 105.;
259 
260  Int_t fgkNVZEROBins = 100;
261  Float_t kMinVZERO = 0.;
262  Float_t kMaxVZERO = 25000;
263 
264  const Int_t fgkNEPatch = 100;
265  Float_t kMinEPatch = 0.;
266  Float_t kMaxEPatch = 200.;
267 
268  const Int_t fgkNADC = 100;
269  Float_t kMinADC = 0.;
270  Float_t kMaxADC = 1500.;
271 
272  const Int_t fgkNEta = 10;
273  const Int_t fgkNPhi = 10;
274 
275  const Int_t nDim = 8;//cent;V0mult;ptjet1;ptjet2;Epatch;ADCpatch;EtaPatch;PhiPatch
276  const Int_t nBins[nDim] = {fgkNCentBins,fgkNVZEROBins,fgkNEta,fgkNPhi,fgkNEPatch,fgkNADC,fgkNEta,fgkNPhi};
277  const Double_t xmin0[nDim] = {kMinCent,kMinVZERO,-0.7,1.4,kMinEPatch,kMinADC,-0.7,1.4};
278  const Double_t xmax0[nDim] = {kMaxCent,kMaxVZERO,0.7,3.14,kMaxEPatch,kMaxADC, 0.7,3.14};
279  fhnTriggerInfo = new THnSparseF("fhnTriggerInfo", "hnTriggerInfo;cent;V0mult;EtaCMPatch;PhiCMPatch;Epatch;ADCpatch;EtaGeoPatch;PhiGeoPatch",nDim,nBins,xmin0,xmax0);
280 
281  fOutput->Add(fhnTriggerInfo);
292 
293  }//Fill Histograms
294 
295  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
296  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
297  if(!inputHandler) {
298  AliFatal("Input handler needed");
299  return;
300  }
301 
302 
303  UInt_t bitcoded1 = 0;
304  bitcoded1 = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<10 | 1<<11 ;
305  fhnJetTrigger = NewTHnSparseDJetTrigger("fhnJetTrigger", bitcoded1);
306 
307  cout << "_______________Created Sparse__________________" << endl;
308  fOutput->Add(fhnJetTrigger);
309 
310  PostData(1, fOutput);
311 
312 }
313 
314 //________________________________________________________
316 {
317  // Initialize the analysis
319 
320  if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
321  if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
322 } // end of ExecOnce
323 
324 //________________________________________________________________________
326 {
327  // check to see if we have any tracks
328  if (!fTracks) return kTRUE;
329  if (!fJets) return kTRUE;
330 
331  UInt_t trig = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
332  if(fFillHists>0){
333  if(trig & AliVEvent::kMB) fHistEventSelection->Fill(1);
334  if(trig & AliVEvent::kINT7) fHistEventSelection->Fill(2);
335  if(trig & AliVEvent::kHighMult) fHistEventSelection->Fill(3);
336  if(trig & AliVEvent::kEMC1) fHistEventSelection->Fill(4);
337  if(trig & AliVEvent::kCINT5) fHistEventSelection->Fill(5);
338  if(trig & AliVEvent::kEMC7) fHistEventSelection->Fill(6);
339  if(trig & AliVEvent::kEMC8) fHistEventSelection->Fill(7);
340  if(trig & AliVEvent::kEMCEJE) fHistEventSelection->Fill(8);
341  if(trig & AliVEvent::kEMCEGA) fHistEventSelection->Fill(9);
342  if(trig & AliVEvent::kCentral) fHistEventSelection->Fill(10);
343  if(trig & AliVEvent::kSemiCentral) fHistEventSelection->Fill(11);
344  if(trig & AliVEvent::kZED) fHistEventSelection->Fill(12);
345  if(trig & AliVEvent::kINT8) fHistEventSelection->Fill(13);
346  if(trig & AliVEvent::kFastOnly) fHistEventSelection->Fill(14);
347  if(trig & AliVEvent::kAnyINT) fHistEventSelection->Fill(15);
348  if(trig & AliVEvent::kAny) fHistEventSelection->Fill(16);
349  fHistEventSelection->Fill(17);
350  }
351  // what kind of event do we have: AOD or ESD?
352  Bool_t useAOD;
353  if (dynamic_cast<AliAODEvent*>(InputEvent())) useAOD = kTRUE;
354  else useAOD = kFALSE;
355 
356  fEventTrigEMCALL1Gamma1 = kFALSE;
357  fEventTrigEMCALL1Gamma2 = kFALSE;
358 
359  // if we have ESD event, set up ESD object
360  if(!useAOD){
361  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
362  if (!fESD) {
363  AliError(Form("ERROR: fESD not available\n"));
364  return kTRUE;
365  }
366  }
367 
368  // if we have AOD event, set up AOD object
369  if(useAOD){
370  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
371  if(!fAOD) {
372  AliError(Form("ERROR: fAOD not available\n"));
373  return kTRUE;
374  }
375  }
376 
377  // get vertex information
378  Double_t fvertex[3]={0,0,0};
379  InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
380 
381  // create pointer to list of input event
382  TList *list = InputEvent()->GetList();
383  if(!list) {
384  AliError(Form("ERROR: list not attached\n"));
385  return kTRUE;
386  }
387 
388  // background density
389  fRhoVal = fRho->GetVal();
390 
391  // initialize TClonesArray pointers to jets and tracks
392  TClonesArray *jets = 0;
393  // get Jets object
394  jets = dynamic_cast<TClonesArray*>(list->FindObject(fJets));
395  if(!jets){
396  AliError(Form("Pointer to jets %s == 0", fJets->GetName()));
397  return kTRUE;
398  } // verify existence of jets
399 
400  // get number of jets and tracks
401  const Int_t Njets = jets->GetEntries();
402  if(Njets<1) return kTRUE;
403 
404  event++;
405  //cout<<"Event #: "<<event<<" Number of Clusters: "<<fCaloClustersCont->GetNClusters()<<" Number of Tracks: "<<fTracksCont->GetNParticles()<<" Number Of Jets: "<< Njets <<endl;
406 
407  if (fCaloClustersCont) {
408  fCaloClustersCont->ResetCurrentID();
409  AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster();
410  while(cluster) {
411  TLorentzVector nPart;
412  cluster->GetMomentum(nPart, fVertex);
413  if(fFillHists>0) fHistEventClusSpect->Fill(cluster->E());
415 
416  }
417  }
418 
419  // Start Jet Analysis
420  // initialize jet parameters
421  Int_t ijethi=-1;
422  Double_t highestjetpt=0.0;
423  Double_t GAeta = 0.0;
424  Double_t GAphi = 0.0;
425  Double_t GAetaCM = -999.;
426  Double_t GAphiCM = -999.;
427  Double_t PatchEE = 0.0;
428  Double_t VZEROAmp = (Double_t)(InputEvent()->GetVZEROData()->GetTriggerChargeA() + InputEvent()->GetVZEROData()->GetTriggerChargeC());
429  Double_t PatchEtaMax, PatchEtaMin, PatchPhiMax, PatchPhiMin, PatchAreaE, PatchAreaP;
430  // Get GA Trigger Info
431  if(fTriggerPatchInfo) {
433  else if(fMainPatchType==kEmcalJet){
435  }
436  PatchEE = fMaxPatch->GetPatchE();
437  GAphi = fMaxPatch->GetPhiGeo();
438  GAeta = fMaxPatch->GetEtaGeo();
439  GAetaCM = fMaxPatch->GetEtaCM();
440  GAphiCM = fMaxPatch->GetPhiCM();
441  PatchEtaMax = fMaxPatch->GetEtaMax();
442  PatchEtaMin = fMaxPatch->GetEtaMin();
443  PatchPhiMax = fMaxPatch->GetPhiMax();
444  PatchPhiMin = fMaxPatch->GetPhiMin();
445  PatchAreaE = PatchEtaMax - PatchEtaMin;
446  PatchAreaP = PatchPhiMax - PatchPhiMin;
447 
448  if(fFillHists>0){
449  if (fMaxPatch->IsLevel0() == 1 ) fHistMaxTriggerBitInfo->Fill(1);
450  if (fMaxPatch->IsJetLow() == 1 ) fHistMaxTriggerBitInfo->Fill(2);
451  if (fMaxPatch->IsJetHigh() == 1 ) fHistMaxTriggerBitInfo->Fill(3);
452  if (fMaxPatch->IsGammaLow() == 1 ) fHistMaxTriggerBitInfo->Fill(4);
453  if (fMaxPatch->IsGammaHigh() == 1 ) fHistMaxTriggerBitInfo->Fill(5);
454  if (fMaxPatch->IsMainTrigger() == 1 ) fHistMaxTriggerBitInfo->Fill(6);
455  if (fMaxPatch->IsJetLowSimple() == 1 ) fHistMaxTriggerBitInfo->Fill(7);
456  if (fMaxPatch->IsJetHighSimple() == 1 ) fHistMaxTriggerBitInfo->Fill(8);
457  if (fMaxPatch->IsGammaLowSimple() == 1 ) fHistMaxTriggerBitInfo->Fill(9);
458  if (fMaxPatch->IsGammaHighSimple() == 1 ) fHistMaxTriggerBitInfo->Fill(10);
459  if (fMaxPatch->IsMainTriggerSimple() == 1 ) fHistMaxTriggerBitInfo->Fill(11);
460  if (fMaxPatch->IsOfflineSimple() == 1 ) fHistMaxTriggerBitInfo->Fill(12);
461  if (fMaxPatch->IsRecalcJet() == 1 ) fHistMaxTriggerBitInfo->Fill(13);
462 
463  if (fMaxPatch->IsRecalcGamma() == 1 ){
464  fHistMaxTriggerBitInfo->Fill(14);
465  fHistRecalcGASize->Fill(PatchAreaE/0.014, PatchAreaP/0.014);
466  fHistRecalcGAEnergy->Fill(PatchEE);
467  }
468  }
469 
470  Double_t var[8] = {fCent, VZEROAmp, GAetaCM, GAphiCM, fMaxPatch->GetPatchE(), (Double_t)fMaxPatch->GetADCAmp(), fMaxPatch->GetEtaGeo(), fMaxPatch->GetPhiGeo()
471  };
472  if(fFillHists>0) fhnTriggerInfo->Fill(var);
473  if(fFillHists>0) fHistdEtaPatchvdPhiPatchCMtoGeo->Fill(GAetaCM - GAeta, GAphiCM - GAphi);
474 
475  }
476  // **********************************************************************
477  // JET LOOP
478  // **********************************************************************
479  // loop over jets in the event and make appropriate cuts
480  for (Int_t iJets = 0; iJets < Njets; ++iJets) {
481  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
482  if (!jet) continue;
483 
484  // phi of jet, constrained to 1.6 < Phi < 2.94
485  float jetphi = jet->Phi(); // phi of jet
486  // apply jet cuts
487  if(!AcceptMyJet(jet)) continue;
488  Int_t JetClusters = jet->GetNumberOfClusters();
489  Int_t JetTracks = jet -> GetNumberOfTracks();
490  // Initializations and Calculations
491  Double_t jetptraw = jet->Pt(); // raw pT of jet
492  Double_t jetPt = -500; // initialize corr jet pt LOCAL
493  Double_t jetarea = -500; // initialize jet area
494  jetarea = jet->Area(); // jet area
495  jetPt = jet->Pt() - jetarea*fRhoVal; // semi-corrected pT of jet from GLOBAL rho value
496 
497  if(jet->Pt() > fJetHIpt) {
498 
499 
500  Float_t jeteta = -999.;
501  Float_t jetphi = -999.;
502  jeteta = jet->Eta();
503  jetphi = jet->Phi();
504 
505  Double_t deta = 999;
506  Double_t dphi = 999;
507 
508  //******************************Cluster Matched
509  //**************************************************************
510  for(int iCluster = 0; iCluster <= fCaloClustersCont->GetNClusters(); iCluster++){
511  //Get closest track to cluster to track matching!!!!!
512  //AliVCluster *cluster = fCaloClustersCont->GetNextAcceptedCluster(iCluster);
513  AliVCluster *clusMatch = fCaloClustersCont->GetCluster(iCluster);
514  if(!clusMatch) continue;
515  if(! IsJetCluster(jet, iCluster, kFALSE)) continue;
516 
517  Double_t mClusterE = clusMatch->E();
518  Float_t pos_mc[3];
519  clusMatch->GetPosition(pos_mc); // Get cluster position
520  TVector3 mcp(pos_mc);
521  Double_t maxMatchedPatchEta = (fMatchDist * 0.014);
522  Double_t maxMatchedPatchPhi = (fMatchDist * 0.014);
523  Double_t mtchPosition = maxMatchedPatchPhi * maxMatchedPatchPhi + maxMatchedPatchEta * maxMatchedPatchEta;
524  Double_t rGApatch = TMath::Sqrt(mtchPosition);
525  Double_t rJetCluster = -999.;
526 
527  Double_t PatchClusterdiff = (mcp.Phi() - fMaxPatch->GetPhiGeo())*(mcp.Phi() - fMaxPatch->GetPhiGeo()) + (mcp.PseudoRapidity() - fMaxPatch->GetEtaGeo())*(mcp.PseudoRapidity() - fMaxPatch->GetEtaGeo());
528  rJetCluster = TMath::Sqrt(PatchClusterdiff);
529 
530  //Matched Trigger Jet
531  if ( rJetCluster <= rGApatch ){
532 
533  Double_t dEtaPatch = 0.0;
534  Double_t dPhiPatch = 0.0;
535  dEtaPatch = fMaxPatch->GetEtaGeo() - mcp.PseudoRapidity();
536  dPhiPatch = fMaxPatch->GetPhiGeo() - mcp.Phi();
537  if(fFillHists>0) fHistClusEvPatchE->Fill(mClusterE,fMaxPatch->GetPatchE());
538  if(fFillHists>0) fHistdEtaPatchvdPhiPatch->Fill(dEtaPatch,dPhiPatch);
539  if(fFillHists>0) fHistRawJetPtvPatchE->Fill(jetptraw,fMaxPatch->GetPatchE());
540  Double_t JetTrig[18] = {jetptraw, jetPt, jeteta, jetphi, jet->E(), fMaxPatch->GetPatchE(), fMaxPatch->GetEtaGeo(), fMaxPatch->GetPhiGeo(), (Double_t)fMaxPatch->GetADCAmp(), mClusterE, mcp.PseudoRapidity(), mcp.Phi() };
541  fhnJetTrigger->Fill(JetTrig);
542  } // matching cluster of a jet to trigger patch
543  } // cluster for
544  } // highest pt jet cut
545  } // LOOP over JETS in event
546 
547  return kTRUE;
548 
549 }
550 //________________________________________________________________________
552 {
553  cout<<"###########################"<<endl;
554  cout<<"#### Task Finished ####"<<endl;
555  cout<<"###########################"<<endl;
556  cout<<"###########################"<<endl;
557 } // end of terminate
558 
559 
560 //________________________________________________________________________
562 
563  //Find main trigger
564  if(!fTriggerPatchInfo)
565  return;
566 
567  //number of patches in event
568  Int_t nPatch = fTriggerPatchInfo->GetEntriesFast();
569  //extract main trigger patch
570  Double_t emax = -1.;
571  for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
572  AliEMCALTriggerPatchInfo *patch = (AliEMCALTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
573  if (!patch) continue;
574 
575  if(fFillHists>0){
576  if (patch->IsLevel0() == 1 ) fHistTriggerBitInfo->Fill(1);
577  if (patch->IsJetLow() == 1 ) fHistTriggerBitInfo->Fill(2);
578  if (patch->IsJetHigh() == 1 ) fHistTriggerBitInfo->Fill(3);
579  if (patch->IsGammaLow() == 1 ) fHistTriggerBitInfo->Fill(4);
580  if (patch->IsGammaHigh() == 1 ) fHistTriggerBitInfo->Fill(5);
581  if (patch->IsMainTrigger() == 1 ) fHistTriggerBitInfo->Fill(6);
582  if (patch->IsJetLowSimple() == 1 ) fHistTriggerBitInfo->Fill(7);
583  if (patch->IsJetHighSimple() == 1 ) fHistTriggerBitInfo->Fill(8);
584  if (patch->IsGammaLowSimple() == 1 ) fHistTriggerBitInfo->Fill(9);
585  if (patch->IsGammaHighSimple() == 1 ) fHistTriggerBitInfo->Fill(10);
586  if (patch->IsMainTriggerSimple() == 1 ) fHistTriggerBitInfo->Fill(11);
587  if (patch->IsOfflineSimple() == 1 ) fHistTriggerBitInfo->Fill(12);
588  if (patch->IsRecalcJet() == 1 ) fHistTriggerBitInfo->Fill(13);
589  if (patch->IsRecalcGamma() == 1 ) fHistTriggerBitInfo->Fill(14);
590  }
591  //Force ExtractMainPatch to return Recalc GA info
592  if (patch->IsRecalcGamma() == 1 ){
593  if(patch->GetPatchE()>emax) {
594  fMaxPatch = patch;
595  emax = patch->GetPatchE();
596  }
597  }
598  }
599 }
600 
601 //________________________________________________________________________
603  //applies all jet cuts except pt
604  if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) return 0;
605  if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) return 0;
606  if (jet->Area()<fAreacut) return 0;
607  // prevents 0 area jets from sneaking by when area cut == 0
608  if (jet->Area()==0) return 0;
609  //exclude jets with extremely high pt tracks which are likely misreconstructed
610  if(jet->MaxTrackPt()>100) return 0;
611 
612  //passed all above cuts
613  return 1;
614 }
615 
616 
618 {
619  // generate new THnSparseD JetQA, axes are defined in GetDimParamsJetQA()
620  Int_t count = 0;
621  UInt_t tmp = entries;
622  while(tmp!=0){
623  count++;
624  tmp = tmp &~ -tmp; // clear lowest bit
625  }
626 
627  TString hnTitle(name);
628  const Int_t dim = count;
629  Int_t nbins[dim];
630  Double_t xmin[dim];
631  Double_t xmax[dim];
632 
633  Int_t i=0;
634  Int_t c=0;
635  while(c<dim && i<32){
636  if(entries&(1<<i)){
637 
638  TString label("");
639  GetDimParamsJetTrigger(i, label, nbins[c], xmin[c], xmax[c]);
640  hnTitle += Form(";%s",label.Data());
641  c++;
642  }
643 
644  i++;
645  }
646  hnTitle += ";";
647 
648  return new THnSparseD(name, hnTitle.Data(), dim, nbins, xmin, xmax);
649 } // end of NewTHnSparseF JetTrigger
650 
651 
652 
654 {
655  // stores label and binning of axis for THnSparse
656  const Double_t pi = TMath::Pi();
657 
658  switch(iEntry){
659 
660  case 0:
661  label = "Jet Pt";
662  nbins = 50;
663  xmin = 0.;
664  xmax = 250.;
665  break;
666 
667  case 1:
668  label = "Jet Corr Pt";
669  nbins = 50;
670  xmin = -50.;
671  xmax = 200.;
672  break;
673 
674  case 2:
675  label = "Jet Eta";
676  nbins = 24;
677  xmin = -1.2;
678  xmax = 1.2;
679  break;
680 
681  case 3:
682  label = "Jet Phi";
683  nbins = 72;
684  xmin = 0;
685  xmax = 2*pi;
686  break;
687 
688  case 4:
689  label = "Jet E";
690  nbins = 50;
691  xmin = 0.;
692  xmax = 250.;
693  break;
694 
695  case 5:
696  label = "Trigger Patch E";
697  nbins = 100;
698  xmin = 0.;
699  xmax = 100.;
700  break;
701 
702  case 6:
703  label = "Trigger Patch Eta";
704  nbins = 24;
705  xmin = -1.2;
706  xmax = 1.2;
707  break;
708 
709  case 7:
710  label = "Trigger Patch Phi";
711  nbins = 72;
712  xmin = 0;
713  xmax = 2*pi;
714  break;
715 
716  case 8:
717  label = "Trigger Patch ADC";
718  nbins = 150;
719  xmin = 0;
720  xmax = 1500;
721  break;
722 
723  case 9:
724  label = "Cluster E";
725  nbins = 100;
726  xmin = 0.;
727  xmax = 100.;
728  break;
729 
730  case 10:
731  label = "Cluster Eta";
732  nbins = 24;
733  xmin = -1.2;
734  xmax = 1.2;
735  break;
736 
737  case 11:
738  label = "Cluster Phi";
739  nbins = 72;
740  xmin = 0;
741  xmax = 2*pi;
742  break;
743 
744  } // end of switch
745 } // end of getting dim-params
Double_t Area() const
Definition: AliEmcalJet.h:130
AliEMCALTriggerPatchInfo * fMaxPatch
// AOD Object
double Double_t
Definition: External.C:58
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:27
Double_t Eta() const
Definition: AliEmcalJet.h:121
Double_t Phi() const
Definition: AliEmcalJet.h:117
TCanvas * c
Definition: TestFitELoss.C:172
Double_t E() const
Definition: AliEmcalJet.h:119
virtual THnSparse * NewTHnSparseDJetTrigger(const char *name, UInt_t entries)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Bool_t IsJetCluster(AliEmcalJet *jet, Int_t iclus, Bool_t sorted=kFALSE) const
int Int_t
Definition: External.C:63
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:138
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
AliRhoParameter * fRho
! event rho
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
Double_t MaxTrackPt() const
Definition: AliEmcalJet.h:155
Int_t GetNClusters() const
void ExecOnce()
Perform steps needed to initialize the analysis.
Double_t fCent
!event centrality
AliEMCALTriggerPatchInfo * GetMainTriggerPatch(TriggerCategory triggersel=kTriggerLevel1Jet, Bool_t doSimpleOffline=kFALSE)
Get main trigger match.
AliVCluster * GetCluster(Int_t i) const
TClonesArray * fJets
! jets
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliEmcalList * fOutput
!output list
TClonesArray * fTracks
!tracks
Double_t fVertex[3]
!event vertex
TH1F * fHistEventClusSpect
Plot of Trigger Geo pos - Trigger COM pos.
Bool_t fCreateHisto
whether or not create histograms
void SetMakeGeneralHistograms(Bool_t g)
TClonesArray * fTriggerPatchInfo
!trigger patch info array
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
const Double_t pi
const Int_t nbins
bool Bool_t
Definition: External.C:53
virtual void GetDimParamsJetTrigger(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
Double_t fRhoVal
! event rho value, same for local rho
AliVCluster * GetNextAcceptCluster()
TH2F * fHistdEtaPatchvdPhiPatch
Matched Cluster E vs Patch E.
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.