PWGJE/EMCALJetTasks/AliEmcalPicoTrackInGridMaker.cxx

   1 // $Id$
   2 //
   3 // Class to put collection of tracks into grid of PicoTracks
   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 "AliPicoTrack.h"
  15 #include "AliVTrack.h"
  16 #include "AliEmcalJet.h"
  17 #include "AliParticleContainer.h"
  18 #include "AliJetContainer.h"
  19
  20 #include "AliEmcalPicoTrackInGridMaker.h"
  21
  22 ClassImp(AliEmcalPicoTrackInGridMaker)
  23
  24 //________________________________________________________________________
  25 AliEmcalPicoTrackInGridMaker::AliEmcalPicoTrackInGridMaker() :
  26 AliAnalysisTaskEmcalJet("AliEmcalPicoTrackInGridMaker",kTRUE),
  27   fTracksOutName("PicoTracksInGrid"),
  28   fTracksOut(0),
  29   fL1Slide(0),
  30   fCellSize(0.0145),
  31   fMinCellE(0.15),
  32   fExclLeadingPatch(0),
  33   fPatchSub(3),
  34   fRhoMean(184.),
  35   fNCells(-1),
  36   fNCellsEMCal(-1),
  37   fNCellsDCal(-1),
  38   fMultVsRho(0)
  39 {
  40   // Constructor.
  41
  42   for(Int_t i = 0; i<2; i++) {
  43     fCellGrid[i] = 0;
  44     fMiniPatchGrid[i] = 0;
  45     fActiveAreaMP[i] = 0;
  46   }
  47
  48   fPhiMin[0] = 1.405;
  49   fPhiMax[0] = 1.405+TMath::DegToRad()*110.;
  50   fPhiMin[1] = 4.547;
  51   fPhiMax[1] = 5.71;
  52   fEtaMin[0] = -0.7;
  53   fEtaMax[0] = 0.7;
  54   fEtaMin[1] = -0.7;
  55   fEtaMax[1] = 0.7;
  56
  57   for(Int_t i = 0; i<5; i++) {
  58     fNPatchesEMCal[i] = 0;
  59     fh1RhoEmcal[i] = 0;
  60     fh1RhoDcal[i] = 0;
  61     fPatchEnVsActivityEmcal[i] = 0;
  62     fPatchEnVsActivityDcal[i]  = 0;
  63
  64     for(Int_t j = 0; j<2; j++) {
  65       fPatchGrid[j][i] = 0;
  66       fActiveAreaMPP[j][i] = 0;
  67       fActiveAreaCP[j][i]  = 0;
  68       fPatchECorr[j][i] = 0;
  69       fPatchECorrPar[j][i] = 0;
  70       fPatchERaw[j][i] = 0;
  71       fPatchECorrRho[j][i] = 0;
  72       fPatchECorrECorrRho[j][i] = 0;
  73       fh2JetPtPatchECorr[j][i] = 0;
  74     }
  75     fh2PatchEtaPhiEmcal[i] = 0;
  76     fh2PatchEtaPhiDcal[i]  = 0;
  77   }
  78
  79   for(Int_t i = 0; i<3; i++) {
  80     fh2MedianTypeEmcal[i] = 0;
  81     fh2MedianTypeDcal[i] = 0;
  82     fpMedianTypeEmcal[i] = 0;
  83     fpMedianTypeDcal[i] = 0;
  84   }
  85   SetMakeGeneralHistograms(kTRUE);
  86 }
  87
  88 //________________________________________________________________________
  89 AliEmcalPicoTrackInGridMaker::AliEmcalPicoTrackInGridMaker(const char *name) :
  90   AliAnalysisTaskEmcalJet(name,kTRUE),
  91   fTracksOutName("PicoTracksInGrid"),
  92   fTracksOut(0),
  93   fL1Slide(0),
  94   fCellSize(0.0145),
  95   fMinCellE(0.15),
  96   fExclLeadingPatch(0),
  97   fPatchSub(3),
  98   fRhoMean(184.),
  99   fNCells(-1),
 100   fNCellsEMCal(-1),
 101   fNCellsDCal(-1),
 102   fMultVsRho(0)
 103 {
 104   // Constructor.
 105
 106   for(Int_t i = 0; i<2; i++) {
 107     fCellGrid[i] = 0;
 108     fMiniPatchGrid[i] = 0;
 109     fActiveAreaMP[i] = 0;
 110   }
 111
 112   fPhiMin[0] = 1.405;
 113   fPhiMax[0] = 1.405+TMath::DegToRad()*110.;
 114   fPhiMin[1] = 4.547;
 115   fPhiMax[1] = 5.71;
 116   fEtaMin[0] = -0.7;
 117   fEtaMax[0] = 0.7;
 118   fEtaMin[1] = -0.7;
 119   fEtaMax[1] = 0.7;
 120
 121   for(Int_t i = 0; i<5; i++) {
 122     fNPatchesEMCal[i] = 0;
 123     fh1RhoEmcal[i] = 0;
 124     fh1RhoDcal[i] = 0;
 125     fPatchEnVsActivityEmcal[i] = 0;
 126     fPatchEnVsActivityDcal[i]  = 0;
 127
 128     for(Int_t j = 0; j<2; j++) {
 129       fPatchGrid[j][i] = 0;
 130       fActiveAreaMPP[j][i] = 0;
 131       fActiveAreaCP[j][i]  = 0;
 132       fPatchECorr[j][i] = 0;
 133       fPatchECorrPar[j][i] = 0;
 134       fPatchERaw[j][i] = 0;
 135       fPatchECorrRho[j][i] = 0;
 136       fPatchECorrECorrRho[j][i] = 0;
 137       fh2JetPtPatchECorr[j][i] = 0;
 138     }
 139     fh2PatchEtaPhiEmcal[i] = 0;
 140     fh2PatchEtaPhiDcal[i]  = 0;
 141   }
 142
 143   for(Int_t i = 0; i<3; i++) {
 144     fh2MedianTypeEmcal[i] = 0;
 145     fh2MedianTypeDcal[i] = 0;
 146     fpMedianTypeEmcal[i] = 0;
 147     fpMedianTypeDcal[i] = 0;
 148   }
 149
 150   SetMakeGeneralHistograms(kTRUE);
 151 }
 152
 153 //________________________________________________________________________
 154 AliEmcalPicoTrackInGridMaker::~AliEmcalPicoTrackInGridMaker()
 155 {
 156   // Destructor.
 157 }
 158
 159 //________________________________________________________________________
 160 void AliEmcalPicoTrackInGridMaker::UserCreateOutputObjects()
 161 {
 162   // Create my user objects.
 163   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 164
 165   Int_t nBinsMed = 500;
 166   Double_t minMed = 0.;
 167   Double_t maxMed = 500.;
 168
 169   Int_t nBinsPhiEmcal = 132+64;
 170   Double_t phiMinEmcal = 1.436931 - 32.*fCellSize;
 171   Double_t phiMaxEmcal = 3.292931 + 32.*fCellSize;
 172
 173   Int_t nBinsPhiDcal = 80+64;
 174   Double_t phiMinDcal = 4.664500 - 32.*fCellSize;
 175   Double_t phiMaxDcal = 5.592500 + 32.*fCellSize;
 176
 177   Int_t nBinsEta = 96+64;
 178   Double_t etaMin = -0.696 - 32.*fCellSize;
 179   Double_t etaMax =  0.696 + 32.*fCellSize;
 180
 181   for(Int_t i = 0; i<3; i++) {
 182     fh2MedianTypeEmcal[i] = new TH2F(Form("fh2MedianTypeEmcalAreaType%d",i),Form("fh2MedianTypeEmcalAreaType%d",i),5,0.5,5.5,nBinsMed,minMed,maxMed);
 183     fOutput->Add(fh2MedianTypeEmcal[i]);
 184
 185     fh2MedianTypeDcal[i] = new TH2F(Form("fh2MedianTypeDcalAreaType%d",i),Form("fh2MedianTypeDcalAreaType%d",i),5,0.5,5.5,nBinsMed,minMed,maxMed);
 186     fOutput->Add(fh2MedianTypeDcal[i]);
 187
 188     fpMedianTypeEmcal[i] = new TProfile(Form("fpMedianTypeEmcalAreaType%d",i),Form("fpMedianTypeEmcalAreaType%d",i),5,0.5,5.5,"s");
 189     fOutput->Add(fpMedianTypeEmcal[i]);
 190
 191     fpMedianTypeDcal[i] = new TProfile(Form("fpMedianTypeDcalAreaType%d",i),Form("fpMedianTypeDcalAreaType%d",i),5,0.5,5.5,"s");
 192     fOutput->Add(fpMedianTypeDcal[i]);
 193   }
 194
 195   TString det[2] = {"Emcal","Dcal"};
 196   for(Int_t i = 0; i<5; i++) { //loop over patch types
 197     fh1RhoEmcal[i] = new TH1F(Form("fh1RhoEmcal_%d",i),Form("fh1RhoEmcal_%d",i),500,0.,1000.);
 198     fOutput->Add(fh1RhoEmcal[i]);
 199     fh1RhoDcal[i] = new TH1F(Form("fh1RhoDcal_%d",i),Form("fh1RhoDcal_%d",i),500,0.,1000.);
 200     fOutput->Add(fh1RhoDcal[i]);
 201
 202     fPatchEnVsActivityEmcal[i] = new TH2F(Form("fh2PatchEnVsActivityEmcal_%d",i),Form("fh2PatchEnVsActivityEmcal_%d",i),300,0.,300.,150,-0.5,149.5);
 203     fOutput->Add(fPatchEnVsActivityEmcal[i]);
 204
 205     fPatchEnVsActivityDcal[i] = new TH2F(Form("fh2PatchEnVsActivityDcal_%d",i),Form("fh2PatchEnVsActivityDcal_%d",i),300,0.,300.,150,-0.5,149.5);
 206     fOutput->Add(fPatchEnVsActivityDcal[i]);
 207
 208     for(Int_t j = 0; j<2; j++) {
 209       fPatchECorr[j][i] = new TH1F(Form("fPatchECorr%s_%d",det[j].Data(),i),Form("fPatchECorr%s_%d;#it{E}_{patch}^{corr}",det[j].Data(),i),250,-50.,200.);
 210       fOutput->Add(fPatchECorr[j][i]);
 211
 212       fPatchECorrPar[j][i] = new TH1F(Form("fPatchECorrPar%s_%d",det[j].Data(),i),Form("fPatchECorrPar%s_%d;#it{E}_{patch}^{corr}",det[j].Data(),i),250,-50.,200.);
 213       fOutput->Add(fPatchECorrPar[j][i]);
 214
 215       fPatchERaw[j][i] = new TH1F(Form("fPatchERaw%s_%d",det[j].Data(),i),Form("fPatchERaw%s_%d;#it{E}_{patch}^{corr}",det[j].Data(),i),250,-50.,200.);
 216       fOutput->Add(fPatchERaw[j][i]);
 217
 218       fPatchECorrRho[j][i] = new TH2F(Form("fPatchECorrRho%s_%d",det[j].Data(),i),Form("fPatchECorrRho%s_%d;#it{E}_{patch}^{corr};#rho",det[j].Data(),i),250,-50.,200.,500,0.,500.);
 219       fOutput->Add(fPatchECorrRho[j][i]);
 220
 221       fPatchECorrECorrRho[j][i] = new TH3F(Form("fPatchECorrECorrRho%s_%d",det[j].Data(),i),Form("fPatchECorrECorrRho%s_%d;#it{E}_{patch,det1}^{corr};#it{E}_{patch,det2}^{corr};#rho",det[j].Data(),i),210,-30.,180.,210,-30.,180.,250,0.,250.);
 222       fOutput->Add(fPatchECorrECorrRho[j][i]);
 223
 224       fh2JetPtPatchECorr[j][i] = new TH2F(Form("fh2JetPtPatchECorr%s_%d",det[j].Data(),i),Form("fh2JetPtPatchECorr%s_%d",det[j].Data(),i),250,-50.,200.,250,-50.,200.);
 225       fOutput->Add(fh2JetPtPatchECorr[j][i]);
 226     }
 227
 228     fh2PatchEtaPhiEmcal[i] = new TH2F(Form("fh2PatchEtaPhiEmcal_%d",i),Form("fh2PatchEtaPhiEmcal_%d;#eta;#phi",i),nBinsEta,etaMin,etaMax,nBinsPhiEmcal,phiMinEmcal,phiMaxEmcal);
 229     fOutput->Add(fh2PatchEtaPhiEmcal[i]);
 230
 231     fh2PatchEtaPhiDcal[i] = new TH2F(Form("fh2PatchEtaPhiDcal_%d",i),Form("fh2PatchEtaPhiDcal_%d;#eta;#phi",i),nBinsEta,etaMin,etaMax,nBinsPhiDcal,phiMinDcal,phiMaxDcal);
 232     fOutput->Add(fh2PatchEtaPhiDcal[i]);
 233   }
 234
 235   fMultVsRho = new TH2F("fMultVsRho","fMultVsRho",3000,0,3000,400,0,400);
 236   fOutput->Add(fMultVsRho);
 237
 238   PostData(1, fOutput);
 239 }
 240
 241 //________________________________________________________________________
 242 Bool_t AliEmcalPicoTrackInGridMaker::Run()
 243 {
 244   // Main loop, called for each event.
 245
 246   Bool_t b = CreateGridCells();
 247   if(!b) return kFALSE;
 248   b = CreateGridMiniPatches();
 249   if(!b) return kFALSE;
 250
 251   //L0 single shower trigger
 252   CreateGridPatches(4,0);
 253   //  return kTRUE;
 254
 255   //L1 triggers: sliding window
 256   CreateGridPatches(4,1);
 257   CreateGridPatches(8,1);
 258   CreateGridPatches(16,1);
 259   CreateGridPatches(32,1);
 260
 261
 262   Double_t medL0 = CalculateMedian(0,0);
 263   fh2MedianTypeEmcal[0]->Fill(0.5,medL0);
 264   fpMedianTypeEmcal[0]->Fill(0.5,medL0);
 265   medL0 = CalculateMedian(0,1);
 266   fh2MedianTypeDcal[0]->Fill(0.5,medL0);
 267   fpMedianTypeDcal[0]->Fill(0.5,medL0);
 268
 269   Double_t medL1[4][2];
 270   for(Int_t i = 0; i<4; i++) { //patches
 271     for(Int_t type = 0; type<2; type++) { //EMCal or DCal
 272       for(Int_t areaT = 0; areaT<1; areaT++) { //areay type (passive vs active)
 273         medL1[i][type] = CalculateMedian(i+1,type,areaT);
 274         if(type==0) {
 275           fh2MedianTypeEmcal[areaT]->Fill((Double_t)(i+1)+0.5,medL1[i][type]);
 276           fpMedianTypeEmcal[areaT]->Fill((Double_t)(i+1)+0.5,medL1[i][type]);
 277         }
 278         if(type==1) {
 279           fh2MedianTypeDcal[areaT]->Fill((Double_t)(i+1)+0.5,medL1[i][type]);
 280           fpMedianTypeDcal[areaT]->Fill((Double_t)(i+1)+0.5,medL1[i][type]);
 281         }
 282       }
 283     }
 284   }
 285
 286   // subtract energy density and store energy distributions of corrected patches in histo
 287   Int_t EleadID[5][2];
 288   Double_t Elead[5][2];
 289   Double_t EleadRaw[5][2];
 290   for(Int_t i = 1; i<5; i++) { //patch types
 291     Int_t stepSize = GetTriggerPatchIdStepSizeNoOverlap(GetPatchDim(i));
 292     for(Int_t type = 0; type<2; type++) {
 293       EleadID[i][type] = -1;
 294       Elead[i][type] = -1e6;
 295       EleadRaw[i][type] = -1e6;
 296       Int_t subType = 1;
 297       if(type==1) subType = 0;
 298       //      for(Int_t j = 0; j<(fPatchGrid[type][i].GetSize()-stepSize+1); j+=stepSize) { //patches
 299       for(Int_t k = 0; k<GetNColTriggerPatches(type,GetPatchDim(i),i); k+=stepSize) {
 300         for(Int_t l = 0; l<GetNRowTriggerPatches(type,GetPatchDim(i),i); l+=stepSize) {
 301           Int_t id = GetTriggerPatchID(k,l,type,GetPatchDim(i),i);
 302           //      if(type==1 && i==4) Printf("id: %d/%d k: %d/%d l: %d/%d",id,fPatchGrid[type][i].GetSize(),k,GetNColTriggerPatches(type,GetPatchDim(i),i),l,GetNRowTriggerPatches(type,GetPatchDim(i),i));
 303           if(fPatchGrid[type][i].At(id)>0.) { //don't do anything with empty patches
 304             Double_t sub = medL1[fPatchSub-1][subType]*GetPatchArea(i);
 305             fPatchECorr[type][i]->Fill(fPatchGrid[type][i].At(id) - sub);
 306             fPatchECorrPar[type][i]->Fill(fPatchGrid[type][i].At(id) - fRhoMean*GetPatchArea(i));
 307
 308             //Bookkeep leading patches
 309             if((fPatchGrid[type][i].At(id)-sub)>Elead[i][type]) {
 310               EleadID[i][type] = id;
 311               Elead[i][type] = fPatchGrid[type][i].At(id)-sub;
 312             }
 313             if(fPatchGrid[type][i].At(id)>EleadRaw[i][type])
 314               EleadRaw[i][type] = fPatchGrid[type][i].At(id);
 315           }
 316         }//cols
 317       }//rows
 318     }//type
 319
 320     AliJetContainer *cont = GetJetContainer(0);
 321     for(Int_t k = 0; k<2; k++) {
 322       Int_t subType = 1;
 323       if(k==1) subType=0;
 324       fPatchECorrRho[k][i]->Fill(Elead[i][k],medL1[fPatchSub-1][subType]);
 325       fPatchECorrECorrRho[k][i]->Fill(Elead[i][k],Elead[i][subType],medL1[fPatchSub-1][subType]);
 326       fPatchERaw[k][i]->Fill(EleadRaw[i][k]);
 327       //Save jet spectra within EMCal and DCal fiducial acceptance
 328       //jet pT vs highest energy patch for preferred trigger patch types
 329       if(cont) {
 330         Double_t r = cont->GetJetRadius();
 331         cont->SetJetEtaLimits(fEtaMin[k]+r,fEtaMax[k]-r);
 332         cont->SetJetPhiLimits(fPhiMin[k]+r,fPhiMax[k]-r);
 333         Double_t rho = cont->GetRhoVal();
 334         AliEmcalJet *jet = NULL;
 335         cont->ResetCurrentID();
 336         while((jet = cont->GetNextAcceptJet())) {
 337           Double_t jetpt = jet->Pt() - rho*jet->Area();
 338           fh2JetPtPatchECorr[k][i]->Fill(jetpt,Elead[i][k]);
 339         }
 340       }
 341     }
 342     Double_t eta = 0.; Double_t phi = 0.;
 343     GetEtaPhiFromTriggerPatchID(EleadID[i][0],0,GetPatchDim(i),1,eta,phi);
 344     fh2PatchEtaPhiEmcal[i]->Fill(eta,phi);
 345
 346     GetEtaPhiFromTriggerPatchID(EleadID[i][1],1,GetPatchDim(i),1,eta,phi);
 347     fh2PatchEtaPhiDcal[i]->Fill(eta,phi);
 348   } //patch types
 349
 350   fMultVsRho->Fill(GetParticleContainer(0)->GetNParticles(),medL1[3][0]);
 351   return kTRUE;
 352 }
 353
 354 //________________________________________________________________________
 355 AliEmcalJet* AliEmcalPicoTrackInGridMaker::GetClosestJet(const Double_t eta, const Double_t phi, const Int_t icont) const {
 356
 357   AliJetContainer *cont = GetJetContainer(icont);
 358   if(!cont) return NULL;
 359
 360   Double_t closest_dr = 1e6;
 361   Int_t closest_id = -1;
 362   AliEmcalJet *jet = NULL;
 363   cont->ResetCurrentID();
 364   while((jet = cont->GetNextAcceptJet())) {
 365     Double_t dphi = jet->Phi() - phi;
 366     Double_t deta = jet->Eta() - eta;
 367     dphi = TVector2::Phi_mpi_pi(dphi);
 368     Double_t dr = TMath::Sqrt ( dphi * dphi + deta * deta );
 369     if(dr<closest_dr) {
 370       closest_dr = dr;
 371       closest_id = cont->GetCurrentID();
 372     }
 373   }
 374   jet = cont->GetJet(closest_id);
 375   return jet;
 376 }
 377
 378 //________________________________________________________________________
 379 Double_t AliEmcalPicoTrackInGridMaker::CalculateSum(const Int_t patchType) const {
 380   //calc total energy of all patches
 381   Int_t n = fPatchGrid[0][patchType].GetSize();
 382   if(n<1) return -1.;
 383
 384   Double_t sum = 0.;
 385   Int_t count = 0;
 386   Int_t stepSize = GetTriggerPatchIdStepSizeNoOverlap(GetPatchDim(patchType));
 387   for(Int_t type = 0; type<2; type++) {
 388     for(Int_t i = 0; i<fPatchGrid[type][patchType].GetSize(); i+=stepSize) {
 389       if(fPatchGrid[type][patchType].At(i)>0.) count++;
 390       sum+=fPatchGrid[type][patchType].At(i);
 391     }
 392   }
 393   return sum;
 394 }
 395
 396 //________________________________________________________________________
 397 Double_t AliEmcalPicoTrackInGridMaker::CalculateMedian(const Int_t patchType, const Int_t type, const Int_t areaType) {
 398   //areaType:
 399   //0: passive area
 400   //1: active area mini patches
 401   //2: active arear cells
 402
 403   Int_t n = fPatchGrid[type][patchType].GetSize();
 404   if(n<1) return -1.;
 405
 406   Int_t level = 0;
 407   if(patchType>0) level = 1;
 408   Int_t dim = GetPatchDim(patchType);
 409   Double_t area = GetPatchArea(patchType);
 410   Int_t stepSize = GetTriggerPatchIdStepSizeNoOverlap(GetPatchDim(patchType),level);
 411   //Printf("patchType: %d dim: %d stepSizeNoOverlap: %d ",patchType,GetPatchDim(patchType),stepSize);
 412
 413   static Double_t arr[999];
 414   Int_t c = 0;
 415
 416   //find patch with highest energy
 417   Int_t imax = -1;
 418   Double_t max = 0.;
 419   for(Int_t i = 0; i<(GetNColTriggerPatches(type,dim,patchType)); i+=stepSize) {
 420     for(Int_t j = 0; j<(GetNRowTriggerPatches(type,dim,patchType)); j+=stepSize) {
 421       Int_t id = GetTriggerPatchID(i,j,type,dim,patchType);
 422       //      Printf("id: %d/%d i: %d/%d j:%d/%d",id,fPatchGrid[type][patchType].GetSize(),i,GetNColTriggerPatches(type,dim,patchType),j,GetNRowTriggerPatches(type,dim,patchType));
 423       if(fPatchGrid[type][patchType].At(id)>max) {
 424         imax = id;
 425         max = fPatchGrid[type][patchType].At(id);
 426       }
 427     }//cols
 428   }//rows
 429
 430   for(Int_t i = 0; i<GetNColTriggerPatches(type,dim,patchType); i+=stepSize) {
 431     for(Int_t j = 0; j<GetNRowTriggerPatches(type,dim,patchType); j+=stepSize) {
 432       Int_t id = GetTriggerPatchID(i,j,type,dim,patchType);
 433       if(fExclLeadingPatch>0 && id==imax) continue;
 434       if(fPatchGrid[type][patchType].At(id)>0.) {
 435         Int_t active = 99;
 436         if(areaType==1) active = fActiveAreaMPP[type][patchType].At(id);
 437         else if(areaType==2) active = fActiveAreaCP[type][patchType].At(id);
 438         if(areaType>0) area = GetPatchAreaActive(id,type,patchType,areaType-1);
 439
 440         if(area>0. && active>1) {
 441           arr[c] = fPatchGrid[type][patchType].At(id)/area;
 442           c++;
 443         }
 444         if(type==0 && areaType==0) {
 445           fh1RhoEmcal[patchType]->Fill(arr[c-1]);
 446           fPatchEnVsActivityEmcal[patchType]->Fill(fPatchGrid[type][patchType].At(id),fActiveAreaCP[type][patchType].At(id));
 447         }
 448         if(type==1 && areaType==0) {
 449           fh1RhoDcal[patchType]->Fill(arr[c-1]);
 450           fPatchEnVsActivityDcal[patchType]->Fill(fPatchGrid[type][patchType].At(id),fActiveAreaCP[type][patchType].At(id));
 451         }
 452       }
 453     }
 454   }
 455   Double_t med = TMath::Median(c,arr);
 456   return med;
 457 }
 458
 459 //________________________________________________________________________
 460 Bool_t AliEmcalPicoTrackInGridMaker::CreateGridCells() {
 461   //create cells from track input
 462   if(!InitCells()) return kFALSE;
 463
 464   AliVParticle *track = NULL;
 465   AliParticleContainer *trackCont = GetParticleContainer(0);
 466   if(!trackCont) return kFALSE;
 467   trackCont->ResetCurrentID();
 468   while((track = trackCont->GetNextAcceptParticle())) {
 469     if(track->Pt()<fMinCellE) continue;
 470     Int_t id = GetGridID(track);
 471     Int_t type = GetCellType(track);
 472     if(id>-1)
 473       fCellGrid[type].AddAt(fCellGrid[type].At(id)+track->Pt(),id);
 474     }
 475   return kTRUE;
 476 }
 477
 478 //________________________________________________________________________
 479 Bool_t AliEmcalPicoTrackInGridMaker::CreateGridMiniPatches() {
 480   //create mini patches (2x2 cells)
 481   if(!InitMiniPatches()) return kFALSE;
 482
 483   //loop over edges of mini patches
 484   Int_t nm = 0; //mini patch number
 485   for(Int_t type = 0; type<2; type++) {
 486     nm = 0;
 487     for(Int_t i = 0; i<(GetNCellsCol(type)-1); i+=2) {
 488       for(Int_t j = 0; j<(GetNCellsRow(type)-1); j+=2) {
 489         //loop over cells in mini patch
 490         for(Int_t k = 0; k<2; k++) {
 491           for(Int_t l = 0; l<2; l++) {
 492             Int_t id = GetGridID(i+k,j+l,type);
 493             fMiniPatchGrid[type].AddAt(fMiniPatchGrid[type].At(nm)+fCellGrid[type].At(id),nm);
 494             if(fCellGrid[type].At(id)>0.)
 495               fActiveAreaMP[type].AddAt(fActiveAreaMP[type].At(nm)+1,nm);
 496           }
 497         }
 498         nm++;
 499       }
 500     }
 501   }
 502   return kTRUE;
 503 }
 504
 505 //________________________________________________________________________
 506 Int_t AliEmcalPicoTrackInGridMaker::GetNRowMiniPatches(const Int_t type) const {
 507   //returns number of rows of mini patches in detector of type X (0: EMCal 1: DCal)
 508   Int_t nRows = TMath::FloorNint(0.5*GetNCellsCol(type));
 509   return nRows;
 510 }
 511
 512 //________________________________________________________________________
 513 Int_t AliEmcalPicoTrackInGridMaker::GetNColMiniPatches(const Int_t type) const {
 514   //returns number of rows of mini patches in detector of type X (0: EMCal 1: DCal)
 515   Int_t nCols = TMath::FloorNint(0.5*GetNCellsRow(type));
 516   return nCols;
 517 }
 518
 519 //________________________________________________________________________
 520 Bool_t AliEmcalPicoTrackInGridMaker::CreateGridPatches(const Int_t dim, const Int_t level) {
 521   //create trigger patches
 522   if(!InitPatches(dim,level)) return kFALSE;
 523
 524   Int_t pt = GetPatchType(dim,level);
 525   Int_t nm = (Int_t)(dim/2.);    //size of trigger patch in number of mini patches
 526   Int_t stepm = (Int_t)(dim/2.); //step size through grid in mini patches
 527   if(level==1 && fL1Slide)  stepm = GetSlidingStepSizeMiniPatches(dim,level);
 528   //loop over edges of mini patches
 529   for(Int_t type = 0; type<2; type++) {
 530     Int_t np = 0; //patch number
 531       for(Int_t j = 0; j<=(GetNColMiniPatches(type)-nm); j+=stepm) {
 532     for(Int_t i = 0; i<=(GetNRowMiniPatches(type)-nm); i+=stepm) {
 533       //      for(Int_t j = 0; j<=(GetNColMiniPatches(type)-nm); j+=stepm) {
 534         //loop over mini patches in patch
 535         for(Int_t k = 0; k<nm; k++) {
 536           for(Int_t l = 0; l<nm; l++) {
 537             Int_t row = i+k;
 538             Int_t col = j+l;
 539             Int_t id = GetMiniPatchID(row,col,type);
 540             fPatchGrid[type][pt].AddAt(fPatchGrid[type][pt].At(np)+fMiniPatchGrid[type].At(id),np);
 541             if(fMiniPatchGrid[type].At(id)>0.) {
 542               fActiveAreaMPP[type][pt].AddAt(fActiveAreaMPP[type][pt].At(np)+1,np);
 543               fActiveAreaCP[type][pt].AddAt(fActiveAreaCP[type][pt].At(np)+fActiveAreaMP[type].At(id),np);
 544             }
 545           }
 546         }
 547         np++;
 548       }
 549     }
 550   }
 551   return kTRUE;
 552 }
 553
 554 //________________________________________________________________________
 555 Int_t AliEmcalPicoTrackInGridMaker::GetTriggerPatchID(const Int_t row, const Int_t col, const Int_t type, const Int_t dim, const Int_t patchType) const {
 556   Int_t id = row*GetNRowTriggerPatches(type,dim,patchType) + col;
 557   return id;
 558 }
 559
 560 //________________________________________________________________________
 561 Int_t AliEmcalPicoTrackInGridMaker::GetMiniPatchID(const Int_t row, const Int_t col, const Int_t type) const {
 562   Int_t id = row*GetNColMiniPatches(type) + col;
 563   return id;
 564 }
 565
 566 //________________________________________________________________________
 567 Int_t AliEmcalPicoTrackInGridMaker::GetCellType(const Double_t eta, const Double_t phi) const {
 568   //cell in EMCal (0) or DCal (1)
 569   for(Int_t i = 0; i<2; i++) {
 570     if(eta>fEtaMin[i] && eta<fEtaMax[i] && phi>fPhiMin[i] && phi<fPhiMax[i]) return i;
 571   }
 572   return -1;
 573 }
 574
 575 //________________________________________________________________________
 576 Int_t AliEmcalPicoTrackInGridMaker::GetGridID(const Int_t row, const Int_t col, const Int_t type) const {
 577   Int_t id = row*GetNCellsRow(type) + col;
 578   return id;
 579 }
 580
 581 //________________________________________________________________________
 582 Int_t AliEmcalPicoTrackInGridMaker::GetGridID(const Double_t eta, const Double_t phi) const {
 583
 584   Int_t type = GetCellType(eta,phi);
 585   if(type<0 || type>1) return -1; //position is not in EMCal or DCal
 586
 587   // grid ID convention:
 588   // upper left corner (min phi, min eta) is first ID
 589   // then walk through grid from upper left to lower right accross the rows in phi
 590
 591   Int_t id = -1;
 592   Double_t etaRel = eta-fEtaMin[type];
 593   Double_t phiRel = phi-fPhiMin[type];
 594   Int_t row = TMath::FloorNint(etaRel/fCellSize);
 595   Int_t col = TMath::FloorNint(phiRel/fCellSize);
 596   id = GetGridID(row,col,type);
 597
 598   if(id>=fNCells) {
 599     Printf("Got too large id %d %d type: %d",id,fNCells,type);
 600     Printf("eta: %f phi: %f",eta,phi);
 601     Printf("etaRel: %f phiRel: %f",etaRel,phiRel);
 602     Printf("row: %d col: %d  ->  %d + %d = %d",row,col,row*GetNCellsRow(type) + col,fNCellsEMCal,id);
 603     Printf("n cells row: %d",GetNCellsRow(type));
 604     Printf("\n");
 605   }
 606   return id;
 607 }
 608
 609 //________________________________________________________________________
 610 void AliEmcalPicoTrackInGridMaker::GetEtaPhiFromGridID(const Int_t id, const Int_t type, Double_t &eta, Double_t &phi) const {
 611   //returns eta phi of cell at lower right edge (lowest eta, lowest phi)
 612   Int_t row = TMath::FloorNint(id/GetNCellsRow(type));
 613   Int_t col = id - row * GetNCellsRow(type);
 614   eta = fEtaMin[type] + row*fCellSize;
 615   phi = fPhiMin[type] + col*fCellSize;
 616   AliDebug(2,Form("id: %d type: %d row: %d col: %d eta: %f phi: %f",id,type,row,col,eta,phi));
 617 }
 618
 619 //________________________________________________________________________
 620 void AliEmcalPicoTrackInGridMaker::GetEtaPhiFromMiniPatchID(const Int_t id, const Int_t type, Double_t &eta, Double_t &phi) const {
 621   //returns eta phi of mini patch at lower right edge (lowest eta, lowest phi)
 622   Int_t row = TMath::FloorNint(id/GetNColMiniPatches(type));
 623   Int_t col = id - row * GetNColMiniPatches(type);
 624   eta = fEtaMin[type] + row*2.*fCellSize;
 625   phi = fPhiMin[type] + col*2.*fCellSize;
 626 }
 627
 628 //________________________________________________________________________
 629 void AliEmcalPicoTrackInGridMaker::GetEtaPhiFromTriggerPatchID(const Int_t id, const Int_t type, const Int_t dim, const Int_t level, Double_t &eta, Double_t &phi) const {
 630   //returns eta phi of mini patch at lower right edge (lowest eta, lowest phi)
 631   Int_t step = GetSlidingStepSizeCells(dim,level); //id: 8/96 k(row): 0/8 l(col): 8/12
 632   //  Int_t offset = dim/step;
 633   Int_t row = TMath::FloorNint(id/GetNColTriggerPatches(type,dim,GetPatchType(dim,level)));
 634   Int_t col = id - row * GetNColTriggerPatches(type,dim,GetPatchType(dim,level));
 635   eta = fEtaMin[type] + row*step*fCellSize;
 636   phi = fPhiMin[type] + col*step*fCellSize;
 637   // if(dim==32) {
 638   //   Printf("dim: %d step: %d offset: %d",dim,step,offset);
 639   //   Printf("dim: %d id: %d row: %d col: %d eta: %f phi: %f",dim,id,row,col,eta,phi);
 640   // }
 641 }
 642
 643 //________________________________________________________________________
 644 Int_t AliEmcalPicoTrackInGridMaker::GetNColTriggerPatches(const Int_t type, const Int_t dim, const Int_t patchType) const {
 645   //returns number of trigger patch columns
 646   Int_t level = 0;
 647   if(patchType>0) level = 1;
 648   Double_t stepmp = (Double_t)GetSlidingStepSizeMiniPatches(dim,level);
 649   Int_t nmp = GetNColMiniPatches(type);
 650   Int_t ntc = TMath::FloorNint(nmp/stepmp);
 651   //  Printf("dim: %d stepmp: %f nmp: %d ntc: %d",dim,stepmp,nmp,ntc);
 652   return ntc;
 653 }
 654
 655 //________________________________________________________________________
 656 Int_t AliEmcalPicoTrackInGridMaker::GetNRowTriggerPatches(const Int_t type, const Int_t dim, const Int_t patchType) const {
 657   //returns number of trigger patch rows
 658   Int_t level = 0;
 659   if(patchType>0) level = 1;
 660   Double_t stepmp = (Double_t)GetSlidingStepSizeMiniPatches(dim,level);
 661   Int_t nmp = GetNRowMiniPatches(type);
 662   Int_t ntr = TMath::FloorNint(nmp/stepmp);
 663   return ntr;
 664 }
 665
 666 //________________________________________________________________________
 667 Int_t AliEmcalPicoTrackInGridMaker::GetNCellsCol(const Int_t type) const {
 668   //returns number of cells in column
 669   Double_t deta = fEtaMax[type] - fEtaMin[type];
 670   Int_t nCellsCol = TMath::FloorNint(deta/fCellSize);
 671   return nCellsCol;
 672 }
 673
 674 //________________________________________________________________________
 675 Int_t AliEmcalPicoTrackInGridMaker::GetNCellsRow(const Int_t type) const {
 676   //returns number of cells in row
 677   Double_t dPhi = fPhiMax[type] - fPhiMin[type];
 678   Int_t nCellsRow = TMath::FloorNint(dPhi/fCellSize);
 679   return nCellsRow;
 680 }
 681
 682 //________________________________________________________________________
 683 Bool_t AliEmcalPicoTrackInGridMaker::InitCells() {
 684   //initialize cells array
 685   CheckEdges();
 686   if(!CheckEdges()) return kFALSE;
 687
 688   //number of cells in EMCal acceptance
 689   Int_t nCellsPhiE = GetNCellsCol(0);
 690   Int_t nCellsEtaE = GetNCellsRow(0);
 691   fNCellsEMCal = nCellsPhiE*nCellsEtaE;
 692
 693   //number of cells in DCal acceptance
 694   Int_t nCellsPhiD = GetNCellsCol(1);
 695   Int_t nCellsEtaD = GetNCellsRow(1);
 696   fNCellsDCal = nCellsPhiD*nCellsEtaD;
 697
 698   //total number of cells
 699   fNCells = fNCellsEMCal + fNCellsDCal;
 700
 701   AliDebug(2,Form("EMCal: %d x %d",nCellsEtaE,nCellsPhiE));
 702   AliDebug(2,Form("DCal: row: %d x col: %d",nCellsEtaD,nCellsPhiD));
 703   AliDebug(2,Form("fNCells: %d fNCellsE: %d fnCellsD: %d",fNCells,fNCellsEMCal,fNCellsDCal));
 704
 705   fCellGrid[0].Set(fNCellsEMCal);
 706   fCellGrid[1].Set(fNCellsDCal);
 707   fCellGrid[0].Reset(0);
 708   fCellGrid[1].Reset(0);
 709   return kTRUE;
 710 }
 711
 712 //________________________________________________________________________
 713 Bool_t AliEmcalPicoTrackInGridMaker::InitMiniPatches() {
 714   //initialize mini patch array
 715   if(fCellGrid[0].GetSize()<0) return kFALSE;
 716   if(fCellGrid[1].GetSize()<0) return kFALSE;
 717   Double_t conv = 0.25; //dimension of mini patch is 2x2 cells
 718   Int_t nMiniPatches[2];
 719   nMiniPatches[0] = (Int_t)(conv*fNCellsEMCal);
 720   nMiniPatches[1] = (Int_t)(conv*fNCellsDCal);
 721   for(Int_t i = 0; i<2; i++) {
 722     fMiniPatchGrid[i].Set(nMiniPatches[i]);
 723     fMiniPatchGrid[i].Reset(0.);
 724
 725     fActiveAreaMP[i].Set(nMiniPatches[i]);
 726     fActiveAreaMP[i].Reset(0);
 727   }
 728   return kTRUE;
 729 }
 730
 731 //________________________________________________________________________
 732 Int_t AliEmcalPicoTrackInGridMaker::GetNTriggerPatches(const Int_t type, const Int_t dim, const Int_t level) const {
 733   //get number of trigger patches in EMCal or DCal
 734   Double_t dimd = (Double_t)dim;
 735   Double_t conv = 1./(dimd*dimd);
 736   if(level==1) {
 737     Double_t step = (Double_t)GetSlidingStepSizeCells(dim);
 738     conv = 1./(step*step);
 739   }
 740   Int_t nPatches = 0;
 741   if(type==0) nPatches = (Int_t)(conv*fNCellsEMCal);
 742   else if(type==1) nPatches = (Int_t)(conv*fNCellsDCal);
 743   return nPatches;
 744 }
 745
 746 //________________________________________________________________________
 747 Bool_t AliEmcalPicoTrackInGridMaker::InitPatches(const Int_t dim, const Int_t level) {
 748   // dimensions in cell units
 749   // if level==1: sliding window will be applied
 750   // L1 4x4: slide by 2 cells: 1 mini patch
 751   // L1 8x8: slide by 4 cells: 2 mini patches
 752   // L1 16x16: slide by 4 cells: 2 mini patches
 753   // L1 32x32: slide by 8 cells: 4 mini patches
 754
 755   if(fCellGrid[0].GetSize()<0) return kFALSE;
 756   if(fCellGrid[1].GetSize()<0) return kFALSE;
 757
 758   Int_t type = GetPatchType(dim,level);
 759   if(type<0 || type>4) return kFALSE;
 760
 761   Int_t nPatches[2]; //number of trigger patches in EMCal and DCal
 762   for(Int_t i = 0; i<2; i++)
 763     nPatches[i] = GetNTriggerPatches(i,dim,level);
 764   //total number of trigger patches
 765   Int_t nPatchesT = nPatches[0] + nPatches[1];
 766
 767   fNPatchesEMCal[type] = nPatches[0];
 768   AliDebug(2,Form("Create trigger patch of type %d with dim %d and with %d patches EMCAL: %d DCAL: %d",type,dim,nPatchesT,nPatches[0],nPatches[1]));
 769   //Printf("Create trigger patch of type %d with dim %d and with %d patches EMCAL: %d DCAL: %d",type,dim,nPatchesT,nPatches[0],nPatches[1]);
 770   for(Int_t i = 0; i<2; i++) {
 771     fPatchGrid[i][type].Set(nPatches[i]);
 772     fPatchGrid[i][type].Reset(0);
 773     fActiveAreaMPP[i][type].Set(nPatches[i]);
 774     fActiveAreaMPP[i][type].Reset(0);
 775     fActiveAreaCP[i][type].Set(nPatches[i]);
 776     fActiveAreaCP[i][type].Reset(0);
 777   }
 778
 779   return kTRUE;
 780 }
 781
 782 //________________________________________________________________________
 783 Int_t AliEmcalPicoTrackInGridMaker::GetPatchDim(const Int_t ipatch) const {
 784   //returns total area of patch
 785   Int_t ncell = 4;
 786   if(ipatch==0) ncell = 4;
 787   if(ipatch==1) ncell = 4;
 788   if(ipatch==2) ncell = 8;
 789   if(ipatch==3) ncell = 16;
 790   if(ipatch==4) ncell = 32;
 791
 792   return ncell;
 793 }
 794
 795 //________________________________________________________________________
 796 Double_t AliEmcalPicoTrackInGridMaker::GetPatchArea(const Int_t ipatch) const {
 797   //returns total area of patch
 798   Double_t ncell = (Double_t)GetPatchDim(ipatch);
 799   Double_t area = ncell*ncell*fCellSize*fCellSize;
 800   return area;
 801 }
 802
 803 //________________________________________________________________________
 804 Double_t AliEmcalPicoTrackInGridMaker::GetPatchAreaActive(const Int_t id, const Int_t type, const Int_t ipatch, const Int_t atype) const {
 805   //atype = 0 : active area from mini patches
 806   //atype = 1 : active area from cells
 807
 808   Int_t active = 0;
 809   if(atype==0) active = fActiveAreaMPP[type][ipatch].At(id);
 810   else if(atype==1) active = fActiveAreaCP[type][ipatch].At(id);
 811   else return -1;
 812
 813   Double_t fac = 1.;
 814   if(atype==0) fac = 4.;
 815   Double_t area = active*fac*fCellSize*fCellSize;
 816   return area;
 817 }
 818
 819 //________________________________________________________________________
 820 Int_t AliEmcalPicoTrackInGridMaker::GetSlidingStepSizeCells(const Int_t dim, const Int_t level) const {
 821   //get step size for mock-up L1 trigger
 822   if(!fL1Slide || level==0) return dim;
 823
 824   if(dim==4) return 2;
 825   else if(dim==8) return 4;
 826   else if(dim==16) return 4;
 827   else if(dim==32) return 8;
 828   else return -1;
 829 }
 830
 831 //________________________________________________________________________
 832 Int_t AliEmcalPicoTrackInGridMaker::GetSlidingStepSizeMiniPatches(const Int_t dim, const Int_t level) const {
 833   //returns step size in mini patches
 834   Double_t step = (Double_t)GetSlidingStepSizeCells(dim,level);
 835   if(step<0) return step;
 836   Int_t stepMiniPatch = (Int_t)(step/2.);
 837   return stepMiniPatch;
 838 }
 839
 840 //________________________________________________________________________
 841 Int_t AliEmcalPicoTrackInGridMaker::GetTriggerPatchIdStepSizeNoOverlap(const Int_t dim, const Int_t level) const {
 842   //return step for trigger patch id's without overlapping patches
 843   if(!fL1Slide) return 1;
 844
 845   Int_t cellStep = GetSlidingStepSizeCells(dim,level);
 846   Int_t step = TMath::FloorNint((Double_t)(dim/cellStep));
 847   return step;
 848 }
 849
 850 //________________________________________________________________________
 851 Int_t AliEmcalPicoTrackInGridMaker::GetPatchType(const Int_t dim, const Int_t level) const {
 852   //type of trigger patch (size)
 853   Int_t type = -1;
 854   if(level==0 && dim==4) type = 0;
 855   if(level==1) {
 856     if(dim==4)  type = 1;
 857     if(dim==8)  type = 2;
 858     if(dim==16) type = 3;
 859     if(dim==32) type = 4;
 860   }
 861   return type;
 862 }
 863
 864 //________________________________________________________________________
 865 Bool_t AliEmcalPicoTrackInGridMaker::CheckEdges() {
 866   //Check if defined edges of EMCal and DCal make sense
 867   if(fPhiMin[0]<0. || fPhiMax[0]<fPhiMin[0]) {
 868     AliDebug(11,Form("EMCal phi edges not defined %f-%f",fPhiMin[0],fPhiMax[0]));
 869     return kFALSE;
 870   }
 871
 872   if(fPhiMin[1]<0. || fPhiMax[1]<fPhiMin[1]) {
 873     AliDebug(11,Form("DCal phi edges not defined %f-%f",fPhiMin[1],fPhiMax[1]));
 874     return kFALSE;
 875   }
 876
 877   if(fEtaMin[0]<-10. || fEtaMax[0]<fEtaMin[1]) {
 878     AliDebug(11,Form("EMCal eta edges not well defined %f-%f",fEtaMin[0],fEtaMax[0]));
 879     return kFALSE;
 880   }
 881
 882   if(fEtaMin[1]<-10. || fEtaMax[1]<fEtaMin[1]) {
 883     AliDebug(11,Form("DCal eta edges not well defined %f-%f",fEtaMin[1],fEtaMax[1]));
 884     return kFALSE;
 885   }
 886
 887   for(Int_t type = 0; type<2; type++) {
 888     Double_t dphi = fPhiMax[type] - fPhiMin[type];
 889     Double_t nPatchPhi32 = TMath::Floor(dphi/(32.*fCellSize));
 890     Double_t nCellsColExact = nPatchPhi32 * 32. ;
 891
 892     Double_t deta = fEtaMax[type] - fEtaMin[type];
 893     Double_t nPatchEta32 = TMath::Floor(deta/(32.*fCellSize));
 894     Double_t nCellsRowExact = nPatchEta32 * 32. ;
 895
 896     Double_t col_extra = dphi/fCellSize - TMath::Floor(nCellsColExact);
 897     Double_t phi_extra = col_extra*fCellSize;
 898     fPhiMin[type] += phi_extra*0.5;
 899     fPhiMax[type] -= phi_extra*0.5;
 900
 901     Double_t row_extra = deta/fCellSize - TMath::Floor(nCellsRowExact);
 902     Double_t eta_extra = row_extra*fCellSize;
 903     fEtaMin[type] += eta_extra*0.5;
 904     fEtaMax[type] -= eta_extra*0.5;
 905     AliDebug(2,Form("type: %d  exact: col: %f row: %f",type,nCellsColExact,nCellsRowExact));
 906     //    Printf("type: %d  exact: col: %f row: %f",type,nCellsColExact,nCellsRowExact);
 907     //  PrintAcceptance();
 908   }
 909   return kTRUE;
 910 }
 911
 912 //________________________________________________________________________
 913 void AliEmcalPicoTrackInGridMaker::PrintAcceptance() const {
 914   Printf("EMCal");
 915   Printf("phi: %f-%f",fPhiMin[0],fPhiMax[0]);
 916   Printf("eta: %f-%f",fEtaMin[0],fEtaMax[0]);
 917
 918   Printf("DCal");
 919   Printf("phi: %f-%f",fPhiMin[1],fPhiMax[1]);
 920   Printf("eta: %f-%f",fEtaMin[1],fEtaMax[1]);
 921 }