1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 //_________________________________________________________________________
19 // Class for trigger analysis.
20 // Digits are grouped in TRU's (Trigger Units). A TRU consists of 384
21 // modules ordered fNTRUPhi x fNTRUEta. The algorithm searches all possible 2x2
22 // and nxn (n is a multiple of 2) cell combinations per each TRU, adding the
23 // digits amplitude and finding the maximum. If found, look if it is isolated.
24 // Maxima are transformed in ADC time samples. Each time bin is compared to the trigger
25 // threshold until it is larger and then, triggers are set. Thresholds need to be fixed.
26 // Thresholds need to be fixed. Last 2 modules are half size in Phi, I considered
27 // that the number of TRU is maintained for the last modules but decision not taken.
28 // If different, then this must be changed.
31 // //Inside the event loop
32 // AliEMCALTrigger *tr = new AliEMCALTrigger();//Init Trigger
33 // tr->SetL0Threshold(100); //Arbitrary threshold values
34 // tr->SetL1GammaLowPtThreshold(1000);
35 // tr->SetL1GammaMediumPtThreshold(10000);
36 // tr->SetL1GammaHighPtThreshold(20000);
38 // tr->Trigger(); //Execute Trigger
39 // tr->Print(""); //Print results
41 //*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN)
42 //////////////////////////////////////////////////////////////////////////////
46 // --- ROOT system ---
52 // --- ALIROOT system ---
54 #include "AliRunLoader.h"
55 #include "AliTriggerInput.h"
57 #include "AliEMCALLoader.h"
58 #include "AliEMCALDigit.h"
59 #include "AliEMCALTrigger.h"
60 #include "AliEMCALGeometry.h"
61 #include "AliEMCALRawUtils.h"
63 #include "AliCaloConstants.h"
64 #include "AliEMCALRawResponse.h"
68 ClassImp(AliEMCALTrigger)
70 TString AliEMCALTrigger::fgNameOfJetTriggers("EMCALJetTriggerL1");
72 //______________________________________________________________________
73 AliEMCALTrigger::AliEMCALTrigger()
74 : AliTriggerDetector(), fGeom(0),
75 f2x2MaxAmp(-1), f2x2ModulePhi(-1), f2x2ModuleEta(-1),
77 fnxnMaxAmp(-1), fnxnModulePhi(-1), fnxnModuleEta(-1),
79 fADCValuesHighnxn(0),fADCValuesLownxn(0),
80 fADCValuesHigh2x2(0),fADCValuesLow2x2(0),
82 fL0Threshold(100),fL1GammaLowPtThreshold(200),
83 fL1GammaMediumPtThreshold(500), fL1GammaHighPtThreshold(1000),
84 fPatchSize(1), fIsolPatchSize(1),
85 f2x2AmpOutOfPatch(-1), fnxnAmpOutOfPatch(-1),
86 f2x2AmpOutOfPatchThres(100000), fnxnAmpOutOfPatchThres(100000),
87 fIs2x2Isol(kFALSE), fIsnxnIsol(kFALSE),
88 fSimulation(kTRUE), fIsolateInSuperModule(kTRUE), fTimeKey(kFALSE),
89 fAmpTrus(0),fTimeRtrus(0),fAmpSMods(0),
90 fTriggerPosition(6), fTriggerAmplitudes(4),
91 fNJetPatchPhi(3), fNJetPatchEta(3), fNJetThreshold(3), fL1JetThreshold(0), fJetMaxAmp(0),
92 fAmpJetMatrix(0), fJetMatrixE(0), fAmpJetMax(6,1), fVZER0Mult(0.)
95 fADCValuesHighnxn = 0x0; //new Int_t[fTimeBins];
96 fADCValuesLownxn = 0x0; //new Int_t[fTimeBins];
97 fADCValuesHigh2x2 = 0x0; //new Int_t[fTimeBins];
98 fADCValuesLow2x2 = 0x0; //new Int_t[fTimeBins];
101 // Define jet threshold - can not change from outside now
102 // fNJetThreshold = 7; // For MB Pythia suppression
103 fNJetThreshold = 10; // Hijing
104 fL1JetThreshold = new Double_t[fNJetThreshold];
105 if(fNJetThreshold == 7) {
106 fL1JetThreshold[0] = 5./0.0153;
107 fL1JetThreshold[1] = 8./0.0153;
108 fL1JetThreshold[2] = 10./0.0153;
109 fL1JetThreshold[3] = 12./0.0153;
110 fL1JetThreshold[4] = 13./0.0153;
111 fL1JetThreshold[5] = 14./0.0153;
112 fL1JetThreshold[6] = 15./0.0153;
113 } else if(fNJetThreshold == 10) {
114 Double_t thGev[10]={5.,8.,10., 12., 13.,14.,15., 17., 20., 25.};
115 for(Int_t i=0; i<fNJetThreshold; i++) fL1JetThreshold[i] = thGev[i]/0.0153;
117 fL1JetThreshold[0] = 5./0.0153;
118 fL1JetThreshold[1] = 10./0.0153;
119 fL1JetThreshold[2] = 15./0.0153;
120 fL1JetThreshold[3] = 20./0.0153;
121 fL1JetThreshold[4] = 25./0.0153;
126 fInputs.SetName("TriggersInputs");
132 //____________________________________________________________________________
133 AliEMCALTrigger::AliEMCALTrigger(const AliEMCALTrigger & trig)
134 : AliTriggerDetector(trig),
136 f2x2MaxAmp(trig.f2x2MaxAmp),
137 f2x2ModulePhi(trig.f2x2ModulePhi),
138 f2x2ModuleEta(trig.f2x2ModuleEta),
140 fnxnMaxAmp(trig.fnxnMaxAmp),
141 fnxnModulePhi(trig.fnxnModulePhi),
142 fnxnModuleEta(trig.fnxnModuleEta),
144 fADCValuesHighnxn(trig.fADCValuesHighnxn),
145 fADCValuesLownxn(trig.fADCValuesLownxn),
146 fADCValuesHigh2x2(trig.fADCValuesHigh2x2),
147 fADCValuesLow2x2(trig.fADCValuesLow2x2),
148 fDigitsList(trig.fDigitsList),
149 fL0Threshold(trig.fL0Threshold),
150 fL1GammaLowPtThreshold(trig.fL1GammaLowPtThreshold),
151 fL1GammaMediumPtThreshold(trig.fL1GammaMediumPtThreshold),
152 fL1GammaHighPtThreshold(trig.fL1GammaHighPtThreshold),
153 fPatchSize(trig.fPatchSize),
154 fIsolPatchSize(trig.fIsolPatchSize),
155 f2x2AmpOutOfPatch(trig.f2x2AmpOutOfPatch),
156 fnxnAmpOutOfPatch(trig.fnxnAmpOutOfPatch),
157 f2x2AmpOutOfPatchThres(trig.f2x2AmpOutOfPatchThres),
158 fnxnAmpOutOfPatchThres(trig.fnxnAmpOutOfPatchThres),
159 fIs2x2Isol(trig.fIs2x2Isol),
160 fIsnxnIsol(trig.fIsnxnIsol),
161 fSimulation(trig.fSimulation),
162 fIsolateInSuperModule(trig.fIsolateInSuperModule),
163 fTimeKey(trig.fTimeKey),
164 fAmpTrus(trig.fAmpTrus),
165 fTimeRtrus(trig.fTimeRtrus),
166 fAmpSMods(trig.fAmpSMods),
167 fTriggerPosition(trig.fTriggerPosition),
168 fTriggerAmplitudes(trig.fTriggerAmplitudes),
169 fNJetPatchPhi(trig.fNJetPatchPhi),
170 fNJetPatchEta(trig.fNJetPatchEta),
171 fNJetThreshold(trig.fNJetThreshold),
172 fL1JetThreshold(trig.fL1JetThreshold),
173 fJetMaxAmp(trig.fJetMaxAmp),
174 fAmpJetMatrix(trig.fAmpJetMatrix),
175 fJetMatrixE(trig.fJetMatrixE),
176 fAmpJetMax(trig.fAmpJetMax),
177 fVZER0Mult(trig.fVZER0Mult)
182 //____________________________________________________________________________
183 AliEMCALTrigger::~AliEMCALTrigger() {
188 delete [] fADCValuesHighnxn;
189 delete [] fADCValuesLownxn;
190 delete [] fADCValuesHigh2x2;
191 delete [] fADCValuesLow2x2;
193 if(fAmpTrus) {fAmpTrus->Delete(); delete fAmpTrus;}
194 if(fTimeRtrus) {fTimeRtrus->Delete(); delete fTimeRtrus;}
195 if(fAmpSMods) {fAmpSMods->Delete(); delete fAmpSMods;}
196 if(fAmpJetMatrix) delete fAmpJetMatrix;
197 if(fJetMatrixE) delete fJetMatrixE;
198 if(fL1JetThreshold) delete [] fL1JetThreshold;
201 //----------------------------------------------------------------------
202 void AliEMCALTrigger::CreateInputs()
206 // Do not create inputs again!!
207 if( fInputs.GetEntriesFast() > 0 ) return;
209 // Second parameter should be detector name = "EMCAL"
210 TString det("EMCAL"); // Apr 29, 2008
211 fInputs.AddLast( new AliTriggerInput( det+"_L0", det, 0x02) );
212 fInputs.AddLast( new AliTriggerInput( det+"_GammaHPt_L1", det, 0x04 ) );
213 fInputs.AddLast( new AliTriggerInput( det+"_GammaMPt_L1", det, 0x08 ) );
214 fInputs.AddLast( new AliTriggerInput( det+"_GammaLPt_L1", det, 0x016 ) );
215 fInputs.AddLast( new AliTriggerInput( det+"_JetHPt_L1", det, 0x032 ) );
216 fInputs.AddLast( new AliTriggerInput( det+"_JetMPt_L1", det, 0x048 ) );
217 fInputs.AddLast( new AliTriggerInput( det+"_JetLPt_L1", det, 0x064 ) );
219 if(fNJetThreshold<=0) return;
221 UInt_t level = 0x032;
222 for(Int_t i=0; i<fNJetThreshold; i++ ) {
223 TString name(GetNameOfJetTrigger(i));
224 TString title("EMCAL Jet triger L1 :"); // unused now
225 // 0.0153 - hard coded now
226 title += Form("Th %i(%5.1f GeV) :", (Int_t)fL1JetThreshold[i], fL1JetThreshold[i] * 0.0153);
227 title += Form("patch %ix%i~(%3.2f(#phi)x%3.2f(#eta)) ",
228 fNJetPatchPhi, fNJetPatchEta, 0.11*(fNJetPatchPhi), 0.11*(fNJetPatchEta));
229 fInputs.AddLast( new AliTriggerInput(name, det, UChar_t(level)) );
235 //____________________________________________________________________________
236 Bool_t AliEMCALTrigger::IsPatchIsolated(Int_t iPatchType, const TClonesArray * ampmatrixes, const Int_t iSM, const Int_t mtru, const Float_t maxamp, const Int_t maxphi, const Int_t maxeta) {
239 // EMCAL RTU size is 4modules(phi) x 24modules (eta)
240 // So maximum size of patch is 4modules x 4modules (EMCAL L0 trigger).
241 // Calculate if the maximum patch found is isolated, find amplitude around maximum (2x2 or nxn) patch,
242 // inside isolation patch . iPatchType = 0 means calculation for 2x2 patch,
243 // iPatchType = 1 means calculation for nxn patch.
244 // In the next table there is an example of the different options of patch size and isolation patch size:
245 // Patch Size (fPatchSize)
247 // fIsolPatchSize 0 2x2 (not overlap) 4x4 (overlapped)
251 if(!ampmatrixes) return kFALSE;
253 // Get matrix of TRU or Module with maximum amplitude patch.
254 Int_t itru = mtru + iSM * fGeom->GetNTRU(); //number of tru, min 0 max 3*12=36.
255 TMatrixD * ampmatrix = 0x0;
258 static int keyPrint = 0;
259 if(keyPrint) AliDebug(2,Form(" IsPatchIsolated : iSM %i mtru %i itru %i maxphi %i maxeta %i \n", iSM, mtru, itru, maxphi, maxeta));
261 if(fIsolateInSuperModule){ // ?
262 ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(iSM)) ;
263 rowborder = fGeom->GetNPhi();
264 colborder = fGeom->GetNZ();
265 AliDebug(2,"Isolate trigger in Module");
267 ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(itru)) ;
268 rowborder = fGeom->GetNModulesInTRUPhi();
269 colborder = fGeom->GetNModulesInTRUEta();
270 AliDebug(2,"Isolate trigger in TRU");
272 if(iSM>9) rowborder /= 2; // half size in phi
274 if(!ampmatrixes || !ampmatrix){
275 AliError("Could not recover the matrix with the amplitudes");
279 //Define patch modules - what is this ??
280 Int_t isolmodules = fIsolPatchSize*(1+iPatchType);
281 Int_t ipatchmodules = 2*(1+fPatchSize*iPatchType);
282 Int_t minrow = maxphi - isolmodules;
283 Int_t mincol = maxeta - isolmodules;
284 Int_t maxrow = maxphi + isolmodules + ipatchmodules;
285 Int_t maxcol = maxeta + isolmodules + ipatchmodules;
287 minrow = minrow<0?0 :minrow;
288 mincol = mincol<0?0 :mincol;
290 maxrow = maxrow>rowborder?rowborder :maxrow;
291 maxcol = maxcol>colborder?colborder :maxcol;
293 //printf("%s\n",Form("Number of added Isol Modules %d, Patch Size %d",isolmodules, ipatchmodules));
294 //printf("%s\n",Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
295 // AliDebug(2,Form("Number of added Isol Modules %d, Patch Size %d",isolmodules, ipatchmodules));
296 //AliDebug(2,Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
298 //Add amplitudes in all isolation patch
300 for(Int_t irow = minrow ; irow < maxrow; irow ++)
301 for(Int_t icol = mincol ; icol < maxcol ; icol ++)
302 amp += (*ampmatrix)(irow,icol);
304 AliDebug(2,Form("Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
307 // AliError(Form("Bad sum: Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
308 // ampmatrix->Print();
311 amp-=maxamp; //Calculate energy in isolation patch that do not comes from maximum patch.
314 AliDebug(2, Form("Maximum amplitude %f, Out of patch %f",maxamp, amp));
316 //Fill isolation amplitude data member and say if patch is isolated.
317 if(iPatchType == 0){ //2x2 case
318 f2x2AmpOutOfPatch = amp;
319 if(amp < f2x2AmpOutOfPatchThres) b=kTRUE;
320 } else if(iPatchType == 1){ //nxn case
321 fnxnAmpOutOfPatch = amp;
322 if(amp < fnxnAmpOutOfPatchThres) b=kTRUE;
325 if(keyPrint) AliDebug(2,Form(" IsPatchIsolated - OUT \n"));
332 //____________________________________________________________________________
333 void AliEMCALTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus, const Int_t isupermod,TMatrixD &max2, TMatrixD &maxn){
335 //Sums energy of all possible 2x2 (L0) and nxn (L1) modules per each TRU.
336 //Fast signal in the experiment is given by 2x2 modules,
337 //for this reason we loop inside the TRU modules by 2.
339 //Declare and initialize variables
340 Int_t nCellsPhi = fGeom->GetNCellsInTRUPhi();
342 nCellsPhi = nCellsPhi / 2 ; //Half size SM. Not Final.
343 // 12(tow)*2(cell)/1 TRU, cells in Phi in one TRU
344 Int_t nCellsEta = fGeom->GetNCellsInTRUEta();
345 Int_t nTRU = fGeom->GetNTRU();
346 // 24(mod)*2(tower)/3 TRU, cells in Eta in one TRU
347 //Int_t nTRU = geom->GeNTRU();//3 TRU per super module
351 for(Int_t i = 0; i < 4; i++){
352 for(Int_t j = 0; j < nTRU; j++){
358 //Create matrix that will contain 2x2 amplitude sums
359 //used to calculate the nxn sums
360 TMatrixD tru2x2(nCellsPhi/2,nCellsEta/2) ;
361 for(Int_t i = 0; i < nCellsPhi/2; i++)
362 for(Int_t j = 0; j < nCellsEta/2; j++)
365 //Loop over all TRUS in a supermodule
366 for(Int_t itru = 0 + isupermod * nTRU ; itru < (isupermod+1)*nTRU ; itru++) {
367 TMatrixD * amptru = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
368 TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
369 Int_t mtru = itru-isupermod*nTRU ; //Number of TRU in Supermodule
371 //Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
372 for(Int_t irow = 0 ; irow < nCellsPhi; irow += 2){
373 for(Int_t icol = 0 ; icol < nCellsEta ; icol += 2){
374 amp2 = (*amptru)(irow,icol)+(*amptru)(irow+1,icol)+
375 (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
377 //Fill matrix with added 2x2 cells for use in nxn sums
378 tru2x2(irow/2,icol/2) = amp2 ;
379 //Select 2x2 maximum sums to select L0
380 if(amp2 > ampmax2(0,mtru)){
381 ampmax2(0,mtru) = amp2 ;
382 ampmax2(1,mtru) = irow;
383 ampmax2(2,mtru) = icol;
388 //Find most recent time in the selected 2x2 cell
389 ampmax2(3,mtru) = 1 ;
390 Int_t row2 = static_cast <Int_t> (ampmax2(1,mtru));
391 Int_t col2 = static_cast <Int_t> (ampmax2(2,mtru));
392 for(Int_t i = 0; i<2; i++){
393 for(Int_t j = 0; j<2; j++){
394 if((*amptru)(row2+i,col2+j) > 0 && (*timeRtru)(row2+i,col2+j)> 0){
395 if((*timeRtru)(row2+i,col2+j) < ampmax2(3,mtru) )
396 ampmax2(3,mtru) = (*timeRtru)(row2+i,col2+j);
401 //Sliding nxn, add nxn amplitudes (OVERLAP)
403 for(Int_t irow = 0 ; irow < nCellsPhi/2; irow++){
404 for(Int_t icol = 0 ; icol < nCellsEta/2 ; icol++){
406 if( (irow+fPatchSize) < nCellsPhi/2 && (icol+fPatchSize) < nCellsEta/2){//Avoid exit the TRU
407 for(Int_t i = 0 ; i <= fPatchSize ; i++)
408 for(Int_t j = 0 ; j <= fPatchSize ; j++)
409 ampn += tru2x2(irow+i,icol+j);
410 //Select nxn maximum sums to select L1
411 if(ampn > ampmaxn(0,mtru)){
412 ampmaxn(0,mtru) = ampn ;
413 ampmaxn(1,mtru) = irow*2;
414 ampmaxn(2,mtru) = icol*2;
420 //Find most recent time in selected nxn cell
421 ampmaxn(3,mtru) = 1 ;
422 Int_t rown = static_cast <Int_t> (ampmaxn(1,mtru));
423 Int_t coln = static_cast <Int_t> (ampmaxn(2,mtru));
424 for(Int_t i = 0; i<4*fPatchSize; i++){
425 for(Int_t j = 0; j<4*fPatchSize; j++){
426 if( (rown+i) < nCellsPhi && (coln+j) < nCellsEta){//Avoid exit the TRU
427 if((*amptru)(rown+i,coln+j) > 0 && (*timeRtru)(rown+i,coln+j)> 0){
428 if((*timeRtru)(rown+i,coln+j) < ampmaxn(3,mtru) )
429 ampmaxn(3,mtru) = (*timeRtru)(rown+i,coln+j);
436 ampmaxn(0,mtru) = ampmax2(0,mtru);
437 ampmaxn(1,mtru) = ampmax2(1,mtru);
438 ampmaxn(2,mtru) = ampmax2(2,mtru);
439 ampmaxn(3,mtru) = ampmax2(3,mtru);
444 //____________________________________________________________________________
445 void AliEMCALTrigger::MakeSlidingTowers(const TClonesArray * amptrus, const TClonesArray * timeRtrus,
446 const Int_t isupermod,TMatrixD &max2, TMatrixD &maxn){
448 // Output from module (2x2 cells from one module)
449 Int_t nModulesPhi = fGeom->GetNModulesInTRUPhi(); // now 4 modules (3 div in phi)
451 nModulesPhi = nModulesPhi / 2 ; // Half size SM. Not Final.
453 Int_t nModulesEta = fGeom->GetNModulesInTRUEta(); // now 24 modules (no division in eta)
454 Int_t nTRU = fGeom->GetNTRU();
455 static int keyPrint = 0;
456 if(keyPrint) AliDebug(2,Form("MakeSlidingTowers : nTRU %i nModulesPhi %i nModulesEta %i ",
457 nTRU, nModulesPhi, nModulesEta ));
461 for(Int_t i = 0; i < 4; i++){
462 for(Int_t j = 0; j < nTRU; j++){
463 ampmax2(i,j) = ampmaxn(i,j) = -1;
467 // Create matrix that will contain 2x2 amplitude sums
468 // used to calculate the nxn sums
469 TMatrixD tru2x2(nModulesPhi/2,nModulesEta/2);
471 // Loop over all TRUS in a supermodule
472 for(Int_t itru = 0 + isupermod * nTRU ; itru < (isupermod+1)*nTRU ; itru++) {
473 TMatrixD * amptru = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
474 TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
475 Int_t mtru = itru - isupermod*nTRU ; // Number of TRU in Supermodule !!
477 if(!amptru || !timeRtru){
478 AliError("Amplitude or Time TRU matrix not available")
482 // Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
483 for(Int_t irow = 0 ; irow < nModulesPhi; irow +=2){
484 for(Int_t icol = 0 ; icol < nModulesEta ; icol +=2){
485 amp2 = (*amptru)(irow,icol) +(*amptru)(irow+1,icol)+
486 (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
488 //Fill matrix with added 2x2 towers for use in nxn sums
489 tru2x2(irow/2,icol/2) = amp2 ;
490 //Select 2x2 maximum sums to select L0
491 if(amp2 > ampmax2(0,mtru)){
492 ampmax2(0,mtru) = amp2 ;
493 ampmax2(1,mtru) = irow;
494 ampmax2(2,mtru) = icol;
499 ampmax2(3,mtru) = 0.;
501 // Find most recent time in the selected 2x2 towers
502 Int_t row2 = static_cast <Int_t> (ampmax2(1,mtru));
503 Int_t col2 = static_cast <Int_t> (ampmax2(2,mtru));
504 for(Int_t i = 0; i<2; i++){
505 for(Int_t j = 0; j<2; j++){
506 if((*amptru)(row2+i,col2+j) > 0 && (*timeRtru)(row2+i,col2+j)> 0){
507 if((*timeRtru)(row2+i,col2+j) > ampmax2(3,mtru) )
508 ampmax2(3,mtru) = (*timeRtru)(row2+i,col2+j); // max time
514 //Sliding nxn, add nxn amplitudes (OVERLAP)
516 for(Int_t irow = 0 ; irow < nModulesPhi/2; irow++){
517 for(Int_t icol = 0 ; icol < nModulesEta/2; icol++){
519 if( (irow+fPatchSize) < nModulesPhi/2 && (icol+fPatchSize) < nModulesEta/2){ //Avoid exit the TRU
520 for(Int_t i = 0 ; i <= fPatchSize ; i++)
521 for(Int_t j = 0 ; j <= fPatchSize ; j++)
522 ampn += tru2x2(irow+i,icol+j);
523 //Select nxn maximum sums to select L1
524 if(ampn > ampmaxn(0,mtru)){
525 ampmaxn(0,mtru) = ampn ;
526 ampmaxn(1,mtru) = irow;
527 ampmaxn(2,mtru) = icol;
533 ampmaxn(3,mtru) = 0.; // Was 1 , I don't know why
535 //Find most recent time in selected nxn cell
536 Int_t rown = static_cast <Int_t> (ampmaxn(1,mtru));
537 Int_t coln = static_cast <Int_t> (ampmaxn(2,mtru));
538 for(Int_t i = 0; i<4*fPatchSize; i++){
539 for(Int_t j = 0; j<4*fPatchSize; j++){
540 if( (rown+i) < nModulesPhi && (coln+j) < nModulesEta){//Avoid exit the TRU
541 if((*amptru)(rown+i,coln+j) > 0 && (*timeRtru)(rown+i,coln+j)> 0){
542 if((*timeRtru)(rown+i,coln+j) > ampmaxn(3,mtru) )
543 ampmaxn(3,mtru) = (*timeRtru)(rown+i,coln+j); // max time
549 } else { // copy 2x2 to nxn
550 ampmaxn(0,mtru) = ampmax2(0,mtru);
551 ampmaxn(1,mtru) = ampmax2(1,mtru);
552 ampmaxn(2,mtru) = ampmax2(2,mtru);
553 ampmaxn(3,mtru) = ampmax2(3,mtru);
556 if(keyPrint) AliDebug(2,Form(" : MakeSlidingTowers -OUt \n"));
559 //____________________________________________________________________________
560 void AliEMCALTrigger::Print(const Option_t * opt) const
563 //Prints main parameters
567 AliTriggerInput* in = 0x0 ;
568 AliInfo(Form(" fSimulation %i (input option) : #digits %i\n", fSimulation, fDigitsList->GetEntries()));
569 AliInfo(Form(" fTimeKey %i \n ", fTimeKey));
571 AliInfo(Form("\t Maximum Amplitude after Sliding Cell, \n")) ;
572 AliInfo(Form("\t -2x2 cells sum (not overlapped): %10.2f, in Super Module %d\n",
573 f2x2MaxAmp,f2x2SM)) ;
574 AliInfo(Form("\t -2x2 from row %d to row %d and from column %d to column %d\n", f2x2ModulePhi, f2x2ModulePhi+2, f2x2ModuleEta, f2x2ModuleEta+2));
575 AliInfo(Form("\t -2x2 Isolation Patch %d x %d, Amplitude out of 2x2 patch is %f, threshold %f, Isolated? %d \n", 2*fIsolPatchSize+2, 2*fIsolPatchSize+2, f2x2AmpOutOfPatch, f2x2AmpOutOfPatchThres,static_cast<Int_t> (fIs2x2Isol)));
577 AliInfo(Form("\t Patch Size, n x n: %d x %d cells\n",2*(fPatchSize+1), 2*(fPatchSize+1)));
578 AliInfo(Form("\t -nxn cells sum (overlapped) : %10.2f, in Super Module %d\n", fnxnMaxAmp,fnxnSM));
579 AliInfo(Form("\t -nxn from row %d to row %d and from column %d to column %d\n", fnxnModulePhi, fnxnModulePhi+4*fPatchSize, fnxnModuleEta, fnxnModuleEta+4*fPatchSize)) ;
580 AliInfo(Form("\t -nxn Isolation Patch %d x %d, Amplitude out of nxn patch is %f, threshold %f, Isolated? %d \n", 4*fIsolPatchSize+2*(fPatchSize+1),4*fIsolPatchSize+2*(fPatchSize+1) , fnxnAmpOutOfPatch, fnxnAmpOutOfPatchThres,static_cast<Int_t> (fIsnxnIsol) ));
583 AliInfo(Form("\t Isolate in SuperModule? %d\n", fIsolateInSuperModule)) ;
584 AliInfo(Form("\t Threshold for LO %10.2f\n", fL0Threshold));
586 in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_L0" );
588 AliInfo(Form("\t *** EMCAL LO is set ***\n"));
590 AliInfo(Form("\t Gamma Low Pt Threshold for L1 %10.2f\n", fL1GammaLowPtThreshold));
591 in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_GammaLPt_L1" );
593 AliInfo(Form("\t *** EMCAL Gamma Low Pt for L1 is set ***\n"));
595 AliInfo(Form("\t Gamma Medium Pt Threshold for L1 %10.2f\n", fL1GammaMediumPtThreshold));
596 in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_GammaMPt_L1" );
598 AliInfo(Form("\t *** EMCAL Gamma Medium Pt for L1 is set ***\n"));
600 AliInfo(Form("\t Gamma High Pt Threshold for L1 %10.2f\n", fL1GammaHighPtThreshold));
601 in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_GammaHPt_L1" );
603 AliInfo(Form("\t *** EMCAL Gamma High Pt for L1 is set ***\n")) ;
607 //____________________________________________________________________________
608 void AliEMCALTrigger::SetTriggers(const TClonesArray * ampmatrix,const Int_t iSM,
609 const TMatrixD &max2,
610 const TMatrixD &maxn)
612 //Checks the 2x2 and nxn maximum amplitude per each TRU and
613 //compares with the different L0 and L1 triggers thresholds
614 Float_t max2[] = {-1,-1,-1,-1} ;
615 Float_t maxn[] = {-1,-1,-1,-1} ;
619 Int_t nTRU = fGeom->GetNTRU();
621 //Find maximum summed amplitude of all the TRU
623 for(Int_t i = 0 ; i < nTRU ; i++){
624 if(max2[0] < ampmax2(0,i) ){
625 max2[0] = ampmax2(0,i) ; // 2x2 summed max amplitude
626 max2[1] = ampmax2(1,i) ; // corresponding phi position in TRU
627 max2[2] = ampmax2(2,i) ; // corresponding eta position in TRU
628 max2[3] = ampmax2(3,i) ; // corresponding most recent time
631 if(maxn[0] < ampmaxn(0,i) ){
632 maxn[0] = ampmaxn(0,i) ; // nxn summed max amplitude
633 maxn[1] = ampmaxn(1,i) ; // corresponding phi position in TRU
634 maxn[2] = ampmaxn(2,i) ; // corresponding eta position in TRU
635 maxn[3] = ampmaxn(3,i) ; // corresponding most recent time
640 //--------Set max amplitude if larger than in other Super Modules------------
641 Float_t maxtimeR2 = -1 ;
642 Float_t maxtimeRn = -1 ;
643 static AliEMCALRawUtils rawUtil;
644 // Int_t nTimeBins = rawUtil.GetRawFormatTimeBins() ;
645 Int_t nTimeBins = TIMEBINS; //changed by PTH
647 //Set max of 2x2 amplitudes and select L0 trigger
648 if(max2[0] > f2x2MaxAmp ){
649 // if(max2[0] > 5) printf(" L0 : iSM %i: max2[0] %5.0f : max2[3] %5.0f (maxtimeR2) \n",
650 // iSM, max2[0], max2[3]);
651 f2x2MaxAmp = max2[0] ;
653 maxtimeR2 = max2[3] ;
654 fGeom->GetModulePhiEtaIndexInSModuleFromTRUIndex(mtru2,
655 static_cast<Int_t>(max2[1]),
656 static_cast<Int_t>(max2[2]),
657 f2x2ModulePhi,f2x2ModuleEta);
660 if(fIsolateInSuperModule)
661 fIs2x2Isol = IsPatchIsolated(0, ampmatrix, iSM, mtru2, f2x2MaxAmp, f2x2ModulePhi,f2x2ModuleEta) ;
663 fIs2x2Isol = IsPatchIsolated(0, ampmatrix, iSM, mtru2, f2x2MaxAmp, static_cast<Int_t>(max2[1]), static_cast<Int_t>(max2[2])) ;
666 //Transform digit amplitude in Raw Samples
667 if (fADCValuesLow2x2 == 0) {
668 fADCValuesLow2x2 = new Int_t[nTimeBins];
669 fADCValuesHigh2x2 = new Int_t[nTimeBins];
671 //printf(" maxtimeR2 %12.5e (1)\n", maxtimeR2);
672 // rawUtil.RawSampledResponse(maxtimeR2 * AliEMCALRawUtils::GetRawFormatTimeBin(),
673 // f2x2MaxAmp, fADCValuesHigh2x2, fADCValuesLow2x2) ;
675 // rawUtil.RawSampledResponse(maxtimeR2*TIMEBINMAX/TIMEBINS,
676 // f2x2MaxAmp, fADCValuesHigh2x2, fADCValuesLow2x2) ;
678 AliEMCALRawResponse::RawSampledResponse( maxtimeR2*TIMEBINMAX/TIMEBINS,
679 f2x2MaxAmp, fADCValuesHigh2x2, fADCValuesLow2x2) ;
681 // Set Trigger Inputs, compare ADC time bins until threshold is attained
683 for(Int_t i = 0 ; i < nTimeBins ; i++){
684 // printf(" fADCValuesHigh2x2[%i] %i : %i \n", i, fADCValuesHigh2x2[i], fADCValuesLow2x2[i]);
685 if(fADCValuesHigh2x2[i] >= fL0Threshold || fADCValuesLow2x2[i] >= fL0Threshold){
686 SetInput("EMCAL_L0") ;
691 // Nov 5 - no analysis of time information
692 if(f2x2MaxAmp >= fL0Threshold) { // should add the low amp too
693 SetInput("EMCAL_L0");
698 //------------Set max of nxn amplitudes and select L1 trigger---------
699 if(maxn[0] > fnxnMaxAmp ){
700 fnxnMaxAmp = maxn[0] ;
702 maxtimeRn = maxn[3] ;
703 fGeom->GetModulePhiEtaIndexInSModuleFromTRUIndex(mtrun,
704 static_cast<Int_t>(maxn[1]),
705 static_cast<Int_t>(maxn[2]),
706 fnxnModulePhi,fnxnModuleEta) ;
709 if(fIsolateInSuperModule)
710 fIsnxnIsol = IsPatchIsolated(1, ampmatrix, iSM, mtrun, fnxnMaxAmp, fnxnModulePhi, fnxnModuleEta) ;
712 fIsnxnIsol = IsPatchIsolated(1, ampmatrix, iSM, mtrun, fnxnMaxAmp, static_cast<Int_t>(maxn[1]), static_cast<Int_t>(maxn[2])) ;
715 //Transform digit amplitude in Raw Samples
716 if (fADCValuesLownxn == 0) {
717 fADCValuesHighnxn = new Int_t[nTimeBins];
718 fADCValuesLownxn = new Int_t[nTimeBins];
720 // rawUtil.RawSampledResponse(maxtimeRn * AliEMCALRawUtils::GetRawFormatTimeBin(),
721 // fnxnMaxAmp, fADCValuesHighnxn, fADCValuesLownxn) ;
723 //rawUtil.RawSampledResponse(maxtimeRn*TIMEBINMAX/TIMEBINS,
724 // fnxnMaxAmp, fADCValuesHighnxn, fADCValuesLownxn) ;
726 AliEMCALRawResponse::RawSampledResponse (maxtimeRn*TIMEBINMAX/TIMEBINS,
727 fnxnMaxAmp, fADCValuesHighnxn, fADCValuesLownxn) ;
729 //Set Trigger Inputs, compare ADC time bins until threshold is attained
731 for(Int_t i = 0 ; i < nTimeBins ; i++){
732 if(fADCValuesHighnxn[i] >= fL1GammaLowPtThreshold || fADCValuesLownxn[i] >= fL1GammaLowPtThreshold){
733 SetInput("EMCAL_GammaLPt_L1") ;
739 for(Int_t i = 0 ; i < nTimeBins ; i++){
740 if(fADCValuesHighnxn[i] >= fL1GammaMediumPtThreshold || fADCValuesLownxn[i] >= fL1GammaMediumPtThreshold){
741 SetInput("EMCAL_GammaMPt_L1") ;
747 for(Int_t i = 0 ; i < nTimeBins ; i++){
748 if(fADCValuesHighnxn[i] >= fL1GammaHighPtThreshold || fADCValuesLownxn[i] >= fL1GammaHighPtThreshold){
749 SetInput("EMCAL_GammaHPt_L1") ;
754 // Nov 5 - no analysis of time information
755 if(fnxnMaxAmp >= fL1GammaLowPtThreshold) { // should add the low amp too
756 SetInput("EMCAL_GammaLPt_L1") ; //SetL1 Low
758 if(fnxnMaxAmp >= fL1GammaMediumPtThreshold) { // should add the low amp too
759 SetInput("EMCAL_GammaMPt_L1") ; //SetL1 Medium
761 if(fnxnMaxAmp >= fL1GammaHighPtThreshold) { // should add the low amp too
762 SetInput("EMCAL_GammaHPt_L1") ; //SetL1 High
768 //____________________________________________________________________________
769 void AliEMCALTrigger::FillTRU(const TClonesArray * digits, TClonesArray * ampmatrix, TClonesArray * ampmatrixsmod, TClonesArray * timeRmatrix) {
771 // Orders digits ampitudes list in fNTRU TRUs (384 cells) per supermodule.
772 // Each TRU is a TMatrixD, and they are kept in TClonesArrays. The number of
773 // TRU in phi is fNTRUPhi, and the number of TRU in eta is fNTRUEta.
774 // Last 2 modules are half size in Phi, I considered that the number of TRU
775 // is maintained for the last modules but decision not taken. If different,
776 // then this must be changed. Also fill a matrix with all amplitudes in supermodule for isolation studies.
778 // Initilize and declare variables
779 // List of TRU matrices initialized to 0.
780 // printf("<I> AliEMCALTrigger::FillTRU() started : # digits %i\n", digits->GetEntriesFast());
783 // One input per EMCAL module so size of matrix is reduced by 4 (2x2 division case)
785 Int_t nPhi = fGeom->GetNPhi();
786 Int_t nZ = fGeom->GetNZ();
787 Int_t nTRU = fGeom->GetNTRU();
788 // Int_t nTRUPhi = fGeom->GetNTRUPhi();
789 Int_t nModulesPhi = fGeom->GetNModulesInTRUPhi();
790 Int_t nModulesPhi2 = fGeom->GetNModulesInTRUPhi();
791 Int_t nModulesEta = fGeom->GetNModulesInTRUEta();
792 // printf("<I> AliEMCALTrigger::FillTRU() nTRU %i nTRUPhi %i : nModulesPhi %i nModulesEta %i \n",
793 // nTRU, nTRUPhi, nModulesPhi, nModulesEta);
804 // iphim, ietam - module indexes in SM
808 //List of TRU matrices initialized to 0.
809 Int_t nSup = fGeom->GetNumberOfSuperModules();
810 for(Int_t k = 0; k < nTRU*nSup; k++){
811 TMatrixD amptrus(nModulesPhi,nModulesEta) ;
812 TMatrixD timeRtrus(nModulesPhi,nModulesEta) ;
813 // Do we need to initialise? I think TMatrixD does it by itself...
814 for(Int_t i = 0; i < nModulesPhi; i++){
815 for(Int_t j = 0; j < nModulesEta; j++){
817 timeRtrus(i,j) = 0.0;
820 new((*ampmatrix)[k]) TMatrixD(amptrus) ;
821 new((*timeRmatrix)[k]) TMatrixD(timeRtrus) ;
824 // List of Modules matrices initialized to 0.
825 for(Int_t k = 0; k < nSup ; k++){
827 // if(nSup>9) mphi = nPhi/2; // the same size
828 TMatrixD ampsmods( mphi, nZ);
829 for(Int_t i = 0; i < mphi; i++){
830 for(Int_t j = 0; j < nZ; j++){
834 new((*ampmatrixsmod)[k]) TMatrixD(ampsmods) ;
837 AliEMCALDigit * dig ;
839 //Digits loop to fill TRU matrices with amplitudes.
840 for(Int_t idig = 0 ; idig < digits->GetEntriesFast() ; idig++){
842 dig = dynamic_cast<AliEMCALDigit *>(digits->At(idig)) ;
844 amp = Float_t(dig->GetAmplitude()); // Energy of the digit (arbitrary units)
845 id = dig->GetId() ; // Id label of the cell
846 timeR = dig->GetTimeR() ; // Earliest time of the digit
847 if(amp<=0.0) AliDebug(1,Form(" id %i amp %f \n", id, amp));
848 // printf(" FILLTRU : timeR %10.5e time %10.5e : amp %10.5e \n", timeR, dig->GetTime(), amp);
849 // Get eta and phi cell position in supermodule
850 Bool_t bCell = fGeom->GetCellIndex(id, iSupMod, nModule, nIphi, nIeta) ;
852 AliError(Form("%i Wrong cell id number %i ", idig, id)) ;
854 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
855 // iphim, ietam - module indexes in SM
856 fGeom->GetModuleIndexesFromCellIndexesInSModule(iSupMod,iphi,ieta, iphim, ietam, nModule);
858 //printf("iSupMod %i nModule %i iphi %i ieta %i iphim %i ietam %i \n",
859 //iSupMod,nModule, iphi, ieta, iphim, ietam);
861 // Check to which TRU in the supermodule belongs the cell.
862 // Supermodules are divided in a TRU matrix of dimension
863 // (fNTRUPhi,fNTRUEta).
864 // Each TRU is a cell matrix of dimension (nModulesPhi,nModulesEta)
866 // First calculate the row and column in the supermodule
867 // of the TRU to which the cell belongs.
868 Int_t row = iphim / nModulesPhi;
869 Int_t col = ietam / nModulesEta;
870 //Calculate label number of the TRU
871 Int_t itru = fGeom->GetAbsTRUNumberFromNumberInSm(row, col, iSupMod);
873 //Fill TRU matrix with cell values
874 TMatrixD * amptrus = dynamic_cast<TMatrixD *>(ampmatrix->At(itru)) ;
875 TMatrixD * timeRtrus = dynamic_cast<TMatrixD *>(timeRmatrix->At(itru)) ;
877 if(!amptrus || !timeRtrus){
878 AliError("Could not recover the TRU matrix with amplitudes or times");
881 //Calculate row and column of the module inside the TRU with number itru
882 Int_t irow = iphim - row * nModulesPhi;
884 irow = iphim - row * nModulesPhi2; // size of matrix the same
885 Int_t icol = ietam - col * nModulesEta;
887 (*amptrus)(irow,icol) += amp ;
888 if((*timeRtrus)(irow,icol) <0.0 || (*timeRtrus)(irow,icol) <= timeR){ // ??
889 (*timeRtrus)(irow,icol) = timeR ;
892 //printf(" ieta %i iphi %i iSM %i || col %i row %i : itru %i -> amp %f\n",
893 // ieta, iphi, iSupMod, col, row, itru, amp);
894 //####################SUPERMODULE MATRIX ##################
895 TMatrixD * ampsmods = dynamic_cast<TMatrixD *>(ampmatrixsmod->At(iSupMod)) ;
897 AliError("Could not recover the matrix per SM");
900 (*ampsmods)(iphim,ietam) += amp ;
901 // printf(" id %i iphim %i ietam %i SM %i : irow %i icol %i itru %i : amp %6.0f\n",
902 //id, iphim, ietam, iSupMod, irow, icol, itru, amp);
904 else AliError("Could not recover the digit");
907 //printf("<I> AliEMCALTrigger::FillTRU() is ended \n");
910 //____________________________________________________________________________
911 void AliEMCALTrigger::Trigger()
913 //Main Method to select triggers.
914 TH1::AddDirectory(0);
916 AliRunLoader *runLoader = AliRunLoader::Instance();
917 AliEMCALLoader *emcalLoader = 0;
919 emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
922 //Load EMCAL Geometry
923 if (runLoader && runLoader->GetAliRun()){
924 AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(runLoader->GetAliRun()->GetDetector("EMCAL"));
925 if(emcal)fGeom = emcal->GetGeometry();
929 fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
932 AliFatal("Did not get geometry from EMCALLoader");
935 Int_t nSuperModules = fGeom->GetNumberOfSuperModules() ; //12 SM in EMCAL
936 Int_t nTRU = fGeom->GetNTRU(); // 3 TRU per super module
938 //Intialize data members each time the trigger is called in event loop
939 f2x2MaxAmp = -1; f2x2ModulePhi = -1; f2x2ModuleEta = -1;
940 fnxnMaxAmp = -1; fnxnModulePhi = -1; fnxnModuleEta = -1;
942 // Take the digits list if simulation
943 if(fSimulation && runLoader && emcalLoader){ // works than run seperate macros
944 runLoader->LoadDigits("EMCAL");
945 fDigitsList = emcalLoader->Digits() ;
946 runLoader->LoadSDigits("EMCAL");
948 // Digits list should be set by method SetDigitsList(TClonesArray * digits)
950 AliFatal("Digits not found !") ;
952 //Take the digits list
954 // Delete old if unzero
955 if(fAmpTrus) {fAmpTrus->Delete(); delete fAmpTrus;}
956 if(fTimeRtrus) {fTimeRtrus->Delete(); delete fTimeRtrus;}
957 if(fAmpSMods) {fAmpSMods->Delete(); delete fAmpSMods;}
958 // Fill TRU and SM matrix
959 fAmpTrus = new TClonesArray("TMatrixD",nTRU);
960 fAmpTrus->SetName("AmpTrus");
961 fTimeRtrus = new TClonesArray("TMatrixD",nTRU);
962 fTimeRtrus->SetName("TimeRtrus");
963 fAmpSMods = new TClonesArray("TMatrixD",nSuperModules);
964 fAmpSMods->SetName("AmpSMods");
966 FillTRU(fDigitsList, fAmpTrus, fAmpSMods, fTimeRtrus);
968 // Jet stuff - only one case, no freedom here
969 if(fGeom->GetNEtaSubOfTRU() == 6) {
970 if(fAmpJetMatrix) {delete fAmpJetMatrix; fAmpJetMatrix=0;}
971 if(fJetMatrixE) {delete fJetMatrixE; fJetMatrixE=0;}
973 fAmpJetMatrix = new TMatrixD(17,12); // 17-phi(row), 12-eta(col)
974 fJetMatrixE = new TH2F("fJetMatrixE"," E of max patch in (#phi,#eta)",
975 17, 80.*TMath::DegToRad(), (180.+20.*2/3.)*TMath::DegToRad(), 12, -0.7, 0.7);
976 for(Int_t row=0; row<fAmpJetMatrix->GetNrows(); row++) {
977 for(Int_t col=0; col<fAmpJetMatrix->GetNcols(); col++) {
978 (*fAmpJetMatrix)(row,col) = 0.;
981 FillJetMatrixFromSMs(fAmpSMods, fAmpJetMatrix, fGeom);
983 if(!CheckConsistentOfMatrixes()) assert(0);
985 // Do Tower Sliding and select Trigger
986 // Initialize varible that will contain maximum amplitudes and
987 // its corresponding tower position in eta and phi, and time.
988 TMatrixD ampmax2(4,nTRU) ; // 0-max amp, 1-irow, 2-icol, 3-timeR
989 TMatrixD ampmaxn(4,nTRU) ;
991 for(Int_t iSM = 0 ; iSM < nSuperModules ; iSM++) {
992 //Do 2x2 and nxn sums, select maximums.
994 MakeSlidingTowers(fAmpTrus, fTimeRtrus, iSM, ampmax2, ampmaxn);
997 if(fIsolateInSuperModule) // here some discripency between tru and SM
998 SetTriggers(fAmpSMods,iSM,ampmax2,ampmaxn) ;
999 if(!fIsolateInSuperModule)
1000 SetTriggers(fAmpTrus,iSM,ampmax2,ampmaxn) ;
1003 // Do patch sliding and select Jet Trigger
1004 // 0-max amp-meanFromVZERO(if), 1-irow, 2-icol, 3-timeR,
1005 // 4-max amp , 5-meanFromVZERO (Nov 25, 2007)
1007 MakeSlidingPatch((*fAmpJetMatrix), fNJetPatchPhi, fAmpJetMax); // no timing information here
1013 //____________________________________________________________________________
1014 void AliEMCALTrigger::GetTriggerInfo(TArrayF &triggerPosition, TArrayF &triggerAmplitudes) const
1016 // Template - should be defined; Nov 5, 2007
1017 triggerPosition[0] = 0.;
1018 triggerAmplitudes[0] = 0.;
1021 //____________________________________________________________________________
1022 void AliEMCALTrigger::FillJetMatrixFromSMs(TClonesArray *ampmatrixsmod, TMatrixD* jetMat, AliEMCALGeometry *g)
1025 // Fill matrix for jet trigger from SM matrixes of modules
1027 static int keyPrint = 0;
1029 if(ampmatrixsmod==0 || jetMat==0 || g==0) return;
1030 Double_t amp = 0.0, ampSum=0.0;
1032 Int_t nEtaModSum = g->GetNZ() / g->GetNEtaSubOfTRU(); // should be 4
1033 Int_t nPhiModSum = g->GetNPhi() / g->GetNTRUPhi(); // should be 4
1035 if(keyPrint) AliDebug(2,Form("%s",Form(" AliEMCALTrigger::FillJetMatrixFromSMs | nEtaModSum %i : nPhiModSum %i \n", nEtaModSum, nPhiModSum)));
1036 Int_t jrow=0, jcol=0; // indexes of jet matrix
1037 Int_t nEtaSM=0, nPhiSM=0;
1038 for(Int_t iSM=0; iSM<ampmatrixsmod->GetEntries(); iSM++) {
1039 TMatrixD * ampsmods = dynamic_cast<TMatrixD *>(ampmatrixsmod->At(iSM));
1041 if(!ampsmods) return;
1043 Int_t nrow = ampsmods->GetNrows();
1044 Int_t ncol = ampsmods->GetNcols();
1045 //printf("%s",Form(" ######## SM %i : nrow %i : ncol %i ##### \n", iSM, nrow, ncol));
1046 for(Int_t row=0; row<nrow; row++) {
1047 for(Int_t col=0; col<ncol; col++) {
1048 amp = (*ampsmods)(row,col);
1052 if(keyPrint) AliDebug(2,Form("%s",Form(" ** nPhiSm %i : nEtaSM %i : row %2.2i : col %2.2i -> ", nPhiSM, nEtaSM, row, col)));
1053 if(nEtaSM == 0) { // positive Z
1054 jrow = 3*nPhiSM + row/nPhiModSum;
1055 jcol = 6 + col / nEtaModSum;
1056 } else { // negative Z
1057 if(iSM<=9) jrow = 3*nPhiSM + 2 - row/nPhiModSum;
1058 else jrow = 3*nPhiSM + 1 - row/nPhiModSum; // half size
1059 jcol = 5 - col / nEtaModSum;
1061 if(keyPrint) AliDebug(2,Form("%s",Form(" jrow %2.2i : jcol %2.2i : amp %f (jetMat) \n", jrow, jcol, amp)));
1063 (*jetMat)(jrow,jcol) += amp;
1064 ampSum += amp; // For controling
1065 } else if(amp<0.0) {
1066 AliDebug(1,Form(" jrow %2.2i : jcol %2.2i : amp %f (jetMat: amp<0) \n", jrow, jcol, amp));
1072 if(ampSum <= 0.0) AliDebug(1,Form("ampSum %f (<=0.0) ", ampSum));
1075 //____________________________________________________________________________
1076 void AliEMCALTrigger::MakeSlidingPatch(const TMatrixD &jm, const Int_t nPatchSize, TMatrixD &JetMax)
1078 // Sliding patch : nPatchSize x nPatchSize (OVERLAP)
1079 static int keyPrint = 0;
1080 if(keyPrint) AliDebug(2,Form(" AliEMCALTrigger::MakeSlidingPatch() was started \n"));
1081 Double_t ampCur = 0.0, e=0.0;
1082 ampJetMax(0,0) = 0.0;
1083 ampJetMax(3,0) = 0.0; // unused now
1084 ampJetMax(4,0) = ampJetMax(5,0) = 0.0;
1085 for(Int_t row=0; row<fAmpJetMatrix->GetNrows(); row ++) {
1086 for(Int_t col=0; col<fAmpJetMatrix->GetNcols(); col++) {
1088 // check on patch size
1089 if( (row+nPatchSize-1) < fAmpJetMatrix->GetNrows() && (col+nPatchSize-1) < fAmpJetMatrix->GetNcols()){
1090 for(Int_t i = 0 ; i < nPatchSize ; i++) {
1091 for(Int_t j = 0 ; j < nPatchSize ; j++) {
1092 ampCur += jm(row+i, col+j);
1094 } // end cycle on patch
1095 if(ampCur > ampJetMax(0,0)){
1096 ampJetMax(0,0) = ampCur;
1097 ampJetMax(1,0) = row;
1098 ampJetMax(2,0) = col;
1100 } // check on patch size
1103 if(keyPrint) AliDebug(2,Form(" ampJetMax %i row %2i->%2i col %2i->%2i \n", Int_t(ampJetMax(0,0)), Int_t(ampJetMax(1,0)), Int_t(ampJetMax(1,0))+nPatchSize-1, Int_t(ampJetMax(2,0)), Int_t(ampJetMax(2,0))+nPatchSize-1));
1105 Double_t eCorrJetMatrix=0.0;
1106 if(fVZER0Mult > 0.0) {
1107 // Correct patch energy (adc) and jet patch matrix energy
1108 Double_t meanAmpBG = GetMeanEmcalPatchEnergy(Int_t(fVZER0Mult), nPatchSize)/0.0153;
1109 ampJetMax(4,0) = ampJetMax(0,0);
1110 ampJetMax(5,0) = meanAmpBG;
1112 Double_t eCorr = ampJetMax(0,0) - meanAmpBG;
1113 AliDebug(2,Form(" ampJetMax(0,0) %f meanAmpBG %f eCorr %f : ampJetMax(4,0) %f \n",
1114 ampJetMax(0,0), meanAmpBG, eCorr, ampJetMax(5,0)));
1115 ampJetMax(0,0) = eCorr;
1117 eCorrJetMatrix = GetMeanEmcalEnergy(Int_t(fVZER0Mult)) / 208.;
1119 // Fill patch energy matrix
1120 for(int row=Int_t(ampJetMax(1,0)); row<Int_t(ampJetMax(1,0))+nPatchSize; row++) {
1121 for(int col=Int_t(ampJetMax(2,0)); col<Int_t(ampJetMax(2,0))+nPatchSize; col++) {
1122 e = Double_t(jm(row,col)*0.0153); // 0.0153 - hard coded now
1123 if(eCorrJetMatrix > 0.0) { // BG subtraction case
1124 e -= eCorrJetMatrix;
1125 fJetMatrixE->SetBinContent(row+1, col+1, e);
1126 } else if(e > 0.0) {
1127 fJetMatrixE->SetBinContent(row+1, col+1, e);
1131 // PrintJetMatrix();
1132 // Set the jet trigger(s), multiple threshold now, Nov 19,2007
1133 for(Int_t i=0; i<fNJetThreshold; i++ ) {
1134 if(ampJetMax(0,0) >= fL1JetThreshold[i]) {
1135 SetInput(GetNameOfJetTrigger(i));
1140 //____________________________________________________________________________
1141 Double_t AliEMCALTrigger::GetEmcalSumAmp() const
1143 // Return sum of amplidutes from EMCal
1144 // Used calibration coefficeint for transition to energy
1145 return fAmpJetMatrix >0 ?fAmpJetMatrix->Sum() :0.0;
1148 //____________________________________________________________________________
1149 void AliEMCALTrigger::PrintJetMatrix() const
1151 // fAmpJetMatrix : (17,12); // 17-phi(row), 12-eta(col)
1152 if(fAmpJetMatrix == 0) return;
1154 AliInfo(Form("\n #### jetMatrix : (%i,%i) ##### \n ",
1155 fAmpJetMatrix->GetNrows(), fAmpJetMatrix->GetNcols()));
1156 PrintMatrix(*fAmpJetMatrix);
1159 //____________________________________________________________________________
1160 void AliEMCALTrigger::PrintAmpTruMatrix(Int_t ind) const
1162 // Print matrix with TRU patches
1163 TMatrixD * tru = dynamic_cast<TMatrixD *>(fAmpTrus->At(ind));
1164 if(tru == 0) return;
1165 AliInfo(Form("\n #### Amp TRU matrix(%i) : (%i,%i) ##### \n ",
1166 ind, tru->GetNrows(), tru->GetNcols()));
1170 //____________________________________________________________________________
1171 void AliEMCALTrigger::PrintAmpSmMatrix(Int_t ind) const
1173 // Print matrix with SM amplitudes
1174 TMatrixD * sm = dynamic_cast<TMatrixD *>(fAmpSMods->At(ind));
1176 AliInfo(Form("\n #### Amp SM matrix(%i) : (%i,%i) ##### \n ",
1177 ind, sm->GetNrows(), sm->GetNcols()));
1181 //____________________________________________________________________________
1182 void AliEMCALTrigger::PrintMatrix(const TMatrixD &mat) const
1184 //Print matrix object
1185 for(Int_t col=0; col<mat.GetNcols(); col++) AliInfo(Form(" %3i ", col));
1186 AliInfo(Form("\n -- \n"));
1187 for(Int_t row=0; row<mat.GetNrows(); row++) {
1188 AliInfo(Form(" row:%2i ", row));
1189 for(Int_t col=0; col<mat.GetNcols(); col++) {
1190 AliInfo(Form(" %4i", (Int_t)mat(row,col)));
1196 //____________________________________________________________________________
1197 Bool_t AliEMCALTrigger::CheckConsistentOfMatrixes(const Int_t pri)
1199 // Check consitency of matrices
1200 Double_t sumSM = 0.0, smCur=0.0;
1201 Double_t sumTru=0.0, sumTruInSM = 0.0, truSum=0.0;
1202 // Bool_t key = kTRUE;
1204 for(Int_t i=0; i<fAmpSMods->GetEntries(); i++) {
1205 TMatrixD * sm = dynamic_cast<TMatrixD *>(fAmpSMods->At(i));
1211 for(Int_t itru=0; itru<3; itru++) { // Cycle on tru inside SM
1212 Int_t ind = 3*i + itru;
1213 TMatrixD *tru = dynamic_cast<TMatrixD *>(fAmpTrus->At(ind));
1215 truSum = tru->Sum();
1216 sumTruInSM += truSum;
1219 sumTru += sumTruInSM;
1221 if(sumTruInSM != smCur) {
1222 AliDebug(1,Form(" sm %i : smCur %f -> sumTruInSM %f \n", i, smCur, sumTruInSM));
1227 Double_t sumJetMat = fAmpJetMatrix->Sum();
1228 if(pri || TMath::Abs(sumSM-sumTru)>0.0001 || TMath::Abs(sumSM-sumJetMat) > 0.0001)
1229 AliDebug(1,Form(" sumSM %f : sumTru %f : sumJetMat %f \n", sumSM, sumTru, sumJetMat));
1230 if(TMath::Abs(sumSM - sumTru)>0.0001 || TMath::Abs(sumSM-sumJetMat) > 0.0001) return kFALSE;
1234 //____________________________________________________________________________
1235 void AliEMCALTrigger::Browse(TBrowser* b)
1238 if(&fInputs) b->Add(&fInputs);
1239 if(fAmpTrus) b->Add(fAmpTrus);
1240 if(fTimeRtrus) b->Add(fTimeRtrus);
1241 if(fAmpSMods) b->Add(fAmpSMods);
1242 if(fAmpJetMatrix) b->Add(fAmpJetMatrix);
1243 if(fJetMatrixE) b->Add(fJetMatrixE);