AliPhysics  29d4213 (29d4213)
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Macros
AliEmcalPatchFromCellMaker.cxx
Go to the documentation of this file.
1 // $Id$
2 //
3 // Class to put cells into trigger patches
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 "AliEMCALGeometry.h"
15 #include "AliEmcalTriggerPatchInfoAP.h"
16 #include "AliEmcalTriggerBitConfig.h"
17 
19 
21 
22 //________________________________________________________________________
25  fCaloTriggersOutName("EmcalPatches32x32"),
26  fCaloTriggersOut(0),
27  fPatchDim(32),
28  fMinCellE(0.05),
29  fCellTimeMin(485e-9),
30  fCellTimeMax(685e-9),
31  fL1Slide(0),
32  fTriggerBitConfig(0x0),
33  fh3EEtaPhiCell(0),
34  fh2CellEnergyVsTime(0),
35  fh1CellEnergySum(0)
36 {
37  // Constructor.
38  for (Int_t i = 0; i < kPatchCols; i++) {
39  for (Int_t j = 0; j < kPatchRows; j++) {
40  fPatchADCSimple[i][j] = 0.;
41  fPatchESimple[i][j] = 0.;
42  }
43  }
44 
45  SetMakeGeneralHistograms(kTRUE);
46 }
47 
48 //________________________________________________________________________
50  AliAnalysisTaskEmcal(name,kTRUE),
51  fCaloTriggersOutName("EmcalPatches32x32"),
52  fCaloTriggersOut(0),
53  fPatchDim(32),
54  fMinCellE(0.05),
55  fCellTimeMin(485e-9),
56  fCellTimeMax(685e-9),
57  fL1Slide(0),
58  fTriggerBitConfig(0x0),
59  fh3EEtaPhiCell(0),
60  fh2CellEnergyVsTime(0),
61  fh1CellEnergySum(0)
62 {
63  // Constructor.
64  for (Int_t i = 0; i < kPatchCols; i++) {
65  for (Int_t j = 0; j < kPatchRows; j++) {
66  fPatchADCSimple[i][j] = 0.;
67  fPatchESimple[i][j] = 0.;
68  }
69  }
70 
72 }
73 
74 //________________________________________________________________________
76 {
77  // Destructor.
78 }
79 
80 //________________________________________________________________________
82 {
83  // Init the analysis.
84 
86 
87  if (!fInitialized)
88  return;
89 
90  if (!fCaloTriggersOutName.IsNull()) {
91  fCaloTriggersOut = new TClonesArray("AliEmcalTriggerPatchInfo");
93 
94  if (!(InputEvent()->FindListObject(fCaloTriggersOutName))) {
95  InputEvent()->AddObject(fCaloTriggersOut);
96  }
97  else {
98  fInitialized = kFALSE;
99  AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fCaloTriggersOutName.Data()));
100  return;
101  }
102  }
103 
104  if(!fTriggerBitConfig)
105  fTriggerBitConfig = new AliEmcalTriggerBitConfigNew();
106 
107 }
108 
109 //________________________________________________________________________
111 {
112  // Create user output.
113 
115 
116  Int_t fgkNPhiBins = 18*8;
117  Float_t kMinPhi = 0.;
118  Float_t kMaxPhi = 2.*TMath::Pi();
119  Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
120  for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
121 
122  Int_t fgkNEtaBins = 100;
123  Float_t fgkEtaMin = -1.;
124  Float_t fgkEtaMax = 1.;
125  Double_t *binsEta=new Double_t[fgkNEtaBins+1];
126  for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
127 
128  Int_t fgkNTimeBins = 600;
129  Float_t kMinTime = -200.;
130  Float_t kMaxTime = 1000;
131  Double_t *binsTime = new Double_t[fgkNTimeBins+1];
132  for(Int_t i=0; i<=fgkNTimeBins; i++) binsTime[i]=(Double_t)kMinTime + (kMaxTime-kMinTime)/fgkNTimeBins*(Double_t)i ;
133 
134  Double_t enBinEdges[3][2];
135  enBinEdges[0][0] = 1.; //10 bins
136  enBinEdges[0][1] = 0.1;
137  enBinEdges[1][0] = 5.; //8 bins
138  enBinEdges[1][1] = 0.5;
139  enBinEdges[2][0] = 100.;//95 bins
140  enBinEdges[2][1] = 1.;
141 
142  const Float_t enmin1 = 0;
143  const Float_t enmax1 = enBinEdges[0][0];
144  const Float_t enmin2 = enmax1 ;
145  const Float_t enmax2 = enBinEdges[1][0];
146  const Float_t enmin3 = enmax2 ;
147  const Float_t enmax3 = enBinEdges[2][0];//fgkEnMax;
148  const Int_t nbin11 = (int)((enmax1-enmin1)/enBinEdges[0][1]);
149  const Int_t nbin12 = (int)((enmax2-enmin2)/enBinEdges[1][1])+nbin11;
150  const Int_t nbin13 = (int)((enmax3-enmin3)/enBinEdges[2][1])+nbin12;
151 
152  Int_t fgkNEnBins=nbin13;
153  Double_t *binsEn=new Double_t[fgkNEnBins+1];
154  for(Int_t i=0; i<=fgkNEnBins; i++) {
155  if(i<=nbin11) binsEn[i]=(Double_t)enmin1 + (enmax1-enmin1)/nbin11*(Double_t)i ;
156  if(i<=nbin12 && i>nbin11) binsEn[i]=(Double_t)enmin2 + (enmax2-enmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
157  if(i<=nbin13 && i>nbin12) binsEn[i]=(Double_t)enmin3 + (enmax3-enmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
158  }
159 
160  fh3EEtaPhiCell = new TH3F("fh3EEtaPhiCell","fh3EEtaPhiCell;E_{cell};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
161  fOutput->Add(fh3EEtaPhiCell);
162 
163  fh2CellEnergyVsTime = new TH2F("fh2CellEnergyVsTime","fh2CellEnergyVsTime;E_{cell};time",fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
165 
166  fh1CellEnergySum = new TH1F("fh1CellEnergySum","fh1CellEnergySum;E_{cell};time",fgkNEnBins,binsEn);
168 
169  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
170 
171  if(binsEn) delete [] binsEn;
172  if(binsPhi) delete [] binsPhi;
173  if(binsEta) delete [] binsEta;
174  if(binsTime) delete [] binsTime;
175 }
176 
177 //________________________________________________________________________
179 {
180  // Main loop, called for each event.
181 
182  fCaloTriggersOut->Delete();
183 
184  if (!fCaloCells) {
185  AliError(Form("Calo cells container %s not available.", fCaloCellsName.Data()));
186  return kFALSE;
187  }
188 
189  for (Int_t i = 0; i < kPatchCols; i++) {
190  for (Int_t j = 0; j < kPatchRows; j++) {
191  fPatchADCSimple[i][j] = 0.;
192  fPatchESimple[i][j] = 0.;
193  }
194  }
195 
196  if(!FillPatchADCSimple()) {
197  AliError(Form("%s Could not create simple ADC patches",GetName()));
198  return kFALSE;
199  }
200 
202 
203  Double_t sum = 0.;
204  for (Int_t i = 0; i < kPatchCols; i++) {
205  for (Int_t j = 0; j < kPatchRows; j++) {
206  sum+=fPatchESimple[i][j];
207  }
208  }
209 
210  return kTRUE;
211 }
212 
213 //________________________________________________________________________
215 {
216 
217  // fill the array for offline trigger processing
218 
219  // fill the patch ADCs from cells
220  Double_t sum = 0.;
221  Int_t nCell = fCaloCells->GetNumberOfCells();
222  for(Int_t iCell = 0; iCell < nCell; ++iCell) {
223  // get the cell info, based in index in array
224  Short_t cellId = fCaloCells->GetCellNumber(iCell);
225 
226  Double_t cellT = fCaloCells->GetCellTime(cellId);
227  Double_t amp = fCaloCells->GetAmplitude(iCell);
228  fh2CellEnergyVsTime->Fill(amp,cellT*1e9);
229 
230  //timing cuts
231  if(cellT<fCellTimeMin || cellT>fCellTimeMax) continue;
232  //energy cut
233  if(amp<fMinCellE) continue;
234  sum+=amp;
235 
236  // get position
237  Int_t absId=-1;
238  fGeom->GetFastORIndexFromCellIndex(cellId, absId);
239  Int_t globCol=-1, globRow=-1;
240  fGeom->GetPositionInEMCALFromAbsFastORIndex(absId, globCol, globRow);
241  // add
242  fPatchADCSimple[globCol][globRow] += amp/kEMCL1ADCtoGeV;
243  fPatchESimple[globCol][globRow] += amp;
244 
245  TVector3 pos;
246  fGeom->GetGlobal(cellId, pos);
247  TLorentzVector lv(pos,amp);
248  Double_t cellEta = lv.Eta();
249  Double_t cellPhi = lv.Phi();
250  if(cellPhi<0.) cellPhi+=TMath::TwoPi();
251  if(cellPhi>TMath::TwoPi()) cellPhi-=TMath::TwoPi();
252  fh3EEtaPhiCell->Fill(amp,cellEta,cellPhi);
253  }
254  fh1CellEnergySum->Fill(sum);
255 
256  return kTRUE;
257 }
258 
259 //________________________________________________________________________
261 {
262  // Runs a simple offline trigger algorithm.
263  // It creates separate patches with dimension fPatchDim
264 
265  // run the trigger algo, stepping by stepsize (in trigger tower units)
266  Int_t itrig = 0;
267  Int_t patchSize = GetDimFastor();
268  Int_t stepSize = GetSlidingStepSizeFastor();
269  Int_t maxCol = kPatchCols - patchSize;
270  // Int_t maxRow = kPatchRows - patchSize;
271  Int_t maxRow = fGeom->GetNTotalTRU()*2 - patchSize; //apparently this is the total TRU in phi ??
272  // Printf("fGeom->GetNTotalTRU(): %d %d = %d x %d ; %d x %d",fGeom->GetNTRU(),fGeom->GetNTotalTRU(),fGeom->GetNTRUEta(),fGeom->GetNTRUPhi(),fGeom->GetNModulesInTRUEta(),fGeom->GetNModulesInTRUPhi());
273 
274  for (Int_t i = 0; i <= maxCol; i += stepSize) {
275  for (Int_t j = 0; j <= maxRow; j += stepSize) {
276  // get the trigger towers composing the patch
277  Int_t adcAmp = 0;
278  Double_t enAmp = 0.;
279  // window
280  for (Int_t k = 0; k < patchSize; ++k) {
281  for (Int_t l = 0; l < patchSize; ++l) {
282  // add amplitudes
283  adcAmp += (ULong64_t)fPatchADCSimple[i+k][j+l];
284  enAmp += fPatchESimple[i+k][j+l];
285  }
286  }
287 
288  if (adcAmp == 0) {
289  AliDebug(2,"EMCal trigger patch with 0 ADC counts.");
290  continue;
291  }
292 
293  Int_t absId=-1;
294  Int_t cellAbsId[4]={-1,-1,-1,-1};
295 
296  // get low left edge (eta max, phi min)
297  fGeom->GetAbsFastORIndexFromPositionInEMCAL(i, j, absId);
298  // convert to the 4 absId of the cells composing the trigger channel
299  fGeom->GetCellIndexFromFastORIndex(absId, cellAbsId);
300  TVector3 edge1;
301  fGeom->GetGlobal(cellAbsId[0], edge1);
302 
303  // get up right edge (eta min, phi max)
304  fGeom->GetAbsFastORIndexFromPositionInEMCAL(i+patchSize-1, j+patchSize-1, absId);
305  fGeom->GetCellIndexFromFastORIndex(absId, cellAbsId);
306  TVector3 edge2;
307  fGeom->GetGlobal(cellAbsId[3], edge2);
308 
309  // get the center of the patch
310  Int_t offsetCenter = TMath::FloorNint(0.5*patchSize);
311  fGeom->GetAbsFastORIndexFromPositionInEMCAL(i+offsetCenter-1, j+offsetCenter-1, absId);
312  fGeom->GetCellIndexFromFastORIndex(absId, cellAbsId);
313  TVector3 center1;
314  fGeom->GetGlobal(cellAbsId[3], center1);
315 
316  fGeom->GetAbsFastORIndexFromPositionInEMCAL(i+offsetCenter, j+offsetCenter, absId);
317  fGeom->GetCellIndexFromFastORIndex(absId, cellAbsId);
318  TVector3 center2;
319  fGeom->GetGlobal(cellAbsId[0], center2);
320 
321  TVector3 centerGeo(center1);
322  centerGeo += center2;
323  centerGeo *= 0.5;
324 
325  // save the trigger object
326  AliEmcalTriggerPatchInfo *trigger =
327  new ((*fCaloTriggersOut)[itrig]) AliEmcalTriggerPatchInfo();
328  itrig++;
329  trigger->SetTriggerBitConfig(fTriggerBitConfig);
330  trigger->SetCenterGeo(centerGeo, enAmp);
331  trigger->SetEdge1(edge1, enAmp);
332  trigger->SetEdge2(edge2, enAmp);
333  trigger->SetADCAmp(adcAmp);
334  trigger->SetEdgeCell(i*2, j*2); // from triggers to cells
335  }
336  } // trigger algo
337  AliDebug(2,Form("Created %d trigger patches (%d) in this event",itrig,patchSize));
338 
339 }
340 
341 //________________________________________________________________________
343 
344  Int_t dim = TMath::FloorNint((Double_t)(fPatchDim/2.));
345  return dim;
346 }
347 
348 //________________________________________________________________________
350 
351  Int_t dim = GetDimFastor();
352  if(!fL1Slide) return dim;
353 
354  if(dim==2) return 2;
355  else if(dim==4) return 4;
356  else if(dim==8) return 4;
357  else if(dim==16) return 8;
358  else return -1;
359 }
360 
361 
362 
363 
AliEmcalTriggerBitConfig * fTriggerBitConfig
Double_t fPatchESimple[kPatchCols][kPatchRows]
TList * fOutput
x-section from pythia header
AliEMCALGeometry * fGeom
whether it's an ESD analysis
TH2F * fh2CellEnergyVsTime
cell E, eta, phi
AliVCaloCells * fCaloCells
clusters
Double_t fPatchADCSimple[kPatchCols][kPatchRows]
trigger array out
ClassImp(AliEmcalPatchFromCellMaker) AliEmcalPatchFromCellMaker
void SetMakeGeneralHistograms(Bool_t g)
TH1F * fh1CellEnergySum
emcal cell energy vs time