classes for trigger studies
authormverweij <marta.verweij@cern.ch>
Mon, 27 Oct 2014 20:30:48 +0000 (21:30 +0100)
committermverweij <marta.verweij@cern.ch>
Mon, 27 Oct 2014 21:15:36 +0000 (22:15 +0100)
fix some warnings

move class to correct lib (PatchFromCellMaker to PWG/EMCAL)

add event selection

13 files changed:
PWG/CMakelibPWGEMCAL.pkg
PWG/EMCAL/AliEmcalPatchFromCellMaker.cxx [new file with mode: 0644]
PWG/EMCAL/AliEmcalPatchFromCellMaker.h [new file with mode: 0644]
PWG/EMCAL/macros/AddTaskEmcalPatchFromCellMaker.C [new file with mode: 0644]
PWG/PWGEMCALLinkDef.h
PWGJE/CMakelibPWGJEEMCALJetTasks.pkg
PWGJE/EMCALJetTasks/AliEmcalPicoTrackInGridMaker.cxx [new file with mode: 0644]
PWGJE/EMCALJetTasks/AliEmcalPicoTrackInGridMaker.h [new file with mode: 0644]
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalHighMultTrigger.cxx [new file with mode: 0644]
PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalHighMultTrigger.h [new file with mode: 0644]
PWGJE/EMCALJetTasks/macros/AddTaskEmcalHighMultTrigger.C [new file with mode: 0644]
PWGJE/EMCALJetTasks/macros/AddTaskGridMaker.C [new file with mode: 0644]
PWGJE/PWGJEEMCALJetTasksLinkDef.h

