AliPhysics  fe039ad (fe039ad)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliEmcalPicoTrackInGridMaker.cxx
Go to the documentation of this file.
1 // $Id$
2 //
3 // Class to put collection of tracks into grid of PicoTracks
4 //
5 // Author: M. Verweij
6 
7 #include <TClonesArray.h>
8 #include <TRandom3.h>
9 #include <TProfile.h>
10 #include <TH3F.h>
11 
12 #include "AliAnalysisManager.h"
13 #include "AliLog.h"
14 #include "AliPicoTrack.h"
15 #include "AliVTrack.h"
16 #include "AliEmcalJet.h"
17 #include "AliParticleContainer.h"
18 #include "AliJetContainer.h"
19 
21 
23 
24 //________________________________________________________________________
27  fTracksOutName("PicoTracksInGrid"),
28  fTracksOut(0),
29  fL1Slide(0),
30  fCellSize(0.0145),
31  fMinCellE(0.15),
32  fExclLeadingPatch(0),
33  fPatchSub(3),
34  fRhoMean(184.),
35  fNCells(-1),
36  fNCellsEMCal(-1),
37  fNCellsDCal(-1),
38  fMultVsRho(0)
39 {
40  // Constructor.
41 
42  for(Int_t i = 0; i<2; i++) {
43  fCellGrid[i] = 0;
44  fMiniPatchGrid[i] = 0;
45  fActiveAreaMP[i] = 0;
46  }
47 
48  fPhiMin[0] = 1.405;
49  fPhiMax[0] = 1.405+TMath::DegToRad()*110.;
50  fPhiMin[1] = 4.547;
51  fPhiMax[1] = 5.71;
52  fEtaMin[0] = -0.7;
53  fEtaMax[0] = 0.7;
54  fEtaMin[1] = -0.7;
55  fEtaMax[1] = 0.7;
56 
57  for(Int_t i = 0; i<5; i++) {
58  fNPatchesEMCal[i] = 0;
59  fh1RhoEmcal[i] = 0;
60  fh1RhoDcal[i] = 0;
61  fPatchEnVsActivityEmcal[i] = 0;
62  fPatchEnVsActivityDcal[i] = 0;
63 
64  for(Int_t j = 0; j<2; j++) {
65  fPatchGrid[j][i] = 0;
66  fActiveAreaMPP[j][i] = 0;
67  fActiveAreaCP[j][i] = 0;
68  fPatchECorr[j][i] = 0;
69  fPatchECorrPar[j][i] = 0;
70  fPatchERaw[j][i] = 0;
71  fPatchECorrRho[j][i] = 0;
72  fPatchECorrECorrRho[j][i] = 0;
73  fh2JetPtPatchECorr[j][i] = 0;
74  }
75  fh2PatchEtaPhiEmcal[i] = 0;
76  fh2PatchEtaPhiDcal[i] = 0;
77  }
78 
79  for(Int_t i = 0; i<3; i++) {
80  fh2MedianTypeEmcal[i] = 0;
81  fh2MedianTypeDcal[i] = 0;
82  fpMedianTypeEmcal[i] = 0;
83  fpMedianTypeDcal[i] = 0;
84  }
85  SetMakeGeneralHistograms(kTRUE);
86 }
87 
88 //________________________________________________________________________
90  AliAnalysisTaskEmcalJet(name,kTRUE),
91  fTracksOutName("PicoTracksInGrid"),
92  fTracksOut(0),
93  fL1Slide(0),
94  fCellSize(0.0145),
95  fMinCellE(0.15),
96  fExclLeadingPatch(0),
97  fPatchSub(3),
98  fRhoMean(184.),
99  fNCells(-1),
100  fNCellsEMCal(-1),
101  fNCellsDCal(-1),
102  fMultVsRho(0)
103 {
104  // Constructor.
105 
106  for(Int_t i = 0; i<2; i++) {
107  fCellGrid[i] = 0;
108  fMiniPatchGrid[i] = 0;
109  fActiveAreaMP[i] = 0;
110  }
111 
112  fPhiMin[0] = 1.405;
113  fPhiMax[0] = 1.405+TMath::DegToRad()*110.;
114  fPhiMin[1] = 4.547;
115  fPhiMax[1] = 5.71;
116  fEtaMin[0] = -0.7;
117  fEtaMax[0] = 0.7;
118  fEtaMin[1] = -0.7;
119  fEtaMax[1] = 0.7;
120 
121  for(Int_t i = 0; i<5; i++) {
122  fNPatchesEMCal[i] = 0;
123  fh1RhoEmcal[i] = 0;
124  fh1RhoDcal[i] = 0;
126  fPatchEnVsActivityDcal[i] = 0;
127 
128  for(Int_t j = 0; j<2; j++) {
129  fPatchGrid[j][i] = 0;
130  fActiveAreaMPP[j][i] = 0;
131  fActiveAreaCP[j][i] = 0;
132  fPatchECorr[j][i] = 0;
133  fPatchECorrPar[j][i] = 0;
134  fPatchERaw[j][i] = 0;
135  fPatchECorrRho[j][i] = 0;
136  fPatchECorrECorrRho[j][i] = 0;
137  fh2JetPtPatchECorr[j][i] = 0;
138  }
139  fh2PatchEtaPhiEmcal[i] = 0;
140  fh2PatchEtaPhiDcal[i] = 0;
141  }
142 
143  for(Int_t i = 0; i<3; i++) {
144  fh2MedianTypeEmcal[i] = 0;
145  fh2MedianTypeDcal[i] = 0;
146  fpMedianTypeEmcal[i] = 0;
147  fpMedianTypeDcal[i] = 0;
148  }
149 
151 }
152 
153 //________________________________________________________________________
155 {
156  // Destructor.
157 }
158 
159 //________________________________________________________________________
161 {
162  // Create my user objects.
164 
165  Int_t nBinsMed = 500;
166  Double_t minMed = 0.;
167  Double_t maxMed = 500.;
168 
169  Int_t nBinsPhiEmcal = 132+64;
170  Double_t phiMinEmcal = 1.436931 - 32.*fCellSize;
171  Double_t phiMaxEmcal = 3.292931 + 32.*fCellSize;
172 
173  Int_t nBinsPhiDcal = 80+64;
174  Double_t phiMinDcal = 4.664500 - 32.*fCellSize;
175  Double_t phiMaxDcal = 5.592500 + 32.*fCellSize;
176 
177  Int_t nBinsEta = 96+64;
178  Double_t etaMin = -0.696 - 32.*fCellSize;
179  Double_t etaMax = 0.696 + 32.*fCellSize;
180 
181  for(Int_t i = 0; i<3; i++) {
182  fh2MedianTypeEmcal[i] = new TH2F(Form("fh2MedianTypeEmcalAreaType%d",i),Form("fh2MedianTypeEmcalAreaType%d",i),5,0.5,5.5,nBinsMed,minMed,maxMed);
183  fOutput->Add(fh2MedianTypeEmcal[i]);
184 
185  fh2MedianTypeDcal[i] = new TH2F(Form("fh2MedianTypeDcalAreaType%d",i),Form("fh2MedianTypeDcalAreaType%d",i),5,0.5,5.5,nBinsMed,minMed,maxMed);
186  fOutput->Add(fh2MedianTypeDcal[i]);
187 
188  fpMedianTypeEmcal[i] = new TProfile(Form("fpMedianTypeEmcalAreaType%d",i),Form("fpMedianTypeEmcalAreaType%d",i),5,0.5,5.5,"s");
189  fOutput->Add(fpMedianTypeEmcal[i]);
190 
191  fpMedianTypeDcal[i] = new TProfile(Form("fpMedianTypeDcalAreaType%d",i),Form("fpMedianTypeDcalAreaType%d",i),5,0.5,5.5,"s");
192  fOutput->Add(fpMedianTypeDcal[i]);
193  }
194 
195  TString det[2] = {"Emcal","Dcal"};
196  for(Int_t i = 0; i<5; i++) { //loop over patch types
197  fh1RhoEmcal[i] = new TH1F(Form("fh1RhoEmcal_%d",i),Form("fh1RhoEmcal_%d",i),500,0.,1000.);
198  fOutput->Add(fh1RhoEmcal[i]);
199  fh1RhoDcal[i] = new TH1F(Form("fh1RhoDcal_%d",i),Form("fh1RhoDcal_%d",i),500,0.,1000.);
200  fOutput->Add(fh1RhoDcal[i]);
201 
202  fPatchEnVsActivityEmcal[i] = new TH2F(Form("fh2PatchEnVsActivityEmcal_%d",i),Form("fh2PatchEnVsActivityEmcal_%d",i),300,0.,300.,150,-0.5,149.5);
204 
205  fPatchEnVsActivityDcal[i] = new TH2F(Form("fh2PatchEnVsActivityDcal_%d",i),Form("fh2PatchEnVsActivityDcal_%d",i),300,0.,300.,150,-0.5,149.5);
207 
208  for(Int_t j = 0; j<2; j++) {
209  fPatchECorr[j][i] = new TH1F(Form("fPatchECorr%s_%d",det[j].Data(),i),Form("fPatchECorr%s_%d;#it{E}_{patch}^{corr}",det[j].Data(),i),250,-50.,200.);
210  fOutput->Add(fPatchECorr[j][i]);
211 
212  fPatchECorrPar[j][i] = new TH1F(Form("fPatchECorrPar%s_%d",det[j].Data(),i),Form("fPatchECorrPar%s_%d;#it{E}_{patch}^{corr}",det[j].Data(),i),250,-50.,200.);
213  fOutput->Add(fPatchECorrPar[j][i]);
214 
215  fPatchERaw[j][i] = new TH1F(Form("fPatchERaw%s_%d",det[j].Data(),i),Form("fPatchERaw%s_%d;#it{E}_{patch}^{corr}",det[j].Data(),i),250,-50.,200.);
216  fOutput->Add(fPatchERaw[j][i]);
217 
218  fPatchECorrRho[j][i] = new TH2F(Form("fPatchECorrRho%s_%d",det[j].Data(),i),Form("fPatchECorrRho%s_%d;#it{E}_{patch}^{corr};#rho",det[j].Data(),i),250,-50.,200.,500,0.,500.);
219  fOutput->Add(fPatchECorrRho[j][i]);
220 
221  fPatchECorrECorrRho[j][i] = new TH3F(Form("fPatchECorrECorrRho%s_%d",det[j].Data(),i),Form("fPatchECorrECorrRho%s_%d;#it{E}_{patch,det1}^{corr};#it{E}_{patch,det2}^{corr};#rho",det[j].Data(),i),210,-30.,180.,210,-30.,180.,250,0.,250.);
222  fOutput->Add(fPatchECorrECorrRho[j][i]);
223 
224  fh2JetPtPatchECorr[j][i] = new TH2F(Form("fh2JetPtPatchECorr%s_%d",det[j].Data(),i),Form("fh2JetPtPatchECorr%s_%d",det[j].Data(),i),250,-50.,200.,250,-50.,200.);
225  fOutput->Add(fh2JetPtPatchECorr[j][i]);
226  }
227 
228  fh2PatchEtaPhiEmcal[i] = new TH2F(Form("fh2PatchEtaPhiEmcal_%d",i),Form("fh2PatchEtaPhiEmcal_%d;#eta;#phi",i),nBinsEta,etaMin,etaMax,nBinsPhiEmcal,phiMinEmcal,phiMaxEmcal);
229  fOutput->Add(fh2PatchEtaPhiEmcal[i]);
230 
231  fh2PatchEtaPhiDcal[i] = new TH2F(Form("fh2PatchEtaPhiDcal_%d",i),Form("fh2PatchEtaPhiDcal_%d;#eta;#phi",i),nBinsEta,etaMin,etaMax,nBinsPhiDcal,phiMinDcal,phiMaxDcal);
232  fOutput->Add(fh2PatchEtaPhiDcal[i]);
233  }
234 
235  fMultVsRho = new TH2F("fMultVsRho","fMultVsRho",3000,0,3000,400,0,400);
236  fOutput->Add(fMultVsRho);
237 
238  PostData(1, fOutput);
239 }
240 
241 //________________________________________________________________________
243 {
244  // Main loop, called for each event.
245 
246  Bool_t b = CreateGridCells();
247  if(!b) return kFALSE;
248  b = CreateGridMiniPatches();
249  if(!b) return kFALSE;
250 
251  //L0 single shower trigger
252  CreateGridPatches(4,0);
253  // return kTRUE;
254 
255  //L1 triggers: sliding window
256  CreateGridPatches(4,1);
257  CreateGridPatches(8,1);
258  CreateGridPatches(16,1);
259  CreateGridPatches(32,1);
260 
261 
262  Double_t medL0 = CalculateMedian(0,0);
263  fh2MedianTypeEmcal[0]->Fill(0.5,medL0);
264  fpMedianTypeEmcal[0]->Fill(0.5,medL0);
265  medL0 = CalculateMedian(0,1);
266  fh2MedianTypeDcal[0]->Fill(0.5,medL0);
267  fpMedianTypeDcal[0]->Fill(0.5,medL0);
268 
269  Double_t medL1[4][2];
270  for(Int_t i = 0; i<4; i++) { //patches
271  for(Int_t type = 0; type<2; type++) { //EMCal or DCal
272  for(Int_t areaT = 0; areaT<1; areaT++) { //areay type (passive vs active)
273  medL1[i][type] = CalculateMedian(i+1,type,areaT);
274  if(type==0) {
275  fh2MedianTypeEmcal[areaT]->Fill((Double_t)(i+1)+0.5,medL1[i][type]);
276  fpMedianTypeEmcal[areaT]->Fill((Double_t)(i+1)+0.5,medL1[i][type]);
277  }
278  if(type==1) {
279  fh2MedianTypeDcal[areaT]->Fill((Double_t)(i+1)+0.5,medL1[i][type]);
280  fpMedianTypeDcal[areaT]->Fill((Double_t)(i+1)+0.5,medL1[i][type]);
281  }
282  }
283  }
284  }
285 
286  // subtract energy density and store energy distributions of corrected patches in histo
287  Int_t EleadID[5][2];
288  Double_t Elead[5][2];
289  Double_t EleadRaw[5][2];
290  for(Int_t i = 1; i<5; i++) { //patch types
292  for(Int_t type = 0; type<2; type++) {
293  EleadID[i][type] = -1;
294  Elead[i][type] = -1e6;
295  EleadRaw[i][type] = -1e6;
296  Int_t subType = 1;
297  if(type==1) subType = 0;
298  // for(Int_t j = 0; j<(fPatchGrid[type][i].GetSize()-stepSize+1); j+=stepSize) { //patches
299  for(Int_t k = 0; k<GetNColTriggerPatches(type,GetPatchDim(i),i); k+=stepSize) {
300  for(Int_t l = 0; l<GetNRowTriggerPatches(type,GetPatchDim(i),i); l+=stepSize) {
301  Int_t id = GetTriggerPatchID(k,l,type,GetPatchDim(i),i);
302  // if(type==1 && i==4) Printf("id: %d/%d k: %d/%d l: %d/%d",id,fPatchGrid[type][i].GetSize(),k,GetNColTriggerPatches(type,GetPatchDim(i),i),l,GetNRowTriggerPatches(type,GetPatchDim(i),i));
303  if(fPatchGrid[type][i].At(id)>0.) { //don't do anything with empty patches
304  Double_t sub = medL1[fPatchSub-1][subType]*GetPatchArea(i);
305  fPatchECorr[type][i]->Fill(fPatchGrid[type][i].At(id) - sub);
306  fPatchECorrPar[type][i]->Fill(fPatchGrid[type][i].At(id) - fRhoMean*GetPatchArea(i));
307 
308  //Bookkeep leading patches
309  if((fPatchGrid[type][i].At(id)-sub)>Elead[i][type]) {
310  EleadID[i][type] = id;
311  Elead[i][type] = fPatchGrid[type][i].At(id)-sub;
312  }
313  if(fPatchGrid[type][i].At(id)>EleadRaw[i][type])
314  EleadRaw[i][type] = fPatchGrid[type][i].At(id);
315  }
316  }//cols
317  }//rows
318  }//type
319 
320  AliJetContainer *cont = GetJetContainer(0);
321  for(Int_t k = 0; k<2; k++) {
322  Int_t subType = 1;
323  if(k==1) subType=0;
324  fPatchECorrRho[k][i]->Fill(Elead[i][k],medL1[fPatchSub-1][subType]);
325  fPatchECorrECorrRho[k][i]->Fill(Elead[i][k],Elead[i][subType],medL1[fPatchSub-1][subType]);
326  fPatchERaw[k][i]->Fill(EleadRaw[i][k]);
327  //Save jet spectra within EMCal and DCal fiducial acceptance
328  //jet pT vs highest energy patch for preferred trigger patch types
329  if(cont) {
330  Double_t r = cont->GetJetRadius();
331  cont->SetJetEtaLimits(fEtaMin[k]+r,fEtaMax[k]-r);
332  cont->SetJetPhiLimits(fPhiMin[k]+r,fPhiMax[k]-r);
333  Double_t rho = cont->GetRhoVal();
334  AliEmcalJet *jet = NULL;
335  cont->ResetCurrentID();
336  while((jet = cont->GetNextAcceptJet())) {
337  Double_t jetpt = jet->Pt() - rho*jet->Area();
338  fh2JetPtPatchECorr[k][i]->Fill(jetpt,Elead[i][k]);
339  }
340  }
341  }
342  Double_t eta = 0.; Double_t phi = 0.;
343  GetEtaPhiFromTriggerPatchID(EleadID[i][0],0,GetPatchDim(i),1,eta,phi);
344  fh2PatchEtaPhiEmcal[i]->Fill(eta,phi);
345 
346  GetEtaPhiFromTriggerPatchID(EleadID[i][1],1,GetPatchDim(i),1,eta,phi);
347  fh2PatchEtaPhiDcal[i]->Fill(eta,phi);
348  } //patch types
349 
350  fMultVsRho->Fill(GetParticleContainer(0)->GetNParticles(),medL1[3][0]);
351  return kTRUE;
352 }
353 
354 //________________________________________________________________________
356 
357  AliJetContainer *cont = GetJetContainer(icont);
358  if(!cont) return NULL;
359 
360  Double_t closest_dr = 1e6;
361  Int_t closest_id = -1;
362  AliEmcalJet *jet = NULL;
363  cont->ResetCurrentID();
364  while((jet = cont->GetNextAcceptJet())) {
365  Double_t dphi = jet->Phi() - phi;
366  Double_t deta = jet->Eta() - eta;
367  dphi = TVector2::Phi_mpi_pi(dphi);
368  Double_t dr = TMath::Sqrt ( dphi * dphi + deta * deta );
369  if(dr<closest_dr) {
370  closest_dr = dr;
371  closest_id = cont->GetCurrentID();
372  }
373  }
374  jet = cont->GetJet(closest_id);
375  return jet;
376 }
377 
378 //________________________________________________________________________
380  //calc total energy of all patches
381  Int_t n = fPatchGrid[0][patchType].GetSize();
382  if(n<1) return -1.;
383 
384  Double_t sum = 0.;
385  Int_t count = 0;
387  for(Int_t type = 0; type<2; type++) {
388  for(Int_t i = 0; i<fPatchGrid[type][patchType].GetSize(); i+=stepSize) {
389  if(fPatchGrid[type][patchType].At(i)>0.) count++;
390  sum+=fPatchGrid[type][patchType].At(i);
391  }
392  }
393  return sum;
394 }
395 
396 //________________________________________________________________________
397 Double_t AliEmcalPicoTrackInGridMaker::CalculateMedian(const Int_t patchType, const Int_t type, const Int_t areaType) {
398  //areaType:
399  //0: passive area
400  //1: active area mini patches
401  //2: active arear cells
402 
403  Int_t n = fPatchGrid[type][patchType].GetSize();
404  if(n<1) return -1.;
405 
406  Int_t level = 0;
407  if(patchType>0) level = 1;
408  Int_t dim = GetPatchDim(patchType);
409  Double_t area = GetPatchArea(patchType);
410  Int_t stepSize = GetTriggerPatchIdStepSizeNoOverlap(GetPatchDim(patchType),level);
411  //Printf("patchType: %d dim: %d stepSizeNoOverlap: %d ",patchType,GetPatchDim(patchType),stepSize);
412 
413  static Double_t arr[999];
414  Int_t c = 0;
415 
416  //find patch with highest energy
417  Int_t imax = -1;
418  Double_t max = 0.;
419  for(Int_t i = 0; i<(GetNColTriggerPatches(type,dim,patchType)); i+=stepSize) {
420  for(Int_t j = 0; j<(GetNRowTriggerPatches(type,dim,patchType)); j+=stepSize) {
421  Int_t id = GetTriggerPatchID(i,j,type,dim,patchType);
422  // Printf("id: %d/%d i: %d/%d j:%d/%d",id,fPatchGrid[type][patchType].GetSize(),i,GetNColTriggerPatches(type,dim,patchType),j,GetNRowTriggerPatches(type,dim,patchType));
423  if(fPatchGrid[type][patchType].At(id)>max) {
424  imax = id;
425  max = fPatchGrid[type][patchType].At(id);
426  }
427  }//cols
428  }//rows
429 
430  for(Int_t i = 0; i<GetNColTriggerPatches(type,dim,patchType); i+=stepSize) {
431  for(Int_t j = 0; j<GetNRowTriggerPatches(type,dim,patchType); j+=stepSize) {
432  Int_t id = GetTriggerPatchID(i,j,type,dim,patchType);
433  if(fExclLeadingPatch>0 && id==imax) continue;
434  if(fPatchGrid[type][patchType].At(id)>0.) {
435  Int_t active = 99;
436  if(areaType==1) active = fActiveAreaMPP[type][patchType].At(id);
437  else if(areaType==2) active = fActiveAreaCP[type][patchType].At(id);
438  if(areaType>0) area = GetPatchAreaActive(id,type,patchType,areaType-1);
439 
440  if(area>0. && active>1) {
441  arr[c] = fPatchGrid[type][patchType].At(id)/area;
442  c++;
443  }
444  if(type==0 && areaType==0) {
445  fh1RhoEmcal[patchType]->Fill(arr[c-1]);
446  fPatchEnVsActivityEmcal[patchType]->Fill(fPatchGrid[type][patchType].At(id),fActiveAreaCP[type][patchType].At(id));
447  }
448  if(type==1 && areaType==0) {
449  fh1RhoDcal[patchType]->Fill(arr[c-1]);
450  fPatchEnVsActivityDcal[patchType]->Fill(fPatchGrid[type][patchType].At(id),fActiveAreaCP[type][patchType].At(id));
451  }
452  }
453  }
454  }
455  Double_t med = TMath::Median(c,arr);
456  return med;
457 }
458 
459 //________________________________________________________________________
461  //create cells from track input
462  if(!InitCells()) return kFALSE;
463 
464  AliVParticle *track = NULL;
466  if(!trackCont) return kFALSE;
467  trackCont->ResetCurrentID();
468  while((track = trackCont->GetNextAcceptParticle())) {
469  if(track->Pt()<fMinCellE) continue;
470  Int_t id = GetGridID(track);
471  Int_t type = GetCellType(track);
472  if(id>-1)
473  fCellGrid[type].AddAt(fCellGrid[type].At(id)+track->Pt(),id);
474  }
475  return kTRUE;
476 }
477 
478 //________________________________________________________________________
480  //create mini patches (2x2 cells)
481  if(!InitMiniPatches()) return kFALSE;
482 
483  //loop over edges of mini patches
484  Int_t nm = 0; //mini patch number
485  for(Int_t type = 0; type<2; type++) {
486  nm = 0;
487  for(Int_t i = 0; i<(GetNCellsCol(type)-1); i+=2) {
488  for(Int_t j = 0; j<(GetNCellsRow(type)-1); j+=2) {
489  //loop over cells in mini patch
490  for(Int_t k = 0; k<2; k++) {
491  for(Int_t l = 0; l<2; l++) {
492  Int_t id = GetGridID(i+k,j+l,type);
493  fMiniPatchGrid[type].AddAt(fMiniPatchGrid[type].At(nm)+fCellGrid[type].At(id),nm);
494  if(fCellGrid[type].At(id)>0.)
495  fActiveAreaMP[type].AddAt(fActiveAreaMP[type].At(nm)+1,nm);
496  }
497  }
498  nm++;
499  }
500  }
501  }
502  return kTRUE;
503 }
504 
505 //________________________________________________________________________
507  //returns number of rows of mini patches in detector of type X (0: EMCal 1: DCal)
508  Int_t nRows = TMath::FloorNint(0.5*GetNCellsCol(type));
509  return nRows;
510 }
511 
512 //________________________________________________________________________
514  //returns number of rows of mini patches in detector of type X (0: EMCal 1: DCal)
515  Int_t nCols = TMath::FloorNint(0.5*GetNCellsRow(type));
516  return nCols;
517 }
518 
519 //________________________________________________________________________
521  //create trigger patches
522  if(!InitPatches(dim,level)) return kFALSE;
523 
524  Int_t pt = GetPatchType(dim,level);
525  if(pt<0) return kFALSE;
526  Int_t nm = (Int_t)(dim/2.); //size of trigger patch in number of mini patches
527  Int_t stepm = (Int_t)(dim/2.); //step size through grid in mini patches
528  if(level==1 && fL1Slide) stepm = GetSlidingStepSizeMiniPatches(dim,level);
529  //loop over edges of mini patches
530  for(Int_t type = 0; type<2; type++) {
531  Int_t np = 0; //patch number
532  for(Int_t j = 0; j<=(GetNColMiniPatches(type)-nm); j+=stepm) {
533  for(Int_t i = 0; i<=(GetNRowMiniPatches(type)-nm); i+=stepm) {
534  // for(Int_t j = 0; j<=(GetNColMiniPatches(type)-nm); j+=stepm) {
535  //loop over mini patches in patch
536  for(Int_t k = 0; k<nm; k++) {
537  for(Int_t l = 0; l<nm; l++) {
538  Int_t row = i+k;
539  Int_t col = j+l;
540  Int_t id = GetMiniPatchID(row,col,type);
541  fPatchGrid[type][pt].AddAt(fPatchGrid[type][pt].At(np)+fMiniPatchGrid[type].At(id),np);
542  if(fMiniPatchGrid[type].At(id)>0.) {
543  fActiveAreaMPP[type][pt].AddAt(fActiveAreaMPP[type][pt].At(np)+1,np);
544  fActiveAreaCP[type][pt].AddAt(fActiveAreaCP[type][pt].At(np)+fActiveAreaMP[type].At(id),np);
545  }
546  }
547  }
548  np++;
549  }
550  }
551  }
552  return kTRUE;
553 }
554 
555 //________________________________________________________________________
556 Int_t AliEmcalPicoTrackInGridMaker::GetTriggerPatchID(const Int_t row, const Int_t col, const Int_t type, const Int_t dim, const Int_t patchType) const {
557  Int_t id = row*GetNRowTriggerPatches(type,dim,patchType) + col;
558  return id;
559 }
560 
561 //________________________________________________________________________
562 Int_t AliEmcalPicoTrackInGridMaker::GetMiniPatchID(const Int_t row, const Int_t col, const Int_t type) const {
563  Int_t id = row*GetNColMiniPatches(type) + col;
564  return id;
565 }
566 
567 //________________________________________________________________________
569  //cell in EMCal (0) or DCal (1)
570  for(Int_t i = 0; i<2; i++) {
571  if(eta>fEtaMin[i] && eta<fEtaMax[i] && phi>fPhiMin[i] && phi<fPhiMax[i]) return i;
572  }
573  return -1;
574 }
575 
576 //________________________________________________________________________
577 Int_t AliEmcalPicoTrackInGridMaker::GetGridID(const Int_t row, const Int_t col, const Int_t type) const {
578  Int_t id = row*GetNCellsRow(type) + col;
579  return id;
580 }
581 
582 //________________________________________________________________________
584 
585  Int_t type = GetCellType(eta,phi);
586  if(type<0 || type>1) return -1; //position is not in EMCal or DCal
587 
588  // grid ID convention:
589  // upper left corner (min phi, min eta) is first ID
590  // then walk through grid from upper left to lower right accross the rows in phi
591 
592  Int_t id = -1;
593  Double_t etaRel = eta-fEtaMin[type];
594  Double_t phiRel = phi-fPhiMin[type];
595  Int_t row = TMath::FloorNint(etaRel/fCellSize);
596  Int_t col = TMath::FloorNint(phiRel/fCellSize);
597  id = GetGridID(row,col,type);
598 
599  if(id>=fNCells) {
600  Printf("Got too large id %d %d type: %d",id,fNCells,type);
601  Printf("eta: %f phi: %f",eta,phi);
602  Printf("etaRel: %f phiRel: %f",etaRel,phiRel);
603  Printf("row: %d col: %d -> %d + %d = %d",row,col,row*GetNCellsRow(type) + col,fNCellsEMCal,id);
604  Printf("n cells row: %d",GetNCellsRow(type));
605  Printf("\n");
606  }
607  return id;
608 }
609 
610 //________________________________________________________________________
611 void AliEmcalPicoTrackInGridMaker::GetEtaPhiFromGridID(const Int_t id, const Int_t type, Double_t &eta, Double_t &phi) const {
612  //returns eta phi of cell at lower right edge (lowest eta, lowest phi)
613  Int_t row = TMath::FloorNint(id/GetNCellsRow(type));
614  Int_t col = id - row * GetNCellsRow(type);
615  eta = fEtaMin[type] + row*fCellSize;
616  phi = fPhiMin[type] + col*fCellSize;
617  AliDebug(2,Form("id: %d type: %d row: %d col: %d eta: %f phi: %f",id,type,row,col,eta,phi));
618 }
619 
620 //________________________________________________________________________
622  //returns eta phi of mini patch at lower right edge (lowest eta, lowest phi)
623  Int_t row = TMath::FloorNint(id/GetNColMiniPatches(type));
624  Int_t col = id - row * GetNColMiniPatches(type);
625  eta = fEtaMin[type] + row*2.*fCellSize;
626  phi = fPhiMin[type] + col*2.*fCellSize;
627 }
628 
629 //________________________________________________________________________
630 void AliEmcalPicoTrackInGridMaker::GetEtaPhiFromTriggerPatchID(const Int_t id, const Int_t type, const Int_t dim, const Int_t level, Double_t &eta, Double_t &phi) const {
631  //returns eta phi of mini patch at lower right edge (lowest eta, lowest phi)
632  Int_t step = GetSlidingStepSizeCells(dim,level); //id: 8/96 k(row): 0/8 l(col): 8/12
633  // Int_t offset = dim/step;
634  Int_t row = TMath::FloorNint(id/GetNColTriggerPatches(type,dim,GetPatchType(dim,level)));
635  Int_t col = id - row * GetNColTriggerPatches(type,dim,GetPatchType(dim,level));
636  eta = fEtaMin[type] + row*step*fCellSize;
637  phi = fPhiMin[type] + col*step*fCellSize;
638  // if(dim==32) {
639  // Printf("dim: %d step: %d offset: %d",dim,step,offset);
640  // Printf("dim: %d id: %d row: %d col: %d eta: %f phi: %f",dim,id,row,col,eta,phi);
641  // }
642 }
643 
644 //________________________________________________________________________
645 Int_t AliEmcalPicoTrackInGridMaker::GetNColTriggerPatches(const Int_t type, const Int_t dim, const Int_t patchType) const {
646  //returns number of trigger patch columns
647  Int_t level = 0;
648  if(patchType>0) level = 1;
649  Double_t stepmp = (Double_t)GetSlidingStepSizeMiniPatches(dim,level);
650  Int_t nmp = GetNColMiniPatches(type);
651  Int_t ntc = TMath::FloorNint(nmp/stepmp);
652  // Printf("dim: %d stepmp: %f nmp: %d ntc: %d",dim,stepmp,nmp,ntc);
653  return ntc;
654 }
655 
656 //________________________________________________________________________
657 Int_t AliEmcalPicoTrackInGridMaker::GetNRowTriggerPatches(const Int_t type, const Int_t dim, const Int_t patchType) const {
658  //returns number of trigger patch rows
659  Int_t level = 0;
660  if(patchType>0) level = 1;
661  Double_t stepmp = (Double_t)GetSlidingStepSizeMiniPatches(dim,level);
662  Int_t nmp = GetNRowMiniPatches(type);
663  Int_t ntr = TMath::FloorNint(nmp/stepmp);
664  return ntr;
665 }
666 
667 //________________________________________________________________________
669  //returns number of cells in column
670  Double_t deta = fEtaMax[type] - fEtaMin[type];
671  Int_t nCellsCol = TMath::FloorNint(deta/fCellSize);
672  return nCellsCol;
673 }
674 
675 //________________________________________________________________________
677  //returns number of cells in row
678  Double_t dPhi = fPhiMax[type] - fPhiMin[type];
679  Int_t nCellsRow = TMath::FloorNint(dPhi/fCellSize);
680  return nCellsRow;
681 }
682 
683 //________________________________________________________________________
685  //initialize cells array
686  CheckEdges();
687  if(!CheckEdges()) return kFALSE;
688 
689  //number of cells in EMCal acceptance
690  Int_t nCellsPhiE = GetNCellsCol(0);
691  Int_t nCellsEtaE = GetNCellsRow(0);
692  fNCellsEMCal = nCellsPhiE*nCellsEtaE;
693 
694  //number of cells in DCal acceptance
695  Int_t nCellsPhiD = GetNCellsCol(1);
696  Int_t nCellsEtaD = GetNCellsRow(1);
697  fNCellsDCal = nCellsPhiD*nCellsEtaD;
698 
699  //total number of cells
701 
702  AliDebug(2,Form("EMCal: %d x %d",nCellsEtaE,nCellsPhiE));
703  AliDebug(2,Form("DCal: row: %d x col: %d",nCellsEtaD,nCellsPhiD));
704  AliDebug(2,Form("fNCells: %d fNCellsE: %d fnCellsD: %d",fNCells,fNCellsEMCal,fNCellsDCal));
705 
706  fCellGrid[0].Set(fNCellsEMCal);
707  fCellGrid[1].Set(fNCellsDCal);
708  fCellGrid[0].Reset(0);
709  fCellGrid[1].Reset(0);
710  return kTRUE;
711 }
712 
713 //________________________________________________________________________
715  //initialize mini patch array
716  if(fCellGrid[0].GetSize()<0) return kFALSE;
717  if(fCellGrid[1].GetSize()<0) return kFALSE;
718  Double_t conv = 0.25; //dimension of mini patch is 2x2 cells
719  Int_t nMiniPatches[2];
720  nMiniPatches[0] = (Int_t)(conv*fNCellsEMCal);
721  nMiniPatches[1] = (Int_t)(conv*fNCellsDCal);
722  for(Int_t i = 0; i<2; i++) {
723  fMiniPatchGrid[i].Set(nMiniPatches[i]);
724  fMiniPatchGrid[i].Reset(0.);
725 
726  fActiveAreaMP[i].Set(nMiniPatches[i]);
727  fActiveAreaMP[i].Reset(0);
728  }
729  return kTRUE;
730 }
731 
732 //________________________________________________________________________
733 Int_t AliEmcalPicoTrackInGridMaker::GetNTriggerPatches(const Int_t type, const Int_t dim, const Int_t level) const {
734  //get number of trigger patches in EMCal or DCal
735  Double_t dimd = (Double_t)dim;
736  Double_t conv = 1./(dimd*dimd);
737  if(level==1) {
739  conv = 1./(step*step);
740  }
741  Int_t nPatches = 0;
742  if(type==0) nPatches = (Int_t)(conv*fNCellsEMCal);
743  else if(type==1) nPatches = (Int_t)(conv*fNCellsDCal);
744  return nPatches;
745 }
746 
747 //________________________________________________________________________
749  // dimensions in cell units
750  // if level==1: sliding window will be applied
751  // L1 4x4: slide by 2 cells: 1 mini patch
752  // L1 8x8: slide by 4 cells: 2 mini patches
753  // L1 16x16: slide by 4 cells: 2 mini patches
754  // L1 32x32: slide by 8 cells: 4 mini patches
755 
756  if(fCellGrid[0].GetSize()<0) return kFALSE;
757  if(fCellGrid[1].GetSize()<0) return kFALSE;
758 
759  Int_t type = GetPatchType(dim,level);
760  if(type<0 || type>4) return kFALSE;
761 
762  Int_t nPatches[2]; //number of trigger patches in EMCal and DCal
763  for(Int_t i = 0; i<2; i++)
764  nPatches[i] = GetNTriggerPatches(i,dim,level);
765  //total number of trigger patches
766  Int_t nPatchesT = nPatches[0] + nPatches[1];
767 
768  fNPatchesEMCal[type] = nPatches[0];
769  AliDebug(2,Form("Create trigger patch of type %d with dim %d and with %d patches EMCAL: %d DCAL: %d",type,dim,nPatchesT,nPatches[0],nPatches[1]));
770  //Printf("Create trigger patch of type %d with dim %d and with %d patches EMCAL: %d DCAL: %d",type,dim,nPatchesT,nPatches[0],nPatches[1]);
771  for(Int_t i = 0; i<2; i++) {
772  fPatchGrid[i][type].Set(nPatches[i]);
773  fPatchGrid[i][type].Reset(0);
774  fActiveAreaMPP[i][type].Set(nPatches[i]);
775  fActiveAreaMPP[i][type].Reset(0);
776  fActiveAreaCP[i][type].Set(nPatches[i]);
777  fActiveAreaCP[i][type].Reset(0);
778  }
779 
780  return kTRUE;
781 }
782 
783 //________________________________________________________________________
785  //returns total area of patch
786  Int_t ncell = 4;
787  if(ipatch==0) ncell = 4;
788  if(ipatch==1) ncell = 4;
789  if(ipatch==2) ncell = 8;
790  if(ipatch==3) ncell = 16;
791  if(ipatch==4) ncell = 32;
792 
793  return ncell;
794 }
795 
796 //________________________________________________________________________
798  //returns total area of patch
799  Double_t ncell = (Double_t)GetPatchDim(ipatch);
800  Double_t area = ncell*ncell*fCellSize*fCellSize;
801  return area;
802 }
803 
804 //________________________________________________________________________
805 Double_t AliEmcalPicoTrackInGridMaker::GetPatchAreaActive(const Int_t id, const Int_t type, const Int_t ipatch, const Int_t atype) const {
806  //atype = 0 : active area from mini patches
807  //atype = 1 : active area from cells
808 
809  Int_t active = 0;
810  if(atype==0) active = fActiveAreaMPP[type][ipatch].At(id);
811  else if(atype==1) active = fActiveAreaCP[type][ipatch].At(id);
812  else return -1;
813 
814  Double_t fac = 1.;
815  if(atype==0) fac = 4.;
816  Double_t area = active*fac*fCellSize*fCellSize;
817  return area;
818 }
819 
820 //________________________________________________________________________
822  //get step size for mock-up L1 trigger
823  if(!fL1Slide || level==0) return dim;
824 
825  if(dim==4) return 2;
826  else if(dim==8) return 4;
827  else if(dim==16) return 4;
828  else if(dim==32) return 8;
829  else return -1;
830 }
831 
832 //________________________________________________________________________
834  //returns step size in mini patches
835  Double_t step = (Double_t)GetSlidingStepSizeCells(dim,level);
836  if(step<0) return step;
837  Int_t stepMiniPatch = (Int_t)(step/2.);
838  return stepMiniPatch;
839 }
840 
841 //________________________________________________________________________
843  //return step for trigger patch id's without overlapping patches
844  if(!fL1Slide) return 1;
845 
846  Int_t cellStep = GetSlidingStepSizeCells(dim,level);
847  Int_t step = TMath::FloorNint((Double_t)(dim/cellStep));
848  return step;
849 }
850 
851 //________________________________________________________________________
853  //type of trigger patch (size)
854  Int_t type = -1;
855  if(level==0 && dim==4) type = 0;
856  if(level==1) {
857  if(dim==4) type = 1;
858  if(dim==8) type = 2;
859  if(dim==16) type = 3;
860  if(dim==32) type = 4;
861  }
862  return type;
863 }
864 
865 //________________________________________________________________________
867  //Check if defined edges of EMCal and DCal make sense
868  if(fPhiMin[0]<0. || fPhiMax[0]<fPhiMin[0]) {
869  AliDebug(11,Form("EMCal phi edges not defined %f-%f",fPhiMin[0],fPhiMax[0]));
870  return kFALSE;
871  }
872 
873  if(fPhiMin[1]<0. || fPhiMax[1]<fPhiMin[1]) {
874  AliDebug(11,Form("DCal phi edges not defined %f-%f",fPhiMin[1],fPhiMax[1]));
875  return kFALSE;
876  }
877 
878  if(fEtaMin[0]<-10. || fEtaMax[0]<fEtaMin[1]) {
879  AliDebug(11,Form("EMCal eta edges not well defined %f-%f",fEtaMin[0],fEtaMax[0]));
880  return kFALSE;
881  }
882 
883  if(fEtaMin[1]<-10. || fEtaMax[1]<fEtaMin[1]) {
884  AliDebug(11,Form("DCal eta edges not well defined %f-%f",fEtaMin[1],fEtaMax[1]));
885  return kFALSE;
886  }
887 
888  for(Int_t type = 0; type<2; type++) {
889  Double_t dphi = fPhiMax[type] - fPhiMin[type];
890  Double_t nPatchPhi32 = TMath::Floor(dphi/(32.*fCellSize));
891  Double_t nCellsColExact = nPatchPhi32 * 32. ;
892 
893  Double_t deta = fEtaMax[type] - fEtaMin[type];
894  Double_t nPatchEta32 = TMath::Floor(deta/(32.*fCellSize));
895  Double_t nCellsRowExact = nPatchEta32 * 32. ;
896 
897  Double_t col_extra = dphi/fCellSize - TMath::Floor(nCellsColExact);
898  Double_t phi_extra = col_extra*fCellSize;
899  fPhiMin[type] += phi_extra*0.5;
900  fPhiMax[type] -= phi_extra*0.5;
901 
902  Double_t row_extra = deta/fCellSize - TMath::Floor(nCellsRowExact);
903  Double_t eta_extra = row_extra*fCellSize;
904  fEtaMin[type] += eta_extra*0.5;
905  fEtaMax[type] -= eta_extra*0.5;
906  AliDebug(2,Form("type: %d exact: col: %f row: %f",type,nCellsColExact,nCellsRowExact));
907  // Printf("type: %d exact: col: %f row: %f",type,nCellsColExact,nCellsRowExact);
908  // PrintAcceptance();
909  }
910  return kTRUE;
911 }
912 
913 //________________________________________________________________________
915  Printf("EMCal");
916  Printf("phi: %f-%f",fPhiMin[0],fPhiMax[0]);
917  Printf("eta: %f-%f",fEtaMin[0],fEtaMax[0]);
918 
919  Printf("DCal");
920  Printf("phi: %f-%f",fPhiMin[1],fPhiMax[1]);
921  Printf("eta: %f-%f",fEtaMin[1],fEtaMax[1]);
922 }
void GetEtaPhiFromMiniPatchID(const Int_t id, const Int_t type, Double_t &eta, Double_t &phi) const
Double_t Area() const
Definition: AliEmcalJet.h:130
virtual AliVParticle * GetNextAcceptParticle()
Bool_t InitPatches(const Int_t dim, const Int_t level)
Double_t GetRhoVal() const
TH1F * fPatchECorrPar[2][5]
corrected patch energy for EMCal and DCal
double Double_t
Definition: External.C:58
Definition: External.C:260
Definition: External.C:236
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:121
Double_t Phi() const
Definition: AliEmcalJet.h:117
void GetEtaPhiFromGridID(const Int_t id, const Int_t type, Double_t &eta, Double_t &phi) const
Int_t GetTriggerPatchID(const Int_t row, const Int_t col, const Int_t type, const Int_t dim, const Int_t patchType) const
TProfile * fpMedianTypeDcal[3]
median vs patch type for 3 types of area calculation
Bool_t CreateGridPatches(const Int_t dim, const Int_t level)
TCanvas * c
Definition: TestFitELoss.C:172
Int_t GetMiniPatchID(const Int_t row, const Int_t col, const Int_t type) const
Int_t GetPatchDim(const Int_t ipatch) const
Int_t GetNColTriggerPatches(const Int_t type, const Int_t dim, const Int_t patchType) const
TH1F * fh1RhoEmcal[5]
median vs patch type for 3 types of area calculation
Container for particles within the EMCAL framework.
TH2F * fPatchEnVsActivityDcal[5]
patch energy vs active cells
Int_t GetNRowTriggerPatches(const Int_t type, const Int_t dim, const Int_t patchType) const
Int_t GetSlidingStepSizeMiniPatches(const Int_t dim, const Int_t level=1) const
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Int_t GetNParticles(Int_t i=0) const
Get number of particles in container attached to this task with index i.
TH2F * fMultVsRho
jet pt vs leading patch energy
int Int_t
Definition: External.C:63
void SetJetPhiLimits(Float_t min, Float_t max)
TH2F * fPatchEnVsActivityEmcal[5]
rho distributions for passive area
Double_t CalculateMedian(const Int_t patchType, const Int_t type, const Int_t areaType=0)
Int_t GetNCellsCol(const Int_t type) const
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.
Int_t GetNCellsRow(const Int_t type) const
Int_t GetNTriggerPatches(const Int_t type, const Int_t dim, const Int_t level) const
AliEmcalJet * GetClosestJet(const Double_t eta, const Double_t phi, const Int_t icont=0) const
TH1F * fPatchECorr[2][5]
patch energy vs active cells
Int_t GetSlidingStepSizeCells(const Int_t dim, const Int_t level=1) const
AliEmcalJet * GetNextAcceptJet()
TH2F * fh2PatchEtaPhiEmcal[5]
Ecorr,det1 vs Ecorr,det2 vs rho,det2 opposite side for dijet in acceptance like events.
Int_t GetNRowMiniPatches(const Int_t type) const
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Double_t Pt() const
Definition: AliEmcalJet.h:109
TH3F * fPatchECorrECorrRho[2][5]
corrected patch energy vs rho opposite side
Float_t GetJetRadius() const
AliEmcalList * fOutput
!output list
Int_t GetCellType(const Double_t eta, const Double_t phi) const
TH2F * fh2MedianTypeDcal[3]
median vs patch type for 3 types of area calculation
Double_t GetPatchAreaActive(const Int_t id, const Int_t type, const Int_t ipatch, const Int_t atype) const
void SetMakeGeneralHistograms(Bool_t g)
TH1F * fPatchERaw[2][5]
corrected patch energy with inclusive mean rho for EMCal and DCal
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
Double_t GetPatchArea(const Int_t ipatch) const
TProfile * fpMedianTypeEmcal[3]
median vs patch type for 3 types of area calculation
void UserCreateOutputObjects()
Main initialization function on the worker.
TH1F * fh1RhoDcal[5]
rho distributions for passive area
bool Bool_t
Definition: External.C:53
Int_t GetGridID(const AliVParticle *vp) const
TH2F * fh2JetPtPatchECorr[2][5]
patch positions in DCal
void SetJetEtaLimits(Float_t min, Float_t max)
TH2F * fPatchECorrRho[2][5]
uncorrected patch energy for EMCal and DCal
Int_t GetTriggerPatchIdStepSizeNoOverlap(const Int_t dim, const Int_t level=1) const
void GetEtaPhiFromTriggerPatchID(const Int_t id, const Int_t type, const Int_t dim, const Int_t level, Double_t &eta, Double_t &phi) const
Int_t GetNColMiniPatches(const Int_t type) const
Double_t CalculateSum(const Int_t patchType) const
Container for jet within the EMCAL jet framework.
TH2F * fh2PatchEtaPhiDcal[5]
patch positions in EMCal
Int_t GetPatchType(const Int_t dim, const Int_t level) const
AliEmcalJet * GetJet(Int_t i) const