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 "AliVCaloCells.h"
  20
  21 ClassImp(AliEmcalTriggerMaker)
  22
  23 using namespace std;
  24
  25 //________________________________________________________________________
  26 AliEmcalTriggerMaker::AliEmcalTriggerMaker() :
  27   AliAnalysisTaskEmcal("AliEmcalTriggerMaker",kFALSE),
  28   fCaloTriggersOutName("EmcalTriggers"),
  29   fCaloTriggerSetupOutName("EmcalTriggersSetup"),
  30   fCaloTriggersOut(0),
  31   fCaloTriggerSetupOut(0)
  32 {
  33   // Constructor.
  34 }
  35
  36 //________________________________________________________________________
  37 AliEmcalTriggerMaker::AliEmcalTriggerMaker(const char *name) :
  38   AliAnalysisTaskEmcal(name,kFALSE),
  39   fCaloTriggersOutName("EmcalTriggers"),
  40   fCaloTriggerSetupOutName("EmcalTriggersSetup"),
  41   fCaloTriggersOut(0),
  42   fCaloTriggerSetupOut(0)
  43 {
  44   // Constructor.
  45
  46 }
  47
  48 //________________________________________________________________________
  49 AliEmcalTriggerMaker::~AliEmcalTriggerMaker()
  50 {
  51   // Destructor.
  52 }
  53
  54 //________________________________________________________________________
  55 void AliEmcalTriggerMaker::ExecOnce()
  56 {
  57   // Init the analysis.
  58
  59   AliAnalysisTaskEmcal::ExecOnce();
  60
  61   if (!fInitialized)
  62     return;
  63
  64   if (!fCaloTriggersOutName.IsNull()) {
  65     fCaloTriggersOut = new TClonesArray("AliEmcalTriggerPatchInfo");
  66     fCaloTriggersOut->SetName(fCaloTriggersOutName);
  67
  68     if (!(InputEvent()->FindListObject(fCaloTriggersOutName))) {
  69       InputEvent()->AddObject(fCaloTriggersOut);
  70     }
  71     else {
  72       fInitialized = kFALSE;
  73       AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fCaloTriggersOutName.Data()));
  74       return;
  75     }
  76   }
  77
  78   if (!fCaloTriggerSetupOutName.IsNull()) {
  79     fCaloTriggerSetupOut = new AliEmcalTriggerSetupInfo();
  80     fCaloTriggerSetupOut->SetName(fCaloTriggerSetupOutName);
  81
  82     if (!(InputEvent()->FindListObject(fCaloTriggerSetupOutName))) {
  83       InputEvent()->AddObject(fCaloTriggerSetupOut);
  84     }
  85     else {
  86       fInitialized = kFALSE;
  87       AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fCaloTriggerSetupOutName.Data()));
  88       return;
  89     }
  90   }
  91
  92 }
  93
  94 //________________________________________________________________________
  95 Bool_t AliEmcalTriggerMaker::Run()
  96 {
  97   // Create the emcal particles
  98
  99   Int_t globCol, globRow, tBits, cellAbsId[4];
 100   Int_t absId, adcAmp;
 101   Int_t i, j, k, iMain, cmCol, cmRow, cmiCellCol, cmiCellRow;
 102   Int_t jetTrigger, iTriggers;
 103         Int_t patchADC[48][64];
 104   Double_t amp, ca, eMain, cmiCol, cmiRow;
 105
 106   TVector3 centerGeo, center1, center2, centerMass, edge1, edge2, vertex;
 107
 108   AliEmcalTriggerPatchInfo *trigger;
 109
 110   // delete patch array, clear setup object
 111   fCaloTriggersOut->Delete();
 112   fCaloTriggerSetupOut->Clean();
 113
 114   if( !fCaloTriggers ){
 115     AliError(Form("Calo triggers container %s not available.", fCaloTriggersName.Data()));
 116     return kTRUE;
 117   }
 118   if( !fCaloCells ){
 119     AliError(Form("Calo cells container %s not available.", fCaloCellsName.Data()));
 120     return kTRUE;
 121   }
 122
 123   // do not process, if sooner than 11h period
 124   if( InputEvent()->GetRunNumber() < 167693 )
 125     return kTRUE;
 126
 127   // do not process any MC, since no MC was generated with correct
 128   // EMCal trigger L1 jet trigger simulation, yet
 129   // productions will be enabled, once some correct once are produced
 130   if( MCEvent() != 0 )
 131     return kTRUE;
 132
 133   // must reset before usage, or the class will fail
 134   fCaloTriggers->Reset();
 135
 136   // first run over the patch array to compose a map of 2x2 patch energies
 137   // which is then needed to construct the full patch ADC energy
 138   // class is not empty
 139   if( fCaloTriggers->GetEntries() > 0 ){
 140
 141                 // zero the array
 142                 for( i = 0; i < 48; i++ )
 143                         for( j = 0; j < 64; j++ )
 144                                 patchADC[i][j] = 0;
 145
 146     // go throuth the trigger channels
 147     while( fCaloTriggers->Next() ){
 148       // get position in global 2x2 tower coordinates
 149       // A0 left bottom (0,0)
 150       fCaloTriggers->GetPosition( globCol, globRow );
 151
 152       // for some strange reason some ADC amps are initialized in reconstruction
 153       // as -1, neglect those :\\ wth
 154       fCaloTriggers->GetL1TimeSum( adcAmp );
 155                         if( adcAmp > -1 )
 156                                 patchADC[globCol][globRow] = adcAmp;
 157                 } // patches
 158         } // array not empty
 159
 160   // reset for re-run
 161   fCaloTriggers->Reset();
 162
 163   // dig out common data (thresholds)
 164   // 0 - jet high, 1 - gamma high, 2 - jet low, 3 - gamma low
 165   fCaloTriggerSetupOut->SetThresholds( fCaloTriggers->GetL1Threshold( 0 ),
 166                                        fCaloTriggers->GetL1Threshold( 1 ),
 167                                        fCaloTriggers->GetL1Threshold( 2 ),
 168                                        fCaloTriggers->GetL1Threshold( 3 ));
 169
 170   // class is not empty
 171   if( fCaloTriggers->GetEntries() > 0 ){
 172
 173     iTriggers = 0;
 174     iMain = -1;
 175     eMain = -1;
 176
 177     // save primary vertex in vector
 178     vertex.SetXYZ( fVertex[0], fVertex[1], fVertex[2] );
 179
 180     // go throuth the trigger channels
 181     while( fCaloTriggers->Next() ){
 182
 183       // check if jet trigger low or high
 184       fCaloTriggers->GetTriggerBits( tBits );
 185
 186       jetTrigger = 0;
 187       if(( tBits >> ( kTriggerTypeEnd + kL1JetLow )) & 1 )
 188         jetTrigger = 1;
 189       if(( tBits >> ( kTriggerTypeEnd + kL1JetHigh )) & 1)
 190         jetTrigger = jetTrigger | 2;
 191
 192       if( jetTrigger == 0 )
 193         continue;
 194
 195       // get position in global 2x2 tower coordinates
 196       // A0 left bottom (0,0)
 197       fCaloTriggers->GetPosition( globCol, globRow );
 198
 199       // get the absolute trigger ID
 200       fGeom->GetAbsFastORIndexFromPositionInEMCAL( globCol, globRow, absId );
 201       // convert to the 4 absId of the cells composing the trigger channel
 202       fGeom->GetCellIndexFromFastORIndex( absId, cellAbsId );
 203
 204       // get low left edge (eta max, phi min)
 205       fGeom->GetGlobal( cellAbsId[0], edge1 );
 206
 207       // sum the available energy in the 32/32 window of cells
 208       // step over trigger channels and get all the corresponding cells
 209       // make CM
 210       amp = 0;
 211       cmiCol = 0;
 212       cmiRow = 0;
 213       adcAmp = 0;
 214       for( i = 0; i < 16; i++ ){
 215         for( j = 0; j < 16; j++ ){
 216           // get the 4 cells composing the trigger channel
 217           fGeom->GetAbsFastORIndexFromPositionInEMCAL( globCol+i, globRow+j, absId );
 218           fGeom->GetCellIndexFromFastORIndex( absId, cellAbsId );
 219           // add amplitudes and find patch edges
 220           for( k = 0; k < 4; k++ ){
 221             ca = fCaloCells->GetCellAmplitude( cellAbsId[k] );
 222             amp += ca;
 223             cmiCol += ca*(Double_t)i;
 224             cmiRow += ca*(Double_t)j;
 225           }
 226           // add the STU ADCs in the patch
 227           adcAmp += patchADC[globCol+i][globRow+j];
 228         }
 229       } // 32x32 cell window
 230       if( amp == 0 ){
 231         AliDebug(2,"EMCal trigger patch with 0 energy.");
 232         continue;
 233       }
 234
 235       // get the CM and patch index
 236       cmiCol /= amp;
 237       cmiRow /= amp;
 238       cmCol = globCol + (Int_t)cmiCol;
 239       cmRow = globRow + (Int_t)cmiRow;
 240
 241       // get the patch and corresponding cells
 242       fGeom->GetAbsFastORIndexFromPositionInEMCAL( cmCol, cmRow, absId );
 243       fGeom->GetCellIndexFromFastORIndex( absId, cellAbsId );
 244
 245       // find which out of the 4 cells is closest to CM and get it's position
 246       cmiCellCol = TMath::Nint( cmiCol * 2. );
 247       cmiCellRow = TMath::Nint( cmiRow * 2. );
 248       fGeom->GetGlobal( cellAbsId[(cmiCellRow%2)*2 + cmiCellCol%2], centerMass );
 249
 250       // get up right edge (eta min, phi max)
 251       // get the absolute trigger ID
 252       fGeom->GetAbsFastORIndexFromPositionInEMCAL( globCol+15, globRow+15, absId );
 253       // convert to the 4 absId of the cells composing the trigger channel
 254       fGeom->GetCellIndexFromFastORIndex( absId, cellAbsId );
 255
 256       fGeom->GetGlobal( cellAbsId[3], edge2 );
 257
 258       // get the geometrical center as an average of two diagonally
 259       // adjacent patches in the center
 260       // picking two diagonally closest cells from the patches
 261       fGeom->GetAbsFastORIndexFromPositionInEMCAL( globCol+7, globRow+7, absId );
 262       fGeom->GetCellIndexFromFastORIndex( absId, cellAbsId );
 263       fGeom->GetGlobal( cellAbsId[3], center1 );
 264
 265       fGeom->GetAbsFastORIndexFromPositionInEMCAL( globCol+8, globRow+8, absId );
 266       fGeom->GetCellIndexFromFastORIndex( absId, cellAbsId );
 267       fGeom->GetGlobal( cellAbsId[0], center2 );
 268
 269       centerGeo = center1;
 270       centerGeo += center2;
 271       centerGeo *= 0.5;
 272
 273       // relate all to primary vertex
 274       centerGeo -= vertex;
 275       centerMass -= vertex;
 276       edge1 -= vertex;
 277       edge2 -= vertex;
 278
 279       // save the trigger object
 280       new ((*fCaloTriggersOut)[iTriggers])AliEmcalTriggerPatchInfo();
 281       trigger = (AliEmcalTriggerPatchInfo*)fCaloTriggersOut->At( iTriggers );
 282       iTriggers++;
 283
 284       trigger->SetCenterGeo( centerGeo, amp );
 285       trigger->SetCenterMass( centerMass, amp );
 286       trigger->SetEdge1( edge1, amp );
 287       trigger->SetEdge2( edge2, amp );
 288       trigger->SetADCAmp( adcAmp );
 289       trigger->SetTriggerBits( tBits );
 290       trigger->SetEdgeCell( globCol*2, globRow*2 ); // from triggers to cells
 291
 292       // check if more energetic than others for main patch marking
 293       if( eMain < amp ){
 294         eMain = amp;
 295         iMain = iTriggers - 1;
 296       }
 297
 298 //       cout << " pi:" << trigger->GetPhiMin() << " px:" << trigger->GetPhiMax();
 299 //       cout << " pg:" << trigger->GetPhiGeo() << " " << (trigger->GetPhiMin()+trigger->GetPhiMax()) / 2.;
 300 //       cout << " pc:" << trigger->GetPhiCM();
 301 //       cout << " ei:" << trigger->GetEtaMin() << " ex:" << trigger->GetEtaMax();
 302 //       cout << " eg:" << trigger->GetEtaGeo() << " " << (trigger->GetEtaMin()+trigger->GetEtaMax()) / 2.;
 303 //       cout << " ec:" << trigger->GetEtaCM();
 304 //       cout << " e:" << trigger->GetPatchE();
 305 //       cout << " jl:" << trigger->IsJetLow() << " jh:" << trigger->IsJetHigh() << endl;
 306
 307     } // triggers
 308
 309     // mark the most energetic patch as main
 310     if( iMain > -1 ){
 311       trigger = (AliEmcalTriggerPatchInfo*)fCaloTriggersOut->At( iMain );
 312       tBits = trigger->GetTriggerBits();
 313       // main trigger flag
 314       tBits = tBits | ( 1 << 24 );
 315       trigger->SetTriggerBits( tBits );
 316     }
 317
 318   } // there are some triggers
 319
 320   return kTRUE;
 321 }