index 2dd3b27..a83339a 100644 (file)
@@ -42,6 +42,7 @@ set ( SRCS
  EMCAL/AliEmcalMCTrackSelector.cxx
  EMCAL/AliEmcalParticle.cxx
  EMCAL/AliEmcalParticleMaker.cxx
+ EMCAL/AliEmcalPatchFromCellMaker.cxx
  EMCAL/AliEmcalPhysicsSelection.cxx
  EMCAL/AliEmcalPhysicsSelectionTask.cxx
  EMCAL/AliEmcalPicoTrackMaker.cxx
diff --git a/PWG/EMCAL/AliEmcalPatchFromCellMaker.cxx b/PWG/EMCAL/AliEmcalPatchFromCellMaker.cxx
new file mode 100644 (file)
index 0000000..ad5f9e4
--- /dev/null
@@ -0,0 +1,354 @@
+// $Id$
+//
+// Class to put cells into trigger patches
+//
+// Author: M. Verweij
+
+#include <TClonesArray.h>
+#include <TRandom3.h>
+#include <TProfile.h>
+#include <TH3F.h>
+
+#include "AliAnalysisManager.h"
+#include "AliLog.h"
+#include "AliEMCALGeometry.h"
+#include "AliEmcalTriggerPatchInfo.h"
+
+#include "AliEmcalPatchFromCellMaker.h"
+
+ClassImp(AliEmcalPatchFromCellMaker)
+
+//________________________________________________________________________
+AliEmcalPatchFromCellMaker::AliEmcalPatchFromCellMaker() : 
+  AliAnalysisTaskEmcal("AliEmcalPatchFromCellMaker",kTRUE),
+  fCaloTriggersOutName("EmcalPatches32x32"),
+  fCaloTriggersOut(0),
+  fPatchDim(32),
+  fMinCellE(0.15),
+  fCellTimeMin(485e-9),
+  fCellTimeMax(685e-9),
+  fL1Slide(0),
+  fh3EEtaPhiCell(0),
+  fh2CellEnergyVsTime(0),
+  fh1CellEnergySum(0)
+{
+  // Constructor.
+ for (Int_t i = 0; i < kPatchCols; i++) {
+    for (Int_t j = 0; j < kPatchRows; j++) {
+      fPatchADCSimple[i][j] = 0.;
+      fPatchESimple[i][j] = 0.;
+    }
+ }
+
+  SetMakeGeneralHistograms(kTRUE);
+}
+
+//________________________________________________________________________
+AliEmcalPatchFromCellMaker::AliEmcalPatchFromCellMaker(const char *name) : 
+  AliAnalysisTaskEmcal(name,kTRUE),
+  fCaloTriggersOutName("EmcalPatches32x32"),
+  fCaloTriggersOut(0),
+  fPatchDim(32),
+  fMinCellE(0.15),
+  fCellTimeMin(485e-9),
+  fCellTimeMax(685e-9),
+  fL1Slide(0),
+  fh3EEtaPhiCell(0),
+  fh2CellEnergyVsTime(0),
+  fh1CellEnergySum(0)
+{
+  // Constructor.
+ for (Int_t i = 0; i < kPatchCols; i++) {
+    for (Int_t j = 0; j < kPatchRows; j++) {
+      fPatchADCSimple[i][j] = 0.;
+      fPatchESimple[i][j] = 0.;
+    }
+ }
+
+  SetMakeGeneralHistograms(kTRUE);
+}
+
+//________________________________________________________________________
+AliEmcalPatchFromCellMaker::~AliEmcalPatchFromCellMaker()
+{
+  // Destructor.
+}
+
+//________________________________________________________________________
+void AliEmcalPatchFromCellMaker::ExecOnce() 
+{
+  // Init the analysis.
+
+  AliAnalysisTaskEmcal::ExecOnce();
+
+  if (!fInitialized)
+    return;
+
+  if (!fCaloTriggersOutName.IsNull()) {
+    fCaloTriggersOut = new TClonesArray("AliEmcalTriggerPatchInfo");
+    fCaloTriggersOut->SetName(fCaloTriggersOutName);
+
+    if (!(InputEvent()->FindListObject(fCaloTriggersOutName))) {
+      InputEvent()->AddObject(fCaloTriggersOut);
+    }
+    else {
+      fInitialized = kFALSE;
+      AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fCaloTriggersOutName.Data()));
+      return;
+    }
+  }
+
+}
+
+//________________________________________________________________________
+void AliEmcalPatchFromCellMaker::UserCreateOutputObjects()
+{
+  // Create user output.
+
+  AliAnalysisTaskEmcal::UserCreateOutputObjects();
+
+  Int_t fgkNPhiBins = 18*8;
+  Float_t kMinPhi   = 0.;
+  Float_t kMaxPhi   = 2.*TMath::Pi();
+  Double_t *binsPhi = new Double_t[fgkNPhiBins+1];
+  for(Int_t i=0; i<=fgkNPhiBins; i++) binsPhi[i]=(Double_t)kMinPhi + (kMaxPhi-kMinPhi)/fgkNPhiBins*(Double_t)i ;
+
+  Int_t fgkNEtaBins = 100;
+  Float_t fgkEtaMin = -1.;
+  Float_t fgkEtaMax =  1.;
+  Double_t *binsEta=new Double_t[fgkNEtaBins+1];
+  for(Int_t i=0; i<=fgkNEtaBins; i++) binsEta[i]=(Double_t)fgkEtaMin + (fgkEtaMax-fgkEtaMin)/fgkNEtaBins*(Double_t)i ;
+
+  Int_t fgkNTimeBins = 600;
+  Float_t kMinTime   = -200.;
+  Float_t kMaxTime   = 1000;
+  Double_t *binsTime = new Double_t[fgkNTimeBins+1];
+  for(Int_t i=0; i<=fgkNTimeBins; i++) binsTime[i]=(Double_t)kMinTime + (kMaxTime-kMinTime)/fgkNTimeBins*(Double_t)i ;
+
+  Double_t enBinEdges[3][2];
+  enBinEdges[0][0] = 1.; //10 bins
+  enBinEdges[0][1] = 0.1;
+  enBinEdges[1][0] = 5.; //8 bins
+  enBinEdges[1][1] = 0.5;
+  enBinEdges[2][0] = 100.;//95 bins
+  enBinEdges[2][1] = 1.;
+
+  const Float_t enmin1 =  0;
+  const Float_t enmax1 =  enBinEdges[0][0];
+  const Float_t enmin2 =  enmax1 ;
+  const Float_t enmax2 =  enBinEdges[1][0];
+  const Float_t enmin3 =  enmax2 ;
+  const Float_t enmax3 =  enBinEdges[2][0];//fgkEnMax;
+  const Int_t nbin11 = (int)((enmax1-enmin1)/enBinEdges[0][1]);
+  const Int_t nbin12 = (int)((enmax2-enmin2)/enBinEdges[1][1])+nbin11;
+  const Int_t nbin13 = (int)((enmax3-enmin3)/enBinEdges[2][1])+nbin12;
+
+  Int_t fgkNEnBins=nbin13;
+  Double_t *binsEn=new Double_t[fgkNEnBins+1];
+  for(Int_t i=0; i<=fgkNEnBins; i++) {
+    if(i<=nbin11) binsEn[i]=(Double_t)enmin1 + (enmax1-enmin1)/nbin11*(Double_t)i ;
+    if(i<=nbin12 && i>nbin11) binsEn[i]=(Double_t)enmin2 + (enmax2-enmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
+    if(i<=nbin13 && i>nbin12) binsEn[i]=(Double_t)enmin3 + (enmax3-enmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
+  }
+
+  fh3EEtaPhiCell = new TH3F("fh3EEtaPhiCell","fh3EEtaPhiCell;E_{cell};#eta;#phi",fgkNEnBins,binsEn,fgkNEtaBins,binsEta,fgkNPhiBins,binsPhi);
+  fOutput->Add(fh3EEtaPhiCell);
+
+  fh2CellEnergyVsTime = new TH2F("fh2CellEnergyVsTime","fh2CellEnergyVsTime;E_{cell};time",fgkNEnBins,binsEn,fgkNTimeBins,binsTime);
+  fOutput->Add(fh2CellEnergyVsTime);
+
+  fh1CellEnergySum = new TH1F("fh1CellEnergySum","fh1CellEnergySum;E_{cell};time",fgkNEnBins,binsEn);
+  fOutput->Add(fh1CellEnergySum);
+
+  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
+
+  if(binsEn)                delete [] binsEn;
+  if(binsPhi)               delete [] binsPhi;
+  if(binsEta)               delete [] binsEta;
+  if(binsTime)              delete [] binsTime;
+}
+
+//________________________________________________________________________
+Bool_t AliEmcalPatchFromCellMaker::Run() 
+{
+  // Main loop, called for each event.
+  
+  fCaloTriggersOut->Delete();
+
+  if (!fCaloCells) {
+    AliError(Form("Calo cells container %s not available.", fCaloCellsName.Data()));
+    return kFALSE;
+  }
+
+ for (Int_t i = 0; i < kPatchCols; i++) {
+    for (Int_t j = 0; j < kPatchRows; j++) {
+      fPatchADCSimple[i][j] = 0.;
+      fPatchESimple[i][j] = 0.;
+    }
+ }
+
+  if(!FillPatchADCSimple()) {
+    AliError(Form("%s Could not create simple ADC patches",GetName()));
+    return kFALSE;
+  }
+
+  RunSimpleOfflineTrigger();
+
+  Double_t sum = 0.;
+  for (Int_t i = 0; i < kPatchCols; i++) {
+    for (Int_t j = 0; j < kPatchRows; j++) {
+      sum+=fPatchESimple[i][j];
+    }
+  }
+
+  return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliEmcalPatchFromCellMaker::FillPatchADCSimple()
+{
+
+  // fill the array for offline trigger processing
+
+  // fill the patch ADCs from cells
+  Double_t sum = 0.;
+  Int_t nCell = fCaloCells->GetNumberOfCells();
+  for(Int_t iCell = 0; iCell < nCell; ++iCell) {
+    // get the cell info, based in index in array
+    Short_t cellId = fCaloCells->GetCellNumber(iCell);
+
+    Double_t cellT = fCaloCells->GetCellTime(cellId); 
+    Double_t amp = fCaloCells->GetAmplitude(iCell);
+    fh2CellEnergyVsTime->Fill(amp,cellT*1e9);
+
+    //timing cuts
+    if(cellT<fCellTimeMin || cellT>fCellTimeMax) continue;
+    //energy cut
+    if(amp<fMinCellE) continue;
+    sum+=amp;
+
+    // get position
+    Int_t absId=-1;
+    fGeom->GetFastORIndexFromCellIndex(cellId, absId);
+    Int_t globCol=-1, globRow=-1;
+    fGeom->GetPositionInEMCALFromAbsFastORIndex(absId, globCol, globRow);
+    // add
+    fPatchADCSimple[globCol][globRow] += amp/kEMCL1ADCtoGeV;
+    fPatchESimple[globCol][globRow] += amp;
+
+    TVector3 pos;
+    fGeom->GetGlobal(cellId, pos);
+    TLorentzVector lv(pos,amp);
+    Double_t cellEta = lv.Eta();
+    Double_t cellPhi = lv.Phi();
+    if(cellPhi<0.) cellPhi+=TMath::TwoPi();
+    if(cellPhi>TMath::TwoPi()) cellPhi-=TMath::TwoPi();
+    fh3EEtaPhiCell->Fill(amp,cellEta,cellPhi); 
+  }
+  fh1CellEnergySum->Fill(sum);
+
+  return kTRUE;
+}
+
+//________________________________________________________________________
+void AliEmcalPatchFromCellMaker::RunSimpleOfflineTrigger()
+{
+  // Runs a simple offline trigger algorithm.
+  // It creates separate patches with dimension fPatchDim
+
+  // run the trigger algo, stepping by stepsize (in trigger tower units)
+  Int_t itrig = 0;
+  Int_t patchSize = GetDimFastor();
+  Int_t stepSize = GetSlidingStepSizeFastor();
+  Int_t maxCol = kPatchCols - patchSize;
+  Int_t maxRow = kPatchRows - patchSize;
+
+  for (Int_t i = 0; i <= maxCol; i += stepSize) {
+    for (Int_t j = 0; j <= maxRow; j += stepSize) {
+      // get the trigger towers composing the patch
+      Int_t   adcAmp = 0;
+      Double_t enAmp = 0.;
+      // window
+      for (Int_t k = 0; k < patchSize; ++k) {
+        for (Int_t l = 0; l < patchSize; ++l) {
+         // add amplitudes
+          adcAmp += (ULong64_t)fPatchADCSimple[i+k][j+l];
+         enAmp += fPatchESimple[i+k][j+l];
+       }
+      }
+
+      if (adcAmp == 0) {
+       AliDebug(2,"EMCal trigger patch with 0 ADC counts.");
+       continue;
+      }
+
+      Int_t absId=-1;
+      Int_t cellAbsId[4]={-1,-1,-1,-1};
+
+      // get low left edge (eta max, phi min)
+      fGeom->GetAbsFastORIndexFromPositionInEMCAL(i, j, absId);
+      // convert to the 4 absId of the cells composing the trigger channel
+      fGeom->GetCellIndexFromFastORIndex(absId, cellAbsId);
+      TVector3 edge1;
+      fGeom->GetGlobal(cellAbsId[0], edge1);
+
+      // get up right edge (eta min, phi max)
+      fGeom->GetAbsFastORIndexFromPositionInEMCAL(i+patchSize-1, j+patchSize-1, absId);
+      fGeom->GetCellIndexFromFastORIndex(absId, cellAbsId);
+      TVector3 edge2;
+      fGeom->GetGlobal(cellAbsId[3], edge2);
+
+      // get the center of the patch
+      Int_t offsetCenter = TMath::FloorNint(0.5*patchSize);
+      fGeom->GetAbsFastORIndexFromPositionInEMCAL(i+offsetCenter-1, j+offsetCenter-1, absId);
+      fGeom->GetCellIndexFromFastORIndex(absId, cellAbsId);
+      TVector3 center1;
+      fGeom->GetGlobal(cellAbsId[3], center1);
+
+      fGeom->GetAbsFastORIndexFromPositionInEMCAL(i+offsetCenter, j+offsetCenter, absId);
+      fGeom->GetCellIndexFromFastORIndex(absId, cellAbsId);
+      TVector3 center2;
+      fGeom->GetGlobal(cellAbsId[0], center2);
+
+      TVector3 centerGeo(center1);
+      centerGeo += center2;
+      centerGeo *= 0.5;
+
+      // save the trigger object
+      AliEmcalTriggerPatchInfo *trigger =
+       new ((*fCaloTriggersOut)[itrig]) AliEmcalTriggerPatchInfo();
+      itrig++;
+      trigger->SetCenterGeo(centerGeo, enAmp);
+      trigger->SetEdge1(edge1, enAmp);
+      trigger->SetEdge2(edge2, enAmp);
+      trigger->SetADCAmp(adcAmp);
+      trigger->SetEdgeCell(i*2, j*2); // from triggers to cells
+    }
+  } // trigger algo
+  AliDebug(2,Form("Created %d trigger patches (%d) in this event",itrig,patchSize));
+
+}
+
+//________________________________________________________________________
+Int_t AliEmcalPatchFromCellMaker::GetDimFastor() const {
+
+  Int_t dim = TMath::FloorNint((Double_t)(fPatchDim/2.));
+  return dim;
+}
+
+//________________________________________________________________________
+Int_t AliEmcalPatchFromCellMaker::GetSlidingStepSizeFastor() const {
+
+  Int_t dim = GetDimFastor();
+  if(!fL1Slide) return dim;
+
+  if(dim==2) return 2;
+  else if(dim==4) return 4;
+  else if(dim==8) return 4;
+  else if(dim==16) return 8;
+  else return -1;
+}
+
+
+
+
diff --git a/PWG/EMCAL/AliEmcalPatchFromCellMaker.h b/PWG/EMCAL/AliEmcalPatchFromCellMaker.h
new file mode 100644 (file)
index 0000000..33c7127
--- /dev/null
@@ -0,0 +1,64 @@
+#ifndef ALIEMCALPATCHFROMCELLMAKER_H
+#define ALIEMCALPATCHFROMCELLMAKER_H
+
+class TClonesArray;
+class TH3F;
+
+#include "AliEMCALTriggerTypes.h"
+#include "AliAnalysisTaskEmcal.h"
+
+class AliEmcalPatchFromCellMaker : public AliAnalysisTaskEmcal {
+ public:
+  AliEmcalPatchFromCellMaker();
+  AliEmcalPatchFromCellMaker(const char *name);
+  virtual ~AliEmcalPatchFromCellMaker();
+
+  void               SetCaloTriggersOutName(const char *name)          { fCaloTriggersOutName = name; }
+  void               SetPatchDimension(Int_t i)                        { fPatchDim            = i;    }
+  void               SetMinCellE(Double_t e)                           { fMinCellE            = e;    }
+  void               SetCellTimeCuts(Double_t min, Double_t max)       { fCellTimeMin = min; fCellTimeMax = max; }
+  void               ActivateSlidingPatch(Bool_t b)                    { fL1Slide             = b;    }
+
+ protected:
+  enum{
+    kPatchCols = 48,
+    kPatchRows = 64
+  };
+
+  void               ExecOnce();
+  Bool_t             Run();
+  // Bool_t             FillHistograms();
+  void               UserCreateOutputObjects();
+
+  Bool_t             FillPatchADCSimple();
+  void               RunSimpleOfflineTrigger();
+
+  //Getters
+  Int_t              GetPatchDimension() const                         { return fPatchDim;  }
+  Double_t           GetPatchArea() const                              { return (Double_t)(fPatchDim*fPatchDim)*0.014*0.014; }
+  Int_t              GetDimFastor() const;
+  Int_t              GetSlidingStepSizeFastor() const;
+
+  TString            fCaloTriggersOutName;  // name of output patch array
+  TClonesArray      *fCaloTriggersOut;      //!trigger array out
+      
+  Double_t           fPatchADCSimple[kPatchCols][kPatchRows];   // patch map for simple offline trigger
+  Double_t           fPatchESimple[kPatchCols][kPatchRows];     // patch map for simple offline trigger
+
+  Int_t              fPatchDim;             // dimension of patch in #cells
+  Double_t           fMinCellE;             // minimum cell energy
+  Double_t           fCellTimeMin;          // minimum time cell
+  Double_t           fCellTimeMax;          // maximum time cell
+  Bool_t             fL1Slide;              // sliding window on
+
+ private:
+  TH3F     *fh3EEtaPhiCell;                    //! cell E, eta, phi
+  TH2F     *fh2CellEnergyVsTime;               //! emcal cell energy vs time
+  TH1F     *fh1CellEnergySum;                  //! sum of energy in all emcal cells
+
+  AliEmcalPatchFromCellMaker(const AliEmcalPatchFromCellMaker&);            // not implemented
+  AliEmcalPatchFromCellMaker &operator=(const AliEmcalPatchFromCellMaker&); // not implemented
+
+  ClassDef(AliEmcalPatchFromCellMaker, 1); // Task to make PicoTracks in a grid corresponding to EMCAL/DCAL acceptance
+};
+#endif
diff --git a/PWG/EMCAL/macros/AddTaskEmcalPatchFromCellMaker.C b/PWG/EMCAL/macros/AddTaskEmcalPatchFromCellMaker.C
new file mode 100644 (file)
index 0000000..c7637cf
--- /dev/null
@@ -0,0 +1,71 @@
+// $Id$
+
+AliEmcalPatchFromCellMaker* AddTaskEmcalPatchFromCellMaker(
+  Int_t       patchDim            = 32,
+  Double_t    cellMinE            = 0.15,
+  Bool_t      bSlideL1            = kFALSE,
+  const char *patchOutName        = "EmcalPatches",
+  const char *cellsName           = 0,
+  const char *taskName            = "PatchFromCellMaker")
+{  
+  // Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr)
+  {
+    ::Error("AddTaskEmcalPatchFromCellMaker", "No analysis manager to connect to.");
+    return NULL;
+  }  
+  
+  // Check the analysis type using the event handlers connected to the analysis manager.
+  //==============================================================================
+  AliVEventHandler *evhand = mgr->GetInputEventHandler();
+  if (!evhand) {
+    ::Error("AddTaskEmcalPatchFromCellMaker", "This task requires an input event handler");
+    return NULL;
+  }
+
+  TString strCellsName(cellsName);
+  if(strCellsName.IsNull()) {
+    if (evhand->InheritsFrom("AliESDInputHandler")) {
+      strCellsName = "EMCALCells";
+      ::Info("AddTaskEmcalPatchFromCellMaker", Form( "ESD analysis, cellsName = \"%s\"", strCellsName.Data() ));
+    }
+    else {
+      strCellsName = "emcalCells";
+      ::Info("AddTaskEmcalPatchFromCellMaker", Form( "AOD analysis, cellsName = \"%s\"", strCellsName.Data() ));
+    }
+  }
+
+  TString strPatchOutName = Form("%s%dx%d",patchOutName,patchDim,patchDim);
+  TString name = Form("%s_%s",taskName,strPatchOutName.Data());
+   //-------------------------------------------------------
+  // Init the task and do settings
+  //-------------------------------------------------------
+
+  AliEmcalPatchFromCellMaker *eTask = new AliEmcalPatchFromCellMaker(name);
+  eTask->SetCaloTriggersOutName(strPatchOutName.Data());
+  eTask->SetCaloCellsName(strCellsName.Data());
+  eTask->SetPatchDimension(patchDim);
+  eTask->SetMinCellE(cellMinE);
+  eTask->ActivateSlidingPatch(bSlideL1);
+
+  //-------------------------------------------------------
+  // Final settings, pass to manager and set the containers
+  //-------------------------------------------------------
+  mgr->AddTask(eTask);
+  
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput1  = mgr->GetCommonInputContainer();
+  TString contname(name);
+  contname += "_histos";
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(),
+                                                           TList::Class(),AliAnalysisManager::kOutputContainer,
+                                                           Form("%s", AliAnalysisManager::GetCommonFileName()));
+
+  mgr->ConnectInput  (eTask, 0,  cinput1 );
+  mgr->ConnectOutput (eTask, 1, coutput1 );
+
+  return eTask;
+}
index 9b4c490..090a182 100644 (file)
@@ -19,6 +19,7 @@
 #pragma link C++ class AliEmcalMCTrackSelector+;
 #pragma link C++ class AliEmcalParticle+;
 #pragma link C++ class AliEmcalParticleMaker+;
+#pragma link C++ class AliEmcalPatchFromCellMaker+;
 #pragma link C++ class AliEmcalPhysicsSelection+;
 #pragma link C++ class AliEmcalPhysicsSelectionTask+;
 #pragma link C++ class AliEmcalPicoTrackMaker+;
index 9e685a0..920c1f4 100644 (file)
@@ -37,6 +37,7 @@ set ( SRCS
  EMCALJetTasks/AliAnalysisTaskLocalRho.cxx
  EMCALJetTasks/AliAnalysisTaskScale.cxx
  EMCALJetTasks/AliEmcalJet.cxx
+ EMCALJetTasks/AliEmcalPicoTrackInGridMaker.cxx
  EMCALJetTasks/AliJetConstituentTagCopier.cxx
  EMCALJetTasks/AliJetContainer.cxx
  EMCALJetTasks/AliJetEmbeddingFromGenTask.cxx
@@ -61,6 +62,7 @@ set ( SRCS
  EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalDiJetBase.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalDiJetAna.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalDiJetResponse.cxx
+ EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalHighMultTrigger.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetHMEC.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetHadCorQA.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetHadEPpid.cxx
diff --git a/PWGJE/EMCALJetTasks/AliEmcalPicoTrackInGridMaker.cxx b/PWGJE/EMCALJetTasks/AliEmcalPicoTrackInGridMaker.cxx
new file mode 100644 (file)
index 0000000..db909a9
--- /dev/null
@@ -0,0 +1,813 @@
+// $Id$
+//
+// Class to put collection of tracks into grid of PicoTracks
+//
+// Author: M. Verweij
+
+#include <TClonesArray.h>
+#include <TRandom3.h>
+#include <TProfile.h>
+#include <TH3F.h>
+#include "AliAODEvent.h"
+#include "AliAODTrack.h"
+#include "AliAnalysisManager.h"
+#include "AliESDtrack.h"
+#include "AliESDtrackCuts.h"
+#include "AliLog.h"
+#include "AliPicoTrack.h"
+#include "AliVTrack.h"
+#include "AliEmcalJet.h"
+#include "AliEmcalPicoTrackInGridMaker.h"
+
+ClassImp(AliEmcalPicoTrackInGridMaker)
+
+//________________________________________________________________________
+AliEmcalPicoTrackInGridMaker::AliEmcalPicoTrackInGridMaker() : 
+AliAnalysisTaskSE("AliEmcalPicoTrackInGridMaker"),
+  fTracksOutName("PicoTracksInGrid"),
+  fTracksInName("tracks"),
+  fTracksIn(0),
+  fTracksOut(0),
+  fL1Slide(0),
+  fCellSize(0.0145),
+  fMinCellE(0.15),
+  fExclLeadingPatch(0),
+  fPatchSub(3),
+  fRhoMean(184.),
+  fNCells(-1),
+  fNCellsEMCal(-1),
+  fNCellsDCal(-1),
+  fCellGrid(),
+  fMiniPatchGrid(),
+  fActiveAreaMP(),
+  fMultVsRho(0),
+  fHistList(0)
+{
+  // Constructor.
+
+  fPhiMin[0] = 1.405;
+  fPhiMax[0] = 1.405+TMath::DegToRad()*110.;
+  fPhiMin[1] = 4.547;
+  fPhiMax[1] = 5.71;
+  fEtaMin[0] = -0.7;
+  fEtaMax[0] = 0.7;
+  fEtaMin[1] = -0.7;
+  fEtaMax[1] = 0.7;
+
+  for(Int_t i = 0; i<5; i++) {
+    fPatchGrid[i] = 0;
+    fNPatchesEMCal[i] = 0;
+    fActiveAreaMPP[i] = 0;
+    fActiveAreaCP[i]  = 0;
+
+    fh1RhoEmcal[i] = 0;
+    fh1RhoDcal[i] = 0;
+    fPatchEnVsActivityEmcal[i] = 0;
+    fPatchEnVsActivityDcal[i]  = 0;
+
+    for(Int_t j = 0; j<2; j++) {
+      fPatchECorr[j][i] = 0;
+      fPatchECorrPar[j][i] = 0;
+      fPatchECorrRho[j][i] = 0;
+      fPatchECorrRhoDijet[j][i] = 0;
+      fPatchECorrECorrRho[j][i] = 0;
+    }
+  }
+
+  for(Int_t i = 0; i<3; i++) {
+    fh2MedianTypeEmcal[i] = 0;
+    fh2MedianTypeDcal[i] = 0;
+    fpMedianTypeEmcal[i] = 0;
+    fpMedianTypeDcal[i] = 0;
+  }
+
+}
+
+//________________________________________________________________________
+AliEmcalPicoTrackInGridMaker::AliEmcalPicoTrackInGridMaker(const char *name) : 
+  AliAnalysisTaskSE(name),
+  fTracksOutName("PicoTracksInGrid"),
+  fTracksInName("tracks"),
+  fTracksIn(0),
+  fTracksOut(0),
+  fL1Slide(0),
+  fCellSize(0.0145),
+  fMinCellE(0.15),
+  fExclLeadingPatch(0),
+  fPatchSub(3),
+  fRhoMean(184.),
+  fNCells(-1),
+  fNCellsEMCal(-1),
+  fNCellsDCal(-1),
+  fCellGrid(),
+  fMiniPatchGrid(),
+  fActiveAreaMP(),
+  fMultVsRho(0),
+  fHistList(0)
+{
+  // Constructor.
+
+  fPhiMin[0] = 1.405;
+  fPhiMax[0] = 1.405+TMath::DegToRad()*110.;
+  fPhiMin[1] = 4.547;
+  fPhiMax[1] = 5.71;
+  fEtaMin[0] = -0.7;
+  fEtaMax[0] = 0.7;
+  fEtaMin[1] = -0.7;
+  fEtaMax[1] = 0.7;
+
+  for(Int_t i = 0; i<5; i++) {
+    fPatchGrid[i] = 0;
+    fNPatchesEMCal[i] = 0;
+    fActiveAreaMPP[i] = 0;
+    fActiveAreaCP[i]  = 0;
+
+    fh1RhoEmcal[i] = 0;
+    fh1RhoDcal[i] = 0;
+
+    fPatchEnVsActivityEmcal[i] = 0;
+    fPatchEnVsActivityDcal[i]  = 0;
+
+    for(Int_t j = 0; j<2; j++) {
+      fPatchECorr[j][i] = 0;
+      fPatchECorrPar[j][i] = 0;
+      fPatchECorrRho[j][i] = 0;
+      fPatchECorrRhoDijet[j][i] = 0;
+      fPatchECorrECorrRho[j][i] = 0;
+    }
+  }
+
+  for(Int_t i = 0; i<3; i++) {
+    fh2MedianTypeEmcal[i] = 0;
+    fh2MedianTypeDcal[i] = 0;
+    fpMedianTypeEmcal[i] = 0;
+    fpMedianTypeDcal[i] = 0;
+  }
+
+  // Output slot #1 write into a TList
+  DefineOutput(1, TList::Class());
+}
+
+//________________________________________________________________________
+AliEmcalPicoTrackInGridMaker::~AliEmcalPicoTrackInGridMaker()
+{
+  // Destructor.
+}
+
+//________________________________________________________________________
+void AliEmcalPicoTrackInGridMaker::UserCreateOutputObjects()
+{
+  // Create my user objects.
+
+  OpenFile(1);
+  fHistList = new TList();
+  fHistList->SetOwner(kTRUE);
+
+  Int_t nBinsMed = 500;
+  Double_t minMed = 0.;
+  Double_t maxMed = 500.;
+
+  // Int_t nBinsMult = 150;
+  // Double_t minMult = 0.;
+  // Double_t maxMult = 4000.;
+
+  for(Int_t i = 0; i<3; i++) {
+    fh2MedianTypeEmcal[i] = new TH2F(Form("fh2MedianTypeEmcalAreaType%d",i),Form("fh2MedianTypeEmcalAreaType%d",i),5,0.5,5.5,nBinsMed,minMed,maxMed);
+    fHistList->Add(fh2MedianTypeEmcal[i]);
+
+    fh2MedianTypeDcal[i] = new TH2F(Form("fh2MedianTypeDcalAreaType%d",i),Form("fh2MedianTypeDcalAreaType%d",i),5,0.5,5.5,nBinsMed,minMed,maxMed);
+    fHistList->Add(fh2MedianTypeDcal[i]);
+
+    fpMedianTypeEmcal[i] = new TProfile(Form("fpMedianTypeEmcalAreaType%d",i),Form("fpMedianTypeEmcalAreaType%d",i),5,0.5,5.5,"s");
+    fHistList->Add(fpMedianTypeEmcal[i]);
+
+    fpMedianTypeDcal[i] = new TProfile(Form("fpMedianTypeDcalAreaType%d",i),Form("fpMedianTypeDcalAreaType%d",i),5,0.5,5.5,"s");
+    fHistList->Add(fpMedianTypeDcal[i]);
+  }
+
+  TString det[2] = {"Emcal","Dcal"};
+  for(Int_t i = 0; i<5; i++) { //loop over patch types
+    fh1RhoEmcal[i] = new TH1F(Form("fh1RhoEmcal_%d",i),Form("fh1RhoEmcal_%d",i),500,0.,1000.);
+    fHistList->Add(fh1RhoEmcal[i]);
+    fh1RhoDcal[i] = new TH1F(Form("fh1RhoDcal_%d",i),Form("fh1RhoDcal_%d",i),500,0.,1000.);
+    fHistList->Add(fh1RhoDcal[i]);
+
+    fPatchEnVsActivityEmcal[i] = new TH2F(Form("fh2PatchEnVsActivityEmcal_%d",i),Form("fh2PatchEnVsActivityEmcal_%d",i),300,0.,300.,150,-0.5,149.5);
+    fHistList->Add(fPatchEnVsActivityEmcal[i]);
+
+    fPatchEnVsActivityDcal[i] = new TH2F(Form("fh2PatchEnVsActivityDcal_%d",i),Form("fh2PatchEnVsActivityDcal_%d",i),300,0.,300.,150,-0.5,149.5);
+    fHistList->Add(fPatchEnVsActivityDcal[i]);
+
+    for(Int_t j = 0; j<2; j++) {
+      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.);
+      fHistList->Add(fPatchECorr[j][i]);
+
+      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.);
+      fHistList->Add(fPatchECorrPar[j][i]);  
+
+      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.);
+      fHistList->Add(fPatchECorrRho[j][i]);  
+
+      fPatchECorrRhoDijet[j][i] = new TH2F(Form("fPatchECorrRhoDijet%s_%d",det[j].Data(),i),Form("fPatchECorrRhoDijet%s_%d;#it{E}_{patch}^{corr};#rho",det[j].Data(),i),250,-50.,200.,500,0.,500.);
+      fHistList->Add(fPatchECorrRhoDijet[j][i]);
+
+      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.);
+      fHistList->Add(fPatchECorrECorrRho[j][i]);
+    }
+  }
+
+  fMultVsRho = new TH2F("fMultVsRho","fMultVsRho",3000,0,3000,400,0,400);
+  fHistList->Add(fMultVsRho);
+
+  PostData(1, fHistList);
+}
+
+//________________________________________________________________________
+void AliEmcalPicoTrackInGridMaker::UserExec(Option_t *) 
+{
+  // Main loop, called for each event.
+
+  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
+  if (!am) {
+    AliError("Manager zero, returning");
+    return;
+  }
+
+  // retrieve tracks from input.
+  if (!fTracksIn) { 
+    fTracksIn = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksInName));
+    if (!fTracksIn) {
+      AliError(Form("Could not retrieve tracks %s!", fTracksInName.Data())); 
+      return;
+    }
+    if (!fTracksIn->GetClass()->GetBaseClass("AliVParticle")) {
+      AliError(Form("%s: Collection %s does not contain AliVParticle objects!", GetName(), fTracksInName.Data())); 
+      return;
+    }
+  }
+
+  Bool_t b = CreateGridCells();
+  if(!b) return;
+  b = CreateGridMiniPatches();
+  if(!b) return;
+
+  //L0 single shower trigger
+  CreateGridPatches(4,0);
+
+  //L1 triggers: sliding window
+  CreateGridPatches(4,1);
+  CreateGridPatches(8,1);
+  CreateGridPatches(16,1);
+  CreateGridPatches(32,1);
+
+  Double_t medL0 = CalculateMedian(0,0);
+  fh2MedianTypeEmcal[0]->Fill(0.5,medL0);
+  fpMedianTypeEmcal[0]->Fill(0.5,medL0);
+  medL0 = CalculateMedian(0,1);
+  fh2MedianTypeDcal[0]->Fill(0.5,medL0);
+  fpMedianTypeDcal[0]->Fill(0.5,medL0);
+
+  Double_t medL1[4][2];
+  for(Int_t i = 0; i<4; i++) { //patches
+    for(Int_t type = 0; type<2; type++) { //EMCal or DCal
+      for(Int_t areaT = 0; areaT<1; areaT++) { //areay type (passive vs active)
+       medL1[i][type] = CalculateMedian(i+1,type,areaT);
+       if(type==0) {
+         fh2MedianTypeEmcal[areaT]->Fill((Double_t)(i+1)+0.5,medL1[i][type]);
+         fpMedianTypeEmcal[areaT]->Fill((Double_t)(i+1)+0.5,medL1[i][type]);
+       }
+       if(type==1) {
+         fh2MedianTypeDcal[areaT]->Fill((Double_t)(i+1)+0.5,medL1[i][type]);
+         fpMedianTypeDcal[areaT]->Fill((Double_t)(i+1)+0.5,medL1[i][type]);
+       }
+      }
+    }
+  }
+
+  // subtract energy density and store energy distributions of corrected patches in histo
+  for(Int_t i = 1; i<5; i++) { //patch types
+    //    Int_t EleadID[2] = {-1,-1};
+    Double_t Elead[2] = {0.,0.};
+    for(Int_t j = 0; j<fPatchGrid[i].GetSize(); j++) { //patches
+      if(fPatchGrid[i].At(j)>0.) { //don't do anything with empty patches
+       Int_t type = 0; //EMCal
+       Int_t subType = 1;
+       if(j>=fNPatchesEMCal[i]) {
+         type = 1; //DCal
+         subType = 0;
+       }
+       Double_t sub = medL1[fPatchSub-1][subType]*GetPatchArea(i);
+       fPatchECorr[type][i]->Fill(fPatchGrid[i].At(j) - sub);
+       fPatchECorrPar[type][i]->Fill(fPatchGrid[i].At(j) - fRhoMean*GetPatchArea(i));
+       
+       //Bookkeep leading patches
+       if((fPatchGrid[i].At(j)-sub)>Elead[type]) {
+         //      EleadID[type] = j;
+         Elead[type] = fPatchGrid[i].At(j)-sub;
+       }
+      }
+    }
+
+    for(Int_t k = 0; k<2; k++) {
+      Int_t subType = 1;
+      if(k==1) subType=0;
+      fPatchECorrRho[k][i]->Fill(Elead[k],medL1[fPatchSub-1][subType]);
+      fPatchECorrECorrRho[k][i]->Fill(Elead[k],Elead[subType],medL1[fPatchSub-1][subType]);
+      if(Elead[subType]>30.) {
+       fPatchECorrRhoDijet[k][i]->Fill(Elead[k],medL1[fPatchSub-1][subType]);
+      }
+    }
+  }
+
+  fMultVsRho->Fill(fTracksIn->GetEntriesFast(),medL1[3][0]);
+}
+
+//________________________________________________________________________
+Double_t AliEmcalPicoTrackInGridMaker::CalculateSum(Int_t patchType) {
+
+  Int_t n = fPatchGrid[patchType].GetSize();
+  if(n<1) return -1.;
+
+  Double_t sum = 0.;
+  Int_t count = 0;
+  for(Int_t i = 0; i<fPatchGrid[patchType].GetSize(); i++) {
+    if(fPatchGrid[patchType].At(i)>0.) count++;
+    sum+=fPatchGrid[patchType].At(i);
+  }
+  //  Double_t avg = -1.;
+  //  if(count>0) avg = sum/((Double_t)count);
+  return sum;
+}
+
+//________________________________________________________________________
+Double_t AliEmcalPicoTrackInGridMaker::CalculateMedian(Int_t patchType, Int_t type, Int_t areaType) {
+  //areaType:
+  //0: passive area
+  //1: active area mini patches
+  //2: active arear cells
+
+  Int_t n = fPatchGrid[patchType].GetSize();
+  if(n<1) return -1.;
+  //  Double_t *arr = fPatchGrid[patchType].GetArray();
+  
+  static Double_t arr[999];
+  Int_t c = 0;
+  Int_t start = 0;
+  Int_t last = n;
+  if(type==0) {
+    start = 0;
+    last = fNPatchesEMCal[patchType];
+  } else if(type==1) {
+    start = fNPatchesEMCal[patchType];
+    last = n;
+  }
+
+  Double_t area = GetPatchArea(patchType);
+
+  //find patch with highest energy
+  Int_t imax = -1;
+  Double_t max = 0.;
+  for(Int_t i = start; i<last; i++) {
+    if(fPatchGrid[patchType].At(i)>max) {
+      imax = i;
+      max = fPatchGrid[patchType].At(i);
+    }
+  }
+  for(Int_t i = start; i<last; i++) {
+    if(fExclLeadingPatch>0 && i==imax) continue;
+    if(fPatchGrid[patchType].At(i)>0.) {
+      Int_t active = 99;
+      if(areaType==1) active = fActiveAreaMPP[patchType].At(i);
+      else if(areaType==2) active = fActiveAreaCP[patchType].At(i);
+      if(areaType>0) area = GetPatchAreaActive(i,patchType,areaType-1);
+
+      if(area>0. && active>1) {
+       arr[c] = fPatchGrid[patchType].At(i)/area;
+       c++;
+      }
+      if(type==0 && areaType==0) {
+       fh1RhoEmcal[patchType]->Fill(arr[c-1]);
+       fPatchEnVsActivityEmcal[patchType]->Fill(fPatchGrid[patchType].At(i),fActiveAreaCP[patchType].At(i));
+       if(fActiveAreaCP[patchType].At(i)==0) Printf("WARNING activity: %d while En=%f active from MP: %d",fActiveAreaCP[patchType].At(i),fPatchGrid[patchType].At(i),fActiveAreaMPP[patchType].At(i));
+      }
+      if(type==1 && areaType==0) {
+       fh1RhoDcal[patchType]->Fill(arr[c-1]);
+       fPatchEnVsActivityDcal[patchType]->Fill(fPatchGrid[patchType].At(i),fActiveAreaCP[patchType].At(i));
+      }
+    }
+  }
+  Double_t med = TMath::Median(c,arr);
+
+  return med;
+}
+
+//________________________________________________________________________
+Bool_t AliEmcalPicoTrackInGridMaker::CreateGridCells() {
+
+  if(!InitCells()) return kFALSE;
+
+  // loop over tracks
+  const Int_t ntr = fTracksIn->GetEntriesFast();
+  for (Int_t itr = 0; itr < ntr; ++itr) {
+
+    AliVParticle *track = static_cast<AliVParticle*>(fTracksIn->At(itr));
+    if (!track)
+      continue;
+
+    if(track->Pt()<fMinCellE) continue;
+
+    Int_t id = GetGridID(track);
+    if(id>-1) {
+      //      Double_t old = fCellGrid.At(id);
+      fCellGrid.AddAt(fCellGrid.At(id)+track->Pt(),id);
+    }
+  }
+
+  return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliEmcalPicoTrackInGridMaker::CreateGridMiniPatches() {
+
+  if(!InitMiniPatches()) return kFALSE;
+
+  //Int_t nMiniPatches = fMiniPatchGrid.GetSize();
+  //Int_t nMiniPatchesE = (Int_t)(0.25*fNCellsEMCal);
+  
+  Int_t type = 0;
+  // Int_t nMiniPatchesRow = TMath::FloorNint(0.5*GetNCellsRow(type));
+  //loop over edges of mini patches
+  Int_t nm = 0; //mini patch number
+  for(Int_t i = 0; i<(GetNCellsCol(type)-1); i+=2) {
+    for(Int_t j = 0; j<(GetNCellsRow(type)-1); j+=2) {
+      //loop over cells in mini patch
+      for(Int_t k = 0; k<2; k++) {
+       for(Int_t l = 0; l<2; l++) {
+         Int_t row = i+k;
+         Int_t col = j+l;
+         Int_t id = GetGridID(row,col,type);
+         fMiniPatchGrid.AddAt(fMiniPatchGrid.At(nm)+fCellGrid.At(id),nm);
+         if(fCellGrid.At(id)>0.)
+           fActiveAreaMP.AddAt(fActiveAreaMP.At(nm)+1,nm);
+       }
+      }
+      nm++;
+    }
+  }
+
+  //DCal mini patches
+  type = 1;
+  for(Int_t i = 0; i<(GetNCellsCol(type)-1); i+=2) {
+    for(Int_t j = 0; j<(GetNCellsRow(type)-1); j+=2) {
+      //loop over cells in mini patch
+      for(Int_t k = 0; k<2; k++) {
+       for(Int_t l = 0; l<2; l++) {
+         Int_t row = i+k;
+         Int_t col = j+l;
+         Int_t id = GetGridID(row,col,type);
+         fMiniPatchGrid.AddAt(fMiniPatchGrid.At(nm)+fCellGrid.At(id),nm);
+         if(fCellGrid.At(id)>0.) 
+           fActiveAreaMP.AddAt(fActiveAreaMP.At(nm)+1,nm);
+       }
+      }
+      nm++;
+    }
+  }
+
+  return kTRUE;
+}
+
+//________________________________________________________________________
+Int_t AliEmcalPicoTrackInGridMaker::GetNRowMiniPatches(Int_t type) {
+  //returns number of rows of mini patches in detector of type X (0: EMCal 1: DCal)
+  Int_t nRows = TMath::FloorNint(0.5*GetNCellsCol(type));
+  return nRows;
+}
+
+//________________________________________________________________________
+Int_t AliEmcalPicoTrackInGridMaker::GetNColMiniPatches(Int_t type) {
+  //returns number of rows of mini patches in detector of type X (0: EMCal 1: DCal)
+  Int_t nCols = TMath::FloorNint(0.5*GetNCellsRow(type));
+  return nCols;
+}
+
+//________________________________________________________________________
+Bool_t AliEmcalPicoTrackInGridMaker::CreateGridPatches(Int_t dim, Int_t level) {
+
+  if(!InitPatches(dim,level)) return kFALSE;
+
+  Int_t pt = GetPatchType(dim,level);
+  Int_t nm = (Int_t)(dim/2.);
+  if(level==1 && fL1Slide)  nm = GetSlidingStepSizeMiniPatches(dim);
+  //loop over edges of mini patches
+  Int_t np = 0; //patch number
+  for(Int_t type = 0; type<2; type++) {
+    for(Int_t i = 0; i<=(GetNColMiniPatches(type)-nm); i+=nm) {
+      for(Int_t j = 0; j<=(GetNRowMiniPatches(type)-nm); j+=nm) {
+       //loop over mini patches in patch
+       for(Int_t k = 0; k<nm; k++) {
+         for(Int_t l = 0; l<nm; l++) {
+           Int_t row = i+k;
+           Int_t col = j+l;
+           Int_t id = GetMiniPatchID(row,col,type);
+           fPatchGrid[pt].AddAt(fPatchGrid[pt].At(np)+fMiniPatchGrid.At(id),np);
+           if(fMiniPatchGrid.At(id)>0.) {
+             fActiveAreaMPP[pt].AddAt(fActiveAreaMPP[pt].At(np)+1,np);
+             fActiveAreaCP[pt].AddAt(fActiveAreaCP[pt].At(np)+fActiveAreaMP.At(id),np);
+           }
+         }
+       }
+       //Sanity check
+       if(fActiveAreaMPP[pt].At(np)>1 && fActiveAreaCP[pt].At(np)==0) {
+         Printf("Activity MP: %d Cells: %d",fActiveAreaMPP[pt].At(np),fActiveAreaCP[pt].At(np));
+         Printf("pt: %d np: %d en: %f",pt,np,fPatchGrid[pt].At(np));
+       }
+       np++;
+      }
+    }
+  }
+
+  return kTRUE;
+}
+
+//________________________________________________________________________
+Int_t AliEmcalPicoTrackInGridMaker::GetMiniPatchID(Int_t row, Int_t col, Int_t type) {
+  Int_t id = row*GetNColMiniPatches(type) + col;
+  //if in DCal move ID to after all EMCal IDs
+  if(type==1) {
+    Int_t ne = GetNColMiniPatches(0) * GetNRowMiniPatches(0);
+    id+=ne;
+  }
+  return id;
+}
+
+//________________________________________________________________________
+Int_t AliEmcalPicoTrackInGridMaker::GetCellType(Double_t eta, Double_t phi) {
+
+  for(Int_t i = 0; i<2; i++) {
+    if(eta>fEtaMin[i] && eta<fEtaMax[i] && phi>fPhiMin[i] && phi<fPhiMax[i]) return i;
+  }
+  return -1;
+}
+
+//________________________________________________________________________
+Int_t AliEmcalPicoTrackInGridMaker::GetGridID(Int_t row, Int_t col, Int_t type) {
+  Int_t id = row*GetNCellsRow(type) + col;
+  //if in DCal move ID to after all EMCal IDs
+  if(type==1) id+=fNCellsEMCal;
+
+  return id;
+}
+
+//________________________________________________________________________
+Int_t AliEmcalPicoTrackInGridMaker::GetGridID(Double_t eta, Double_t phi) {
+  
+  Int_t type = GetCellType(eta,phi);
+  if(type<0 || type>1) return -1; //position is not in EMCal or DCal
+
+  // grid ID convention:
+  // upper left corner (min phi, min eta) is first ID
+  // then walk through grid from upper left to lower right accross the rows in phi
+
+  Int_t id = -1;
+  Double_t etaRel = eta-fEtaMin[type];
+  Double_t phiRel = phi-fPhiMin[type];
+  Int_t row = TMath::FloorNint(etaRel/fCellSize);
+  Int_t col = TMath::FloorNint(phiRel/fCellSize);
+  id = GetGridID(row,col,type);
+
+  if(id>=fNCells) {
+    Printf("Got too large id %d %d type: %d",id,fNCells,type);
+    Printf("eta: %f phi: %f",eta,phi);
+    Printf("etaRel: %f phiRel: %f",etaRel,phiRel);
+    Printf("row: %d col: %d  ->  %d + %d = %d",row,col,row*GetNCellsRow(type) + col,fNCellsEMCal,id);
+    Printf("n cells row: %d",GetNCellsRow(type));
+    Printf("\n");
+  }
+
+  return id;
+}
+
+//________________________________________________________________________
+Int_t AliEmcalPicoTrackInGridMaker::GetNCellsCol(Int_t type) {
+
+  Double_t deta = fEtaMax[type] - fEtaMin[type];
+  Int_t nCellsCol = TMath::FloorNint(deta/fCellSize);
+  return nCellsCol;
+}
+
+//________________________________________________________________________
+Int_t AliEmcalPicoTrackInGridMaker::GetNCellsRow(Int_t type) {
+
+  Double_t dPhi = fPhiMax[type] - fPhiMin[type];
+  Int_t nCellsRow = TMath::FloorNint(dPhi/fCellSize);
+  return nCellsRow;
+}
+
+//________________________________________________________________________
+Bool_t AliEmcalPicoTrackInGridMaker::InitCells() {
+
+  CheckEdges();
+  if(!CheckEdges()) return kFALSE;
+
+  //number of cells in EMCal acceptance
+  Int_t nCellsPhiE = GetNCellsCol(0);
+  Int_t nCellsEtaE = GetNCellsRow(0);
+  fNCellsEMCal = nCellsPhiE*nCellsEtaE;
+  
+  //number of cells in DCal acceptance
+  Int_t nCellsPhiD = GetNCellsCol(1);
+  Int_t nCellsEtaD = GetNCellsRow(1);
+  fNCellsDCal = nCellsPhiD*nCellsEtaD;
+
+  //total number of cells
+  fNCells = fNCellsEMCal + fNCellsDCal;
+
+  /*
+    Printf("EMCal: %d x %d",nCellsEtaE,nCellsPhiE);
+    Printf("DCal: row: %d x col: %d",nCellsEtaD,nCellsPhiD);
+    Printf("fNCells: %d fNCellsE: %d fnCellsD: %d",fNCells,fNCellsEMCal,fNCellsDCal);
+  */
+  fCellGrid.Set(fNCells);
+  fCellGrid.Reset(0);
+  return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliEmcalPicoTrackInGridMaker::InitMiniPatches() {
+  
+  if(fCellGrid.GetSize()<0) return kFALSE;
+  Double_t conv = 0.25; //dimension of mini patch is 2x2 cells
+  Int_t nMiniPatchesE = (Int_t)(conv*fNCellsEMCal);
+  Int_t nMiniPatchesD = (Int_t)(conv*fNCellsDCal);
+
+  //total number of mini patches
+  Int_t nMiniPatches = nMiniPatchesE + nMiniPatchesD;
+
+  fMiniPatchGrid.Set(nMiniPatches);
+  fMiniPatchGrid.Reset(0.);
+
+  fActiveAreaMP.Set(nMiniPatches);
+  fActiveAreaMP.Reset(0);
+  return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliEmcalPicoTrackInGridMaker::InitPatches(Int_t dim, Int_t level) {
+  // dimensions in cell units
+  // if level==1: sliding window will be applied
+  // L1 4x4: slide by 2 cells: 1 mini patch
+  // L1 8x8: slide by 4 cells: 2 mini patches
+  // L1 16x16: slide by 4 cells: 2 mini patches
+  // L1 32x32: slide by 8 cells: 4 mini patches
+  
+  if(fCellGrid.GetSize()<0) return kFALSE;
+
+  Double_t dimd = (Double_t)dim;
+  Double_t conv = 1./(dimd*dimd);
+  if(level==1) {
+    Double_t step = (Double_t)GetSlidingStepSizeCells(dim);
+    conv = 1./(step*step);
+  }
+  Int_t nPatchesE = (Int_t)(conv*fNCellsEMCal);
+  Int_t nPatchesD = (Int_t)(conv*fNCellsDCal);
+
+  //total number of mini patches
+  Int_t nPatches = nPatchesE + nPatchesD;
+
+  Int_t type = GetPatchType(dim,level);
+  if(type<0 || type>4) return kFALSE;
+  fNPatchesEMCal[type] = nPatchesE;
+  AliDebug(2,Form("Create trigger patch of type %d with dim %d and with %d patches EMCAL: %d DCAL: %d",type,dim,nPatches,nPatchesE,nPatchesD));
+  fPatchGrid[type].Set(nPatches);
+  fPatchGrid[type].Reset(0);
+  fActiveAreaMPP[type].Set(nPatches);
+  fActiveAreaMPP[type].Reset(0);
+  fActiveAreaCP[type].Set(nPatches);
+  fActiveAreaCP[type].Reset(0);
+
+  return kTRUE;
+}
+//________________________________________________________________________
+Double_t AliEmcalPicoTrackInGridMaker::GetPatchArea(Int_t ipatch) {
+
+  Double_t ncell = 4.;
+  if(ipatch==0) ncell = 4.;
+  if(ipatch==1) ncell = 4.;
+  if(ipatch==2) ncell = 8.;
+  if(ipatch==3) ncell = 16.;
+  if(ipatch==4) ncell = 32.;
+  Double_t area = ncell*ncell*fCellSize*fCellSize;
+  return area;
+}
+
+//________________________________________________________________________
+Double_t AliEmcalPicoTrackInGridMaker::GetPatchAreaActive(Int_t id, Int_t ipatch, Int_t type) {
+  //type = 0 : active area from mini patches
+  //type = 1 : active area from cells
+
+  Int_t active = 0;
+  if(type==0) active = fActiveAreaMPP[ipatch].At(id);
+  else if(type==1) active = fActiveAreaCP[ipatch].At(id);
+  else return -1;
+
+  Double_t fac = 1.;
+  if(type==0) fac = 4.;
+  Double_t area = active*fac*fCellSize*fCellSize;
+   return area;
+}
+
+//________________________________________________________________________
+Int_t AliEmcalPicoTrackInGridMaker::GetSlidingStepSizeCells(Int_t dim) {
+
+  if(!fL1Slide) return dim;
+
+  if(dim==4) return 2;
+  else if(dim==8) return 4;
+  else if(dim==16) return 4;
+  else if(dim==32) return 8;
+  else return -1;
+}
+
+//________________________________________________________________________
+Int_t AliEmcalPicoTrackInGridMaker::GetSlidingStepSizeMiniPatches(Int_t dim) {
+  //returns step size in mini patches
+  Double_t step = (Double_t)GetSlidingStepSizeCells(dim);
+  if(step<0) return step;
+  Int_t stepMiniPatch = (Int_t)(step/2.);
+  return stepMiniPatch;
+}
+
+//________________________________________________________________________
+Int_t AliEmcalPicoTrackInGridMaker::GetPatchType(Int_t dim, Int_t level) {
+  Int_t type = -1;
+  if(level==0 && dim==4) type = 0;
+  if(level==1) {
+    if(dim==4)  type = 1;
+    if(dim==8)  type = 2;
+    if(dim==16) type = 3;
+    if(dim==32) type = 4;
+  }
+  return type;
+}
+
+//________________________________________________________________________
+Bool_t AliEmcalPicoTrackInGridMaker::CheckEdges() {
+  //
+  //Check if defined edges of EMCal and DCal make sense
+  //
+  if(fPhiMin[0]<0. || fPhiMax[0]<fPhiMin[0]) {
+    AliDebug(11,Form("EMCal phi edges not defined %f-%f",fPhiMin[0],fPhiMax[0]));
+    return kFALSE;
+  }
+
+  if(fPhiMin[1]<0. || fPhiMax[1]<fPhiMin[1]) {
+    AliDebug(11,Form("DCal phi edges not defined %f-%f",fPhiMin[1],fPhiMax[1]));
+    return kFALSE;
+  }
+
+  if(fEtaMin[0]<-10. || fEtaMax[0]<fEtaMin[1]) {
+    AliDebug(11,Form("EMCal eta edges not well defined %f-%f",fEtaMin[0],fEtaMax[0]));
+    return kFALSE;
+  }
+
+  if(fEtaMin[1]<-10. || fEtaMax[1]<fEtaMin[1]) {
+    AliDebug(11,Form("DCal eta edges not well defined %f-%f",fEtaMin[1],fEtaMax[1]));
+    return kFALSE;
+  }
+
+  for(Int_t type = 0; type<2; type++) {
+    Double_t dphi = fPhiMax[type] - fPhiMin[type];
+    Double_t nCellsColExact = dphi/fCellSize;
+    
+    Double_t deta = fEtaMax[type] - fEtaMin[type];
+    Double_t nCellsRowExact = deta/fCellSize;
+
+    Double_t col_extra = nCellsColExact - TMath::Floor(nCellsColExact);
+    Double_t phi_extra = col_extra*fCellSize;
+    fPhiMin[type] += phi_extra*0.5;
+    fPhiMax[type] -= phi_extra*0.5;
+
+    Double_t row_extra = nCellsRowExact - TMath::Floor(nCellsRowExact);
+    Double_t eta_extra = row_extra*fCellSize;
+    fEtaMin[type] += eta_extra*0.5;
+    fEtaMax[type] -= eta_extra*0.5;
+    AliDebug(2,Form("type: %d  exact: col: %f row: %f",type,nCellsColExact,nCellsRowExact));
+    //PrintAcceptance();
+  }
+  return kTRUE;
+}
+
+//________________________________________________________________________
+void AliEmcalPicoTrackInGridMaker::PrintAcceptance() {
+
+  Printf("EMCal");
+  Printf("phi: %f-%f",fPhiMin[0],fPhiMax[0]);
+  Printf("eta: %f-%f",fEtaMin[0],fEtaMax[0]);
+
+  Printf("DCal");
+  Printf("phi: %f-%f",fPhiMin[1],fPhiMax[1]);
+  Printf("eta: %f-%f",fEtaMin[1],fEtaMax[1]);
+
+}
diff --git a/PWGJE/EMCALJetTasks/AliEmcalPicoTrackInGridMaker.h b/PWGJE/EMCALJetTasks/AliEmcalPicoTrackInGridMaker.h
new file mode 100644 (file)
index 0000000..dc755a1
--- /dev/null
@@ -0,0 +1,123 @@
+#ifndef ALIEMCALPICOTRACKINGRIDMAKER_H
+#define ALIEMCALPICOTRACKINGRIDMAKER_H
+
+class TClonesArray;
+class AliVTrack;
+class AliVParticle;
+class TProfile;
+class TH3F;
+
+#include "AliAnalysisTaskSE.h"
+
+class AliEmcalPicoTrackInGridMaker : public AliAnalysisTaskSE {
+ public:
+  AliEmcalPicoTrackInGridMaker();
+  AliEmcalPicoTrackInGridMaker(const char *name);
+  virtual ~AliEmcalPicoTrackInGridMaker();
+
+  void               SetTracksInName(const char *name)                 { fTracksInName      = name; }
+  void               SetTracksOutName(const char *name)                { fTracksOutName     = name; }
+
+  void               SetCellSize(Double_t a)                           { fCellSize          = a;    }
+  void               SetMinCellE(Double_t e)                           { fMinCellE          = e;    }
+
+  void               SetEMCalAcceptance(Double_t phiMin, Double_t phiMax, Double_t etaMin, Double_t etaMax) { fPhiMin[0]=phiMin; fPhiMax[0]=phiMax; fEtaMin[0]=etaMin; fEtaMax[0]=etaMax; }
+  void               SetDCalAcceptance(Double_t phiMin, Double_t phiMax, Double_t etaMin, Double_t etaMax) { fPhiMin[1]=phiMin; fPhiMax[1]=phiMax; fEtaMin[1]=etaMin; fEtaMax[1]=etaMax; }
+
+  void               SetExcludeLeadingPatch(Int_t i)                   { fExclLeadingPatch = i; }
+
+  void               SetPatchTypeForSubtraction(Int_t i)               { fPatchSub = i; }
+  void               SetPatchTypeForSubtraction(Int_t dim, Int_t lev)  { fPatchSub = GetPatchType(dim,lev); }
+  void               SetMeanRho(Double_t r)                            { fRhoMean = r; }
+
+ protected:
+  void               UserCreateOutputObjects();
+  void               UserExec(Option_t *option);
+
+  Bool_t             InitCells();
+  Bool_t             CreateGridCells();
+  Bool_t             CheckEdges();
+  void               PrintAcceptance();
+
+  Bool_t             InitMiniPatches();
+  Bool_t             CreateGridMiniPatches();
+  Bool_t             CreateGridPatches(Int_t dim, Int_t level);
+
+  Bool_t             InitPatches(Int_t dim, Int_t level); //give dimension in cell units
+
+  Double_t           CalculateMedian(Int_t patchType, Int_t type, Int_t areaType = 0);
+  Double_t           CalculateSum(Int_t patchType);
+
+  //Getters
+  Int_t              GetGridID(AliVParticle *vp) {return GetGridID(vp->Eta(),vp->Phi());}
+  Int_t              GetGridID(Double_t eta, Double_t phi);
+  Int_t              GetGridID(Int_t row, Int_t col, Int_t type);
+  Int_t              GetCellType(Double_t eta, Double_t phi);
+  Int_t              GetNCellsRow(Int_t type);
+  Int_t              GetNCellsCol(Int_t type);
+
+  Int_t              GetNRowMiniPatches(Int_t type);
+  Int_t              GetNColMiniPatches(Int_t type);
+  Int_t              GetMiniPatchID(Int_t row, Int_t col, Int_t type);
+
+  Int_t              GetPatchType(Int_t dim, Int_t level);
+  Int_t              GetSlidingStepSizeCells(Int_t dim);
+  Int_t              GetSlidingStepSizeMiniPatches(Int_t dim);
+
+  Double_t           GetPatchArea(Int_t ipatch);
+  Double_t           GetPatchAreaActive(Int_t id, Int_t ipatch, Int_t type);
+
+  TString            fTracksOutName;        // name of output track array
+  TString            fTracksInName;         // name of input jet array
+  TClonesArray      *fTracksIn;             //!jet array in
+  TClonesArray      *fTracksOut;            //!track array out
+  Bool_t             fL1Slide;              // sliding window on
+
+  Double_t           fPhiMin[2];            // min phi of EMCal an DCal
+  Double_t           fPhiMax[2];            // max phi of EMCal an DCa
+  Double_t           fEtaMin[2];            // min eta of EMCal an DCal
+  Double_t           fEtaMax[2];            // max eta of EMCal an DCal
+  Double_t           fCellSize;             // size of cell (equal in eta and phi)
+  Double_t           fMinCellE;             // minimum cell energy
+  Int_t              fExclLeadingPatch;     // exclude leading patch from median calculation
+
+  Int_t              fPatchSub;             // patch type to use for subtraction
+  Double_t           fRhoMean;              // mean rho
+
+  Int_t              fNCells;               // total number of cells
+  Int_t              fNCellsEMCal;          // total number of EMCal cells
+  Int_t              fNCellsDCal;           // total number of DCal cells
+  TArrayD            fCellGrid;             // grid of cells in EMCal and DCal
+  TArrayD            fMiniPatchGrid;        // grid of mini patches in EMCal and DCal
+  TArrayI            fActiveAreaMP;         // active area for each mini patch
+  TArrayD            fPatchGrid[5];         // grid of trigger patches: 4x4 L0, 4x4 L1, 8x8 L1, 16x16 L1, 32x32 L1
+  TArrayI            fActiveAreaMPP[5];     // active area in mini patches for each trigger patch
+  TArrayI            fActiveAreaCP[5];      // active area in cells for each trigger patch
+  Int_t              fNPatchesEMCal[5];     // number of patches in EMCal
+
+ private:
+  AliEmcalPicoTrackInGridMaker(const AliEmcalPicoTrackInGridMaker&);            // not implemented
+  AliEmcalPicoTrackInGridMaker &operator=(const AliEmcalPicoTrackInGridMaker&); // not implemented
+
+  TH2F              *fh2MedianTypeEmcal[3]; //! median vs patch type for 3 types of area calculation
+  TH2F              *fh2MedianTypeDcal[3];  //! median vs patch type for 3 types of area calculation
+  TProfile          *fpMedianTypeEmcal[3];  //! median vs patch type for 3 types of area calculation
+  TProfile          *fpMedianTypeDcal[3];   //! median vs patch type for 3 types of area calculation
+  TH1F              *fh1RhoEmcal[5];        //! rho distributions for passive area
+  TH1F              *fh1RhoDcal[5];         //! rho distributions for passive area
+  TH2F              *fPatchEnVsActivityEmcal[5]; //! patch energy vs active cells
+  TH2F              *fPatchEnVsActivityDcal[5];  //! patch energy vs active cells
+
+  TH1F              *fPatchECorr[2][5];     //! corrected patch energy for EMCal and DCal
+  TH1F              *fPatchECorrPar[2][5];  //! corrected patch energy with inclusive mean rho for EMCal and DCal
+  TH2F              *fPatchECorrRho[2][5];  //! corrected patch energy vs rho opposite side
+  TH2F              *fPatchECorrRhoDijet[2][5]; //! corrected patch energy vs rho opposite side
+  TH3F              *fPatchECorrECorrRho[2][5]; //! Ecorr,det1 vs Ecorr,det2 vs rho,det2 opposite side for dijet in acceptance like events
+
+  TH2F              *fMultVsRho;            //! track multiplicity vs rho from EMCal
+
+  TList             *fHistList;             //! List of Histograms
+
+  ClassDef(AliEmcalPicoTrackInGridMaker, 1); // Task to make PicoTracks in a grid corresponding to EMCAL/DCAL acceptance
+};
+#endif
diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalHighMultTrigger.cxx b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalHighMultTrigger.cxx
new file mode 100644 (file)
index 0000000..10bdbfb
--- /dev/null
@@ -0,0 +1,314 @@
+// $Id$
+//
+// High multiplicity pp trigger analysis task.
+//
+// Author: M. Verweij
+
+#include <TClonesArray.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TList.h>
+#include <TLorentzVector.h>
+
+#include "AliLog.h"
+#include "AliEmcalTriggerPatchInfo.h"
+#include "AliParticleContainer.h"
+#include "AliVVZERO.h"
+
+#include "AliAnalysisTaskEmcalHighMultTrigger.h"
+
+ClassImp(AliAnalysisTaskEmcalHighMultTrigger)
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalHighMultTrigger::AliAnalysisTaskEmcalHighMultTrigger() : 
+  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalHighMultTrigger", kTRUE),
+  fNExLP(1),
+  fNAccPatches(-1),
+  fMedianEnergy(0),
+  fMedianEnergyExLP(0),
+  fSumEnergy(0),
+  fSumEnergyExLP(0),
+  fHistPatchEtaPhi(0),
+  fHistEnergyMedian(0),
+  fHistEnergyMedianExLP(0),
+  fHistEnergySum(0),
+  fHistEnergySumExLP(0),
+  fHistTracks(0),
+  fHistTracklets(0),
+  fHistV0MultSum(0),
+  fHistTracksTracklets(0),
+  fHistTracksV0MultSum(0)
+{
+  // Default constructor.
+
+  const Int_t nMultEst = 3;
+  for(Int_t i = 0; i<nMultEst; i++) {
+    fHistEnergyMedianEst[i] = 0;
+    fHistEnergyMedianExLPEst[i] = 0;
+    fHistEnergySumEst[i] = 0;
+    fHistEnergySumExLPEst[i] = 0;
+    fHistEnergySumAvgEst[i] = 0;
+    fHistEnergySumAvgExLPEst[i] = 0;
+  }
+
+  SetMakeGeneralHistograms(kTRUE);
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalHighMultTrigger::AliAnalysisTaskEmcalHighMultTrigger(const char *name) : 
+  AliAnalysisTaskEmcalJet(name, kTRUE),
+  fNExLP(1),
+  fNAccPatches(-1),
+  fMedianEnergy(0),
+  fMedianEnergyExLP(0),
+  fSumEnergy(0),
+  fSumEnergyExLP(0),
+  fHistPatchEtaPhi(0),
+  fHistEnergyMedian(0),
+  fHistEnergyMedianExLP(0),
+  fHistEnergySum(0),
+  fHistEnergySumExLP(0),
+  fHistTracks(0),
+  fHistTracklets(0),
+  fHistV0MultSum(0),
+  fHistTracksTracklets(0),
+  fHistTracksV0MultSum(0)
+{
+  // Standard constructor.
+
+  const Int_t nMultEst = 3;
+  for(Int_t i = 0; i<nMultEst; i++) {
+    fHistEnergyMedianEst[i] = 0;
+    fHistEnergyMedianExLPEst[i] = 0;
+    fHistEnergySumEst[i] = 0;
+    fHistEnergySumExLPEst[i] = 0;
+    fHistEnergySumAvgEst[i] = 0;
+    fHistEnergySumAvgExLPEst[i] = 0;
+  }
+
+  SetMakeGeneralHistograms(kTRUE);
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalHighMultTrigger::~AliAnalysisTaskEmcalHighMultTrigger()
+{
+  // Destructor.
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcalHighMultTrigger::UserCreateOutputObjects()
+{
+  // Create user output.
+
+  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
+
+  Int_t nE = 500;
+  Double_t minE = 0.;
+  Double_t maxE = 50.;
+
+  TString histName = "";
+  TString histTitle = "";
+
+  histName = Form("fHistPatchEtaPhi");
+  histTitle = Form("%s;#eta;#phi",histName.Data());
+  fHistPatchEtaPhi = new TH2F(histName.Data(),histTitle.Data(),100,-1.,1.,18*8,0.,TMath::TwoPi());
+  fOutput->Add(fHistPatchEtaPhi);
+
+  histName = Form("fHistEnergyMedian");
+  histTitle = Form("%s;med[#it{E}]",histName.Data());
+  fHistEnergyMedian = new TH1F(histName.Data(),histTitle.Data(),nE,minE,maxE/4.);
+  fOutput->Add(fHistEnergyMedian);
+
+  histName = Form("fHistEnergyMedianExLP");
+  histTitle = Form("%s;med[#it{E}]",histName.Data());
+  fHistEnergyMedianExLP = new TH1F(histName.Data(),histTitle.Data(),nE,minE,maxE/4.);
+  fOutput->Add(fHistEnergyMedianExLP);
+
+  histName = Form("fHistEnergySum");
+  histTitle = Form("%s;#sum[#it{E}]",histName.Data());
+  fHistEnergySum = new TH1F(histName.Data(),histTitle.Data(),nE,minE,maxE);
+  fOutput->Add(fHistEnergySum);
+
+  histName = Form("fHistEnergySumExLP");
+  histTitle = Form("%s;#sum[#it{E}]",histName.Data());
+  fHistEnergySumExLP = new TH1F(histName.Data(),histTitle.Data(),nE,minE,maxE);
+  fOutput->Add(fHistEnergySumExLP);
+
+  histName = Form("fHistTracks");
+  histTitle = Form("%s;#it{N}_{tracks}",histName.Data());
+  fHistTracks = new TH1F(histName.Data(),histTitle.Data(),300,0.,300.);
+  fOutput->Add(fHistTracks);
+
+  histName = Form("fHistTracklets");
+  histTitle = Form("%s;#it{N}_{tracklets}",histName.Data());
+  fHistTracklets = new TH1F(histName.Data(),histTitle.Data(),300,0.,300.);
+  fOutput->Add(fHistTracklets);
+
+  histName = Form("fHistV0MultSum");
+  histTitle = Form("%s;mult[V0A+V0C]",histName.Data());
+  fHistV0MultSum = new TH1F(histName.Data(),histTitle.Data(),500,0.,500.);
+  fOutput->Add(fHistV0MultSum);
+
+  const Int_t nMultEst = 3;
+  Int_t nBinsMultEst[nMultEst] = {300,300,500};
+  Double_t multEstMax[nMultEst] = {300.,300.,500.};
+  TString strMultEst[nMultEst] = {"Tracks","Tracklets","V0MultSum"};
+  for(Int_t i = 0; i<nMultEst; i++) {
+    histName = Form("fHistEnergyMedianEst%s",strMultEst[i].Data());
+    histTitle = Form("%s;med[#it{E}];%s",histName.Data(),strMultEst[i].Data());
+    fHistEnergyMedianEst[i] = new TH2F(histName.Data(),histTitle.Data(),nE,minE,maxE/4.,nBinsMultEst[i],0.,multEstMax[i]);
+    fOutput->Add(fHistEnergyMedianEst[i]);
+
+    histName = Form("fHistEnergyMedianExLPEst%s",strMultEst[i].Data());
+    histTitle = Form("%s;med[#it{E}];%s",histName.Data(),strMultEst[i].Data());
+    fHistEnergyMedianExLPEst[i] = new TH2F(histName.Data(),histTitle.Data(),nE,minE,maxE/4.,nBinsMultEst[i],0.,multEstMax[i]);
+    fOutput->Add(fHistEnergyMedianExLPEst[i]);
+
+    histName = Form("fHistEnergySumEst%s",strMultEst[i].Data());
+    histTitle = Form("%s;#sum[#it{E}];%s",histName.Data(),strMultEst[i].Data());
+    fHistEnergySumEst[i] = new TH2F(histName.Data(),histTitle.Data(),nE,minE,maxE,nBinsMultEst[i],0.,multEstMax[i]);
+    fOutput->Add(fHistEnergySumEst[i]);
+
+    histName = Form("fHistEnergySumExLPEst%s",strMultEst[i].Data());
+    histTitle = Form("%s;#sum[#it{E}];%s",histName.Data(),strMultEst[i].Data());
+    fHistEnergySumExLPEst[i] = new TH2F(histName.Data(),histTitle.Data(),nE,minE,maxE,nBinsMultEst[i],0.,multEstMax[i]);
+    fOutput->Add(fHistEnergySumExLPEst[i]);
+
+    histName = Form("fHistEnergySumAvgEst%s",strMultEst[i].Data());
+    histTitle = Form("%s;#sum[#it{E}];%s",histName.Data(),strMultEst[i].Data());
+    fHistEnergySumAvgEst[i] = new TH2F(histName.Data(),histTitle.Data(),nE,minE,maxE/4.,nBinsMultEst[i],0.,multEstMax[i]);
+    fOutput->Add(fHistEnergySumAvgEst[i]);
+
+    histName = Form("fHistEnergySumAvgExLPEst%s",strMultEst[i].Data());
+    histTitle = Form("%s;#sum[#it{E}];%s",histName.Data(),strMultEst[i].Data());
+    fHistEnergySumAvgExLPEst[i] = new TH2F(histName.Data(),histTitle.Data(),nE,minE,maxE/4.,nBinsMultEst[i],0.,multEstMax[i]);
+    fOutput->Add(fHistEnergySumAvgExLPEst[i]);
+  }
+
+  histName = Form("fHistTracksTracklets");
+  histTitle = Form("%s;%s;%s",histName.Data(),strMultEst[0].Data(),strMultEst[1].Data());
+  fHistTracksTracklets = new TH2F(histName.Data(),histTitle.Data(),nBinsMultEst[0],0.,multEstMax[0],nBinsMultEst[1],0.,multEstMax[1]);
+  fOutput->Add(fHistTracksTracklets);
+
+  histName = Form("fHistTracksV0MultSum");
+  histTitle = Form("%s;%s;%s",histName.Data(),strMultEst[0].Data(),strMultEst[2].Data());
+  fHistTracksV0MultSum = new TH2F(histName.Data(),histTitle.Data(),nBinsMultEst[0],0.,multEstMax[0],nBinsMultEst[2],0.,multEstMax[2]);
+  fOutput->Add(fHistTracksV0MultSum);
+
+  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalHighMultTrigger::FillHistograms()
+{
+  // Fill histograms.
+
+  fHistEnergyMedian->Fill(fMedianEnergy);
+  fHistEnergyMedianExLP->Fill(fMedianEnergyExLP);
+  fHistEnergySum->Fill(fSumEnergy);
+  fHistEnergySumExLP->Fill(fSumEnergyExLP);
+
+
+  //Multiplicity estimators
+  Int_t nTracks   = -1;
+  if (GetParticleContainer(0)) nTracks = GetParticleContainer(0)->GetNAcceptedParticles();
+  Int_t nTracklets = InputEvent()->GetMultiplicity()->GetNumberOfTracklets();
+
+  AliVVZERO* vV0 = InputEvent()->GetVZEROData();
+  Float_t multV0A=vV0->GetMTotV0A();
+  Float_t multV0C=vV0->GetMTotV0C();
+
+  fHistTracks->Fill(nTracks);
+  fHistTracklets->Fill(nTracklets);
+  fHistV0MultSum->Fill(multV0A+multV0C);
+
+  const Int_t nMultEst = 3;
+  Float_t multEst[nMultEst] = {(Float_t)(nTracks),(Float_t)(nTracklets),(Float_t)(multV0A+multV0C)};
+  for(Int_t i = 0; i<nMultEst; i++) {
+    fHistEnergyMedianEst[i]->Fill(fMedianEnergy,multEst[i]);
+    fHistEnergyMedianExLPEst[i]->Fill(fMedianEnergyExLP,multEst[i]);
+    fHistEnergySumEst[i]->Fill(fSumEnergy,multEst[i]);
+    fHistEnergySumExLPEst[i]->Fill(fSumEnergyExLP,multEst[i]);
+    if(fNAccPatches>0)
+      fHistEnergySumAvgEst[i]->Fill(fSumEnergy/((Double_t)fNAccPatches),multEst[i]);
+    if(fNAccPatches>1)
+      fHistEnergySumAvgExLPEst[i]->Fill(fSumEnergyExLP/((Double_t)fNAccPatches-1.),multEst[i]);
+  }
+
+  fHistTracksTracklets->Fill(multEst[0],multEst[1]);
+  fHistTracksV0MultSum->Fill(multEst[0],multEst[2]);
+
+  return kTRUE;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcalHighMultTrigger::ExecOnce() {
+
+  AliAnalysisTaskEmcalJet::ExecOnce();
+
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalHighMultTrigger::Run()
+{
+  // Run analysis code here, if needed. It will be executed before FillHistograms().
+
+  if(!fTriggerPatchInfo) {
+    AliDebug(11,Form("Couldn't find patch info object %s",fCaloTriggerPatchInfoName.Data()));
+    return kFALSE;
+  }
+
+  Int_t nPatch = fTriggerPatchInfo->GetEntriesFast();
+
+  //sort patches
+  Double_t ptarr[999] = {0};
+  Int_t indexes[999] = {0};
+  Int_t iacc = 0;
+  for(Int_t i = 0; i<nPatch; i++) {
+    AliEmcalTriggerPatchInfo *patch = dynamic_cast<AliEmcalTriggerPatchInfo*>(fTriggerPatchInfo->At(i));
+    if(!patch) continue;
+    if(patch->GetPatchE()>0.) {
+      ptarr[iacc] = patch->GetPatchE();
+      iacc++;
+    }
+    fHistPatchEtaPhi->Fill(patch->GetEtaMin(),patch->GetPhiMin());
+  }
+  
+  TMath::Sort(nPatch,ptarr,indexes);
+  Double_t ptarrSort[999];
+  // Double_t indexesSort[999];
+  Double_t ptarrSortExLP[999];
+  //  Double_t indexesSortExLP[999];
+  for(Int_t i = 0; i<iacc; i++) {
+    ptarrSort[i] = ptarr[indexes[i]];
+    //  indexesSort[i] = indexes[i];
+    if(i>=fNExLP) {
+      ptarrSortExLP[i] = ptarr[indexes[i]];
+      //  indexesSortExLP[i] = indexes[i];
+    }
+  }
+
+  //Calculate sum energy 
+  fSumEnergy = 0.;
+  fSumEnergyExLP = 0.;
+  for(Int_t i = 0; i<iacc; i++) {
+    fSumEnergy+= ptarr[indexes[i]];
+    if(i>=fNExLP)
+      fSumEnergyExLP+= ptarr[indexes[i]];
+  }
+
+  //calculate median of all and excluding leading patch(es)
+  fMedianEnergy = TMath::Median(iacc,ptarrSort);
+  fMedianEnergyExLP = TMath::Median(iacc-fNExLP,ptarrSortExLP);
+  fNAccPatches = iacc;
+
+   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcalHighMultTrigger::Terminate(Option_t *) 
+{
+  // Called once at the end of the analysis.
+}
+
diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalHighMultTrigger.h b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalHighMultTrigger.h
new file mode 100644 (file)
index 0000000..de973e3
--- /dev/null
@@ -0,0 +1,65 @@
+#ifndef ALIANALYSISTASKEMCALHIGHMULTTRIGGER_H
+#define ALIANALYSISTASKEMCALHIGHMULTTRIGGER_H
+
+// $Id$
+
+class TH1;
+class TH2;
+class TH3;
+
+#include "AliAnalysisTaskEmcalJet.h"
+
+class AliAnalysisTaskEmcalHighMultTrigger : public AliAnalysisTaskEmcalJet {
+ public:
+
+  AliAnalysisTaskEmcalHighMultTrigger();
+  AliAnalysisTaskEmcalHighMultTrigger(const char *name);
+  virtual ~AliAnalysisTaskEmcalHighMultTrigger();
+
+  void                        UserCreateOutputObjects();
+  void                        Terminate(Option_t *option);
+
+  //Setters
+  void                        SetNExcludeLeadingPatches(Int_t n)  {fNExLP = n; }
+
+ protected:
+  void                        ExecOnce();
+  Bool_t                      FillHistograms()   ;
+  Bool_t                      Run()              ;
+  
+ private:
+  Int_t                       fNExLP;                  //nr of leading patched to exclude from estimate
+  Int_t                       fNAccPatches;            //nr of accepted patches
+  Double_t                    fMedianEnergy;           //median event energy
+  Double_t                    fMedianEnergyExLP;       //median event energy
+  Double_t                    fSumEnergy;              //summed energy
+  Double_t                    fSumEnergyExLP;          //summed energy
+
+
+  //Histograms
+  TH2F                       *fHistPatchEtaPhi;         //!
+  TH1F                       *fHistEnergyMedian;        //!
+  TH1F                       *fHistEnergyMedianExLP;    //!
+  TH1F                       *fHistEnergySum;           //!
+  TH1F                       *fHistEnergySumExLP;       //!
+
+  TH1F                       *fHistTracks;              //!
+  TH1F                       *fHistTracklets;           //!
+  TH1F                       *fHistV0MultSum;           //!
+
+  TH2F                       *fHistEnergyMedianEst[3];        //!
+  TH2F                       *fHistEnergyMedianExLPEst[3];    //!
+  TH2F                       *fHistEnergySumEst[3];           //!
+  TH2F                       *fHistEnergySumExLPEst[3];       //!
+  TH2F                       *fHistEnergySumAvgEst[3];        //!
+  TH2F                       *fHistEnergySumAvgExLPEst[3];    //!
+
+  TH2F                       *fHistTracksTracklets;           //!
+  TH2F                       *fHistTracksV0MultSum;           //!
+
+  AliAnalysisTaskEmcalHighMultTrigger(const AliAnalysisTaskEmcalHighMultTrigger&);            // not implemented
+  AliAnalysisTaskEmcalHighMultTrigger &operator=(const AliAnalysisTaskEmcalHighMultTrigger&); // not implemented
+
+  ClassDef(AliAnalysisTaskEmcalHighMultTrigger, 1) // high multiplicity pp trigger analysis task
+};
+#endif
diff --git a/PWGJE/EMCALJetTasks/macros/AddTaskEmcalHighMultTrigger.C b/PWGJE/EMCALJetTasks/macros/AddTaskEmcalHighMultTrigger.C
new file mode 100644 (file)
index 0000000..07e0e7f
--- /dev/null
@@ -0,0 +1,63 @@
+// $Id$
+
+AliAnalysisTaskEmcalHighMultTrigger* AddTaskEmcalHighMultTrigger(
+  const char *ntracks            = "Tracks",
+  const char *nPatches           = "EmcalPatches32x32",
+  Int_t       nCentBins          = 1,
+  const char *taskname           = "HighMultTrigger"
+)
+{
+  // Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr)
+  {
+    ::Error("AddTaskEmcalHighMultTrigger", "No analysis manager to connect to.");
+    return NULL;
+  }
+
+  // Check the analysis type using the event handlers connected to the analysis manager.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler())
+  {
+    ::Error("AddTaskEmcalHighMultTrigger", "This task requires an input event handler");
+    return NULL;
+  }
+
+  //-------------------------------------------------------
+  // Init the task and do settings
+  //-------------------------------------------------------
+
+  TString name = Form("%s_%s",taskname,nPatches);
+
+  Printf("name: %s",name.Data());
+
+  AliAnalysisTaskEmcalHighMultTrigger* task = new AliAnalysisTaskEmcalHighMultTrigger(name);
+  task->SetCentRange(0.,100.);
+  task->SetNCentBins(nCentBins);
+  task->SetCaloTriggerPatchInfoName(nPatches);
+  task->SetVzRange(-10.,10.);
+  task->SetUseAliAnaUtils(kTRUE);
+
+  AliParticleContainer *trackCont  = task->AddParticleContainer(ntracks);
+  trackCont->SetClassName("AliVTrack");
+
+  //-------------------------------------------------------
+  // Final settings, pass to manager and set the containers
+  //-------------------------------------------------------
+
+  mgr->AddTask(task);
+
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput1  = mgr->GetCommonInputContainer()  ;
+  TString contname(name);
+  contname += "_histos";
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(),
+                                                           TList::Class(),AliAnalysisManager::kOutputContainer,
+                                                           Form("%s", AliAnalysisManager::GetCommonFileName()));
+  mgr->ConnectInput  (task, 0,  cinput1 );
+  mgr->ConnectOutput (task, 1, coutput1 );
+
+  return task;
+}
+
diff --git a/PWGJE/EMCALJetTasks/macros/AddTaskGridMaker.C b/PWGJE/EMCALJetTasks/macros/AddTaskGridMaker.C
new file mode 100644 (file)
index 0000000..7e949b0
--- /dev/null
@@ -0,0 +1,52 @@
+// $Id$
+
+AliEmcalPicoTrackInGridMaker* AddTaskGridMaker(
+                                              const char *inname       = "PicoTracks",
+                                              const char *taskName     = "AliEmcalPicoTrackInGridMaker"
+)
+{  
+  // Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr)
+  {
+    ::Error("AddTaskEmcalPicoTrackFromJetMaker", "No analysis manager to connect to.");
+    return NULL;
+  }  
+  
+  // Check the analysis type using the event handlers connected to the analysis manager.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler())
+  {
+    ::Error("AddTaskEmcalPicoTrackFromJetMaker", "This task requires an input event handler");
+    return NULL;
+  }
+
+  TString wagonName = Form("%s",taskName);
+
+  
+  //-------------------------------------------------------
+  // Init the task and do settings
+  //-------------------------------------------------------
+
+  AliEmcalPicoTrackInGridMaker *eTask = new AliEmcalPicoTrackInGridMaker(taskName);
+  eTask->SetTracksInName(inname);
+
+  //-------------------------------------------------------
+  // Final settings, pass to manager and set the containers
+  //-------------------------------------------------------
+  mgr->AddTask(eTask);
+  
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput1  = mgr->GetCommonInputContainer();
+  mgr->ConnectInput(eTask, 0, cinput1 );
+
+  //Connect output
+  TString contName(wagonName);
+  TString outputfile = Form("%s",AliAnalysisManager::GetCommonFileName());
+  Printf("outputfile: %s",outputfile.Data());
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contName.Data(), TList::Class(),AliAnalysisManager::kOutputContainer,outputfile);
+  mgr->ConnectOutput(eTask,1,coutput1);
+  
+  return eTask;
+}
index 4d5579c..8e273fd 100644 (file)
@@ -16,6 +16,7 @@
 #pragma link C++ class AliAnalysisTaskDeltaPt+;
 #pragma link C++ class AliAnalysisTaskScale+;
 #pragma link C++ class AliEmcalJet+;
+#pragma link C++ class AliEmcalPicoTrackInGridMaker+;
 #pragma link C++ class AliJetContainer+;
 #pragma link C++ class AliJetEmbeddingTask+;
 #pragma link C++ class AliJetEmbeddingFromGenTask+;
@@ -42,6 +43,7 @@
 #pragma link C++ class AliAnalysisTaskEmcalDiJetBase+;
 #pragma link C++ class AliAnalysisTaskEmcalDiJetAna+;
 #pragma link C++ class AliAnalysisTaskEmcalDiJetResponse+;
+#pragma link C++ class AliAnalysisTaskEmcalHighMultTrigger+;
 #pragma link C++ class AliAnalysisTaskEmcalJetHMEC+;
 #pragma link C++ class AliAnalysisTaskEmcalJetHadCorQA+;
 #pragma link C++ class AliAnalysisTaskEmcalJetHadEPpid+;