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"
64 ClassImp(AliEMCALTrigger)
66 TString AliEMCALTrigger::fgNameOfJetTriggers("EMCALJetTriggerL1");
68 //______________________________________________________________________
69 AliEMCALTrigger::AliEMCALTrigger()
70 : AliTriggerDetector(), fGeom(0),
71 f2x2MaxAmp(-1), f2x2ModulePhi(-1), f2x2ModuleEta(-1),
73 fnxnMaxAmp(-1), fnxnModulePhi(-1), fnxnModuleEta(-1),
75 fADCValuesHighnxn(0),fADCValuesLownxn(0),
76 fADCValuesHigh2x2(0),fADCValuesLow2x2(0),
78 fL0Threshold(100),fL1GammaLowPtThreshold(200),
79 fL1GammaMediumPtThreshold(500), fL1GammaHighPtThreshold(1000),
80 fPatchSize(1), fIsolPatchSize(1),
81 f2x2AmpOutOfPatch(-1), fnxnAmpOutOfPatch(-1),
82 f2x2AmpOutOfPatchThres(100000), fnxnAmpOutOfPatchThres(100000),
83 fIs2x2Isol(kFALSE), fIsnxnIsol(kFALSE),
84 fSimulation(kTRUE), fIsolateInSuperModule(kTRUE), fTimeKey(kFALSE),
85 fAmpTrus(0),fTimeRtrus(0),fAmpSMods(0),
86 fTriggerPosition(6), fTriggerAmplitudes(4),
87 fNJetPatchPhi(3), fNJetPatchEta(3), fNJetThreshold(3), fL1JetThreshold(0), fJetMaxAmp(0),
88 fAmpJetMatrix(0), fJetMatrixE(0), fAmpJetMax(6,1), fVZER0Mult(0.)
91 fADCValuesHighnxn = 0x0; //new Int_t[fTimeBins];
92 fADCValuesLownxn = 0x0; //new Int_t[fTimeBins];
93 fADCValuesHigh2x2 = 0x0; //new Int_t[fTimeBins];
94 fADCValuesLow2x2 = 0x0; //new Int_t[fTimeBins];
97 // Define jet threshold - can not change from outside now
98 // fNJetThreshold = 7; // For MB Pythia suppression
99 fNJetThreshold = 10; // Hijing
100 fL1JetThreshold = new Double_t[fNJetThreshold];
101 if(fNJetThreshold == 7) {
102 fL1JetThreshold[0] = 5./0.0153;
103 fL1JetThreshold[1] = 8./0.0153;
104 fL1JetThreshold[2] = 10./0.0153;
105 fL1JetThreshold[3] = 12./0.0153;
106 fL1JetThreshold[4] = 13./0.0153;
107 fL1JetThreshold[5] = 14./0.0153;
108 fL1JetThreshold[6] = 15./0.0153;
109 } else if(fNJetThreshold == 10) {
110 Double_t thGev[10]={5.,8.,10., 12., 13.,14.,15., 17., 20., 25.};
111 for(Int_t i=0; i<fNJetThreshold; i++) fL1JetThreshold[i] = thGev[i]/0.0153;
113 fL1JetThreshold[0] = 5./0.0153;
114 fL1JetThreshold[1] = 10./0.0153;
115 fL1JetThreshold[2] = 15./0.0153;
116 fL1JetThreshold[3] = 20./0.0153;
117 fL1JetThreshold[4] = 25./0.0153;
122 fInputs.SetName("TriggersInputs");
128 //____________________________________________________________________________
129 AliEMCALTrigger::AliEMCALTrigger(const AliEMCALTrigger & trig)
130 : AliTriggerDetector(trig),
132 f2x2MaxAmp(trig.f2x2MaxAmp),
133 f2x2ModulePhi(trig.f2x2ModulePhi),
134 f2x2ModuleEta(trig.f2x2ModuleEta),
136 fnxnMaxAmp(trig.fnxnMaxAmp),
137 fnxnModulePhi(trig.fnxnModulePhi),
138 fnxnModuleEta(trig.fnxnModuleEta),
140 fADCValuesHighnxn(trig.fADCValuesHighnxn),
141 fADCValuesLownxn(trig.fADCValuesLownxn),
142 fADCValuesHigh2x2(trig.fADCValuesHigh2x2),
143 fADCValuesLow2x2(trig.fADCValuesLow2x2),
144 fDigitsList(trig.fDigitsList),
145 fL0Threshold(trig.fL0Threshold),
146 fL1GammaLowPtThreshold(trig.fL1GammaLowPtThreshold),
147 fL1GammaMediumPtThreshold(trig.fL1GammaMediumPtThreshold),
148 fL1GammaHighPtThreshold(trig.fL1GammaHighPtThreshold),
149 fPatchSize(trig.fPatchSize),
150 fIsolPatchSize(trig.fIsolPatchSize),
151 f2x2AmpOutOfPatch(trig.f2x2AmpOutOfPatch),
152 fnxnAmpOutOfPatch(trig.fnxnAmpOutOfPatch),
153 f2x2AmpOutOfPatchThres(trig.f2x2AmpOutOfPatchThres),
154 fnxnAmpOutOfPatchThres(trig.fnxnAmpOutOfPatchThres),
155 fIs2x2Isol(trig.fIs2x2Isol),
156 fIsnxnIsol(trig.fIsnxnIsol),
157 fSimulation(trig.fSimulation),
158 fIsolateInSuperModule(trig.fIsolateInSuperModule),
159 fTimeKey(trig.fTimeKey),
160 fAmpTrus(trig.fAmpTrus),
161 fTimeRtrus(trig.fTimeRtrus),
162 fAmpSMods(trig.fAmpSMods),
163 fTriggerPosition(trig.fTriggerPosition),
164 fTriggerAmplitudes(trig.fTriggerAmplitudes),
165 fNJetPatchPhi(trig.fNJetPatchPhi),
166 fNJetPatchEta(trig.fNJetPatchEta),
167 fNJetThreshold(trig.fNJetThreshold),
168 fL1JetThreshold(trig.fL1JetThreshold),
169 fJetMaxAmp(trig.fJetMaxAmp),
170 fAmpJetMatrix(trig.fAmpJetMatrix),
171 fJetMatrixE(trig.fJetMatrixE),
172 fAmpJetMax(trig.fAmpJetMax),
173 fVZER0Mult(trig.fVZER0Mult)
178 //____________________________________________________________________________
179 AliEMCALTrigger::~AliEMCALTrigger() {
184 delete [] fADCValuesHighnxn;
185 delete [] fADCValuesLownxn;
186 delete [] fADCValuesHigh2x2;
187 delete [] fADCValuesLow2x2;
189 if(fAmpTrus) {fAmpTrus->Delete(); delete fAmpTrus;}
190 if(fTimeRtrus) {fTimeRtrus->Delete(); delete fTimeRtrus;}
191 if(fAmpSMods) {fAmpSMods->Delete(); delete fAmpSMods;}
192 if(fAmpJetMatrix) delete fAmpJetMatrix;
193 if(fJetMatrixE) delete fJetMatrixE;
194 if(fL1JetThreshold) delete [] fL1JetThreshold;
197 //----------------------------------------------------------------------
198 void AliEMCALTrigger::CreateInputs()
202 // Do not create inputs again!!
203 if( fInputs.GetEntriesFast() > 0 ) return;
205 // Second parameter should be detector name = "EMCAL"
206 TString det("EMCAL"); // Apr 29, 2008
207 fInputs.AddLast( new AliTriggerInput( det+"_L0", det, 0x02) );
208 fInputs.AddLast( new AliTriggerInput( det+"_GammaHPt_L1", det, 0x04 ) );
209 fInputs.AddLast( new AliTriggerInput( det+"_GammaMPt_L1", det, 0x08 ) );
210 fInputs.AddLast( new AliTriggerInput( det+"_GammaLPt_L1", det, 0x016 ) );
211 fInputs.AddLast( new AliTriggerInput( det+"_JetHPt_L1", det, 0x032 ) );
212 fInputs.AddLast( new AliTriggerInput( det+"_JetMPt_L1", det, 0x048 ) );
213 fInputs.AddLast( new AliTriggerInput( det+"_JetLPt_L1", det, 0x064 ) );
215 if(fNJetThreshold<=0) return;
217 UInt_t level = 0x032;
218 for(Int_t i=0; i<fNJetThreshold; i++ ) {
219 TString name(GetNameOfJetTrigger(i));
220 TString title("EMCAL Jet triger L1 :"); // unused now
221 // 0.0153 - hard coded now
222 title += Form("Th %i(%5.1f GeV) :", (Int_t)fL1JetThreshold[i], fL1JetThreshold[i] * 0.0153);
223 title += Form("patch %ix%i~(%3.2f(#phi)x%3.2f(#eta)) ",
224 fNJetPatchPhi, fNJetPatchEta, 0.11*(fNJetPatchPhi), 0.11*(fNJetPatchEta));
225 fInputs.AddLast( new AliTriggerInput(name, det, UChar_t(level)) );
231 //____________________________________________________________________________
232 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) {
235 // EMCAL RTU size is 4modules(phi) x 24modules (eta)
236 // So maximum size of patch is 4modules x 4modules (EMCAL L0 trigger).
237 // Calculate if the maximum patch found is isolated, find amplitude around maximum (2x2 or nxn) patch,
238 // inside isolation patch . iPatchType = 0 means calculation for 2x2 patch,
239 // iPatchType = 1 means calculation for nxn patch.
240 // In the next table there is an example of the different options of patch size and isolation patch size:
241 // Patch Size (fPatchSize)
243 // fIsolPatchSize 0 2x2 (not overlap) 4x4 (overlapped)
247 if(!ampmatrixes) return kFALSE;
249 // Get matrix of TRU or Module with maximum amplitude patch.
250 Int_t itru = mtru + iSM * fGeom->GetNTRU(); //number of tru, min 0 max 3*12=36.
251 TMatrixD * ampmatrix = 0x0;
254 static int keyPrint = 0;
255 if(keyPrint) AliDebug(2,Form(" IsPatchIsolated : iSM %i mtru %i itru %i maxphi %i maxeta %i \n", iSM, mtru, itru, maxphi, maxeta));
257 if(fIsolateInSuperModule){ // ?
258 ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(iSM)) ;
259 rowborder = fGeom->GetNPhi();
260 colborder = fGeom->GetNZ();
261 AliDebug(2,"Isolate trigger in Module");
263 ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(itru)) ;
264 rowborder = fGeom->GetNModulesInTRUPhi();
265 colborder = fGeom->GetNModulesInTRUEta();
266 AliDebug(2,"Isolate trigger in TRU");
268 if(iSM>9) rowborder /= 2; // half size in phi
270 if(!ampmatrixes || !ampmatrix){
271 AliError("Could not recover the matrix with the amplitudes");
275 //Define patch modules - what is this ??
276 Int_t isolmodules = fIsolPatchSize*(1+iPatchType);
277 Int_t ipatchmodules = 2*(1+fPatchSize*iPatchType);
278 Int_t minrow = maxphi - isolmodules;
279 Int_t mincol = maxeta - isolmodules;
280 Int_t maxrow = maxphi + isolmodules + ipatchmodules;
281 Int_t maxcol = maxeta + isolmodules + ipatchmodules;
283 minrow = minrow<0?0 :minrow;
284 mincol = mincol<0?0 :mincol;
286 maxrow = maxrow>rowborder?rowborder :maxrow;
287 maxcol = maxcol>colborder?colborder :maxcol;
289 //printf("%s\n",Form("Number of added Isol Modules %d, Patch Size %d",isolmodules, ipatchmodules));
290 //printf("%s\n",Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
291 // AliDebug(2,Form("Number of added Isol Modules %d, Patch Size %d",isolmodules, ipatchmodules));
292 //AliDebug(2,Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
294 //Add amplitudes in all isolation patch
296 for(Int_t irow = minrow ; irow < maxrow; irow ++)
297 for(Int_t icol = mincol ; icol < maxcol ; icol ++)
298 amp += (*ampmatrix)(irow,icol);
300 AliDebug(2,Form("Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
303 // AliError(Form("Bad sum: Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
304 // ampmatrix->Print();
307 amp-=maxamp; //Calculate energy in isolation patch that do not comes from maximum patch.
310 AliDebug(2, Form("Maximum amplitude %f, Out of patch %f",maxamp, amp));
312 //Fill isolation amplitude data member and say if patch is isolated.
313 if(iPatchType == 0){ //2x2 case
314 f2x2AmpOutOfPatch = amp;
315 if(amp < f2x2AmpOutOfPatchThres) b=kTRUE;
316 } else if(iPatchType == 1){ //nxn case
317 fnxnAmpOutOfPatch = amp;
318 if(amp < fnxnAmpOutOfPatchThres) b=kTRUE;
321 if(keyPrint) AliDebug(2,Form(" IsPatchIsolated - OUT \n"));
328 //____________________________________________________________________________
329 void AliEMCALTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus, const Int_t isupermod,TMatrixD &max2, TMatrixD &maxn){
331 //Sums energy of all possible 2x2 (L0) and nxn (L1) modules per each TRU.
332 //Fast signal in the experiment is given by 2x2 modules,
333 //for this reason we loop inside the TRU modules by 2.
335 //Declare and initialize variables
336 Int_t nCellsPhi = fGeom->GetNCellsInTRUPhi();
338 nCellsPhi = nCellsPhi / 2 ; //Half size SM. Not Final.
339 // 12(tow)*2(cell)/1 TRU, cells in Phi in one TRU
340 Int_t nCellsEta = fGeom->GetNCellsInTRUEta();
341 Int_t nTRU = fGeom->GetNTRU();
342 // 24(mod)*2(tower)/3 TRU, cells in Eta in one TRU
343 //Int_t nTRU = geom->GeNTRU();//3 TRU per super module
347 for(Int_t i = 0; i < 4; i++){
348 for(Int_t j = 0; j < nTRU; j++){
354 //Create matrix that will contain 2x2 amplitude sums
355 //used to calculate the nxn sums
356 TMatrixD tru2x2(nCellsPhi/2,nCellsEta/2) ;
357 for(Int_t i = 0; i < nCellsPhi/2; i++)
358 for(Int_t j = 0; j < nCellsEta/2; j++)
361 //Loop over all TRUS in a supermodule
362 for(Int_t itru = 0 + isupermod * nTRU ; itru < (isupermod+1)*nTRU ; itru++) {
363 TMatrixD * amptru = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
364 TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
365 Int_t mtru = itru-isupermod*nTRU ; //Number of TRU in Supermodule
367 //Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
368 for(Int_t irow = 0 ; irow < nCellsPhi; irow += 2){
369 for(Int_t icol = 0 ; icol < nCellsEta ; icol += 2){
370 amp2 = (*amptru)(irow,icol)+(*amptru)(irow+1,icol)+
371 (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
373 //Fill matrix with added 2x2 cells for use in nxn sums
374 tru2x2(irow/2,icol/2) = amp2 ;
375 //Select 2x2 maximum sums to select L0
376 if(amp2 > ampmax2(0,mtru)){
377 ampmax2(0,mtru) = amp2 ;
378 ampmax2(1,mtru) = irow;
379 ampmax2(2,mtru) = icol;
384 //Find most recent time in the selected 2x2 cell
385 ampmax2(3,mtru) = 1 ;
386 Int_t row2 = static_cast <Int_t> (ampmax2(1,mtru));
387 Int_t col2 = static_cast <Int_t> (ampmax2(2,mtru));
388 for(Int_t i = 0; i<2; i++){
389 for(Int_t j = 0; j<2; j++){
390 if((*amptru)(row2+i,col2+j) > 0 && (*timeRtru)(row2+i,col2+j)> 0){
391 if((*timeRtru)(row2+i,col2+j) < ampmax2(3,mtru) )
392 ampmax2(3,mtru) = (*timeRtru)(row2+i,col2+j);
397 //Sliding nxn, add nxn amplitudes (OVERLAP)
399 for(Int_t irow = 0 ; irow < nCellsPhi/2; irow++){
400 for(Int_t icol = 0 ; icol < nCellsEta/2 ; icol++){
402 if( (irow+fPatchSize) < nCellsPhi/2 && (icol+fPatchSize) < nCellsEta/2){//Avoid exit the TRU
403 for(Int_t i = 0 ; i <= fPatchSize ; i++)
404 for(Int_t j = 0 ; j <= fPatchSize ; j++)
405 ampn += tru2x2(irow+i,icol+j);
406 //Select nxn maximum sums to select L1
407 if(ampn > ampmaxn(0,mtru)){
408 ampmaxn(0,mtru) = ampn ;
409 ampmaxn(1,mtru) = irow*2;
410 ampmaxn(2,mtru) = icol*2;
416 //Find most recent time in selected nxn cell
417 ampmaxn(3,mtru) = 1 ;
418 Int_t rown = static_cast <Int_t> (ampmaxn(1,mtru));
419 Int_t coln = static_cast <Int_t> (ampmaxn(2,mtru));
420 for(Int_t i = 0; i<4*fPatchSize; i++){
421 for(Int_t j = 0; j<4*fPatchSize; j++){
422 if( (rown+i) < nCellsPhi && (coln+j) < nCellsEta){//Avoid exit the TRU
423 if((*amptru)(rown+i,coln+j) > 0 && (*timeRtru)(rown+i,coln+j)> 0){
424 if((*timeRtru)(rown+i,coln+j) < ampmaxn(3,mtru) )
425 ampmaxn(3,mtru) = (*timeRtru)(rown+i,coln+j);
432 ampmaxn(0,mtru) = ampmax2(0,mtru);
433 ampmaxn(1,mtru) = ampmax2(1,mtru);
434 ampmaxn(2,mtru) = ampmax2(2,mtru);
435 ampmaxn(3,mtru) = ampmax2(3,mtru);
440 //____________________________________________________________________________
441 void AliEMCALTrigger::MakeSlidingTowers(const TClonesArray * amptrus, const TClonesArray * timeRtrus,
442 const Int_t isupermod,TMatrixD &max2, TMatrixD &maxn){
444 // Output from module (2x2 cells from one module)
445 Int_t nModulesPhi = fGeom->GetNModulesInTRUPhi(); // now 4 modules (3 div in phi)
447 nModulesPhi = nModulesPhi / 2 ; // Half size SM. Not Final.
449 Int_t nModulesEta = fGeom->GetNModulesInTRUEta(); // now 24 modules (no division in eta)
450 Int_t nTRU = fGeom->GetNTRU();
451 static int keyPrint = 0;
452 if(keyPrint) AliDebug(2,Form("MakeSlidingTowers : nTRU %i nModulesPhi %i nModulesEta %i ",
453 nTRU, nModulesPhi, nModulesEta ));
457 for(Int_t i = 0; i < 4; i++){
458 for(Int_t j = 0; j < nTRU; j++){
459 ampmax2(i,j) = ampmaxn(i,j) = -1;
463 // Create matrix that will contain 2x2 amplitude sums
464 // used to calculate the nxn sums
465 TMatrixD tru2x2(nModulesPhi/2,nModulesEta/2);
467 // Loop over all TRUS in a supermodule
468 for(Int_t itru = 0 + isupermod * nTRU ; itru < (isupermod+1)*nTRU ; itru++) {
469 TMatrixD * amptru = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
470 TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
471 Int_t mtru = itru - isupermod*nTRU ; // Number of TRU in Supermodule !!
473 if(!amptru || !timeRtru){
474 AliError("Amplitude or Time TRU matrix not available")
478 // Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
479 for(Int_t irow = 0 ; irow < nModulesPhi; irow +=2){
480 for(Int_t icol = 0 ; icol < nModulesEta ; icol +=2){
481 amp2 = (*amptru)(irow,icol) +(*amptru)(irow+1,icol)+
482 (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
484 //Fill matrix with added 2x2 towers for use in nxn sums
485 tru2x2(irow/2,icol/2) = amp2 ;
486 //Select 2x2 maximum sums to select L0
487 if(amp2 > ampmax2(0,mtru)){
488 ampmax2(0,mtru) = amp2 ;
489 ampmax2(1,mtru) = irow;
490 ampmax2(2,mtru) = icol;
495 ampmax2(3,mtru) = 0.;
497 // Find most recent time in the selected 2x2 towers
498 Int_t row2 = static_cast <Int_t> (ampmax2(1,mtru));
499 Int_t col2 = static_cast <Int_t> (ampmax2(2,mtru));
500 for(Int_t i = 0; i<2; i++){
501 for(Int_t j = 0; j<2; j++){
502 if((*amptru)(row2+i,col2+j) > 0 && (*timeRtru)(row2+i,col2+j)> 0){
503 if((*timeRtru)(row2+i,col2+j) > ampmax2(3,mtru) )
504 ampmax2(3,mtru) = (*timeRtru)(row2+i,col2+j); // max time
510 //Sliding nxn, add nxn amplitudes (OVERLAP)
512 for(Int_t irow = 0 ; irow < nModulesPhi/2; irow++){
513 for(Int_t icol = 0 ; icol < nModulesEta/2; icol++){
515 if( (irow+fPatchSize) < nModulesPhi/2 && (icol+fPatchSize) < nModulesEta/2){ //Avoid exit the TRU
516 for(Int_t i = 0 ; i <= fPatchSize ; i++)
517 for(Int_t j = 0 ; j <= fPatchSize ; j++)
518 ampn += tru2x2(irow+i,icol+j);
519 //Select nxn maximum sums to select L1
520 if(ampn > ampmaxn(0,mtru)){
521 ampmaxn(0,mtru) = ampn ;
522 ampmaxn(1,mtru) = irow;
523 ampmaxn(2,mtru) = icol;
529 ampmaxn(3,mtru) = 0.; // Was 1 , I don't know why
531 //Find most recent time in selected nxn cell
532 Int_t rown = static_cast <Int_t> (ampmaxn(1,mtru));
533 Int_t coln = static_cast <Int_t> (ampmaxn(2,mtru));
534 for(Int_t i = 0; i<4*fPatchSize; i++){
535 for(Int_t j = 0; j<4*fPatchSize; j++){
536 if( (rown+i) < nModulesPhi && (coln+j) < nModulesEta){//Avoid exit the TRU
537 if((*amptru)(rown+i,coln+j) > 0 && (*timeRtru)(rown+i,coln+j)> 0){
538 if((*timeRtru)(rown+i,coln+j) > ampmaxn(3,mtru) )
539 ampmaxn(3,mtru) = (*timeRtru)(rown+i,coln+j); // max time
545 } else { // copy 2x2 to nxn
546 ampmaxn(0,mtru) = ampmax2(0,mtru);
547 ampmaxn(1,mtru) = ampmax2(1,mtru);
548 ampmaxn(2,mtru) = ampmax2(2,mtru);
549 ampmaxn(3,mtru) = ampmax2(3,mtru);
552 if(keyPrint) AliDebug(2,Form(" : MakeSlidingTowers -OUt \n"));
555 //____________________________________________________________________________
556 void AliEMCALTrigger::Print(const Option_t * opt) const
559 //Prints main parameters
563 AliTriggerInput* in = 0x0 ;
564 AliInfo(Form(" fSimulation %i (input option) : #digits %i\n", fSimulation, fDigitsList->GetEntries()));
565 AliInfo(Form(" fTimeKey %i \n ", fTimeKey));
567 AliInfo(Form("\t Maximum Amplitude after Sliding Cell, \n")) ;
568 AliInfo(Form("\t -2x2 cells sum (not overlapped): %10.2f, in Super Module %d\n",
569 f2x2MaxAmp,f2x2SM)) ;
570 AliInfo(Form("\t -2x2 from row %d to row %d and from column %d to column %d\n", f2x2ModulePhi, f2x2ModulePhi+2, f2x2ModuleEta, f2x2ModuleEta+2));
571 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)));
573 AliInfo(Form("\t Patch Size, n x n: %d x %d cells\n",2*(fPatchSize+1), 2*(fPatchSize+1)));
574 AliInfo(Form("\t -nxn cells sum (overlapped) : %10.2f, in Super Module %d\n", fnxnMaxAmp,fnxnSM));
575 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)) ;
576 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) ));
579 AliInfo(Form("\t Isolate in SuperModule? %d\n", fIsolateInSuperModule)) ;
580 AliInfo(Form("\t Threshold for LO %10.2f\n", fL0Threshold));
582 in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_L0" );
584 AliInfo(Form("\t *** EMCAL LO is set ***\n"));
586 AliInfo(Form("\t Gamma Low Pt Threshold for L1 %10.2f\n", fL1GammaLowPtThreshold));
587 in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_GammaLPt_L1" );
589 AliInfo(Form("\t *** EMCAL Gamma Low Pt for L1 is set ***\n"));
591 AliInfo(Form("\t Gamma Medium Pt Threshold for L1 %10.2f\n", fL1GammaMediumPtThreshold));
592 in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_GammaMPt_L1" );
594 AliInfo(Form("\t *** EMCAL Gamma Medium Pt for L1 is set ***\n"));
596 AliInfo(Form("\t Gamma High Pt Threshold for L1 %10.2f\n", fL1GammaHighPtThreshold));
597 in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_GammaHPt_L1" );
599 AliInfo(Form("\t *** EMCAL Gamma High Pt for L1 is set ***\n")) ;
603 //____________________________________________________________________________
604 void AliEMCALTrigger::SetTriggers(const TClonesArray * ampmatrix,const Int_t iSM,
605 const TMatrixD &max2,
606 const TMatrixD &maxn)
608 //Checks the 2x2 and nxn maximum amplitude per each TRU and
609 //compares with the different L0 and L1 triggers thresholds
610 Float_t max2[] = {-1,-1,-1,-1} ;
611 Float_t maxn[] = {-1,-1,-1,-1} ;
615 Int_t nTRU = fGeom->GetNTRU();
617 //Find maximum summed amplitude of all the TRU
619 for(Int_t i = 0 ; i < nTRU ; i++){
620 if(max2[0] < ampmax2(0,i) ){
621 max2[0] = ampmax2(0,i) ; // 2x2 summed max amplitude
622 max2[1] = ampmax2(1,i) ; // corresponding phi position in TRU
623 max2[2] = ampmax2(2,i) ; // corresponding eta position in TRU
624 max2[3] = ampmax2(3,i) ; // corresponding most recent time
627 if(maxn[0] < ampmaxn(0,i) ){
628 maxn[0] = ampmaxn(0,i) ; // nxn summed max amplitude
629 maxn[1] = ampmaxn(1,i) ; // corresponding phi position in TRU
630 maxn[2] = ampmaxn(2,i) ; // corresponding eta position in TRU
631 maxn[3] = ampmaxn(3,i) ; // corresponding most recent time
636 //--------Set max amplitude if larger than in other Super Modules------------
637 Float_t maxtimeR2 = -1 ;
638 Float_t maxtimeRn = -1 ;
639 static AliEMCALRawUtils rawUtil;
640 Int_t nTimeBins = rawUtil.GetRawFormatTimeBins() ;
642 //Set max of 2x2 amplitudes and select L0 trigger
643 if(max2[0] > f2x2MaxAmp ){
644 // if(max2[0] > 5) printf(" L0 : iSM %i: max2[0] %5.0f : max2[3] %5.0f (maxtimeR2) \n",
645 // iSM, max2[0], max2[3]);
646 f2x2MaxAmp = max2[0] ;
648 maxtimeR2 = max2[3] ;
649 fGeom->GetModulePhiEtaIndexInSModuleFromTRUIndex(mtru2,
650 static_cast<Int_t>(max2[1]),
651 static_cast<Int_t>(max2[2]),
652 f2x2ModulePhi,f2x2ModuleEta);
655 if(fIsolateInSuperModule)
656 fIs2x2Isol = IsPatchIsolated(0, ampmatrix, iSM, mtru2, f2x2MaxAmp, f2x2ModulePhi,f2x2ModuleEta) ;
658 fIs2x2Isol = IsPatchIsolated(0, ampmatrix, iSM, mtru2, f2x2MaxAmp, static_cast<Int_t>(max2[1]), static_cast<Int_t>(max2[2])) ;
661 //Transform digit amplitude in Raw Samples
662 if (fADCValuesLow2x2 == 0) {
663 fADCValuesLow2x2 = new Int_t[nTimeBins];
664 fADCValuesHigh2x2 = new Int_t[nTimeBins];
666 //printf(" maxtimeR2 %12.5e (1)\n", maxtimeR2);
667 rawUtil.RawSampledResponse(maxtimeR2 * AliEMCALRawUtils::GetRawFormatTimeBin(),
668 f2x2MaxAmp, fADCValuesHigh2x2, fADCValuesLow2x2) ;
670 // Set Trigger Inputs, compare ADC time bins until threshold is attained
672 for(Int_t i = 0 ; i < nTimeBins ; i++){
673 // printf(" fADCValuesHigh2x2[%i] %i : %i \n", i, fADCValuesHigh2x2[i], fADCValuesLow2x2[i]);
674 if(fADCValuesHigh2x2[i] >= fL0Threshold || fADCValuesLow2x2[i] >= fL0Threshold){
675 SetInput("EMCAL_L0") ;
680 // Nov 5 - no analysis of time information
681 if(f2x2MaxAmp >= fL0Threshold) { // should add the low amp too
682 SetInput("EMCAL_L0");
687 //------------Set max of nxn amplitudes and select L1 trigger---------
688 if(maxn[0] > fnxnMaxAmp ){
689 fnxnMaxAmp = maxn[0] ;
691 maxtimeRn = maxn[3] ;
692 fGeom->GetModulePhiEtaIndexInSModuleFromTRUIndex(mtrun,
693 static_cast<Int_t>(maxn[1]),
694 static_cast<Int_t>(maxn[2]),
695 fnxnModulePhi,fnxnModuleEta) ;
698 if(fIsolateInSuperModule)
699 fIsnxnIsol = IsPatchIsolated(1, ampmatrix, iSM, mtrun, fnxnMaxAmp, fnxnModulePhi, fnxnModuleEta) ;
701 fIsnxnIsol = IsPatchIsolated(1, ampmatrix, iSM, mtrun, fnxnMaxAmp, static_cast<Int_t>(maxn[1]), static_cast<Int_t>(maxn[2])) ;
704 //Transform digit amplitude in Raw Samples
705 if (fADCValuesLownxn == 0) {
706 fADCValuesHighnxn = new Int_t[nTimeBins];
707 fADCValuesLownxn = new Int_t[nTimeBins];
709 rawUtil.RawSampledResponse(maxtimeRn * AliEMCALRawUtils::GetRawFormatTimeBin(),
710 fnxnMaxAmp, fADCValuesHighnxn, fADCValuesLownxn) ;
712 //Set Trigger Inputs, compare ADC time bins until threshold is attained
714 for(Int_t i = 0 ; i < nTimeBins ; i++){
715 if(fADCValuesHighnxn[i] >= fL1GammaLowPtThreshold || fADCValuesLownxn[i] >= fL1GammaLowPtThreshold){
716 SetInput("EMCAL_GammaLPt_L1") ;
722 for(Int_t i = 0 ; i < nTimeBins ; i++){
723 if(fADCValuesHighnxn[i] >= fL1GammaMediumPtThreshold || fADCValuesLownxn[i] >= fL1GammaMediumPtThreshold){
724 SetInput("EMCAL_GammaMPt_L1") ;
730 for(Int_t i = 0 ; i < nTimeBins ; i++){
731 if(fADCValuesHighnxn[i] >= fL1GammaHighPtThreshold || fADCValuesLownxn[i] >= fL1GammaHighPtThreshold){
732 SetInput("EMCAL_GammaHPt_L1") ;
737 // Nov 5 - no analysis of time information
738 if(fnxnMaxAmp >= fL1GammaLowPtThreshold) { // should add the low amp too
739 SetInput("EMCAL_GammaLPt_L1") ; //SetL1 Low
741 if(fnxnMaxAmp >= fL1GammaMediumPtThreshold) { // should add the low amp too
742 SetInput("EMCAL_GammaMPt_L1") ; //SetL1 Medium
744 if(fnxnMaxAmp >= fL1GammaHighPtThreshold) { // should add the low amp too
745 SetInput("EMCAL_GammaHPt_L1") ; //SetL1 High
751 //____________________________________________________________________________
752 void AliEMCALTrigger::FillTRU(const TClonesArray * digits, TClonesArray * ampmatrix, TClonesArray * ampmatrixsmod, TClonesArray * timeRmatrix) {
754 // Orders digits ampitudes list in fNTRU TRUs (384 cells) per supermodule.
755 // Each TRU is a TMatrixD, and they are kept in TClonesArrays. The number of
756 // TRU in phi is fNTRUPhi, and the number of TRU in eta is fNTRUEta.
757 // Last 2 modules are half size in Phi, I considered that the number of TRU
758 // is maintained for the last modules but decision not taken. If different,
759 // then this must be changed. Also fill a matrix with all amplitudes in supermodule for isolation studies.
761 // Initilize and declare variables
762 // List of TRU matrices initialized to 0.
763 // printf("<I> AliEMCALTrigger::FillTRU() started : # digits %i\n", digits->GetEntriesFast());
766 // One input per EMCAL module so size of matrix is reduced by 4 (2x2 division case)
768 Int_t nPhi = fGeom->GetNPhi();
769 Int_t nZ = fGeom->GetNZ();
770 Int_t nTRU = fGeom->GetNTRU();
771 // Int_t nTRUPhi = fGeom->GetNTRUPhi();
772 Int_t nModulesPhi = fGeom->GetNModulesInTRUPhi();
773 Int_t nModulesPhi2 = fGeom->GetNModulesInTRUPhi();
774 Int_t nModulesEta = fGeom->GetNModulesInTRUEta();
775 // printf("<I> AliEMCALTrigger::FillTRU() nTRU %i nTRUPhi %i : nModulesPhi %i nModulesEta %i \n",
776 // nTRU, nTRUPhi, nModulesPhi, nModulesEta);
787 // iphim, ietam - module indexes in SM
791 //List of TRU matrices initialized to 0.
792 Int_t nSup = fGeom->GetNumberOfSuperModules();
793 for(Int_t k = 0; k < nTRU*nSup; k++){
794 TMatrixD amptrus(nModulesPhi,nModulesEta) ;
795 TMatrixD timeRtrus(nModulesPhi,nModulesEta) ;
796 // Do we need to initialise? I think TMatrixD does it by itself...
797 for(Int_t i = 0; i < nModulesPhi; i++){
798 for(Int_t j = 0; j < nModulesEta; j++){
800 timeRtrus(i,j) = 0.0;
803 new((*ampmatrix)[k]) TMatrixD(amptrus) ;
804 new((*timeRmatrix)[k]) TMatrixD(timeRtrus) ;
807 // List of Modules matrices initialized to 0.
808 for(Int_t k = 0; k < nSup ; k++){
810 // if(nSup>9) mphi = nPhi/2; // the same size
811 TMatrixD ampsmods( mphi, nZ);
812 for(Int_t i = 0; i < mphi; i++){
813 for(Int_t j = 0; j < nZ; j++){
817 new((*ampmatrixsmod)[k]) TMatrixD(ampsmods) ;
820 AliEMCALDigit * dig ;
822 //Digits loop to fill TRU matrices with amplitudes.
823 for(Int_t idig = 0 ; idig < digits->GetEntriesFast() ; idig++){
825 dig = dynamic_cast<AliEMCALDigit *>(digits->At(idig)) ;
827 amp = Float_t(dig->GetAmplitude()); // Energy of the digit (arbitrary units)
828 id = dig->GetId() ; // Id label of the cell
829 timeR = dig->GetTimeR() ; // Earliest time of the digit
830 if(amp<=0.0) AliDebug(1,Form(" id %i amp %f \n", id, amp));
831 // printf(" FILLTRU : timeR %10.5e time %10.5e : amp %10.5e \n", timeR, dig->GetTime(), amp);
832 // Get eta and phi cell position in supermodule
833 Bool_t bCell = fGeom->GetCellIndex(id, iSupMod, nModule, nIphi, nIeta) ;
835 AliError(Form("%i Wrong cell id number %i ", idig, id)) ;
837 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
838 // iphim, ietam - module indexes in SM
839 fGeom->GetModuleIndexesFromCellIndexesInSModule(iSupMod,iphi,ieta, iphim, ietam, nModule);
841 //printf("iSupMod %i nModule %i iphi %i ieta %i iphim %i ietam %i \n",
842 //iSupMod,nModule, iphi, ieta, iphim, ietam);
844 // Check to which TRU in the supermodule belongs the cell.
845 // Supermodules are divided in a TRU matrix of dimension
846 // (fNTRUPhi,fNTRUEta).
847 // Each TRU is a cell matrix of dimension (nModulesPhi,nModulesEta)
849 // First calculate the row and column in the supermodule
850 // of the TRU to which the cell belongs.
851 Int_t row = iphim / nModulesPhi;
852 Int_t col = ietam / nModulesEta;
853 //Calculate label number of the TRU
854 Int_t itru = fGeom->GetAbsTRUNumberFromNumberInSm(row, col, iSupMod);
856 //Fill TRU matrix with cell values
857 TMatrixD * amptrus = dynamic_cast<TMatrixD *>(ampmatrix->At(itru)) ;
858 TMatrixD * timeRtrus = dynamic_cast<TMatrixD *>(timeRmatrix->At(itru)) ;
860 if(!amptrus || !timeRtrus){
861 AliError("Could not recover the TRU matrix with amplitudes or times");
864 //Calculate row and column of the module inside the TRU with number itru
865 Int_t irow = iphim - row * nModulesPhi;
867 irow = iphim - row * nModulesPhi2; // size of matrix the same
868 Int_t icol = ietam - col * nModulesEta;
870 (*amptrus)(irow,icol) += amp ;
871 if((*timeRtrus)(irow,icol) <0.0 || (*timeRtrus)(irow,icol) <= timeR){ // ??
872 (*timeRtrus)(irow,icol) = timeR ;
875 //printf(" ieta %i iphi %i iSM %i || col %i row %i : itru %i -> amp %f\n",
876 // ieta, iphi, iSupMod, col, row, itru, amp);
877 //####################SUPERMODULE MATRIX ##################
878 TMatrixD * ampsmods = dynamic_cast<TMatrixD *>(ampmatrixsmod->At(iSupMod)) ;
880 AliError("Could not recover the matrix per SM");
883 (*ampsmods)(iphim,ietam) += amp ;
884 // printf(" id %i iphim %i ietam %i SM %i : irow %i icol %i itru %i : amp %6.0f\n",
885 //id, iphim, ietam, iSupMod, irow, icol, itru, amp);
887 else AliError("Could not recover the digit");
890 //printf("<I> AliEMCALTrigger::FillTRU() is ended \n");
893 //____________________________________________________________________________
894 void AliEMCALTrigger::Trigger()
896 //Main Method to select triggers.
897 TH1::AddDirectory(0);
899 AliRunLoader *runLoader = AliRunLoader::Instance();
900 AliEMCALLoader *emcalLoader = 0;
902 emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
905 //Load EMCAL Geometry
906 if (runLoader && runLoader->GetAliRun()){
907 AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(runLoader->GetAliRun()->GetDetector("EMCAL"));
908 if(emcal)fGeom = emcal->GetGeometry();
912 fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
915 AliFatal("Did not get geometry from EMCALLoader");
918 Int_t nSuperModules = fGeom->GetNumberOfSuperModules() ; //12 SM in EMCAL
919 Int_t nTRU = fGeom->GetNTRU(); // 3 TRU per super module
921 //Intialize data members each time the trigger is called in event loop
922 f2x2MaxAmp = -1; f2x2ModulePhi = -1; f2x2ModuleEta = -1;
923 fnxnMaxAmp = -1; fnxnModulePhi = -1; fnxnModuleEta = -1;
925 // Take the digits list if simulation
926 if(fSimulation && runLoader && emcalLoader){ // works than run seperate macros
927 runLoader->LoadDigits("EMCAL");
928 fDigitsList = emcalLoader->Digits() ;
929 runLoader->LoadSDigits("EMCAL");
931 // Digits list should be set by method SetDigitsList(TClonesArray * digits)
933 AliFatal("Digits not found !") ;
935 //Take the digits list
937 // Delete old if unzero
938 if(fAmpTrus) {fAmpTrus->Delete(); delete fAmpTrus;}
939 if(fTimeRtrus) {fTimeRtrus->Delete(); delete fTimeRtrus;}
940 if(fAmpSMods) {fAmpSMods->Delete(); delete fAmpSMods;}
941 // Fill TRU and SM matrix
942 fAmpTrus = new TClonesArray("TMatrixD",nTRU);
943 fAmpTrus->SetName("AmpTrus");
944 fTimeRtrus = new TClonesArray("TMatrixD",nTRU);
945 fTimeRtrus->SetName("TimeRtrus");
946 fAmpSMods = new TClonesArray("TMatrixD",nSuperModules);
947 fAmpSMods->SetName("AmpSMods");
949 FillTRU(fDigitsList, fAmpTrus, fAmpSMods, fTimeRtrus);
951 // Jet stuff - only one case, no freedom here
952 if(fGeom->GetNEtaSubOfTRU() == 6) {
953 if(fAmpJetMatrix) {delete fAmpJetMatrix; fAmpJetMatrix=0;}
954 if(fJetMatrixE) {delete fJetMatrixE; fJetMatrixE=0;}
956 fAmpJetMatrix = new TMatrixD(17,12); // 17-phi(row), 12-eta(col)
957 fJetMatrixE = new TH2F("fJetMatrixE"," E of max patch in (#phi,#eta)",
958 17, 80.*TMath::DegToRad(), (180.+20.*2/3.)*TMath::DegToRad(), 12, -0.7, 0.7);
959 for(Int_t row=0; row<fAmpJetMatrix->GetNrows(); row++) {
960 for(Int_t col=0; col<fAmpJetMatrix->GetNcols(); col++) {
961 (*fAmpJetMatrix)(row,col) = 0.;
964 FillJetMatrixFromSMs(fAmpSMods, fAmpJetMatrix, fGeom);
966 if(!CheckConsistentOfMatrixes()) assert(0);
968 // Do Tower Sliding and select Trigger
969 // Initialize varible that will contain maximum amplitudes and
970 // its corresponding tower position in eta and phi, and time.
971 TMatrixD ampmax2(4,nTRU) ; // 0-max amp, 1-irow, 2-icol, 3-timeR
972 TMatrixD ampmaxn(4,nTRU) ;
974 for(Int_t iSM = 0 ; iSM < nSuperModules ; iSM++) {
975 //Do 2x2 and nxn sums, select maximums.
977 MakeSlidingTowers(fAmpTrus, fTimeRtrus, iSM, ampmax2, ampmaxn);
980 if(fIsolateInSuperModule) // here some discripency between tru and SM
981 SetTriggers(fAmpSMods,iSM,ampmax2,ampmaxn) ;
982 if(!fIsolateInSuperModule)
983 SetTriggers(fAmpTrus,iSM,ampmax2,ampmaxn) ;
986 // Do patch sliding and select Jet Trigger
987 // 0-max amp-meanFromVZERO(if), 1-irow, 2-icol, 3-timeR,
988 // 4-max amp , 5-meanFromVZERO (Nov 25, 2007)
990 MakeSlidingPatch((*fAmpJetMatrix), fNJetPatchPhi, fAmpJetMax); // no timing information here
996 //____________________________________________________________________________
997 void AliEMCALTrigger::GetTriggerInfo(TArrayF &triggerPosition, TArrayF &triggerAmplitudes) const
999 // Template - should be defined; Nov 5, 2007
1000 triggerPosition[0] = 0.;
1001 triggerAmplitudes[0] = 0.;
1004 //____________________________________________________________________________
1005 void AliEMCALTrigger::FillJetMatrixFromSMs(TClonesArray *ampmatrixsmod, TMatrixD* jetMat, AliEMCALGeometry *g)
1008 // Fill matrix for jet trigger from SM matrixes of modules
1010 static int keyPrint = 0;
1012 if(ampmatrixsmod==0 || jetMat==0 || g==0) return;
1013 Double_t amp = 0.0, ampSum=0.0;
1015 Int_t nEtaModSum = g->GetNZ() / g->GetNEtaSubOfTRU(); // should be 4
1016 Int_t nPhiModSum = g->GetNPhi() / g->GetNTRUPhi(); // should be 4
1018 if(keyPrint) AliDebug(2,Form("%s",Form(" AliEMCALTrigger::FillJetMatrixFromSMs | nEtaModSum %i : nPhiModSum %i \n", nEtaModSum, nPhiModSum)));
1019 Int_t jrow=0, jcol=0; // indexes of jet matrix
1020 Int_t nEtaSM=0, nPhiSM=0;
1021 for(Int_t iSM=0; iSM<ampmatrixsmod->GetEntries(); iSM++) {
1022 TMatrixD * ampsmods = dynamic_cast<TMatrixD *>(ampmatrixsmod->At(iSM));
1024 if(!ampsmods) return;
1026 Int_t nrow = ampsmods->GetNrows();
1027 Int_t ncol = ampsmods->GetNcols();
1028 //printf("%s",Form(" ######## SM %i : nrow %i : ncol %i ##### \n", iSM, nrow, ncol));
1029 for(Int_t row=0; row<nrow; row++) {
1030 for(Int_t col=0; col<ncol; col++) {
1031 amp = (*ampsmods)(row,col);
1035 if(keyPrint) AliDebug(2,Form("%s",Form(" ** nPhiSm %i : nEtaSM %i : row %2.2i : col %2.2i -> ", nPhiSM, nEtaSM, row, col)));
1036 if(nEtaSM == 0) { // positive Z
1037 jrow = 3*nPhiSM + row/nPhiModSum;
1038 jcol = 6 + col / nEtaModSum;
1039 } else { // negative Z
1040 if(iSM<=9) jrow = 3*nPhiSM + 2 - row/nPhiModSum;
1041 else jrow = 3*nPhiSM + 1 - row/nPhiModSum; // half size
1042 jcol = 5 - col / nEtaModSum;
1044 if(keyPrint) AliDebug(2,Form("%s",Form(" jrow %2.2i : jcol %2.2i : amp %f (jetMat) \n", jrow, jcol, amp)));
1046 (*jetMat)(jrow,jcol) += amp;
1047 ampSum += amp; // For controling
1048 } else if(amp<0.0) {
1049 AliDebug(1,Form(" jrow %2.2i : jcol %2.2i : amp %f (jetMat: amp<0) \n", jrow, jcol, amp));
1055 if(ampSum <= 0.0) AliDebug(1,Form("ampSum %f (<=0.0) ", ampSum));
1058 //____________________________________________________________________________
1059 void AliEMCALTrigger::MakeSlidingPatch(const TMatrixD &jm, const Int_t nPatchSize, TMatrixD &JetMax)
1061 // Sliding patch : nPatchSize x nPatchSize (OVERLAP)
1062 static int keyPrint = 0;
1063 if(keyPrint) AliDebug(2,Form(" AliEMCALTrigger::MakeSlidingPatch() was started \n"));
1064 Double_t ampCur = 0.0, e=0.0;
1065 ampJetMax(0,0) = 0.0;
1066 ampJetMax(3,0) = 0.0; // unused now
1067 ampJetMax(4,0) = ampJetMax(5,0) = 0.0;
1068 for(Int_t row=0; row<fAmpJetMatrix->GetNrows(); row ++) {
1069 for(Int_t col=0; col<fAmpJetMatrix->GetNcols(); col++) {
1071 // check on patch size
1072 if( (row+nPatchSize-1) < fAmpJetMatrix->GetNrows() && (col+nPatchSize-1) < fAmpJetMatrix->GetNcols()){
1073 for(Int_t i = 0 ; i < nPatchSize ; i++) {
1074 for(Int_t j = 0 ; j < nPatchSize ; j++) {
1075 ampCur += jm(row+i, col+j);
1077 } // end cycle on patch
1078 if(ampCur > ampJetMax(0,0)){
1079 ampJetMax(0,0) = ampCur;
1080 ampJetMax(1,0) = row;
1081 ampJetMax(2,0) = col;
1083 } // check on patch size
1086 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));
1088 Double_t eCorrJetMatrix=0.0;
1089 if(fVZER0Mult > 0.0) {
1090 // Correct patch energy (adc) and jet patch matrix energy
1091 Double_t meanAmpBG = GetMeanEmcalPatchEnergy(Int_t(fVZER0Mult), nPatchSize)/0.0153;
1092 ampJetMax(4,0) = ampJetMax(0,0);
1093 ampJetMax(5,0) = meanAmpBG;
1095 Double_t eCorr = ampJetMax(0,0) - meanAmpBG;
1096 AliDebug(2,Form(" ampJetMax(0,0) %f meanAmpBG %f eCorr %f : ampJetMax(4,0) %f \n",
1097 ampJetMax(0,0), meanAmpBG, eCorr, ampJetMax(5,0)));
1098 ampJetMax(0,0) = eCorr;
1100 eCorrJetMatrix = GetMeanEmcalEnergy(Int_t(fVZER0Mult)) / 208.;
1102 // Fill patch energy matrix
1103 for(int row=Int_t(ampJetMax(1,0)); row<Int_t(ampJetMax(1,0))+nPatchSize; row++) {
1104 for(int col=Int_t(ampJetMax(2,0)); col<Int_t(ampJetMax(2,0))+nPatchSize; col++) {
1105 e = Double_t(jm(row,col)*0.0153); // 0.0153 - hard coded now
1106 if(eCorrJetMatrix > 0.0) { // BG subtraction case
1107 e -= eCorrJetMatrix;
1108 fJetMatrixE->SetBinContent(row+1, col+1, e);
1109 } else if(e > 0.0) {
1110 fJetMatrixE->SetBinContent(row+1, col+1, e);
1114 // PrintJetMatrix();
1115 // Set the jet trigger(s), multiple threshold now, Nov 19,2007
1116 for(Int_t i=0; i<fNJetThreshold; i++ ) {
1117 if(ampJetMax(0,0) >= fL1JetThreshold[i]) {
1118 SetInput(GetNameOfJetTrigger(i));
1123 //____________________________________________________________________________
1124 Double_t AliEMCALTrigger::GetEmcalSumAmp() const
1126 // Return sum of amplidutes from EMCal
1127 // Used calibration coefficeint for transition to energy
1128 return fAmpJetMatrix >0 ?fAmpJetMatrix->Sum() :0.0;
1131 //____________________________________________________________________________
1132 void AliEMCALTrigger::PrintJetMatrix() const
1134 // fAmpJetMatrix : (17,12); // 17-phi(row), 12-eta(col)
1135 if(fAmpJetMatrix == 0) return;
1137 AliInfo(Form("\n #### jetMatrix : (%i,%i) ##### \n ",
1138 fAmpJetMatrix->GetNrows(), fAmpJetMatrix->GetNcols()));
1139 PrintMatrix(*fAmpJetMatrix);
1142 //____________________________________________________________________________
1143 void AliEMCALTrigger::PrintAmpTruMatrix(Int_t ind) const
1145 // Print matrix with TRU patches
1146 TMatrixD * tru = dynamic_cast<TMatrixD *>(fAmpTrus->At(ind));
1147 if(tru == 0) return;
1148 AliInfo(Form("\n #### Amp TRU matrix(%i) : (%i,%i) ##### \n ",
1149 ind, tru->GetNrows(), tru->GetNcols()));
1153 //____________________________________________________________________________
1154 void AliEMCALTrigger::PrintAmpSmMatrix(Int_t ind) const
1156 // Print matrix with SM amplitudes
1157 TMatrixD * sm = dynamic_cast<TMatrixD *>(fAmpSMods->At(ind));
1159 AliInfo(Form("\n #### Amp SM matrix(%i) : (%i,%i) ##### \n ",
1160 ind, sm->GetNrows(), sm->GetNcols()));
1164 //____________________________________________________________________________
1165 void AliEMCALTrigger::PrintMatrix(const TMatrixD &mat) const
1167 //Print matrix object
1168 for(Int_t col=0; col<mat.GetNcols(); col++) AliInfo(Form(" %3i ", col));
1169 AliInfo(Form("\n -- \n"));
1170 for(Int_t row=0; row<mat.GetNrows(); row++) {
1171 AliInfo(Form(" row:%2i ", row));
1172 for(Int_t col=0; col<mat.GetNcols(); col++) {
1173 AliInfo(Form(" %4i", (Int_t)mat(row,col)));
1179 //____________________________________________________________________________
1180 Bool_t AliEMCALTrigger::CheckConsistentOfMatrixes(const Int_t pri)
1182 // Check consitency of matrices
1183 Double_t sumSM = 0.0, smCur=0.0;
1184 Double_t sumTru=0.0, sumTruInSM = 0.0, truSum=0.0;
1185 // Bool_t key = kTRUE;
1187 for(Int_t i=0; i<fAmpSMods->GetEntries(); i++) {
1188 TMatrixD * sm = dynamic_cast<TMatrixD *>(fAmpSMods->At(i));
1194 for(Int_t itru=0; itru<3; itru++) { // Cycle on tru inside SM
1195 Int_t ind = 3*i + itru;
1196 TMatrixD *tru = dynamic_cast<TMatrixD *>(fAmpTrus->At(ind));
1198 truSum = tru->Sum();
1199 sumTruInSM += truSum;
1202 sumTru += sumTruInSM;
1204 if(sumTruInSM != smCur) {
1205 AliDebug(1,Form(" sm %i : smCur %f -> sumTruInSM %f \n", i, smCur, sumTruInSM));
1210 Double_t sumJetMat = fAmpJetMatrix->Sum();
1211 if(pri || TMath::Abs(sumSM-sumTru)>0.0001 || TMath::Abs(sumSM-sumJetMat) > 0.0001)
1212 AliDebug(1,Form(" sumSM %f : sumTru %f : sumJetMat %f \n", sumSM, sumTru, sumJetMat));
1213 if(TMath::Abs(sumSM - sumTru)>0.0001 || TMath::Abs(sumSM-sumJetMat) > 0.0001) return kFALSE;
1217 //____________________________________________________________________________
1218 void AliEMCALTrigger::Browse(TBrowser* b)
1221 if(&fInputs) b->Add(&fInputs);
1222 if(fAmpTrus) b->Add(fAmpTrus);
1223 if(fTimeRtrus) b->Add(fTimeRtrus);
1224 if(fAmpSMods) b->Add(fAmpSMods);
1225 if(fAmpJetMatrix) b->Add(fAmpJetMatrix);
1226 if(fJetMatrixE) b->Add(fJetMatrixE);