PWG/EMCAL/AliEmcalTriggerMaker.cxx

   1 // $Id$
   2 //
   3 // Class to make emcal particles in AOD/ESD events.
   4 //
   5 // Author: J.Kral
   6
   7 #include <TClonesArray.h>
   8 #include <iostream>
   9
  10 #include "AliLog.h"
  11 #include "AliEmcalTriggerPatchInfo.h"
  12 #include "AliEmcalTriggerSetupInfo.h"
  13
  14 #include "AliEmcalTriggerMaker.h"
  15 #include "AliEMCALTriggerTypes.h"
  16 #include "AliEMCALGeometry.h"
  17
  18 #include "AliVCaloTrigger.h"
  19 #include "AliAODCaloTrigger.h"
  20 #include "AliVCaloCells.h"
  21 #include "AliVVZERO.h"
  22
  23 ClassImp(AliEmcalTriggerMaker)
  24
  25 using namespace std;
  26
  27 //________________________________________________________________________
  28 AliEmcalTriggerMaker::AliEmcalTriggerMaker() :
  29   AliAnalysisTaskEmcal("AliEmcalTriggerMaker",kFALSE),
  30   fCaloTriggersOutName("EmcalTriggers"),
  31   fCaloTriggerSetupOutName("EmcalTriggersSetup"),
  32   fV0InName("AliAODVZERO"),
  33   fCaloTriggersOut(0),
  34   fCaloTriggerSetupOut(0),
  35   fSimpleOfflineTriggers(0),
  36   fV0(0)
  37 {
  38   // Constructor.
  39   for( int i = 0; i < 4; i++ )
  40     for( int j = 0; j < 3; j++ )
  41       fThresholdConstants[i][j] = 0;
  42 }
  43
  44 //________________________________________________________________________
  45 AliEmcalTriggerMaker::AliEmcalTriggerMaker(const char *name) :
  46   AliAnalysisTaskEmcal(name,kFALSE),
  47   fCaloTriggersOutName("EmcalTriggers"),
  48   fCaloTriggerSetupOutName("EmcalTriggersSetup"),
  49   fV0InName("AliAODVZERO"),
  50   fCaloTriggersOut(0),
  51   fCaloTriggerSetupOut(0),
  52   fSimpleOfflineTriggers(0),
  53   fV0(0)
  54 {
  55   // Constructor.
  56   for( int i = 0; i < 4; i++ )
  57     for( int j = 0; j < 3; j++ )
  58       fThresholdConstants[i][j] = 0;
  59 }
  60
  61 //________________________________________________________________________
  62 AliEmcalTriggerMaker::~AliEmcalTriggerMaker()
  63 {
  64   // Destructor.
  65 }
  66
  67 //________________________________________________________________________
  68 void AliEmcalTriggerMaker::ExecOnce()
  69 {
  70   // Init the analysis.
  71
  72   AliAnalysisTaskEmcal::ExecOnce();
  73
  74   if (!fInitialized)
  75     return;
  76
  77   if (!fCaloTriggersOutName.IsNull()) {
  78     fCaloTriggersOut = new TClonesArray("AliEmcalTriggerPatchInfo");
  79     fCaloTriggersOut->SetName(fCaloTriggersOutName);
  80
  81     if (!(InputEvent()->FindListObject(fCaloTriggersOutName))) {
  82       InputEvent()->AddObject(fCaloTriggersOut);
  83     }
  84     else {
  85       fInitialized = kFALSE;
  86       AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fCaloTriggersOutName.Data()));
  87       return;
  88     }
  89   }
  90
  91   if (!fCaloTriggerSetupOutName.IsNull()) {
  92     fCaloTriggerSetupOut = new AliEmcalTriggerSetupInfo();
  93     fCaloTriggerSetupOut->SetName(fCaloTriggerSetupOutName);
  94
  95     if (!(InputEvent()->FindListObject(fCaloTriggerSetupOutName))) {
  96       InputEvent()->AddObject(fCaloTriggerSetupOut);
  97     }
  98     else {
  99       fInitialized = kFALSE;
 100       AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fCaloTriggerSetupOutName.Data()));
 101       return;
 102     }
 103   }
 104
 105   if ( ! fV0InName.IsNull()) {
 106     fV0 = (AliVVZERO*)InputEvent()->FindListObject( fV0InName );
 107   }
 108
 109   // container for simple offline trigger processing
 110   fSimpleOfflineTriggers = new AliAODCaloTrigger();
 111   fSimpleOfflineTriggers->Allocate( 0 );
 112 }
 113
 114 //________________________________________________________________________
 115 Bool_t AliEmcalTriggerMaker::Run()
 116 {
 117   // Create the emcal particles
 118
 119   Short_t cellId;
 120   Int_t globCol, globRow, tBits, cellAbsId[4], v0[2];
 121   Int_t absId, adcAmp;
 122   Int_t i, j, k, iMain, iMainSimple, cmCol, cmRow, cmiCellCol, cmiCellRow, nCell, iCell;
 123   Int_t jetTrigger, iTriggers;
 124         Int_t patchADC[48][64];
 125   Double_t amp, ca, eMain, eMainSimple, cmiCol, cmiRow;
 126   ULong64_t v0S, thresh;
 127   Bool_t isOfflineSimple;
 128
 129   TVector3 centerGeo, center1, center2, centerMass, edge1, edge2, vertex;
 130
 131   AliEmcalTriggerPatchInfo *trigger;
 132
 133   // delete patch array, clear setup object
 134   fCaloTriggersOut->Delete();
 135   fCaloTriggerSetupOut->Clean();
 136
 137   if( !fCaloTriggers ){
 138     AliError(Form("Calo triggers container %s not available.", fCaloTriggersName.Data()));
 139     return kTRUE;
 140   }
 141   if( !fCaloCells ){
 142     AliError(Form("Calo cells container %s not available.", fCaloCellsName.Data()));
 143     return kTRUE;
 144   }
 145   if( !fCaloCells ){
 146     AliError(Form("V0 container %s not available.", fV0InName.Data()));
 147     return kTRUE;
 148   }
 149
 150   // must reset before usage, or the class will fail
 151   fCaloTriggers->Reset();
 152   for (i=0; i<2; i++) {
 153     fEGA[i] = 0;
 154     fEJE[i] = 0;
 155   }
 156   // first run over the patch array to compose a map of 2x2 patch energies
 157   // which is then needed to construct the full patch ADC energy
 158   // class is not empty
 159   Int_t isMC = 0;
 160   if (MCEvent()) isMC = 1;
 161   Int_t offSet = (1 - isMC) * kTriggerTypeEnd;
 162
 163   if( fCaloTriggers->GetEntries() > 0 ){
 164
 165     // zero the array
 166     for( i = 0; i < 48; i++ )
 167       for( j = 0; j < 64; j++ )
 168         patchADC[i][j] = 0;
 169
 170     // go throuth the trigger channels
 171     while( fCaloTriggers->Next() ){
 172       // get position in global 2x2 tower coordinates
 173       // A0 left bottom (0,0)
 174       fCaloTriggers->GetPosition( globCol, globRow );
 175
 176       // for some strange reason some ADC amps are initialized in reconstruction
 177       // as -1, neglect those :\\ wth
 178       fCaloTriggers->GetL1TimeSum( adcAmp );
 179                         if( adcAmp > -1 )
 180                                 patchADC[globCol][globRow] = adcAmp;
 181
 182       fCaloTriggers->GetTriggerBits( tBits );
 183
 184       if (tBits) {
 185         if ((tBits >> (offSet + kL1GammaHigh)) & 1 ) fEGA[0] = 1;
 186         if ((tBits >> (offSet + kL1GammaLow )) & 1 ) fEGA[1] = 1;
 187         if ((tBits >> (offSet + kL1JetHigh  )) & 1 ) fEJE[0] = 1;
 188         if ((tBits >> (offSet + kL1JetLow   )) & 1 ) fEJE[1] = 1;
 189       }
 190     } // patches
 191   } // array not empty
 192
 193   // fill the array for offline trigger processing
 194   // using calibrated cell energies
 195   for( i = 0; i < 48; i++ )
 196     for( j = 0; j < 64; j++ )
 197       fPatchADCSimple[i][j] = 0;
 198
 199   // fill the patch ADCs from cells
 200   nCell = fCaloCells->GetNumberOfCells();
 201
 202   for( iCell = 0; iCell < nCell; iCell++ ){
 203     // get the cell info, based in index in array
 204     cellId = fCaloCells->GetCellNumber( iCell );
 205     amp = fCaloCells->GetAmplitude( iCell );
 206
 207     // get position
 208     fGeom->GetFastORIndexFromCellIndex( cellId, absId );
 209     fGeom->GetPositionInEMCALFromAbsFastORIndex( absId, globCol, globRow );
 210
 211     // add
 212     fPatchADCSimple[globCol][globRow] += amp/kEMCL1ADCtoGeV;
 213   }
 214
 215   // dig out common data (thresholds)
 216   // 0 - jet high, 1 - gamma high, 2 - jet low, 3 - gamma low
 217   fCaloTriggerSetupOut->SetThresholds( fCaloTriggers->GetL1Threshold( 0 ),
 218                                        fCaloTriggers->GetL1Threshold( 1 ),
 219                                        fCaloTriggers->GetL1Threshold( 2 ),
 220                                        fCaloTriggers->GetL1Threshold( 3 ));
 221
 222   // get the V0 value and compute and set the offline thresholds
 223   // get V0, compute thresholds and save them as global parameters
 224   v0[0] = fV0->GetTriggerChargeA();
 225   v0[1] = fV0->GetTriggerChargeC();
 226   v0S = v0[0] + v0[1];
 227
 228   fSimpleOfflineTriggers->SetL1V0( v0 );
 229
 230   for( i = 0; i < 4; i++ ){
 231     // A*V0^2/2^32+B*V0/2^16+C
 232     thresh = ( ((ULong64_t)fThresholdConstants[i][0]) * v0S * v0S ) >> 32;
 233     thresh += ( ((ULong64_t)fThresholdConstants[i][1]) * v0S ) >> 16;
 234     thresh += ((ULong64_t)fThresholdConstants[i][2]);
 235     fSimpleOfflineTriggers->SetL1Threshold( i, thresh );
 236   }
 237
 238   // save the thresholds in output object
 239   fCaloTriggerSetupOut->SetThresholdsSimple( fSimpleOfflineTriggers->GetL1Threshold( 0 ),
 240                                        fSimpleOfflineTriggers->GetL1Threshold( 1 ),
 241                                        fSimpleOfflineTriggers->GetL1Threshold( 2 ),
 242                                        fSimpleOfflineTriggers->GetL1Threshold( 3 ));
 243
 244   // run the trigger
 245   RunSimpleOfflineTrigger();
 246
 247   // reset for re-run
 248   fCaloTriggers->Reset();
 249   fSimpleOfflineTriggers->Reset();
 250
 251   // class is not empty
 252   if( fCaloTriggers->GetEntries() > 0 ||  fSimpleOfflineTriggers->GetEntries() > 0 ){
 253
 254     iTriggers = 0;
 255     iMain = -1;
 256     eMain = -1;
 257     iMainSimple = -1;
 258     eMainSimple = -1;
 259
 260     // save primary vertex in vector
 261     vertex.SetXYZ( fVertex[0], fVertex[1], fVertex[2] );
 262
 263     // go throuth the trigger channels, real first, then offline
 264     while( NextTrigger( isOfflineSimple ) ){
 265       // check if jet trigger low or high
 266       if( ! isOfflineSimple )
 267         fCaloTriggers->GetTriggerBits( tBits );
 268       else
 269         fSimpleOfflineTriggers->GetTriggerBits( tBits );
 270
 271       jetTrigger = 0;
 272       if(( tBits >> ( offSet + kL1JetLow )) & 1 )
 273         jetTrigger = 1;
 274       if(( tBits >> ( offSet + kL1JetHigh )) & 1)
 275         jetTrigger = jetTrigger | 2;
 276
 277       if( jetTrigger == 0 )
 278         continue;
 279
 280       // get position in global 2x2 tower coordinates
 281       // A0 left bottom (0,0)
 282       if( ! isOfflineSimple )
 283         fCaloTriggers->GetPosition( globCol, globRow );
 284       else
 285         fSimpleOfflineTriggers->GetPosition( globCol, globRow );
 286
 287       // get the absolute trigger ID
 288       fGeom->GetAbsFastORIndexFromPositionInEMCAL( globCol, globRow, absId );
 289       // convert to the 4 absId of the cells composing the trigger channel
 290       fGeom->GetCellIndexFromFastORIndex( absId, cellAbsId );
 291
 292       // get low left edge (eta max, phi min)
 293       fGeom->GetGlobal( cellAbsId[0], edge1 );
 294
 295       // sum the available energy in the 32/32 window of cells
 296       // step over trigger channels and get all the corresponding cells
 297       // make CM
 298       amp = 0;
 299       cmiCol = 0;
 300       cmiRow = 0;
 301       adcAmp = 0;
 302       for( i = 0; i < 16; i++ ){
 303         for( j = 0; j < 16; j++ ){
 304           // get the 4 cells composing the trigger channel
 305           fGeom->GetAbsFastORIndexFromPositionInEMCAL( globCol+i, globRow+j, absId );
 306           fGeom->GetCellIndexFromFastORIndex( absId, cellAbsId );
 307           // add amplitudes and find patch edges
 308           for( k = 0; k < 4; k++ ){
 309             ca = fCaloCells->GetCellAmplitude( cellAbsId[k] );
 310             amp += ca;
 311             cmiCol += ca*(Double_t)i;
 312             cmiRow += ca*(Double_t)j;
 313           }
 314           // add the STU ADCs in the patch
 315           if( ! isOfflineSimple )
 316           adcAmp += patchADC[globCol+i][globRow+j];
 317           else
 318             adcAmp += fPatchADCSimple[globCol+i][globRow+j];
 319         }
 320       } // 32x32 cell window
 321       if( amp == 0 ){
 322         AliDebug(2,"EMCal trigger patch with 0 energy.");
 323         continue;
 324       }
 325
 326       // get the CM and patch index
 327       cmiCol /= amp;
 328       cmiRow /= amp;
 329       cmCol = globCol + (Int_t)cmiCol;
 330       cmRow = globRow + (Int_t)cmiRow;
 331
 332       // get the patch and corresponding cells
 333       fGeom->GetAbsFastORIndexFromPositionInEMCAL( cmCol, cmRow, absId );
 334       fGeom->GetCellIndexFromFastORIndex( absId, cellAbsId );
 335
 336       // find which out of the 4 cells is closest to CM and get it's position
 337       cmiCellCol = TMath::Nint( cmiCol * 2. );
 338       cmiCellRow = TMath::Nint( cmiRow * 2. );
 339       fGeom->GetGlobal( cellAbsId[(cmiCellRow%2)*2 + cmiCellCol%2], centerMass );
 340
 341       // get up right edge (eta min, phi max)
 342       // get the absolute trigger ID
 343       fGeom->GetAbsFastORIndexFromPositionInEMCAL( globCol+15, globRow+15, absId );
 344       // convert to the 4 absId of the cells composing the trigger channel
 345       fGeom->GetCellIndexFromFastORIndex( absId, cellAbsId );
 346
 347       fGeom->GetGlobal( cellAbsId[3], edge2 );
 348
 349       // get the geometrical center as an average of two diagonally
 350       // adjacent patches in the center
 351       // picking two diagonally closest cells from the patches
 352       fGeom->GetAbsFastORIndexFromPositionInEMCAL( globCol+7, globRow+7, absId );
 353       fGeom->GetCellIndexFromFastORIndex( absId, cellAbsId );
 354       fGeom->GetGlobal( cellAbsId[3], center1 );
 355
 356       fGeom->GetAbsFastORIndexFromPositionInEMCAL( globCol+8, globRow+8, absId );
 357       fGeom->GetCellIndexFromFastORIndex( absId, cellAbsId );
 358       fGeom->GetGlobal( cellAbsId[0], center2 );
 359
 360       centerGeo = center1;
 361       centerGeo += center2;
 362       centerGeo *= 0.5;
 363
 364       // relate all to primary vertex
 365       centerGeo -= vertex;
 366       centerMass -= vertex;
 367       edge1 -= vertex;
 368       edge2 -= vertex;
 369
 370       // save the trigger object
 371       new ((*fCaloTriggersOut)[iTriggers])AliEmcalTriggerPatchInfo();
 372       trigger = (AliEmcalTriggerPatchInfo*)fCaloTriggersOut->At( iTriggers );
 373       iTriggers++;
 374
 375       trigger->SetCenterGeo( centerGeo, amp );
 376       trigger->SetCenterMass( centerMass, amp );
 377       trigger->SetEdge1( edge1, amp );
 378       trigger->SetEdge2( edge2, amp );
 379       trigger->SetADCAmp( adcAmp );
 380       trigger->SetTriggerBits( tBits );
 381       trigger->SetEdgeCell( globCol*2, globRow*2 ); // from triggers to cells
 382       trigger->SetOffSet(offSet);
 383
 384       // check if more energetic than others for main patch marking
 385       if( ! isOfflineSimple && eMain < amp ){
 386         eMain = amp;
 387         iMain = iTriggers - 1;
 388       }
 389       if( isOfflineSimple && eMainSimple < amp ){
 390         eMainSimple = amp;
 391         iMainSimple = iTriggers - 1;
 392       }
 393     } // triggers
 394
 395     // mark the most energetic patch as main
 396     // for real and also simple offline
 397     if( iMain > -1 ){
 398       trigger = (AliEmcalTriggerPatchInfo*)fCaloTriggersOut->At( iMain );
 399       tBits = trigger->GetTriggerBits();
 400       // main trigger flag
 401       tBits = tBits | ( 1 << 24 );
 402       trigger->SetTriggerBits( tBits );
 403     }
 404     if( iMainSimple > -1 ){
 405       trigger = (AliEmcalTriggerPatchInfo*)fCaloTriggersOut->At( iMainSimple );
 406       tBits = trigger->GetTriggerBits();
 407       // main trigger flag
 408       tBits = tBits | ( 1 << 24 );
 409       trigger->SetTriggerBits( tBits );
 410     }
 411
 412   } // there are some triggers
 413
 414   return kTRUE;
 415 }
 416
 417 //________________________________________________________________________
 418 void AliEmcalTriggerMaker::RunSimpleOfflineTrigger()
 419 {
 420   // runs simple offline trigger algorithm
 421
 422   Int_t i, j, k, l, tBits, tSum;
 423   TArrayI tBitsArray, rowArray, colArray;
 424
 425   Int_t isMC = 0;
 426   if (MCEvent()) isMC = 1;
 427   Int_t offSet = (1 - isMC) * kTriggerTypeEnd;
 428
 429   // 0 thresholds = no processing
 430   if( fCaloTriggerSetupOut->GetThresholdJetLowSimple() == 0 &&
 431     fCaloTriggerSetupOut->GetThresholdJetHighSimple() == 0 )
 432     return;
 433
 434   // run the trigger algo, stepping by 8 towers (= 4 trigger channels)
 435   for( i = 0; i < 32; i += 4 ){
 436     for( j = 0; j < 48; j += 4 ){
 437
 438       tSum = 0;
 439       tBits = 0;
 440
 441       // window
 442       for( k = 0; k < 16; k++ )
 443         for( l = 0; l < 16; l++ )
 444           tSum += (ULong64_t)fPatchADCSimple[i+k][j+l];
 445
 446       // check thresholds
 447       if( tSum > fCaloTriggerSetupOut->GetThresholdJetLowSimple() )
 448         tBits = tBits | ( 1 << ( offSet + kL1JetLow ));
 449       if( tSum > fCaloTriggerSetupOut->GetThresholdJetHighSimple() )
 450         tBits = tBits | ( 1 << ( offSet + kL1JetHigh ));
 451
 452       // add trigger values
 453       if( tBits != 0 ){
 454         // add offline bit
 455         tBits = tBits | ( 1 << 25 );
 456
 457         tBitsArray.Set( tBitsArray.GetSize() + 1 );
 458         colArray.Set( colArray.GetSize() + 1 );
 459         rowArray.Set( rowArray.GetSize() + 1 );
 460
 461         tBitsArray[tBitsArray.GetSize()-1] = tBits;
 462         colArray[colArray.GetSize()-1] = i;
 463         rowArray[rowArray.GetSize()-1] = j;
 464       }
 465     }
 466   } // trigger algo
 467
 468   // save in object
 469   fSimpleOfflineTriggers->DeAllocate();
 470   fSimpleOfflineTriggers->Allocate( tBitsArray.GetSize() );
 471
 472   for( i = 0; i < tBitsArray.GetSize(); i++ ){
 473     fSimpleOfflineTriggers->Add( colArray[i], rowArray[i], 0, 0, 0, 0, 0, tBitsArray[i] );
 474   }
 475
 476 }
 477
 478 //________________________________________________________________________
 479 Bool_t AliEmcalTriggerMaker::NextTrigger( Bool_t &isOfflineSimple )
 480 {
 481   // get next trigger
 482
 483   Bool_t loopContinue;
 484
 485   loopContinue = fCaloTriggers->Next();
 486   isOfflineSimple = kFALSE;
 487
 488   if( ! loopContinue ){
 489     loopContinue = fSimpleOfflineTriggers->Next();
 490     isOfflineSimple = kTRUE;
 491   }
 492
 493   return loopContinue;
 494 }