PWG/EMCAL/AliEmcalPatchFromCellMaker.cxx

   1 // $Id$
   2 //
   3 // Class to put cells into trigger patches
   4 //
   5 // Author: M. Verweij
   6
   7 #include <TClonesArray.h>
   8 #include <TRandom3.h>
   9 #include <TProfile.h>
  10 #include <TH3F.h>
  11
  12 #include "AliAnalysisManager.h"
  13 #include "AliLog.h"
  14 #include "AliEMCALGeometry.h"
  15 #include "AliEmcalTriggerPatchInfo.h"
  16
  17 #include "AliEmcalPatchFromCellMaker.h"
  18
  19 ClassImp(AliEmcalPatchFromCellMaker)
  20
  21 //________________________________________________________________________
  22 AliEmcalPatchFromCellMaker::AliEmcalPatchFromCellMaker() :
  23   AliAnalysisTaskEmcal("AliEmcalPatchFromCellMaker",kTRUE),
  24   fCaloTriggersOutName("EmcalPatches32x32"),
  25   fCaloTriggersOut(0),
  26   fPatchDim(32),
  27   fMinCellE(0.05),
  28   fCellTimeMin(485e-9),
  29   fCellTimeMax(685e-9),
  30   fL1Slide(0),
  31   fh3EEtaPhiCell(0),
  32   fh2CellEnergyVsTime(0),
  33   fh1CellEnergySum(0)
  34 {
  35   // Constructor.
  36  for (Int_t i = 0; i < kPatchCols; i++) {
  37     for (Int_t j = 0; j < kPatchRows; j++) {
  38       fPatchADCSimple[i][j] = 0.;
  39       fPatchESimple[i][j] = 0.;
  40     }
  41  }
  42
  43   SetMakeGeneralHistograms(kTRUE);
  44 }
  45
  46 //________________________________________________________________________
  47 AliEmcalPatchFromCellMaker::AliEmcalPatchFromCellMaker(const char *name) :
  48   AliAnalysisTaskEmcal(name,kTRUE),
  49   fCaloTriggersOutName("EmcalPatches32x32"),
  50   fCaloTriggersOut(0),
  51   fPatchDim(32),
  52   fMinCellE(0.05),
  53   fCellTimeMin(485e-9),
  54   fCellTimeMax(685e-9),
  55   fL1Slide(0),
  56   fh3EEtaPhiCell(0),
  57   fh2CellEnergyVsTime(0),
  58   fh1CellEnergySum(0)
  59 {
  60   // Constructor.
  61  for (Int_t i = 0; i < kPatchCols; i++) {
  62     for (Int_t j = 0; j < kPatchRows; j++) {
  63       fPatchADCSimple[i][j] = 0.;
  64       fPatchESimple[i][j] = 0.;
  65     }
  66  }
  67
  68   SetMakeGeneralHistograms(kTRUE);
  69 }
  70
  71 //________________________________________________________________________
  72 AliEmcalPatchFromCellMaker::~AliEmcalPatchFromCellMaker()
  73 {
  74   // Destructor.
  75 }
  76
  77 //________________________________________________________________________
  78 void AliEmcalPatchFromCellMaker::ExecOnce()
  79 {
  80   // Init the analysis.
  81
  82   AliAnalysisTaskEmcal::ExecOnce();
  83
  84   if (!fInitialized)
  85     return;
  86
  87   if (!fCaloTriggersOutName.IsNull()) {
  88     fCaloTriggersOut = new TClonesArray("AliEmcalTriggerPatchInfo");
  89     fCaloTriggersOut->SetName(fCaloTriggersOutName);
  90
  91     if (!(InputEvent()->FindListObject(fCaloTriggersOutName))) {
  92       InputEvent()->AddObject(fCaloTriggersOut);
  93     }
  94     else {
  95       fInitialized = kFALSE;
  96       AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fCaloTriggersOutName.Data()));
  97       return;
  98     }
  99   }
 100
 101 }
 102
 103 //________________________________________________________________________
 104 void AliEmcalPatchFromCellMaker::UserCreateOutputObjects()
 105 {
 106   // Create user output.
 107
 108   AliAnalysisTaskEmcal::UserCreateOutputObjects();
 109
 110   Int_t fgkNPhiBins = 18*8;
 111   Float_t kMinPhi   = 0.;
 112   Float_t kMaxPhi   = 2.*TMath::Pi();
 113   Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
 114   for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
 115
 116   Int_t fgkNEtaBins = 100;
 117   Float_t fgkEtaMin = -1.;
 118   Float_t fgkEtaMax =  1.;
 119   Double_t *binsEta=new Double_t[fgkNEtaBins+1];
 120   for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
 121
 122   Int_t fgkNTimeBins = 600;
 123   Float_t kMinTime   = -200.;
 124   Float_t kMaxTime   = 1000;
 125   Double_t *binsTime = new Double_t[fgkNTimeBins+1];
 126   for(Int_t i=0; i<=fgkNTimeBins; i++) binsTime[i]=(Double_t)kMinTime + (kMaxTime-kMinTime)/fgkNTimeBins*(Double_t)i ;
 127
 128   Double_t enBinEdges[3][2];
 129   enBinEdges[0][0] = 1.; //10 bins
 130   enBinEdges[0][1] = 0.1;
 131   enBinEdges[1][0] = 5.; //8 bins
 132   enBinEdges[1][1] = 0.5;
 133   enBinEdges[2][0] = 100.;//95 bins
 134   enBinEdges[2][1] = 1.;
 135
 136   const Float_t enmin1 =  0;
 137   const Float_t enmax1 =  enBinEdges[0][0];
 138   const Float_t enmin2 =  enmax1 ;
 139   const Float_t enmax2 =  enBinEdges[1][0];
 140   const Float_t enmin3 =  enmax2 ;
 141   const Float_t enmax3 =  enBinEdges[2][0];//fgkEnMax;
 142   const Int_t nbin11 = (int)((enmax1-enmin1)/enBinEdges[0][1]);
 143   const Int_t nbin12 = (int)((enmax2-enmin2)/enBinEdges[1][1])+nbin11;
 144   const Int_t nbin13 = (int)((enmax3-enmin3)/enBinEdges[2][1])+nbin12;
 145
 146   Int_t fgkNEnBins=nbin13;
 147   Double_t *binsEn=new Double_t[fgkNEnBins+1];
 148   for(Int_t i=0; i<=fgkNEnBins; i++) {
 149     if(i<=nbin11) binsEn[i]=(Double_t)enmin1 + (enmax1-enmin1)/nbin11*(Double_t)i ;
 150     if(i<=nbin12 && i>nbin11) binsEn[i]=(Double_t)enmin2 + (enmax2-enmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
 151     if(i<=nbin13 && i>nbin12) binsEn[i]=(Double_t)enmin3 + (enmax3-enmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
 152   }
 153
 154   fh3EEtaPhiCell = new TH3F("fh3EEtaPhiCell","fh3EEtaPhiCell;E_{cell};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
 155   fOutput->Add(fh3EEtaPhiCell);
 156
 157   fh2CellEnergyVsTime = new TH2F("fh2CellEnergyVsTime","fh2CellEnergyVsTime;E_{cell};time",fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
 158   fOutput->Add(fh2CellEnergyVsTime);
 159
 160   fh1CellEnergySum = new TH1F("fh1CellEnergySum","fh1CellEnergySum;E_{cell};time",fgkNEnBins,binsEn);
 161   fOutput->Add(fh1CellEnergySum);
 162
 163   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 164
 165   if(binsEn)                delete [] binsEn;
 166   if(binsPhi)               delete [] binsPhi;
 167   if(binsEta)               delete [] binsEta;
 168   if(binsTime)              delete [] binsTime;
 169 }
 170
 171 //________________________________________________________________________
 172 Bool_t AliEmcalPatchFromCellMaker::Run()
 173 {
 174   // Main loop, called for each event.
 175
 176   fCaloTriggersOut->Delete();
 177
 178   if (!fCaloCells) {
 179     AliError(Form("Calo cells container %s not available.", fCaloCellsName.Data()));
 180     return kFALSE;
 181   }
 182
 183  for (Int_t i = 0; i < kPatchCols; i++) {
 184     for (Int_t j = 0; j < kPatchRows; j++) {
 185       fPatchADCSimple[i][j] = 0.;
 186       fPatchESimple[i][j] = 0.;
 187     }
 188  }
 189
 190   if(!FillPatchADCSimple()) {
 191     AliError(Form("%s Could not create simple ADC patches",GetName()));
 192     return kFALSE;
 193   }
 194
 195   RunSimpleOfflineTrigger();
 196
 197   Double_t sum = 0.;
 198   for (Int_t i = 0; i < kPatchCols; i++) {
 199     for (Int_t j = 0; j < kPatchRows; j++) {
 200       sum+=fPatchESimple[i][j];
 201     }
 202   }
 203
 204   return kTRUE;
 205 }
 206
 207 //________________________________________________________________________
 208 Bool_t AliEmcalPatchFromCellMaker::FillPatchADCSimple()
 209 {
 210
 211   // fill the array for offline trigger processing
 212
 213   // fill the patch ADCs from cells
 214   Double_t sum = 0.;
 215   Int_t nCell = fCaloCells->GetNumberOfCells();
 216   for(Int_t iCell = 0; iCell < nCell; ++iCell) {
 217     // get the cell info, based in index in array
 218     Short_t cellId = fCaloCells->GetCellNumber(iCell);
 219
 220     Double_t cellT = fCaloCells->GetCellTime(cellId);
 221     Double_t amp = fCaloCells->GetAmplitude(iCell);
 222     fh2CellEnergyVsTime->Fill(amp,cellT*1e9);
 223
 224     //timing cuts
 225     if(cellT<fCellTimeMin || cellT>fCellTimeMax) continue;
 226     //energy cut
 227     if(amp<fMinCellE) continue;
 228     sum+=amp;
 229
 230     // get position
 231     Int_t absId=-1;
 232     fGeom->GetFastORIndexFromCellIndex(cellId, absId);
 233     Int_t globCol=-1, globRow=-1;
 234     fGeom->GetPositionInEMCALFromAbsFastORIndex(absId, globCol, globRow);
 235     // add
 236     fPatchADCSimple[globCol][globRow] += amp/kEMCL1ADCtoGeV;
 237     fPatchESimple[globCol][globRow] += amp;
 238
 239     TVector3 pos;
 240     fGeom->GetGlobal(cellId, pos);
 241     TLorentzVector lv(pos,amp);
 242     Double_t cellEta = lv.Eta();
 243     Double_t cellPhi = lv.Phi();
 244     if(cellPhi<0.) cellPhi+=TMath::TwoPi();
 245     if(cellPhi>TMath::TwoPi()) cellPhi-=TMath::TwoPi();
 246     fh3EEtaPhiCell->Fill(amp,cellEta,cellPhi);
 247   }
 248   fh1CellEnergySum->Fill(sum);
 249
 250   return kTRUE;
 251 }
 252
 253 //________________________________________________________________________
 254 void AliEmcalPatchFromCellMaker::RunSimpleOfflineTrigger()
 255 {
 256   // Runs a simple offline trigger algorithm.
 257   // It creates separate patches with dimension fPatchDim
 258
 259   // run the trigger algo, stepping by stepsize (in trigger tower units)
 260   Int_t itrig = 0;
 261   Int_t patchSize = GetDimFastor();
 262   Int_t stepSize = GetSlidingStepSizeFastor();
 263   Int_t maxCol = kPatchCols - patchSize;
 264   Int_t maxRow = kPatchRows - patchSize;
 265
 266   for (Int_t i = 0; i <= maxCol; i += stepSize) {
 267     for (Int_t j = 0; j <= maxRow; j += stepSize) {
 268       // get the trigger towers composing the patch
 269       Int_t   adcAmp = 0;
 270       Double_t enAmp = 0.;
 271       // window
 272       for (Int_t k = 0; k < patchSize; ++k) {
 273         for (Int_t l = 0; l < patchSize; ++l) {
 274           // add amplitudes
 275           adcAmp += (ULong64_t)fPatchADCSimple[i+k][j+l];
 276           enAmp += fPatchESimple[i+k][j+l];
 277         }
 278       }
 279
 280       if (adcAmp == 0) {
 281         AliDebug(2,"EMCal trigger patch with 0 ADC counts.");
 282         continue;
 283       }
 284
 285       Int_t absId=-1;
 286       Int_t cellAbsId[4]={-1,-1,-1,-1};
 287
 288       // get low left edge (eta max, phi min)
 289       fGeom->GetAbsFastORIndexFromPositionInEMCAL(i, j, absId);
 290       // convert to the 4 absId of the cells composing the trigger channel
 291       fGeom->GetCellIndexFromFastORIndex(absId, cellAbsId);
 292       TVector3 edge1;
 293       fGeom->GetGlobal(cellAbsId[0], edge1);
 294
 295       // get up right edge (eta min, phi max)
 296       fGeom->GetAbsFastORIndexFromPositionInEMCAL(i+patchSize-1, j+patchSize-1, absId);
 297       fGeom->GetCellIndexFromFastORIndex(absId, cellAbsId);
 298       TVector3 edge2;
 299       fGeom->GetGlobal(cellAbsId[3], edge2);
 300
 301       // get the center of the patch
 302       Int_t offsetCenter = TMath::FloorNint(0.5*patchSize);
 303       fGeom->GetAbsFastORIndexFromPositionInEMCAL(i+offsetCenter-1, j+offsetCenter-1, absId);
 304       fGeom->GetCellIndexFromFastORIndex(absId, cellAbsId);
 305       TVector3 center1;
 306       fGeom->GetGlobal(cellAbsId[3], center1);
 307
 308       fGeom->GetAbsFastORIndexFromPositionInEMCAL(i+offsetCenter, j+offsetCenter, absId);
 309       fGeom->GetCellIndexFromFastORIndex(absId, cellAbsId);
 310       TVector3 center2;
 311       fGeom->GetGlobal(cellAbsId[0], center2);
 312
 313       TVector3 centerGeo(center1);
 314       centerGeo += center2;
 315       centerGeo *= 0.5;
 316
 317       // save the trigger object
 318       AliEmcalTriggerPatchInfo *trigger =
 319         new ((*fCaloTriggersOut)[itrig]) AliEmcalTriggerPatchInfo();
 320       itrig++;
 321       trigger->SetCenterGeo(centerGeo, enAmp);
 322       trigger->SetEdge1(edge1, enAmp);
 323       trigger->SetEdge2(edge2, enAmp);
 324       trigger->SetADCAmp(adcAmp);
 325       trigger->SetEdgeCell(i*2, j*2); // from triggers to cells
 326     }
 327   } // trigger algo
 328   AliDebug(2,Form("Created %d trigger patches (%d) in this event",itrig,patchSize));
 329
 330 }
 331
 332 //________________________________________________________________________
 333 Int_t AliEmcalPatchFromCellMaker::GetDimFastor() const {
 334
 335   Int_t dim = TMath::FloorNint((Double_t)(fPatchDim/2.));
 336   return dim;
 337 }
 338
 339 //________________________________________________________________________
 340 Int_t AliEmcalPatchFromCellMaker::GetSlidingStepSizeFastor() const {
 341
 342   Int_t dim = GetDimFastor();
 343   if(!fL1Slide) return dim;
 344
 345   if(dim==2) return 2;
 346   else if(dim==4) return 4;
 347   else if(dim==8) return 4;
 348   else if(dim==16) return 8;
 349   else return -1;
 350 }
 351
 352
 353
 354