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)
248 // Get matrix of TRU or Module with maximum amplitude patch.
249 Int_t itru = mtru + iSM * fGeom->GetNTRU(); //number of tru, min 0 max 3*12=36.
250 TMatrixD * ampmatrix = 0x0;
253 static int keyPrint = 0;
254 if(keyPrint) AliDebug(2,Form(" IsPatchIsolated : iSM %i mtru %i itru %i maxphi %i maxeta %i \n", iSM, mtru, itru, maxphi, maxeta));
256 if(fIsolateInSuperModule){ // ?
257 ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(iSM)) ;
258 rowborder = fGeom->GetNPhi();
259 colborder = fGeom->GetNZ();
260 AliDebug(2,"Isolate trigger in Module");
262 ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(itru)) ;
263 rowborder = fGeom->GetNModulesInTRUPhi();
264 colborder = fGeom->GetNModulesInTRUEta();
265 AliDebug(2,"Isolate trigger in TRU");
267 if(iSM>9) rowborder /= 2; // half size in phi
269 //Define patch modules - what is this ??
270 Int_t isolmodules = fIsolPatchSize*(1+iPatchType);
271 Int_t ipatchmodules = 2*(1+fPatchSize*iPatchType);
272 Int_t minrow = maxphi - isolmodules;
273 Int_t mincol = maxeta - isolmodules;
274 Int_t maxrow = maxphi + isolmodules + ipatchmodules;
275 Int_t maxcol = maxeta + isolmodules + ipatchmodules;
277 minrow = minrow<0?0 :minrow;
278 mincol = mincol<0?0 :mincol;
280 maxrow = maxrow>rowborder?rowborder :maxrow;
281 maxcol = maxcol>colborder?colborder :maxcol;
283 //printf("%s\n",Form("Number of added Isol Modules %d, Patch Size %d",isolmodules, ipatchmodules));
284 //printf("%s\n",Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
285 // AliDebug(2,Form("Number of added Isol Modules %d, Patch Size %d",isolmodules, ipatchmodules));
286 //AliDebug(2,Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
288 //Add amplitudes in all isolation patch
290 for(Int_t irow = minrow ; irow < maxrow; irow ++)
291 for(Int_t icol = mincol ; icol < maxcol ; icol ++)
292 amp += (*ampmatrix)(irow,icol);
294 AliDebug(2,Form("Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
297 // AliError(Form("Bad sum: Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
298 // ampmatrix->Print();
301 amp-=maxamp; //Calculate energy in isolation patch that do not comes from maximum patch.
304 AliDebug(2, Form("Maximum amplitude %f, Out of patch %f",maxamp, amp));
306 //Fill isolation amplitude data member and say if patch is isolated.
307 if(iPatchType == 0){ //2x2 case
308 f2x2AmpOutOfPatch = amp;
309 if(amp < f2x2AmpOutOfPatchThres) b=kTRUE;
310 } else if(iPatchType == 1){ //nxn case
311 fnxnAmpOutOfPatch = amp;
312 if(amp < fnxnAmpOutOfPatchThres) b=kTRUE;
315 if(keyPrint) AliDebug(2,Form(" IsPatchIsolated - OUT \n"));
322 //____________________________________________________________________________
323 void AliEMCALTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus, const Int_t isupermod,TMatrixD &max2, TMatrixD &maxn){
325 //Sums energy of all possible 2x2 (L0) and nxn (L1) modules per each TRU.
326 //Fast signal in the experiment is given by 2x2 modules,
327 //for this reason we loop inside the TRU modules by 2.
329 //Declare and initialize variables
330 Int_t nCellsPhi = fGeom->GetNCellsInTRUPhi();
332 nCellsPhi = nCellsPhi / 2 ; //Half size SM. Not Final.
333 // 12(tow)*2(cell)/1 TRU, cells in Phi in one TRU
334 Int_t nCellsEta = fGeom->GetNCellsInTRUEta();
335 Int_t nTRU = fGeom->GetNTRU();
336 // 24(mod)*2(tower)/3 TRU, cells in Eta in one TRU
337 //Int_t nTRU = geom->GeNTRU();//3 TRU per super module
341 for(Int_t i = 0; i < 4; i++){
342 for(Int_t j = 0; j < nTRU; j++){
348 //Create matrix that will contain 2x2 amplitude sums
349 //used to calculate the nxn sums
350 TMatrixD tru2x2(nCellsPhi/2,nCellsEta/2) ;
351 for(Int_t i = 0; i < nCellsPhi/2; i++)
352 for(Int_t j = 0; j < nCellsEta/2; j++)
355 //Loop over all TRUS in a supermodule
356 for(Int_t itru = 0 + isupermod * nTRU ; itru < (isupermod+1)*nTRU ; itru++) {
357 TMatrixD * amptru = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
358 TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
359 Int_t mtru = itru-isupermod*nTRU ; //Number of TRU in Supermodule
361 //Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
362 for(Int_t irow = 0 ; irow < nCellsPhi; irow += 2){
363 for(Int_t icol = 0 ; icol < nCellsEta ; icol += 2){
364 amp2 = (*amptru)(irow,icol)+(*amptru)(irow+1,icol)+
365 (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
367 //Fill matrix with added 2x2 cells for use in nxn sums
368 tru2x2(irow/2,icol/2) = amp2 ;
369 //Select 2x2 maximum sums to select L0
370 if(amp2 > ampmax2(0,mtru)){
371 ampmax2(0,mtru) = amp2 ;
372 ampmax2(1,mtru) = irow;
373 ampmax2(2,mtru) = icol;
378 //Find most recent time in the selected 2x2 cell
379 ampmax2(3,mtru) = 1 ;
380 Int_t row2 = static_cast <Int_t> (ampmax2(1,mtru));
381 Int_t col2 = static_cast <Int_t> (ampmax2(2,mtru));
382 for(Int_t i = 0; i<2; i++){
383 for(Int_t j = 0; j<2; j++){
384 if((*amptru)(row2+i,col2+j) > 0 && (*timeRtru)(row2+i,col2+j)> 0){
385 if((*timeRtru)(row2+i,col2+j) < ampmax2(3,mtru) )
386 ampmax2(3,mtru) = (*timeRtru)(row2+i,col2+j);
391 //Sliding nxn, add nxn amplitudes (OVERLAP)
393 for(Int_t irow = 0 ; irow < nCellsPhi/2; irow++){
394 for(Int_t icol = 0 ; icol < nCellsEta/2 ; icol++){
396 if( (irow+fPatchSize) < nCellsPhi/2 && (icol+fPatchSize) < nCellsEta/2){//Avoid exit the TRU
397 for(Int_t i = 0 ; i <= fPatchSize ; i++)
398 for(Int_t j = 0 ; j <= fPatchSize ; j++)
399 ampn += tru2x2(irow+i,icol+j);
400 //Select nxn maximum sums to select L1
401 if(ampn > ampmaxn(0,mtru)){
402 ampmaxn(0,mtru) = ampn ;
403 ampmaxn(1,mtru) = irow*2;
404 ampmaxn(2,mtru) = icol*2;
410 //Find most recent time in selected nxn cell
411 ampmaxn(3,mtru) = 1 ;
412 Int_t rown = static_cast <Int_t> (ampmaxn(1,mtru));
413 Int_t coln = static_cast <Int_t> (ampmaxn(2,mtru));
414 for(Int_t i = 0; i<4*fPatchSize; i++){
415 for(Int_t j = 0; j<4*fPatchSize; j++){
416 if( (rown+i) < nCellsPhi && (coln+j) < nCellsEta){//Avoid exit the TRU
417 if((*amptru)(rown+i,coln+j) > 0 && (*timeRtru)(rown+i,coln+j)> 0){
418 if((*timeRtru)(rown+i,coln+j) < ampmaxn(3,mtru) )
419 ampmaxn(3,mtru) = (*timeRtru)(rown+i,coln+j);
426 ampmaxn(0,mtru) = ampmax2(0,mtru);
427 ampmaxn(1,mtru) = ampmax2(1,mtru);
428 ampmaxn(2,mtru) = ampmax2(2,mtru);
429 ampmaxn(3,mtru) = ampmax2(3,mtru);
434 //____________________________________________________________________________
435 void AliEMCALTrigger::MakeSlidingTowers(const TClonesArray * amptrus, const TClonesArray * timeRtrus,
436 const Int_t isupermod,TMatrixD &max2, TMatrixD &maxn){
438 // Output from module (2x2 cells from one module)
439 Int_t nModulesPhi = fGeom->GetNModulesInTRUPhi(); // now 4 modules (3 div in phi)
441 nModulesPhi = nModulesPhi / 2 ; // Half size SM. Not Final.
443 Int_t nModulesEta = fGeom->GetNModulesInTRUEta(); // now 24 modules (no division in eta)
444 Int_t nTRU = fGeom->GetNTRU();
445 static int keyPrint = 0;
446 if(keyPrint) AliDebug(2,Form("MakeSlidingTowers : nTRU %i nModulesPhi %i nModulesEta %i ",
447 nTRU, nModulesPhi, nModulesEta ));
451 for(Int_t i = 0; i < 4; i++){
452 for(Int_t j = 0; j < nTRU; j++){
453 ampmax2(i,j) = ampmaxn(i,j) = -1;
457 // Create matrix that will contain 2x2 amplitude sums
458 // used to calculate the nxn sums
459 TMatrixD tru2x2(nModulesPhi/2,nModulesEta/2);
461 // Loop over all TRUS in a supermodule
462 for(Int_t itru = 0 + isupermod * nTRU ; itru < (isupermod+1)*nTRU ; itru++) {
463 TMatrixD * amptru = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
464 TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
465 Int_t mtru = itru - isupermod*nTRU ; // Number of TRU in Supermodule !!
467 // Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
468 for(Int_t irow = 0 ; irow < nModulesPhi; irow +=2){
469 for(Int_t icol = 0 ; icol < nModulesEta ; icol +=2){
470 amp2 = (*amptru)(irow,icol) +(*amptru)(irow+1,icol)+
471 (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
473 //Fill matrix with added 2x2 towers for use in nxn sums
474 tru2x2(irow/2,icol/2) = amp2 ;
475 //Select 2x2 maximum sums to select L0
476 if(amp2 > ampmax2(0,mtru)){
477 ampmax2(0,mtru) = amp2 ;
478 ampmax2(1,mtru) = irow;
479 ampmax2(2,mtru) = icol;
484 ampmax2(3,mtru) = 0.;
486 // Find most recent time in the selected 2x2 towers
487 Int_t row2 = static_cast <Int_t> (ampmax2(1,mtru));
488 Int_t col2 = static_cast <Int_t> (ampmax2(2,mtru));
489 for(Int_t i = 0; i<2; i++){
490 for(Int_t j = 0; j<2; j++){
491 if((*amptru)(row2+i,col2+j) > 0 && (*timeRtru)(row2+i,col2+j)> 0){
492 if((*timeRtru)(row2+i,col2+j) > ampmax2(3,mtru) )
493 ampmax2(3,mtru) = (*timeRtru)(row2+i,col2+j); // max time
499 //Sliding nxn, add nxn amplitudes (OVERLAP)
501 for(Int_t irow = 0 ; irow < nModulesPhi/2; irow++){
502 for(Int_t icol = 0 ; icol < nModulesEta/2; icol++){
504 if( (irow+fPatchSize) < nModulesPhi/2 && (icol+fPatchSize) < nModulesEta/2){ //Avoid exit the TRU
505 for(Int_t i = 0 ; i <= fPatchSize ; i++)
506 for(Int_t j = 0 ; j <= fPatchSize ; j++)
507 ampn += tru2x2(irow+i,icol+j);
508 //Select nxn maximum sums to select L1
509 if(ampn > ampmaxn(0,mtru)){
510 ampmaxn(0,mtru) = ampn ;
511 ampmaxn(1,mtru) = irow;
512 ampmaxn(2,mtru) = icol;
518 ampmaxn(3,mtru) = 0.; // Was 1 , I don't know why
520 //Find most recent time in selected nxn cell
521 Int_t rown = static_cast <Int_t> (ampmaxn(1,mtru));
522 Int_t coln = static_cast <Int_t> (ampmaxn(2,mtru));
523 for(Int_t i = 0; i<4*fPatchSize; i++){
524 for(Int_t j = 0; j<4*fPatchSize; j++){
525 if( (rown+i) < nModulesPhi && (coln+j) < nModulesEta){//Avoid exit the TRU
526 if((*amptru)(rown+i,coln+j) > 0 && (*timeRtru)(rown+i,coln+j)> 0){
527 if((*timeRtru)(rown+i,coln+j) > ampmaxn(3,mtru) )
528 ampmaxn(3,mtru) = (*timeRtru)(rown+i,coln+j); // max time
534 } else { // copy 2x2 to nxn
535 ampmaxn(0,mtru) = ampmax2(0,mtru);
536 ampmaxn(1,mtru) = ampmax2(1,mtru);
537 ampmaxn(2,mtru) = ampmax2(2,mtru);
538 ampmaxn(3,mtru) = ampmax2(3,mtru);
541 if(keyPrint) AliDebug(2,Form(" : MakeSlidingTowers -OUt \n"));
544 //____________________________________________________________________________
545 void AliEMCALTrigger::Print(const Option_t * opt) const
548 //Prints main parameters
552 AliTriggerInput* in = 0x0 ;
553 AliInfo(Form(" fSimulation %i (input option) : #digits %i\n", fSimulation, fDigitsList->GetEntries()));
554 AliInfo(Form(" fTimeKey %i \n ", fTimeKey));
556 AliInfo(Form("\t Maximum Amplitude after Sliding Cell, \n")) ;
557 AliInfo(Form("\t -2x2 cells sum (not overlapped): %10.2f, in Super Module %d\n",
558 f2x2MaxAmp,f2x2SM)) ;
559 AliInfo(Form("\t -2x2 from row %d to row %d and from column %d to column %d\n", f2x2ModulePhi, f2x2ModulePhi+2, f2x2ModuleEta, f2x2ModuleEta+2));
560 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)));
562 AliInfo(Form("\t Patch Size, n x n: %d x %d cells\n",2*(fPatchSize+1), 2*(fPatchSize+1)));
563 AliInfo(Form("\t -nxn cells sum (overlapped) : %10.2f, in Super Module %d\n", fnxnMaxAmp,fnxnSM));
564 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)) ;
565 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) ));
568 AliInfo(Form("\t Isolate in SuperModule? %d\n", fIsolateInSuperModule)) ;
569 AliInfo(Form("\t Threshold for LO %10.2f\n", fL0Threshold));
571 in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_L0" );
573 AliInfo(Form("\t *** EMCAL LO is set ***\n"));
575 AliInfo(Form("\t Gamma Low Pt Threshold for L1 %10.2f\n", fL1GammaLowPtThreshold));
576 in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_GammaLPt_L1" );
578 AliInfo(Form("\t *** EMCAL Gamma Low Pt for L1 is set ***\n"));
580 AliInfo(Form("\t Gamma Medium Pt Threshold for L1 %10.2f\n", fL1GammaMediumPtThreshold));
581 in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_GammaMPt_L1" );
583 AliInfo(Form("\t *** EMCAL Gamma Medium Pt for L1 is set ***\n"));
585 AliInfo(Form("\t Gamma High Pt Threshold for L1 %10.2f\n", fL1GammaHighPtThreshold));
586 in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_GammaHPt_L1" );
588 AliInfo(Form("\t *** EMCAL Gamma High Pt for L1 is set ***\n")) ;
592 //____________________________________________________________________________
593 void AliEMCALTrigger::SetTriggers(const TClonesArray * ampmatrix,const Int_t iSM,
594 const TMatrixD &max2,
595 const TMatrixD &maxn)
597 //Checks the 2x2 and nxn maximum amplitude per each TRU and
598 //compares with the different L0 and L1 triggers thresholds
599 Float_t max2[] = {-1,-1,-1,-1} ;
600 Float_t maxn[] = {-1,-1,-1,-1} ;
604 Int_t nTRU = fGeom->GetNTRU();
606 //Find maximum summed amplitude of all the TRU
608 for(Int_t i = 0 ; i < nTRU ; i++){
609 if(max2[0] < ampmax2(0,i) ){
610 max2[0] = ampmax2(0,i) ; // 2x2 summed max amplitude
611 max2[1] = ampmax2(1,i) ; // corresponding phi position in TRU
612 max2[2] = ampmax2(2,i) ; // corresponding eta position in TRU
613 max2[3] = ampmax2(3,i) ; // corresponding most recent time
616 if(maxn[0] < ampmaxn(0,i) ){
617 maxn[0] = ampmaxn(0,i) ; // nxn summed max amplitude
618 maxn[1] = ampmaxn(1,i) ; // corresponding phi position in TRU
619 maxn[2] = ampmaxn(2,i) ; // corresponding eta position in TRU
620 maxn[3] = ampmaxn(3,i) ; // corresponding most recent time
625 //--------Set max amplitude if larger than in other Super Modules------------
626 Float_t maxtimeR2 = -1 ;
627 Float_t maxtimeRn = -1 ;
628 static AliEMCALRawUtils rawUtil;
629 Int_t nTimeBins = rawUtil.GetRawFormatTimeBins() ;
631 //Set max of 2x2 amplitudes and select L0 trigger
632 if(max2[0] > f2x2MaxAmp ){
633 // if(max2[0] > 5) printf(" L0 : iSM %i: max2[0] %5.0f : max2[3] %5.0f (maxtimeR2) \n",
634 // iSM, max2[0], max2[3]);
635 f2x2MaxAmp = max2[0] ;
637 maxtimeR2 = max2[3] ;
638 fGeom->GetModulePhiEtaIndexInSModuleFromTRUIndex(mtru2,
639 static_cast<Int_t>(max2[1]),
640 static_cast<Int_t>(max2[2]),
641 f2x2ModulePhi,f2x2ModuleEta);
644 if(fIsolateInSuperModule)
645 fIs2x2Isol = IsPatchIsolated(0, ampmatrix, iSM, mtru2, f2x2MaxAmp, f2x2ModulePhi,f2x2ModuleEta) ;
647 fIs2x2Isol = IsPatchIsolated(0, ampmatrix, iSM, mtru2, f2x2MaxAmp, static_cast<Int_t>(max2[1]), static_cast<Int_t>(max2[2])) ;
650 //Transform digit amplitude in Raw Samples
651 if (fADCValuesLow2x2 == 0) {
652 fADCValuesLow2x2 = new Int_t[nTimeBins];
653 fADCValuesHigh2x2 = new Int_t[nTimeBins];
655 //printf(" maxtimeR2 %12.5e (1)\n", maxtimeR2);
656 rawUtil.RawSampledResponse(maxtimeR2 * AliEMCALRawUtils::GetRawFormatTimeBin(),
657 f2x2MaxAmp, fADCValuesHigh2x2, fADCValuesLow2x2) ;
659 // Set Trigger Inputs, compare ADC time bins until threshold is attained
661 for(Int_t i = 0 ; i < nTimeBins ; i++){
662 // printf(" fADCValuesHigh2x2[%i] %i : %i \n", i, fADCValuesHigh2x2[i], fADCValuesLow2x2[i]);
663 if(fADCValuesHigh2x2[i] >= fL0Threshold || fADCValuesLow2x2[i] >= fL0Threshold){
664 SetInput("EMCAL_L0") ;
669 // Nov 5 - no analysis of time information
670 if(f2x2MaxAmp >= fL0Threshold) { // should add the low amp too
671 SetInput("EMCAL_L0");
676 //------------Set max of nxn amplitudes and select L1 trigger---------
677 if(maxn[0] > fnxnMaxAmp ){
678 fnxnMaxAmp = maxn[0] ;
680 maxtimeRn = maxn[3] ;
681 fGeom->GetModulePhiEtaIndexInSModuleFromTRUIndex(mtrun,
682 static_cast<Int_t>(maxn[1]),
683 static_cast<Int_t>(maxn[2]),
684 fnxnModulePhi,fnxnModuleEta) ;
687 if(fIsolateInSuperModule)
688 fIsnxnIsol = IsPatchIsolated(1, ampmatrix, iSM, mtrun, fnxnMaxAmp, fnxnModulePhi, fnxnModuleEta) ;
690 fIsnxnIsol = IsPatchIsolated(1, ampmatrix, iSM, mtrun, fnxnMaxAmp, static_cast<Int_t>(maxn[1]), static_cast<Int_t>(maxn[2])) ;
693 //Transform digit amplitude in Raw Samples
694 if (fADCValuesLownxn == 0) {
695 fADCValuesHighnxn = new Int_t[nTimeBins];
696 fADCValuesLownxn = new Int_t[nTimeBins];
698 rawUtil.RawSampledResponse(maxtimeRn * AliEMCALRawUtils::GetRawFormatTimeBin(),
699 fnxnMaxAmp, fADCValuesHighnxn, fADCValuesLownxn) ;
701 //Set Trigger Inputs, compare ADC time bins until threshold is attained
703 for(Int_t i = 0 ; i < nTimeBins ; i++){
704 if(fADCValuesHighnxn[i] >= fL1GammaLowPtThreshold || fADCValuesLownxn[i] >= fL1GammaLowPtThreshold){
705 SetInput("EMCAL_GammaLPt_L1") ;
711 for(Int_t i = 0 ; i < nTimeBins ; i++){
712 if(fADCValuesHighnxn[i] >= fL1GammaMediumPtThreshold || fADCValuesLownxn[i] >= fL1GammaMediumPtThreshold){
713 SetInput("EMCAL_GammaMPt_L1") ;
719 for(Int_t i = 0 ; i < nTimeBins ; i++){
720 if(fADCValuesHighnxn[i] >= fL1GammaHighPtThreshold || fADCValuesLownxn[i] >= fL1GammaHighPtThreshold){
721 SetInput("EMCAL_GammaHPt_L1") ;
726 // Nov 5 - no analysis of time information
727 if(fnxnMaxAmp >= fL1GammaLowPtThreshold) { // should add the low amp too
728 SetInput("EMCAL_GammaLPt_L1") ; //SetL1 Low
730 if(fnxnMaxAmp >= fL1GammaMediumPtThreshold) { // should add the low amp too
731 SetInput("EMCAL_GammaMPt_L1") ; //SetL1 Medium
733 if(fnxnMaxAmp >= fL1GammaHighPtThreshold) { // should add the low amp too
734 SetInput("EMCAL_GammaHPt_L1") ; //SetL1 High
740 //____________________________________________________________________________
741 void AliEMCALTrigger::FillTRU(const TClonesArray * digits, TClonesArray * ampmatrix, TClonesArray * ampmatrixsmod, TClonesArray * timeRmatrix) {
743 // Orders digits ampitudes list in fNTRU TRUs (384 cells) per supermodule.
744 // Each TRU is a TMatrixD, and they are kept in TClonesArrays. The number of
745 // TRU in phi is fNTRUPhi, and the number of TRU in eta is fNTRUEta.
746 // Last 2 modules are half size in Phi, I considered that the number of TRU
747 // is maintained for the last modules but decision not taken. If different,
748 // then this must be changed. Also fill a matrix with all amplitudes in supermodule for isolation studies.
750 // Initilize and declare variables
751 // List of TRU matrices initialized to 0.
752 // printf("<I> AliEMCALTrigger::FillTRU() started : # digits %i\n", digits->GetEntriesFast());
755 // One input per EMCAL module so size of matrix is reduced by 4 (2x2 division case)
757 Int_t nPhi = fGeom->GetNPhi();
758 Int_t nZ = fGeom->GetNZ();
759 Int_t nTRU = fGeom->GetNTRU();
760 // Int_t nTRUPhi = fGeom->GetNTRUPhi();
761 Int_t nModulesPhi = fGeom->GetNModulesInTRUPhi();
762 Int_t nModulesPhi2 = fGeom->GetNModulesInTRUPhi();
763 Int_t nModulesEta = fGeom->GetNModulesInTRUEta();
764 // printf("<I> AliEMCALTrigger::FillTRU() nTRU %i nTRUPhi %i : nModulesPhi %i nModulesEta %i \n",
765 // nTRU, nTRUPhi, nModulesPhi, nModulesEta);
776 // iphim, ietam - module indexes in SM
780 //List of TRU matrices initialized to 0.
781 Int_t nSup = fGeom->GetNumberOfSuperModules();
782 for(Int_t k = 0; k < nTRU*nSup; k++){
783 TMatrixD amptrus(nModulesPhi,nModulesEta) ;
784 TMatrixD timeRtrus(nModulesPhi,nModulesEta) ;
785 // Do we need to initialise? I think TMatrixD does it by itself...
786 for(Int_t i = 0; i < nModulesPhi; i++){
787 for(Int_t j = 0; j < nModulesEta; j++){
789 timeRtrus(i,j) = 0.0;
792 new((*ampmatrix)[k]) TMatrixD(amptrus) ;
793 new((*timeRmatrix)[k]) TMatrixD(timeRtrus) ;
796 // List of Modules matrices initialized to 0.
797 for(Int_t k = 0; k < nSup ; k++){
799 // if(nSup>9) mphi = nPhi/2; // the same size
800 TMatrixD ampsmods( mphi, nZ);
801 for(Int_t i = 0; i < mphi; i++){
802 for(Int_t j = 0; j < nZ; j++){
806 new((*ampmatrixsmod)[k]) TMatrixD(ampsmods) ;
809 AliEMCALDigit * dig ;
811 //Digits loop to fill TRU matrices with amplitudes.
812 for(Int_t idig = 0 ; idig < digits->GetEntriesFast() ; idig++){
814 dig = dynamic_cast<AliEMCALDigit *>(digits->At(idig)) ;
815 amp = Float_t(dig->GetAmp()); // Energy of the digit (arbitrary units)
816 id = dig->GetId() ; // Id label of the cell
817 timeR = dig->GetTimeR() ; // Earliest time of the digit
818 if(amp<=0.0) AliDebug(1,Form(" id %i amp %f \n", id, amp));
819 // printf(" FILLTRU : timeR %10.5e time %10.5e : amp %10.5e \n", timeR, dig->GetTime(), amp);
820 // Get eta and phi cell position in supermodule
821 Bool_t bCell = fGeom->GetCellIndex(id, iSupMod, nModule, nIphi, nIeta) ;
823 AliError(Form("FillTRU","%i Wrong cell id number %i ", idig, id)) ;
825 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
826 // iphim, ietam - module indexes in SM
827 fGeom->GetModuleIndexesFromCellIndexesInSModule(iSupMod,iphi,ieta, iphim, ietam, nModule);
829 //printf("iSupMod %i nModule %i iphi %i ieta %i iphim %i ietam %i \n",
830 //iSupMod,nModule, iphi, ieta, iphim, ietam);
832 // Check to which TRU in the supermodule belongs the cell.
833 // Supermodules are divided in a TRU matrix of dimension
834 // (fNTRUPhi,fNTRUEta).
835 // Each TRU is a cell matrix of dimension (nModulesPhi,nModulesEta)
837 // First calculate the row and column in the supermodule
838 // of the TRU to which the cell belongs.
839 Int_t row = iphim / nModulesPhi;
840 Int_t col = ietam / nModulesEta;
841 //Calculate label number of the TRU
842 Int_t itru = fGeom->GetAbsTRUNumberFromNumberInSm(row, col, iSupMod);
844 //Fill TRU matrix with cell values
845 TMatrixD * amptrus = dynamic_cast<TMatrixD *>(ampmatrix->At(itru)) ;
846 TMatrixD * timeRtrus = dynamic_cast<TMatrixD *>(timeRmatrix->At(itru)) ;
848 //Calculate row and column of the module inside the TRU with number itru
849 Int_t irow = iphim - row * nModulesPhi;
851 irow = iphim - row * nModulesPhi2; // size of matrix the same
852 Int_t icol = ietam - col * nModulesEta;
854 (*amptrus)(irow,icol) += amp ;
855 if((*timeRtrus)(irow,icol) <0.0 || (*timeRtrus)(irow,icol) <= timeR){ // ??
856 (*timeRtrus)(irow,icol) = timeR ;
858 //printf(" ieta %i iphi %i iSM %i || col %i row %i : itru %i -> amp %f\n",
859 // ieta, iphi, iSupMod, col, row, itru, amp);
860 //####################SUPERMODULE MATRIX ##################
861 TMatrixD * ampsmods = dynamic_cast<TMatrixD *>(ampmatrixsmod->At(iSupMod)) ;
862 (*ampsmods)(iphim,ietam) += amp ;
863 // printf(" id %i iphim %i ietam %i SM %i : irow %i icol %i itru %i : amp %6.0f\n",
864 //id, iphim, ietam, iSupMod, irow, icol, itru, amp);
867 //printf("<I> AliEMCALTrigger::FillTRU() is ended \n");
870 //____________________________________________________________________________
871 void AliEMCALTrigger::Trigger()
873 //Main Method to select triggers.
874 TH1::AddDirectory(0);
876 AliRunLoader *runLoader = AliRunLoader::Instance();
877 AliEMCALLoader *emcalLoader = 0;
879 emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
882 //Load EMCAL Geometry
883 if (runLoader && runLoader->GetAliRun() && runLoader->GetAliRun()->GetDetector("EMCAL"))
884 fGeom = dynamic_cast<AliEMCAL*>(runLoader->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
887 fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
890 AliFatal("Did not get geometry from EMCALLoader");
893 Int_t nSuperModules = fGeom->GetNumberOfSuperModules() ; //12 SM in EMCAL
894 Int_t nTRU = fGeom->GetNTRU(); // 3 TRU per super module
896 //Intialize data members each time the trigger is called in event loop
897 f2x2MaxAmp = -1; f2x2ModulePhi = -1; f2x2ModuleEta = -1;
898 fnxnMaxAmp = -1; fnxnModulePhi = -1; fnxnModuleEta = -1;
900 // Take the digits list if simulation
901 if(fSimulation && runLoader && emcalLoader){ // works than run seperate macros
902 runLoader->LoadDigits("EMCAL");
903 fDigitsList = emcalLoader->Digits() ;
904 runLoader->LoadSDigits("EMCAL");
906 // Digits list should be set by method SetDigitsList(TClonesArray * digits)
908 AliFatal("Digits not found !") ;
910 //Take the digits list
912 // Delete old if unzero
913 if(fAmpTrus) {fAmpTrus->Delete(); delete fAmpTrus;}
914 if(fTimeRtrus) {fTimeRtrus->Delete(); delete fTimeRtrus;}
915 if(fAmpSMods) {fAmpSMods->Delete(); delete fAmpSMods;}
916 // Fill TRU and SM matrix
917 fAmpTrus = new TClonesArray("TMatrixD",nTRU);
918 fAmpTrus->SetName("AmpTrus");
919 fTimeRtrus = new TClonesArray("TMatrixD",nTRU);
920 fTimeRtrus->SetName("TimeRtrus");
921 fAmpSMods = new TClonesArray("TMatrixD",nSuperModules);
922 fAmpSMods->SetName("AmpSMods");
924 FillTRU(fDigitsList, fAmpTrus, fAmpSMods, fTimeRtrus);
926 // Jet stuff - only one case, no freedom here
927 if(fGeom->GetNEtaSubOfTRU() == 6) {
928 if(fAmpJetMatrix) {delete fAmpJetMatrix; fAmpJetMatrix=0;}
929 if(fJetMatrixE) {delete fJetMatrixE; fJetMatrixE=0;}
931 fAmpJetMatrix = new TMatrixD(17,12); // 17-phi(row), 12-eta(col)
932 fJetMatrixE = new TH2F("fJetMatrixE"," E of max patch in (#phi,#eta)",
933 17, 80.*TMath::DegToRad(), (180.+20.*2/3.)*TMath::DegToRad(), 12, -0.7, 0.7);
934 for(Int_t row=0; row<fAmpJetMatrix->GetNrows(); row++) {
935 for(Int_t col=0; col<fAmpJetMatrix->GetNcols(); col++) {
936 (*fAmpJetMatrix)(row,col) = 0.;
939 FillJetMatrixFromSMs(fAmpSMods, fAmpJetMatrix, fGeom);
941 if(!CheckConsistentOfMatrixes()) assert(0);
943 // Do Tower Sliding and select Trigger
944 // Initialize varible that will contain maximum amplitudes and
945 // its corresponding tower position in eta and phi, and time.
946 TMatrixD ampmax2(4,nTRU) ; // 0-max amp, 1-irow, 2-icol, 3-timeR
947 TMatrixD ampmaxn(4,nTRU) ;
949 for(Int_t iSM = 0 ; iSM < nSuperModules ; iSM++) {
950 //Do 2x2 and nxn sums, select maximums.
952 MakeSlidingTowers(fAmpTrus, fTimeRtrus, iSM, ampmax2, ampmaxn);
955 if(fIsolateInSuperModule) // here some discripency between tru and SM
956 SetTriggers(fAmpSMods,iSM,ampmax2,ampmaxn) ;
957 if(!fIsolateInSuperModule)
958 SetTriggers(fAmpTrus,iSM,ampmax2,ampmaxn) ;
961 // Do patch sliding and select Jet Trigger
962 // 0-max amp-meanFromVZERO(if), 1-irow, 2-icol, 3-timeR,
963 // 4-max amp , 5-meanFromVZERO (Nov 25, 2007)
965 MakeSlidingPatch((*fAmpJetMatrix), fNJetPatchPhi, fAmpJetMax); // no timing information here
971 //____________________________________________________________________________
972 void AliEMCALTrigger::GetTriggerInfo(TArrayF &triggerPosition, TArrayF &triggerAmplitudes) const
974 // Template - should be defined; Nov 5, 2007
975 triggerPosition[0] = 0.;
976 triggerAmplitudes[0] = 0.;
979 //____________________________________________________________________________
980 void AliEMCALTrigger::FillJetMatrixFromSMs(TClonesArray *ampmatrixsmod, TMatrixD* jetMat, AliEMCALGeometry *g)
983 // Fill matrix for jet trigger from SM matrixes of modules
985 static int keyPrint = 0;
987 if(ampmatrixsmod==0 || jetMat==0 || g==0) return;
988 Double_t amp = 0.0, ampSum=0.0;
990 Int_t nEtaModSum = g->GetNZ() / g->GetNEtaSubOfTRU(); // should be 4
991 Int_t nPhiModSum = g->GetNPhi() / g->GetNTRUPhi(); // should be 4
993 if(keyPrint) AliDebug(2,Form("%s",Form(" AliEMCALTrigger::FillJetMatrixFromSMs | nEtaModSum %i : nPhiModSum %i \n", nEtaModSum, nPhiModSum)));
994 Int_t jrow=0, jcol=0; // indexes of jet matrix
995 Int_t nEtaSM=0, nPhiSM=0;
996 for(Int_t iSM=0; iSM<ampmatrixsmod->GetEntries(); iSM++) {
997 TMatrixD * ampsmods = dynamic_cast<TMatrixD *>(ampmatrixsmod->At(iSM));
998 Int_t nrow = ampsmods->GetNrows();
999 Int_t ncol = ampsmods->GetNcols();
1000 //printf("%s",Form(" ######## SM %i : nrow %i : ncol %i ##### \n", iSM, nrow, ncol));
1001 for(Int_t row=0; row<nrow; row++) {
1002 for(Int_t col=0; col<ncol; col++) {
1003 amp = (*ampsmods)(row,col);
1007 if(keyPrint) AliDebug(2,Form("%s",Form(" ** nPhiSm %i : nEtaSM %i : row %2.2i : col %2.2i -> ", nPhiSM, nEtaSM, row, col)));
1008 if(nEtaSM == 0) { // positive Z
1009 jrow = 3*nPhiSM + row/nPhiModSum;
1010 jcol = 6 + col / nEtaModSum;
1011 } else { // negative Z
1012 if(iSM<=9) jrow = 3*nPhiSM + 2 - row/nPhiModSum;
1013 else jrow = 3*nPhiSM + 1 - row/nPhiModSum; // half size
1014 jcol = 5 - col / nEtaModSum;
1016 if(keyPrint) AliDebug(2,Form("%s",Form(" jrow %2.2i : jcol %2.2i : amp %f (jetMat) \n", jrow, jcol, amp)));
1018 (*jetMat)(jrow,jcol) += amp;
1019 ampSum += amp; // For controling
1020 } else if(amp<0.0) {
1021 AliDebug(1,Form(" jrow %2.2i : jcol %2.2i : amp %f (jetMat: amp<0) \n", jrow, jcol, amp));
1027 if(ampSum <= 0.0) AliDebug(1,Form("FillJetMatrixFromSMs","ampSum %f (<=0.0) ", ampSum));
1030 //____________________________________________________________________________
1031 void AliEMCALTrigger::MakeSlidingPatch(const TMatrixD &jm, const Int_t nPatchSize, TMatrixD &JetMax)
1033 // Sliding patch : nPatchSize x nPatchSize (OVERLAP)
1034 static int keyPrint = 0;
1035 if(keyPrint) AliDebug(2,Form(" AliEMCALTrigger::MakeSlidingPatch() was started \n"));
1036 Double_t ampCur = 0.0, e=0.0;
1037 ampJetMax(0,0) = 0.0;
1038 ampJetMax(3,0) = 0.0; // unused now
1039 ampJetMax(4,0) = ampJetMax(5,0) = 0.0;
1040 for(Int_t row=0; row<fAmpJetMatrix->GetNrows(); row ++) {
1041 for(Int_t col=0; col<fAmpJetMatrix->GetNcols(); col++) {
1043 // check on patch size
1044 if( (row+nPatchSize-1) < fAmpJetMatrix->GetNrows() && (col+nPatchSize-1) < fAmpJetMatrix->GetNcols()){
1045 for(Int_t i = 0 ; i < nPatchSize ; i++) {
1046 for(Int_t j = 0 ; j < nPatchSize ; j++) {
1047 ampCur += jm(row+i, col+j);
1049 } // end cycle on patch
1050 if(ampCur > ampJetMax(0,0)){
1051 ampJetMax(0,0) = ampCur;
1052 ampJetMax(1,0) = row;
1053 ampJetMax(2,0) = col;
1055 } // check on patch size
1058 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));
1060 Double_t eCorrJetMatrix=0.0;
1061 if(fVZER0Mult > 0.0) {
1062 // Correct patch energy (adc) and jet patch matrix energy
1063 Double_t meanAmpBG = GetMeanEmcalPatchEnergy(Int_t(fVZER0Mult), nPatchSize)/0.0153;
1064 ampJetMax(4,0) = ampJetMax(0,0);
1065 ampJetMax(5,0) = meanAmpBG;
1067 Double_t eCorr = ampJetMax(0,0) - meanAmpBG;
1068 AliDebug(2,Form(" ampJetMax(0,0) %f meanAmpBG %f eCorr %f : ampJetMax(4,0) %f \n",
1069 ampJetMax(0,0), meanAmpBG, eCorr, ampJetMax(5,0)));
1070 ampJetMax(0,0) = eCorr;
1072 eCorrJetMatrix = GetMeanEmcalEnergy(Int_t(fVZER0Mult)) / 208.;
1074 // Fill patch energy matrix
1075 for(int row=Int_t(ampJetMax(1,0)); row<Int_t(ampJetMax(1,0))+nPatchSize; row++) {
1076 for(int col=Int_t(ampJetMax(2,0)); col<Int_t(ampJetMax(2,0))+nPatchSize; col++) {
1077 e = Double_t(jm(row,col)*0.0153); // 0.0153 - hard coded now
1078 if(eCorrJetMatrix > 0.0) { // BG subtraction case
1079 e -= eCorrJetMatrix;
1080 fJetMatrixE->SetBinContent(row+1, col+1, e);
1081 } else if(e > 0.0) {
1082 fJetMatrixE->SetBinContent(row+1, col+1, e);
1086 // PrintJetMatrix();
1087 // Set the jet trigger(s), multiple threshold now, Nov 19,2007
1088 for(Int_t i=0; i<fNJetThreshold; i++ ) {
1089 if(ampJetMax(0,0) >= fL1JetThreshold[i]) {
1090 SetInput(GetNameOfJetTrigger(i));
1095 //____________________________________________________________________________
1096 Double_t AliEMCALTrigger::GetEmcalSumAmp() const
1098 // Return sum of amplidutes from EMCal
1099 // Used calibration coefficeint for transition to energy
1100 return fAmpJetMatrix >0 ?fAmpJetMatrix->Sum() :0.0;
1103 //____________________________________________________________________________
1104 void AliEMCALTrigger::PrintJetMatrix() const
1106 // fAmpJetMatrix : (17,12); // 17-phi(row), 12-eta(col)
1107 if(fAmpJetMatrix == 0) return;
1109 AliInfo(Form("\n #### jetMatrix : (%i,%i) ##### \n ",
1110 fAmpJetMatrix->GetNrows(), fAmpJetMatrix->GetNcols()));
1111 PrintMatrix(*fAmpJetMatrix);
1114 //____________________________________________________________________________
1115 void AliEMCALTrigger::PrintAmpTruMatrix(Int_t ind) const
1117 // Print matrix with TRU patches
1118 TMatrixD * tru = dynamic_cast<TMatrixD *>(fAmpTrus->At(ind));
1119 if(tru == 0) return;
1120 AliInfo(Form("\n #### Amp TRU matrix(%i) : (%i,%i) ##### \n ",
1121 ind, tru->GetNrows(), tru->GetNcols()));
1125 //____________________________________________________________________________
1126 void AliEMCALTrigger::PrintAmpSmMatrix(Int_t ind) const
1128 // Print matrix with SM amplitudes
1129 TMatrixD * sm = dynamic_cast<TMatrixD *>(fAmpSMods->At(ind));
1131 AliInfo(Form("\n #### Amp SM matrix(%i) : (%i,%i) ##### \n ",
1132 ind, sm->GetNrows(), sm->GetNcols()));
1136 //____________________________________________________________________________
1137 void AliEMCALTrigger::PrintMatrix(const TMatrixD &mat) const
1139 //Print matrix object
1140 for(Int_t col=0; col<mat.GetNcols(); col++) AliInfo(Form(" %3i ", col));
1141 AliInfo(Form("\n -- \n"));
1142 for(Int_t row=0; row<mat.GetNrows(); row++) {
1143 AliInfo(Form(" row:%2i ", row));
1144 for(Int_t col=0; col<mat.GetNcols(); col++) {
1145 AliInfo(Form(" %4i", (Int_t)mat(row,col)));
1151 //____________________________________________________________________________
1152 Bool_t AliEMCALTrigger::CheckConsistentOfMatrixes(const Int_t pri)
1154 // Check consitency of matrices
1155 Double_t sumSM = 0.0, smCur=0.0;
1156 Double_t sumTru=0.0, sumTruInSM = 0.0, truSum=0.0;
1157 // Bool_t key = kTRUE;
1159 for(Int_t i=0; i<fAmpSMods->GetEntries(); i++) {
1160 TMatrixD * sm = dynamic_cast<TMatrixD *>(fAmpSMods->At(i));
1166 for(Int_t itru=0; itru<3; itru++) { // Cycle on tru inside SM
1167 Int_t ind = 3*i + itru;
1168 TMatrixD *tru = dynamic_cast<TMatrixD *>(fAmpTrus->At(ind));
1170 truSum = tru->Sum();
1171 sumTruInSM += truSum;
1174 sumTru += sumTruInSM;
1176 if(sumTruInSM != smCur) {
1177 AliDebug(1,Form(" sm %i : smCur %f -> sumTruInSM %f \n", i, smCur, sumTruInSM));
1182 Double_t sumJetMat = fAmpJetMatrix->Sum();
1183 if(pri || TMath::Abs(sumSM-sumTru)>0.0001 || TMath::Abs(sumSM-sumJetMat) > 0.0001)
1184 AliDebug(1,Form(" sumSM %f : sumTru %f : sumJetMat %f \n", sumSM, sumTru, sumJetMat));
1185 if(TMath::Abs(sumSM - sumTru)>0.0001 || TMath::Abs(sumSM-sumJetMat) > 0.0001) return kFALSE;
1189 //____________________________________________________________________________
1190 void AliEMCALTrigger::Browse(TBrowser* b)
1193 if(&fInputs) b->Add(&fInputs);
1194 if(fAmpTrus) b->Add(fAmpTrus);
1195 if(fTimeRtrus) b->Add(fTimeRtrus);
1196 if(fAmpSMods) b->Add(fAmpSMods);
1197 if(fAmpJetMatrix) b->Add(fAmpJetMatrix);
1198 if(fJetMatrixE) b->Add(fJetMatrixE);