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 ClassImp(AliEMCALTrigger)
65 TString AliEMCALTrigger::fgNameOfJetTriggers("EMCALJetTriggerL1");
67 //______________________________________________________________________
68 AliEMCALTrigger::AliEMCALTrigger()
69 : AliTriggerDetector(), fGeom(0),
70 f2x2MaxAmp(-1), f2x2ModulePhi(-1), f2x2ModuleEta(-1),
72 fnxnMaxAmp(-1), fnxnModulePhi(-1), fnxnModuleEta(-1),
74 fADCValuesHighnxn(0),fADCValuesLownxn(0),
75 fADCValuesHigh2x2(0),fADCValuesLow2x2(0),
77 fL0Threshold(100),fL1GammaLowPtThreshold(200),
78 fL1GammaMediumPtThreshold(500), fL1GammaHighPtThreshold(1000),
79 fPatchSize(1), fIsolPatchSize(1),
80 f2x2AmpOutOfPatch(-1), fnxnAmpOutOfPatch(-1),
81 f2x2AmpOutOfPatchThres(100000), fnxnAmpOutOfPatchThres(100000),
82 fIs2x2Isol(kFALSE), fIsnxnIsol(kFALSE),
83 fSimulation(kTRUE), fIsolateInSuperModule(kTRUE), fTimeKey(kFALSE),
84 fAmpTrus(0),fTimeRtrus(0),fAmpSMods(0),
85 fTriggerPosition(6), fTriggerAmplitudes(4),
86 fNJetPatchPhi(3), fNJetPatchEta(3), fNJetThreshold(3), fL1JetThreshold(0), fJetMaxAmp(0),
87 fAmpJetMatrix(0), fJetMatrixE(0), fAmpJetMax(6,1), fVZER0Mult(0.)
90 fADCValuesHighnxn = 0x0; //new Int_t[fTimeBins];
91 fADCValuesLownxn = 0x0; //new Int_t[fTimeBins];
92 fADCValuesHigh2x2 = 0x0; //new Int_t[fTimeBins];
93 fADCValuesLow2x2 = 0x0; //new Int_t[fTimeBins];
96 // Define jet threshold - can not change from outside now
97 // fNJetThreshold = 7; // For MB Pythia suppression
98 fNJetThreshold = 10; // Hijing
99 fL1JetThreshold = new Double_t[fNJetThreshold];
100 if(fNJetThreshold == 7) {
101 fL1JetThreshold[0] = 5./0.0153;
102 fL1JetThreshold[1] = 8./0.0153;
103 fL1JetThreshold[2] = 10./0.0153;
104 fL1JetThreshold[3] = 12./0.0153;
105 fL1JetThreshold[4] = 13./0.0153;
106 fL1JetThreshold[5] = 14./0.0153;
107 fL1JetThreshold[6] = 15./0.0153;
108 } else if(fNJetThreshold == 10) {
109 Double_t thGev[10]={5.,8.,10., 12., 13.,14.,15., 17., 20., 25.};
110 for(Int_t i=0; i<fNJetThreshold; i++) fL1JetThreshold[i] = thGev[i]/0.0153;
112 fL1JetThreshold[0] = 5./0.0153;
113 fL1JetThreshold[1] = 10./0.0153;
114 fL1JetThreshold[2] = 15./0.0153;
115 fL1JetThreshold[3] = 20./0.0153;
116 fL1JetThreshold[4] = 25./0.0153;
121 fInputs.SetName("TriggersInputs");
127 //____________________________________________________________________________
128 AliEMCALTrigger::AliEMCALTrigger(const AliEMCALTrigger & trig)
129 : AliTriggerDetector(trig),
131 f2x2MaxAmp(trig.f2x2MaxAmp),
132 f2x2ModulePhi(trig.f2x2ModulePhi),
133 f2x2ModuleEta(trig.f2x2ModuleEta),
135 fnxnMaxAmp(trig.fnxnMaxAmp),
136 fnxnModulePhi(trig.fnxnModulePhi),
137 fnxnModuleEta(trig.fnxnModuleEta),
139 fADCValuesHighnxn(trig.fADCValuesHighnxn),
140 fADCValuesLownxn(trig.fADCValuesLownxn),
141 fADCValuesHigh2x2(trig.fADCValuesHigh2x2),
142 fADCValuesLow2x2(trig.fADCValuesLow2x2),
143 fDigitsList(trig.fDigitsList),
144 fL0Threshold(trig.fL0Threshold),
145 fL1GammaLowPtThreshold(trig.fL1GammaLowPtThreshold),
146 fL1GammaMediumPtThreshold(trig.fL1GammaMediumPtThreshold),
147 fL1GammaHighPtThreshold(trig.fL1GammaHighPtThreshold),
148 fPatchSize(trig.fPatchSize),
149 fIsolPatchSize(trig.fIsolPatchSize),
150 f2x2AmpOutOfPatch(trig.f2x2AmpOutOfPatch),
151 fnxnAmpOutOfPatch(trig.fnxnAmpOutOfPatch),
152 f2x2AmpOutOfPatchThres(trig.f2x2AmpOutOfPatchThres),
153 fnxnAmpOutOfPatchThres(trig.fnxnAmpOutOfPatchThres),
154 fIs2x2Isol(trig.fIs2x2Isol),
155 fIsnxnIsol(trig.fIsnxnIsol),
156 fSimulation(trig.fSimulation),
157 fIsolateInSuperModule(trig.fIsolateInSuperModule),
158 fTimeKey(trig.fTimeKey),
159 fAmpTrus(trig.fAmpTrus),
160 fTimeRtrus(trig.fTimeRtrus),
161 fAmpSMods(trig.fAmpSMods),
162 fTriggerPosition(trig.fTriggerPosition),
163 fTriggerAmplitudes(trig.fTriggerAmplitudes),
164 fNJetPatchPhi(trig.fNJetPatchPhi),
165 fNJetPatchEta(trig.fNJetPatchEta),
166 fNJetThreshold(trig.fNJetThreshold),
167 fL1JetThreshold(trig.fL1JetThreshold),
168 fJetMaxAmp(trig.fJetMaxAmp),
169 fAmpJetMatrix(trig.fAmpJetMatrix),
170 fJetMatrixE(trig.fJetMatrixE),
171 fAmpJetMax(trig.fAmpJetMax),
172 fVZER0Mult(trig.fVZER0Mult)
177 AliEMCALTrigger::~AliEMCALTrigger() {
179 delete [] fADCValuesHighnxn;
180 delete [] fADCValuesLownxn;
181 delete [] fADCValuesHigh2x2;
182 delete [] fADCValuesLow2x2;
184 if(fAmpTrus) {fAmpTrus->Delete(); delete fAmpTrus;}
185 if(fTimeRtrus) {fTimeRtrus->Delete(); delete fTimeRtrus;}
186 if(fAmpSMods) {fAmpSMods->Delete(); delete fAmpSMods;}
187 if(fAmpJetMatrix) delete fAmpJetMatrix;
188 if(fJetMatrixE) delete fJetMatrixE;
189 if(fL1JetThreshold) delete [] fL1JetThreshold;
192 //----------------------------------------------------------------------
193 void AliEMCALTrigger::CreateInputs()
197 // Do not create inputs again!!
198 if( fInputs.GetEntriesFast() > 0 ) return;
200 // Second parameter should be detector name = "EMCAL"
201 TString det("EMCAL"); // Apr 29, 2008
202 fInputs.AddLast( new AliTriggerInput( det+"_L0", det, 0x02) );
203 fInputs.AddLast( new AliTriggerInput( det+"_GammaHPt_L1", det, 0x04 ) );
204 fInputs.AddLast( new AliTriggerInput( det+"_GammaMPt_L1", det, 0x08 ) );
205 fInputs.AddLast( new AliTriggerInput( det+"_GammaLPt_L1", det, 0x016 ) );
207 if(fNJetThreshold<=0) return;
209 UInt_t level = 0x032;
210 for(Int_t i=0; i<fNJetThreshold; i++ ) {
211 TString name(Form("%s_Th_%2.2i",fgNameOfJetTriggers.Data(),i));
212 TString title("EMCAL Jet triger L1 :"); // unused now
213 // 0.0153 - hard coded now
214 title += Form("Th %i(%5.1f GeV) :", (Int_t)fL1JetThreshold[i], fL1JetThreshold[i] * 0.0153);
215 title += Form("patch %ix%i~(%3.2f(#phi)x%3.2f(#eta)) ",
216 fNJetPatchPhi, fNJetPatchEta, 0.11*(fNJetPatchPhi), 0.11*(fNJetPatchEta));
217 fInputs.AddLast( new AliTriggerInput(name, det, UChar_t(level)) );
223 //____________________________________________________________________________
224 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) {
227 // EMCAL RTU size is 4modules(phi) x 24modules (eta)
228 // So maximum size of patch is 4modules x 4modules (EMCAL L0 trigger).
229 // Calculate if the maximum patch found is isolated, find amplitude around maximum (2x2 or nxn) patch,
230 // inside isolation patch . iPatchType = 0 means calculation for 2x2 patch,
231 // iPatchType = 1 means calculation for nxn patch.
232 // In the next table there is an example of the different options of patch size and isolation patch size:
233 // Patch Size (fPatchSize)
235 // fIsolPatchSize 0 2x2 (not overlap) 4x4 (overlapped)
240 // Get matrix of TRU or Module with maximum amplitude patch.
241 Int_t itru = mtru + iSM * fGeom->GetNTRU(); //number of tru, min 0 max 3*12=36.
242 TMatrixD * ampmatrix = 0x0;
245 static int keyPrint = 0;
246 if(keyPrint) printf(" IsPatchIsolated : iSM %i mtru %i itru %i maxphi %i maxeta %i \n", iSM, mtru, itru, maxphi, maxeta);
248 if(fIsolateInSuperModule){ // ?
249 ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(iSM)) ;
250 rowborder = fGeom->GetNPhi();
251 colborder = fGeom->GetNZ();
252 AliDebug(2,"Isolate trigger in Module");
254 ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(itru)) ;
255 rowborder = fGeom->GetNModulesInTRUPhi();
256 colborder = fGeom->GetNModulesInTRUEta();
257 AliDebug(2,"Isolate trigger in TRU");
259 if(iSM>9) rowborder /= 2; // half size in phi
261 //Define patch modules - what is this ??
262 Int_t isolmodules = fIsolPatchSize*(1+iPatchType);
263 Int_t ipatchmodules = 2*(1+fPatchSize*iPatchType);
264 Int_t minrow = maxphi - isolmodules;
265 Int_t mincol = maxeta - isolmodules;
266 Int_t maxrow = maxphi + isolmodules + ipatchmodules;
267 Int_t maxcol = maxeta + isolmodules + ipatchmodules;
269 minrow = minrow<0?0 :minrow;
270 mincol = mincol<0?0 :mincol;
272 maxrow = maxrow>rowborder?rowborder :maxrow;
273 maxcol = maxcol>colborder?colborder :maxcol;
275 //printf("%s\n",Form("Number of added Isol Modules %d, Patch Size %d",isolmodules, ipatchmodules));
276 //printf("%s\n",Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
277 // AliDebug(2,Form("Number of added Isol Modules %d, Patch Size %d",isolmodules, ipatchmodules));
278 //AliDebug(2,Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
280 //Add amplitudes in all isolation patch
282 for(Int_t irow = minrow ; irow < maxrow; irow ++)
283 for(Int_t icol = mincol ; icol < maxcol ; icol ++)
284 amp += (*ampmatrix)(irow,icol);
286 AliDebug(2,Form("Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
289 // AliError(Form("Bad sum: Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
290 // ampmatrix->Print();
293 amp-=maxamp; //Calculate energy in isolation patch that do not comes from maximum patch.
296 AliDebug(2, Form("Maximum amplitude %f, Out of patch %f",maxamp, amp));
298 //Fill isolation amplitude data member and say if patch is isolated.
299 if(iPatchType == 0){ //2x2 case
300 f2x2AmpOutOfPatch = amp;
301 if(amp < f2x2AmpOutOfPatchThres) b=kTRUE;
302 } else if(iPatchType == 1){ //nxn case
303 fnxnAmpOutOfPatch = amp;
304 if(amp < fnxnAmpOutOfPatchThres) b=kTRUE;
307 if(keyPrint) printf(" IsPatchIsolated - OUT \n");
314 //____________________________________________________________________________
315 void AliEMCALTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus, const Int_t isupermod,TMatrixD &max2, TMatrixD &maxn){
317 //Sums energy of all possible 2x2 (L0) and nxn (L1) modules per each TRU.
318 //Fast signal in the experiment is given by 2x2 modules,
319 //for this reason we loop inside the TRU modules by 2.
321 //Declare and initialize variables
322 Int_t nCellsPhi = fGeom->GetNCellsInTRUPhi();
324 nCellsPhi = nCellsPhi / 2 ; //Half size SM. Not Final.
325 // 12(tow)*2(cell)/1 TRU, cells in Phi in one TRU
326 Int_t nCellsEta = fGeom->GetNCellsInTRUEta();
327 Int_t nTRU = fGeom->GetNTRU();
328 // 24(mod)*2(tower)/3 TRU, cells in Eta in one TRU
329 //Int_t nTRU = geom->GeNTRU();//3 TRU per super module
333 for(Int_t i = 0; i < 4; i++){
334 for(Int_t j = 0; j < nTRU; j++){
340 //Create matrix that will contain 2x2 amplitude sums
341 //used to calculate the nxn sums
342 TMatrixD tru2x2(nCellsPhi/2,nCellsEta/2) ;
343 for(Int_t i = 0; i < nCellsPhi/2; i++)
344 for(Int_t j = 0; j < nCellsEta/2; j++)
347 //Loop over all TRUS in a supermodule
348 for(Int_t itru = 0 + isupermod * nTRU ; itru < (isupermod+1)*nTRU ; itru++) {
349 TMatrixD * amptru = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
350 TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
351 Int_t mtru = itru-isupermod*nTRU ; //Number of TRU in Supermodule
353 //Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
354 for(Int_t irow = 0 ; irow < nCellsPhi; irow += 2){
355 for(Int_t icol = 0 ; icol < nCellsEta ; icol += 2){
356 amp2 = (*amptru)(irow,icol)+(*amptru)(irow+1,icol)+
357 (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
359 //Fill matrix with added 2x2 cells for use in nxn sums
360 tru2x2(irow/2,icol/2) = amp2 ;
361 //Select 2x2 maximum sums to select L0
362 if(amp2 > ampmax2(0,mtru)){
363 ampmax2(0,mtru) = amp2 ;
364 ampmax2(1,mtru) = irow;
365 ampmax2(2,mtru) = icol;
370 //Find most recent time in the selected 2x2 cell
371 ampmax2(3,mtru) = 1 ;
372 Int_t row2 = static_cast <Int_t> (ampmax2(1,mtru));
373 Int_t col2 = static_cast <Int_t> (ampmax2(2,mtru));
374 for(Int_t i = 0; i<2; i++){
375 for(Int_t j = 0; j<2; j++){
376 if((*amptru)(row2+i,col2+j) > 0 && (*timeRtru)(row2+i,col2+j)> 0){
377 if((*timeRtru)(row2+i,col2+j) < ampmax2(3,mtru) )
378 ampmax2(3,mtru) = (*timeRtru)(row2+i,col2+j);
383 //Sliding nxn, add nxn amplitudes (OVERLAP)
385 for(Int_t irow = 0 ; irow < nCellsPhi/2; irow++){
386 for(Int_t icol = 0 ; icol < nCellsEta/2 ; icol++){
388 if( (irow+fPatchSize) < nCellsPhi/2 && (icol+fPatchSize) < nCellsEta/2){//Avoid exit the TRU
389 for(Int_t i = 0 ; i <= fPatchSize ; i++)
390 for(Int_t j = 0 ; j <= fPatchSize ; j++)
391 ampn += tru2x2(irow+i,icol+j);
392 //Select nxn maximum sums to select L1
393 if(ampn > ampmaxn(0,mtru)){
394 ampmaxn(0,mtru) = ampn ;
395 ampmaxn(1,mtru) = irow*2;
396 ampmaxn(2,mtru) = icol*2;
402 //Find most recent time in selected nxn cell
403 ampmaxn(3,mtru) = 1 ;
404 Int_t rown = static_cast <Int_t> (ampmaxn(1,mtru));
405 Int_t coln = static_cast <Int_t> (ampmaxn(2,mtru));
406 for(Int_t i = 0; i<4*fPatchSize; i++){
407 for(Int_t j = 0; j<4*fPatchSize; j++){
408 if( (rown+i) < nCellsPhi && (coln+j) < nCellsEta){//Avoid exit the TRU
409 if((*amptru)(rown+i,coln+j) > 0 && (*timeRtru)(rown+i,coln+j)> 0){
410 if((*timeRtru)(rown+i,coln+j) < ampmaxn(3,mtru) )
411 ampmaxn(3,mtru) = (*timeRtru)(rown+i,coln+j);
418 ampmaxn(0,mtru) = ampmax2(0,mtru);
419 ampmaxn(1,mtru) = ampmax2(1,mtru);
420 ampmaxn(2,mtru) = ampmax2(2,mtru);
421 ampmaxn(3,mtru) = ampmax2(3,mtru);
426 //____________________________________________________________________________
427 void AliEMCALTrigger::MakeSlidingTowers(const TClonesArray * amptrus, const TClonesArray * timeRtrus,
428 const Int_t isupermod,TMatrixD &max2, TMatrixD &maxn){
430 // Output from module (2x2 cells from one module)
431 Int_t nModulesPhi = fGeom->GetNModulesInTRUPhi(); // now 4 modules (3 div in phi)
433 nModulesPhi = nModulesPhi / 2 ; // Half size SM. Not Final.
435 Int_t nModulesEta = fGeom->GetNModulesInTRUEta(); // now 24 modules (no division in eta)
436 Int_t nTRU = fGeom->GetNTRU();
437 static int keyPrint = 0;
438 if(keyPrint) printf("MakeSlidingTowers : nTRU %i nModulesPhi %i nModulesEta %i ",
439 nTRU, nModulesPhi, nModulesEta );
443 for(Int_t i = 0; i < 4; i++){
444 for(Int_t j = 0; j < nTRU; j++){
445 ampmax2(i,j) = ampmaxn(i,j) = -1;
449 // Create matrix that will contain 2x2 amplitude sums
450 // used to calculate the nxn sums
451 TMatrixD tru2x2(nModulesPhi/2,nModulesEta/2);
453 // Loop over all TRUS in a supermodule
454 for(Int_t itru = 0 + isupermod * nTRU ; itru < (isupermod+1)*nTRU ; itru++) {
455 TMatrixD * amptru = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
456 TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
457 Int_t mtru = itru - isupermod*nTRU ; // Number of TRU in Supermodule !!
459 // Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
460 for(Int_t irow = 0 ; irow < nModulesPhi; irow +=2){
461 for(Int_t icol = 0 ; icol < nModulesEta ; icol +=2){
462 amp2 = (*amptru)(irow,icol) +(*amptru)(irow+1,icol)+
463 (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
465 //Fill matrix with added 2x2 towers for use in nxn sums
466 tru2x2(irow/2,icol/2) = amp2 ;
467 //Select 2x2 maximum sums to select L0
468 if(amp2 > ampmax2(0,mtru)){
469 ampmax2(0,mtru) = amp2 ;
470 ampmax2(1,mtru) = irow;
471 ampmax2(2,mtru) = icol;
476 ampmax2(3,mtru) = 0.;
478 // Find most recent time in the selected 2x2 towers
479 Int_t row2 = static_cast <Int_t> (ampmax2(1,mtru));
480 Int_t col2 = static_cast <Int_t> (ampmax2(2,mtru));
481 for(Int_t i = 0; i<2; i++){
482 for(Int_t j = 0; j<2; j++){
483 if((*amptru)(row2+i,col2+j) > 0 && (*timeRtru)(row2+i,col2+j)> 0){
484 if((*timeRtru)(row2+i,col2+j) > ampmax2(3,mtru) )
485 ampmax2(3,mtru) = (*timeRtru)(row2+i,col2+j); // max time
491 //Sliding nxn, add nxn amplitudes (OVERLAP)
493 for(Int_t irow = 0 ; irow < nModulesPhi/2; irow++){
494 for(Int_t icol = 0 ; icol < nModulesEta/2; icol++){
496 if( (irow+fPatchSize) < nModulesPhi/2 && (icol+fPatchSize) < nModulesEta/2){ //Avoid exit the TRU
497 for(Int_t i = 0 ; i <= fPatchSize ; i++)
498 for(Int_t j = 0 ; j <= fPatchSize ; j++)
499 ampn += tru2x2(irow+i,icol+j);
500 //Select nxn maximum sums to select L1
501 if(ampn > ampmaxn(0,mtru)){
502 ampmaxn(0,mtru) = ampn ;
503 ampmaxn(1,mtru) = irow;
504 ampmaxn(2,mtru) = icol;
510 ampmaxn(3,mtru) = 0.; // Was 1 , I don't know why
512 //Find most recent time in selected nxn cell
513 Int_t rown = static_cast <Int_t> (ampmaxn(1,mtru));
514 Int_t coln = static_cast <Int_t> (ampmaxn(2,mtru));
515 for(Int_t i = 0; i<4*fPatchSize; i++){
516 for(Int_t j = 0; j<4*fPatchSize; j++){
517 if( (rown+i) < nModulesPhi && (coln+j) < nModulesEta){//Avoid exit the TRU
518 if((*amptru)(rown+i,coln+j) > 0 && (*timeRtru)(rown+i,coln+j)> 0){
519 if((*timeRtru)(rown+i,coln+j) > ampmaxn(3,mtru) )
520 ampmaxn(3,mtru) = (*timeRtru)(rown+i,coln+j); // max time
526 } else { // copy 2x2 to nxn
527 ampmaxn(0,mtru) = ampmax2(0,mtru);
528 ampmaxn(1,mtru) = ampmax2(1,mtru);
529 ampmaxn(2,mtru) = ampmax2(2,mtru);
530 ampmaxn(3,mtru) = ampmax2(3,mtru);
533 if(keyPrint) printf(" : MakeSlidingTowers -OUt \n");
536 //____________________________________________________________________________
537 void AliEMCALTrigger::Print(const Option_t * opt) const
540 //Prints main parameters
544 AliTriggerInput* in = 0x0 ;
545 printf( " fSimulation %i (input option) : #digits %i\n", fSimulation, fDigitsList->GetEntries());
546 printf( " fTimeKey %i \n ", fTimeKey);
548 printf( " Maximum Amplitude after Sliding Cell, \n") ;
549 printf( " -2x2 cells sum (not overlapped): %10.2f, in Super Module %d\n",
551 printf( " -2x2 from row %d to row %d and from column %d to column %d\n", f2x2ModulePhi, f2x2ModulePhi+2, f2x2ModuleEta, f2x2ModuleEta+2) ;
552 printf( " -2x2 Isolation Patch %d x %d, Amplitude out of 2x2 patch is %f, threshold %f, Isolated? %d \n",
553 2*fIsolPatchSize+2, 2*fIsolPatchSize+2, f2x2AmpOutOfPatch, f2x2AmpOutOfPatchThres,static_cast<Int_t> (fIs2x2Isol)) ;
555 printf( " Patch Size, n x n: %d x %d cells\n",2*(fPatchSize+1), 2*(fPatchSize+1));
556 printf( " -nxn cells sum (overlapped) : %10.2f, in Super Module %d\n",
558 printf( " -nxn from row %d to row %d and from column %d to column %d\n", fnxnModulePhi, fnxnModulePhi+4*fPatchSize, fnxnModuleEta, fnxnModuleEta+4*fPatchSize) ;
559 printf( " -nxn Isolation Patch %d x %d, Amplitude out of nxn patch is %f, threshold %f, Isolated? %d \n",
560 4*fIsolPatchSize+2*(fPatchSize+1),4*fIsolPatchSize+2*(fPatchSize+1) , fnxnAmpOutOfPatch, fnxnAmpOutOfPatchThres,static_cast<Int_t> (fIsnxnIsol) ) ;
563 printf( " Isolate in SuperModule? %d\n",
564 fIsolateInSuperModule) ;
566 printf( " Threshold for LO %10.2f\n",
568 in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_L0" );
570 printf( " *** EMCAL LO is set ***\n") ;
572 printf( " Gamma Low Pt Threshold for L1 %10.2f\n",
573 fL1GammaLowPtThreshold) ;
574 in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_GammaLPt_L1" );
576 printf( " *** EMCAL Gamma Low Pt for L1 is set ***\n") ;
578 printf( " Gamma Medium Pt Threshold for L1 %10.2f\n",
579 fL1GammaMediumPtThreshold) ;
580 in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_GammaMPt_L1" );
582 printf( " *** EMCAL Gamma Medium Pt for L1 is set ***\n") ;
584 printf( " Gamma High Pt Threshold for L1 %10.2f\n",
585 fL1GammaHighPtThreshold) ;
586 in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_GammaHPt_L1" );
588 printf( " *** 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) printf("<I> AliEMCALTrigger::FillTRU : 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 Error("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");
869 //____________________________________________________________________________
870 void AliEMCALTrigger::Trigger()
872 TH1::AddDirectory(0);
873 //Main Method to select triggers.
874 AliRunLoader *runLoader = AliRunLoader::GetRunLoader();
875 AliEMCALLoader *emcalLoader = 0;
877 emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
880 //Load EMCAL Geometry
881 if (runLoader && runLoader->GetAliRun() && runLoader->GetAliRun()->GetDetector("EMCAL"))
882 fGeom = dynamic_cast<AliEMCAL*>(runLoader->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
885 AliFatal("Did not get geometry from EMCALLoader");
888 Int_t nSuperModules = fGeom->GetNumberOfSuperModules() ; //12 SM in EMCAL
889 Int_t nTRU = fGeom->GetNTRU(); // 3 TRU per super module
891 //Intialize data members each time the trigger is called in event loop
892 f2x2MaxAmp = -1; f2x2ModulePhi = -1; f2x2ModuleEta = -1;
893 fnxnMaxAmp = -1; fnxnModulePhi = -1; fnxnModuleEta = -1;
895 // Take the digits list if simulation
896 if(fSimulation && runLoader && emcalLoader){ // works than run seperate macros
897 runLoader->LoadDigits("EMCAL");
898 fDigitsList = emcalLoader->Digits() ;
899 runLoader->LoadSDigits("EMCAL");
901 // Digits list should be set by method SetDigitsList(TClonesArray * digits)
903 AliFatal("Digits not found !") ;
905 //Take the digits list
907 // Delete old if unzero
908 if(fAmpTrus) {fAmpTrus->Delete(); delete fAmpTrus;}
909 if(fTimeRtrus) {fTimeRtrus->Delete(); delete fTimeRtrus;}
910 if(fAmpSMods) {fAmpSMods->Delete(); delete fAmpSMods;}
911 // Fill TRU and SM matrix
912 fAmpTrus = new TClonesArray("TMatrixD",nTRU);
913 fAmpTrus->SetName("AmpTrus");
914 fTimeRtrus = new TClonesArray("TMatrixD",nTRU);
915 fTimeRtrus->SetName("TimeRtrus");
916 fAmpSMods = new TClonesArray("TMatrixD",nSuperModules);
917 fAmpSMods->SetName("AmpSMods");
919 FillTRU(fDigitsList, fAmpTrus, fAmpSMods, fTimeRtrus);
921 // Jet staff - only one case, no fredom here
922 if(fGeom->GetNEtaSubOfTRU() == 6) {
923 if(fAmpJetMatrix) {delete fAmpJetMatrix; fAmpJetMatrix=0;}
924 if(fJetMatrixE) {delete fJetMatrixE; fJetMatrixE=0;}
926 fAmpJetMatrix = new TMatrixD(17,12); // 17-phi(row), 12-eta(col)
927 fJetMatrixE = new TH2F("fJetMatrixE"," E of max patch in (#phi,#eta)",
928 17, 80.*TMath::DegToRad(), (180.+20.*2/3.)*TMath::DegToRad(), 12, -0.7, 0.7);
929 for(Int_t row=0; row<fAmpJetMatrix->GetNrows(); row++) {
930 for(Int_t col=0; col<fAmpJetMatrix->GetNcols(); col++) {
931 (*fAmpJetMatrix)(row,col) = 0.;
934 FillJetMatrixFromSMs(fAmpSMods, fAmpJetMatrix, fGeom);
936 if(!CheckConsistentOfMatrixes()) assert(0);
938 // Do Tower Sliding and select Trigger
939 // Initialize varible that will contain maximum amplitudes and
940 // its corresponding tower position in eta and phi, and time.
941 TMatrixD ampmax2(4,nTRU) ; // 0-max amp, 1-irow, 2-icol, 3-timeR
942 TMatrixD ampmaxn(4,nTRU) ;
944 for(Int_t iSM = 0 ; iSM < nSuperModules ; iSM++) {
945 //Do 2x2 and nxn sums, select maximums.
947 MakeSlidingTowers(fAmpTrus, fTimeRtrus, iSM, ampmax2, ampmaxn);
950 if(fIsolateInSuperModule) // here some discripency between tru and SM
951 SetTriggers(fAmpSMods,iSM,ampmax2,ampmaxn) ;
952 if(!fIsolateInSuperModule)
953 SetTriggers(fAmpTrus,iSM,ampmax2,ampmaxn) ;
956 // Do patch sliding and select Jet Trigger
957 // 0-max amp-meanFromVZERO(if), 1-irow, 2-icol, 3-timeR,
958 // 4-max amp , 5-meanFromVZERO (Nov 25, 2007)
960 MakeSlidingPatch((*fAmpJetMatrix), fNJetPatchPhi, fAmpJetMax); // no timing information here
966 void AliEMCALTrigger::GetTriggerInfo(TArrayF &triggerPosition, TArrayF &triggerAmplitudes)
968 // Template - should be defined; Nov 5, 2007
969 triggerPosition[0] = 0.;
970 triggerAmplitudes[0] = 0.;
973 void AliEMCALTrigger::FillJetMatrixFromSMs(TClonesArray *ampmatrixsmod, TMatrixD* jetMat, AliEMCALGeometry *g)
976 // Fill matrix for jet trigger from SM matrixes of modules
978 static int keyPrint = 0;
980 if(ampmatrixsmod==0 || jetMat==0 || g==0) return;
981 Double_t amp = 0.0, ampSum=0.0;
983 Int_t nEtaModSum = g->GetNZ() / g->GetNEtaSubOfTRU(); // should be 4
984 Int_t nPhiModSum = g->GetNPhi() / g->GetNTRUPhi(); // should be 4
986 if(keyPrint) printf("%s",Form(" AliEMCALTrigger::FillJetMatrixFromSMs | nEtaModSum %i : nPhiModSum %i \n",
987 nEtaModSum, nPhiModSum));
988 Int_t jrow=0, jcol=0; // indexes of jet matrix
989 Int_t nEtaSM=0, nPhiSM=0;
990 for(Int_t iSM=0; iSM<ampmatrixsmod->GetEntries(); iSM++) {
991 TMatrixD * ampsmods = dynamic_cast<TMatrixD *>(ampmatrixsmod->At(iSM));
992 Int_t nrow = ampsmods->GetNrows();
993 Int_t ncol = ampsmods->GetNcols();
994 //printf("%s",Form(" ######## SM %i : nrow %i : ncol %i ##### \n", iSM, nrow, ncol));
995 for(Int_t row=0; row<nrow; row++) {
996 for(Int_t col=0; col<ncol; col++) {
997 amp = (*ampsmods)(row,col);
1001 if(keyPrint) printf("%s",Form(" ** nPhiSm %i : nEtaSM %i : row %2.2i : col %2.2i -> ",
1002 nPhiSM, nEtaSM, row, col));
1003 if(nEtaSM == 0) { // positive Z
1004 jrow = 3*nPhiSM + row/nPhiModSum;
1005 jcol = 6 + col / nEtaModSum;
1006 } else { // negative Z
1007 if(iSM<=9) jrow = 3*nPhiSM + 2 - row/nPhiModSum;
1008 else jrow = 3*nPhiSM + 1 - row/nPhiModSum; // half size
1009 jcol = 5 - col / nEtaModSum;
1011 if(keyPrint) printf("%s",Form(" jrow %2.2i : jcol %2.2i : amp %f (jetMat) \n", jrow, jcol, amp));
1013 (*jetMat)(jrow,jcol) += amp;
1014 ampSum += amp; // For controling
1015 } else if(amp<0.0) {
1016 printf("%s",Form(" jrow %2.2i : jcol %2.2i : amp %f (jetMat: amp<0) \n", jrow, jcol, amp));
1022 if(ampSum <= 0.0) Warning("FillJetMatrixFromSMs","ampSum %f (<=0.0) ", ampSum);
1025 void AliEMCALTrigger::MakeSlidingPatch(const TMatrixD &jm, const Int_t nPatchSize, TMatrixD &JetMax)
1027 // Sliding patch : nPatchSize x nPatchSize (OVERLAP)
1028 static int keyPrint = 0;
1029 if(keyPrint) printf(" AliEMCALTrigger::MakeSlidingPatch() was started \n");
1030 Double_t ampCur = 0.0, e=0.0;
1031 ampJetMax(0,0) = 0.0;
1032 ampJetMax(3,0) = 0.0; // unused now
1033 ampJetMax(4,0) = ampJetMax(5,0) = 0.0;
1034 for(Int_t row=0; row<fAmpJetMatrix->GetNrows(); row ++) {
1035 for(Int_t col=0; col<fAmpJetMatrix->GetNcols(); col++) {
1037 // check on patch size
1038 if( (row+nPatchSize-1) < fAmpJetMatrix->GetNrows() && (col+nPatchSize-1) < fAmpJetMatrix->GetNcols()){
1039 for(Int_t i = 0 ; i < nPatchSize ; i++) {
1040 for(Int_t j = 0 ; j < nPatchSize ; j++) {
1041 ampCur += jm(row+i, col+j);
1043 } // end cycle on patch
1044 if(ampCur > ampJetMax(0,0)){
1045 ampJetMax(0,0) = ampCur;
1046 ampJetMax(1,0) = row;
1047 ampJetMax(2,0) = col;
1049 } // check on patch size
1052 if(keyPrint) printf(" ampJetMax %i row %2i->%2i col %2i->%2i \n",
1053 Int_t(ampJetMax(0,0)), Int_t(ampJetMax(1,0)), Int_t(ampJetMax(1,0))+nPatchSize-1,
1054 Int_t(ampJetMax(2,0)), Int_t(ampJetMax(2,0))+nPatchSize-1);
1056 Double_t eCorrJetMatrix=0.0;
1057 if(fVZER0Mult > 0.0) {
1058 // Correct patch energy (adc) and jet patch matrix energy
1059 Double_t meanAmpBG = GetMeanEmcalPatchEnergy(Int_t(fVZER0Mult), nPatchSize)/0.0153;
1060 ampJetMax(4,0) = ampJetMax(0,0);
1061 ampJetMax(5,0) = meanAmpBG;
1063 Double_t eCorr = ampJetMax(0,0) - meanAmpBG;
1064 printf(" ampJetMax(0,0) %f meanAmpBG %f eCorr %f : ampJetMax(4,0) %f \n",
1065 ampJetMax(0,0), meanAmpBG, eCorr, ampJetMax(5,0));
1066 ampJetMax(0,0) = eCorr;
1068 eCorrJetMatrix = GetMeanEmcalEnergy(Int_t(fVZER0Mult)) / 208.;
1070 // Fill patch energy matrix
1071 for(int row=Int_t(ampJetMax(1,0)); row<Int_t(ampJetMax(1,0))+nPatchSize; row++) {
1072 for(int col=Int_t(ampJetMax(2,0)); col<Int_t(ampJetMax(2,0))+nPatchSize; col++) {
1073 e = Double_t(jm(row,col)*0.0153); // 0.0153 - hard coded now
1074 if(eCorrJetMatrix > 0.0) { // BG subtraction case
1075 e -= eCorrJetMatrix;
1076 fJetMatrixE->SetBinContent(row+1, col+1, e);
1077 } else if(e > 0.0) {
1078 fJetMatrixE->SetBinContent(row+1, col+1, e);
1082 // PrintJetMatrix();
1083 // Set the jet trigger(s), multiple threshold now, Nov 19,2007
1084 for(Int_t i=0; i<fNJetThreshold; i++ ) {
1085 if(ampJetMax(0,0) >= fL1JetThreshold[i]) {
1086 SetInput((Form("%s_Th%2i", fgNameOfJetTriggers.Data(),i)));
1091 Double_t AliEMCALTrigger::GetEmcalSumAmp() const
1093 // Return sum of amplidutes from EMCal
1094 // Used calibration coefficeint for transition to energy
1095 return fAmpJetMatrix >0 ?fAmpJetMatrix->Sum() :0.0;
1099 void AliEMCALTrigger::PrintJetMatrix() const
1101 // fAmpJetMatrix : (17,12); // 17-phi(row), 12-eta(col)
1102 if(fAmpJetMatrix == 0) return;
1104 printf("\n #### jetMatrix : (%i,%i) ##### \n ",
1105 fAmpJetMatrix->GetNrows(), fAmpJetMatrix->GetNcols());
1106 PrintMatrix(*fAmpJetMatrix);
1109 void AliEMCALTrigger::PrintAmpTruMatrix(Int_t ind) const
1111 TMatrixD * tru = dynamic_cast<TMatrixD *>(fAmpTrus->At(ind));
1112 if(tru == 0) return;
1113 printf("\n #### Amp TRU matrix(%i) : (%i,%i) ##### \n ",
1114 ind, tru->GetNrows(), tru->GetNcols());
1118 void AliEMCALTrigger::PrintAmpSmMatrix(Int_t ind) const
1120 TMatrixD * sm = dynamic_cast<TMatrixD *>(fAmpSMods->At(ind));
1122 printf("\n #### Amp SM matrix(%i) : (%i,%i) ##### \n ",
1123 ind, sm->GetNrows(), sm->GetNcols());
1127 void AliEMCALTrigger::PrintMatrix(const TMatrixD &mat) const
1129 for(Int_t col=0; col<mat.GetNcols(); col++) printf(" %3i ", col);
1131 for(Int_t row=0; row<mat.GetNrows(); row++) {
1132 printf(" row:%2i ", row);
1133 for(Int_t col=0; col<mat.GetNcols(); col++) {
1134 printf(" %4i", (Int_t)mat(row,col));
1140 Bool_t AliEMCALTrigger::CheckConsistentOfMatrixes(const Int_t pri)
1142 Double_t sumSM = 0.0, smCur=0.0;
1143 Double_t sumTru=0.0, sumTruInSM = 0.0, truSum=0.0;
1144 // Bool_t key = kTRUE;
1146 for(Int_t i=0; i<fAmpSMods->GetEntries(); i++) {
1147 TMatrixD * sm = dynamic_cast<TMatrixD *>(fAmpSMods->At(i));
1153 for(Int_t itru=0; itru<3; itru++) { // Cycle on tru inside SM
1154 Int_t ind = 3*i + itru;
1155 TMatrixD *tru = dynamic_cast<TMatrixD *>(fAmpTrus->At(ind));
1157 truSum = tru->Sum();
1158 sumTruInSM += truSum;
1161 sumTru += sumTruInSM;
1163 if(sumTruInSM != smCur) {
1164 printf(" sm %i : smCur %f -> sumTruInSM %f \n", i, smCur, sumTruInSM);
1169 Double_t sumJetMat = fAmpJetMatrix->Sum();
1170 if(pri || sumSM != sumTru || sumSM != sumJetMat)
1171 printf(" sumSM %f : sumTru %f : sumJetMat %f \n", sumSM, sumTru, sumJetMat);
1172 if(sumSM != sumTru || sumSM != sumJetMat) return kFALSE;
1177 void AliEMCALTrigger::Browse(TBrowser* b)
1179 if(&fInputs) b->Add(&fInputs);
1180 if(fAmpTrus) b->Add(fAmpTrus);
1181 if(fTimeRtrus) b->Add(fTimeRtrus);
1182 if(fAmpSMods) b->Add(fAmpSMods);
1183 if(fAmpJetMatrix) b->Add(fAmpJetMatrix);
1184 if(fJetMatrixE) b->Add(fJetMatrixE);