PWGJE/EMCALJetTasks/Tracks/AliAnalysisTaskPtEMCalTrigger.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2014, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17  * Analysis task of the pt analysis on EMCal-triggered events
  18  *
  19  *   Author: Markus Fasel
  20  */
  21
  22 #include <map>
  23 #include <cstring>
  24 #include <iostream>
  25 #include <memory>
  26 #include <vector>
  27 #include <string>
  28 #include <sstream>
  29
  30 #include <TClonesArray.h>
  31 #include <TDirectory.h>
  32 #include <TH1.h>
  33 #include <THashList.h>
  34 #include <TList.h>
  35 #include <TKey.h>
  36 #include <TMath.h>
  37 #include <TObjArray.h>
  38 #include <TObjString.h>
  39 #include <TString.h>
  40 #include <TVector2.h>
  41
  42 #include "AliAODEvent.h"
  43 #include "AliESDEvent.h"
  44 #include "AliInputEventHandler.h"
  45 #include "AliMCEvent.h"
  46 #include "AliParticleContainer.h"
  47 #include "AliVCluster.h"
  48 #include "AliVParticle.h"
  49 #include "AliVTrack.h"
  50 #include "AliVVertex.h"
  51
  52 #include "AliClusterContainer.h"
  53 #include "AliEmcalJet.h"
  54 #include "AliEmcalTriggerPatchInfo.h"
  55 #include "AliEmcalPhysicsSelection.h"
  56 #include "AliAnalysisTaskPtEMCalTrigger.h"
  57 #include "AliEMCalHistoContainer.h"
  58 #include "AliEMCalPtTaskVTrackSelection.h"
  59 #include "AliEMCalPtTaskTrackSelectionESD.h"
  60 #include "AliEMCalPtTaskTrackSelectionAOD.h"
  61 #include "AliJetContainer.h"
  62 #include "AliParticleContainer.h"
  63 #include "AliPicoTrack.h"
  64
  65 ClassImp(EMCalTriggerPtAnalysis::AliAnalysisTaskPtEMCalTrigger)
  66
  67 namespace EMCalTriggerPtAnalysis {
  68
  69   /*
  70    * constants
  71    */
  72   const Int_t AliAnalysisTaskPtEMCalTrigger::kNJetRadii = 4;
  73   const Double_t jetRadVals[AliAnalysisTaskPtEMCalTrigger::kNJetRadii] = {0.2, 0.3, 0.4, 0.5};
  74   const Double_t *AliAnalysisTaskPtEMCalTrigger::kJetRadii = jetRadVals;
  75
  76   //______________________________________________________________________________
  77   AliAnalysisTaskPtEMCalTrigger::AliAnalysisTaskPtEMCalTrigger():
  78                                     AliAnalysisTaskEmcalJet(),
  79                                     fHistos(NULL),
  80                                     fListTrackCuts(NULL),
  81                                     fEtaRange(),
  82                                     fPtRange(),
  83                                     fEnergyRange(),
  84                                     fVertexRange(),
  85                                     fJetContainersMC(),
  86                                     fJetContainersData(),
  87                                     fSelectAllTracks(kFALSE),
  88                                     fSwapEta(kFALSE),
  89                                     fUseTriggersFromTriggerMaker(kFALSE)
  90   {
  91     /*
  92      * Dummy constructor, initialising the values with default (NULL) values
  93      */
  94   }
  95
  96   //______________________________________________________________________________
  97   AliAnalysisTaskPtEMCalTrigger::AliAnalysisTaskPtEMCalTrigger(const char *name):
  98                                     AliAnalysisTaskEmcalJet(name, kTRUE),
  99                                     fHistos(NULL),
 100                                     fListTrackCuts(NULL),
 101                                     fEtaRange(),
 102                                     fPtRange(),
 103                         fEnergyRange(),
 104                         fVertexRange(),
 105                         fJetContainersMC(),
 106                         fJetContainersData(),
 107                         fSelectAllTracks(kFALSE),
 108                                     fSwapEta(kFALSE),
 109                                     fUseTriggersFromTriggerMaker(kFALSE)
 110   {
 111     /*
 112      * Main constructor, setting default values for eta and zvertex cut
 113      */
 114
 115     fListTrackCuts = new TList;
 116     fListTrackCuts->SetOwner(false);
 117
 118     // Set default cuts
 119     fEtaRange.SetLimits(-0.8, 0.8);
 120     fPtRange.SetLimits(0.15, 100.);
 121     fEnergyRange.SetLimits(0., 1000.);
 122     fVertexRange.SetLimits(-10., 10.);
 123     SetMakeGeneralHistograms(kTRUE);
 124   }
 125
 126   //______________________________________________________________________________
 127   AliAnalysisTaskPtEMCalTrigger::~AliAnalysisTaskPtEMCalTrigger(){
 128     /*
 129      * Destructor, deleting output
 130      */
 131     //if(fTrackSelection) delete fTrackSelection;
 132     if(fHistos) delete fHistos;
 133     if(fListTrackCuts) delete fListTrackCuts;
 134   }
 135
 136   //______________________________________________________________________________
 137   void AliAnalysisTaskPtEMCalTrigger::UserCreateOutputObjects(){
 138     /*
 139      * Create the list of output objects and define the histograms.
 140      * Also adding the track cuts to the list of histograms.
 141      */
 142     AliAnalysisTaskEmcal::UserCreateOutputObjects();
 143     SetCaloTriggerPatchInfoName("EmcalTriggers");
 144     fHistos = new AliEMCalHistoContainer("PtEMCalTriggerHistograms");
 145     fHistos->ReleaseOwner();
 146
 147     if(fJetContainersMC.GetEntries()){
 148       AliDebug(1,"Jet containers for MC truth:");
 149       TObjString *contname(NULL);
 150       TIter contMCIter(&fJetContainersMC);
 151       while((contname = dynamic_cast<TObjString *>(contMCIter.Next())))
 152         AliDebug(1, Form("Next container: %s", contname->String().Data()));
 153     }
 154     if(fJetContainersData.GetEntries()){
 155       AliDebug(1, "Jet containers for Data:");
 156       TObjString *contname(NULL);
 157       TIter contDataIter(&fJetContainersData);
 158       while((contname = dynamic_cast<TObjString *>(contDataIter.Next())))
 159         AliDebug(1, Form("Next container: %s", contname->String().Data()));
 160     }
 161     if(fJetCollArray.GetEntries()){
 162       AliDebug(1, "Jet containers attached to this task:");
 163       AliJetContainer *cont(NULL);
 164       TIter contIter(&fJetCollArray);
 165       while((cont = dynamic_cast<AliJetContainer *>(contIter.Next())))
 166         AliDebug(1, Form("Container: %s", cont->GetName()));
 167     }
 168
 169     std::map<std::string, std::string> triggerCombinations;
 170     const char *triggernames[12] = {"MinBias", "EMCJHigh", "EMCJLow", "EMCGHigh", "EMCGLow", "NoEMCal", "EMCHighBoth", "EMCHighGammaOnly", "EMCHighJetOnly", "EMCLowBoth", "EMCLowGammaOnly", "EMCLowJetOnly"},
 171         *bitnames[4] = {"CINT7", "EMC7", "kEMCEGA", "kEMCEJE"};
 172     // Define axes for the trigger correlation histogram
 173     const TAxis *triggeraxis[5]; memset(triggeraxis, 0, sizeof(const TAxis *) * 5);
 174     const TAxis *bitaxes[4]; memset(bitaxes, 0, sizeof(TAxis *) * 4);
 175     const char *binlabels[2] = {"OFF", "ON"};
 176     TAxis mytrgaxis[5], mybitaxis[4];
 177     for(int itrg = 0; itrg < 5; ++itrg){
 178       DefineAxis(mytrgaxis[itrg], triggernames[itrg], triggernames[itrg], 2, -0.5, 1.5, binlabels);
 179       triggeraxis[itrg] = mytrgaxis+itrg;
 180       if(itrg < 4){
 181         DefineAxis(mybitaxis[itrg], bitnames[itrg], bitnames[itrg], 2, -0.5, 1.5, binlabels);
 182         bitaxes[itrg] = mybitaxis+itrg;
 183       }
 184     }
 185     // Define names and titles for different triggers in the histogram container
 186     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[0], "min. bias events"));
 187     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[1], "jet-triggered events (high threshold)"));
 188     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[2], "jet-triggered events (low threshold)"));
 189     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[3], "gamma-triggered events (high threshold)"));
 190     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[4], "gamma-triggered events (low threshold)"));
 191     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[5], "non-EMCal-triggered events"));
 192     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[6], "jet and gamma triggered events (high threshold)"));
 193     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[7], "exclusively gamma-triggered events (high threshold)"));
 194     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[8], "exclusively jet-triggered events (high threshold)"));
 195     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[9], "jet and gamma triggered events (low threshold)"));
 196     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[10], "exclusively gamma-triggered events (low threshold)"));
 197     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[11], "exclusively-triggered events (low threshold)"));
 198     // Define axes for the pt histogram
 199     // Dimensions:
 200     // 1. pt
 201     // 2. eta
 202     // 3. phi
 203     // 4. vertex
 204     // 5. pileup (0 = all events, 1 = after pileup rejection)
 205     // 6. track cuts (0 = no cuts; 1 = after std cuts)
 206     TArrayD ptbinning, zvertexBinning, etabinning, pileupaxis(3);
 207     pileupaxis[0] = -0.5; pileupaxis[1] = 0.5; pileupaxis[2] = 1.5;
 208     CreateDefaultPtBinning(ptbinning);
 209     CreateDefaultZVertexBinning(zvertexBinning);
 210     CreateDefaultEtaBinning(etabinning);
 211     TAxis htrackaxes[7];
 212     DefineAxis(htrackaxes[0], "pt", "p_{t} (GeV/c)", ptbinning);
 213     DefineAxis(htrackaxes[1], "eta", "#eta", etabinning);
 214     DefineAxis(htrackaxes[2], "phi", "#phi", 20, 0, 2 * TMath::Pi());
 215     DefineAxis(htrackaxes[3], "zvertex", "z_{V} (cm)", zvertexBinning);
 216     DefineAxis(htrackaxes[4], "pileup", "Pileup rejection", 2, -0.5, 1.5);
 217     DefineAxis(htrackaxes[5], "trackcuts", "Track Cuts", (fListTrackCuts ? fListTrackCuts->GetEntries() : 0) + 1, -0.5, (fListTrackCuts ? fListTrackCuts->GetEntries() : 0) + 0.5);
 218     DefineAxis(htrackaxes[6], "mbtrigger", "Has MB trigger", 2, -0.5, 1.5);
 219     const TAxis *trackaxes[7];
 220     for(int iaxis = 0; iaxis < 7; ++iaxis) trackaxes[iaxis] = htrackaxes + iaxis;
 221     TAxis hclusteraxes[4];
 222     DefineAxis(hclusteraxes[0], "energy", "E (GeV)", ptbinning);
 223     DefineAxis(hclusteraxes[1], "zvertex", "z_{V} (cm)", zvertexBinning);
 224     DefineAxis(hclusteraxes[2], "pileup", "Pileup rejection", 2,1 -0.5, 1.5);
 225     DefineAxis(hclusteraxes[3], "mbtrigger", "Has MB trigger", 2, -0.5, 1.5);
 226     const TAxis *clusteraxes[4];
 227     for(int iaxis = 0; iaxis < 4; ++iaxis) clusteraxes[iaxis] = hclusteraxes + iaxis;
 228     TAxis hpatchenergyaxes[5];
 229     DefineAxis(hpatchenergyaxes[0], "energy", "Patch energy (GeV)", 100, 0., 100);
 230     DefineAxis(hpatchenergyaxes[1], "eta", "#eta", etabinning);
 231     DefineAxis(hpatchenergyaxes[2], "phi", "#phi",  20, 0, 2 * TMath::Pi());
 232     DefineAxis(hpatchenergyaxes[3], "isMain", "Main trigger", 2, -0.5, 1.5);
 233     DefineAxis(hpatchenergyaxes[4],  "emcalgood", "EMCAL good event", 2, -0.5, 1.5);
 234     const TAxis *patchenergyaxes[5];
 235     for(int iaxis = 0; iaxis < 5; ++iaxis) patchenergyaxes[iaxis] = hpatchenergyaxes + iaxis;
 236     TAxis hpatchampaxes[5];
 237     DefineAxis(hpatchampaxes[0], "amplitude", "Patch energy (GeV)", 10000, 0., 10000.);
 238     DefineAxis(hpatchampaxes[1], "eta", "#eta", etabinning);
 239     DefineAxis(hpatchampaxes[2], "phi", "#phi",  20, 0, 2 * TMath::Pi());
 240     DefineAxis(hpatchampaxes[3], "isMain", "Main trigger", 2, -0.5, 1.5);
 241     DefineAxis(hpatchampaxes[4],  "emcalgood", "EMCAL good event", 2, -0.5, 1.5);
 242     const TAxis *patchampaxes[5];
 243     for(int iaxis = 0; iaxis < 5; ++iaxis) patchampaxes[iaxis] = hpatchampaxes + iaxis;
 244     try{
 245       std::string patchnames[] = {"Level0", "JetHigh", "JetLow", "GammaHigh", "GammaLow"};
 246       for(std::string * triggerpatch = patchnames; triggerpatch < patchnames + sizeof(patchnames)/sizeof(std::string); ++triggerpatch){
 247         fHistos->CreateTHnSparse(Form("Energy%s", triggerpatch->c_str()), Form("Patch energy for %s trigger patches", triggerpatch->c_str()), 5, patchenergyaxes, "s");
 248         fHistos->CreateTHnSparse(Form("EnergyRough%s", triggerpatch->c_str()), Form("Rough patch energy for %s trigger patches", triggerpatch->c_str()), 5, patchenergyaxes, "s");
 249         fHistos->CreateTHnSparse(Form("Amplitude%s", triggerpatch->c_str()), Form("Patch amplitude for %s trigger patches", triggerpatch->c_str()), 5, patchampaxes, "s");
 250       }
 251
 252       // Create histogram for MC-truth
 253       fHistos->CreateTHnSparse("hMCtrueParticles", "Particle-based histogram for MC-true particles", 5, trackaxes, "s");
 254       if(fJetCollArray.GetEntries()){
 255         for(int irad = 0; irad < kNJetRadii; irad++){
 256           fHistos->CreateTHnSparse(Form("hMCtrueParticlesRad%02d", int(kJetRadii[irad]*10)),
 257               Form("Particle-based histogram for MC-true particles in Jets with radius %.1f", kJetRadii[irad]*10), 5, trackaxes, "s");
 258         }
 259         // histogram for isolated particles
 260         fHistos->CreateTHnSparse("hMCtrueParticlesIsolated", "Particle-based histogram for isolated MC-true particles", 5, trackaxes, "s");
 261       }
 262       for(std::map<std::string,std::string>::iterator it = triggerCombinations.begin(); it != triggerCombinations.end(); ++it){
 263         const std::string name = it->first, &title = it->second;
 264         // Create event-based histogram
 265         fHistos->CreateTH2(Form("hEventHist%s", name.c_str()), Form("Event-based data for %s events; pileup rejection; z_{V} (cm)", title.c_str()), pileupaxis, zvertexBinning);
 266         // Create track-based histogram
 267         fHistos->CreateTHnSparse(Form("hTrackHist%s", name.c_str()), Form("Track-based data for %s events", title.c_str()), 7, trackaxes, "s");
 268         fHistos->CreateTHnSparse(Form("hTrackInAcceptanceHist%s", name.c_str()), Form("Track-based data for %s events  for tracks matched to EMCal clusters", title.c_str()), 7, trackaxes, "s");
 269         fHistos->CreateTHnSparse(Form("hMCTrackHist%s", name.c_str()), Form("Track-based data for %s events with MC kinematics", title.c_str()), 7, trackaxes, "s");
 270         fHistos->CreateTHnSparse(Form("hMCTrackInAcceptanceHist%s", name.c_str()), Form("Track-based data for %s events with MC kinematics for tracks matched to EMCal clusters", title.c_str()), 7, trackaxes, "s");
 271         // Create cluster-based histogram (Uncalibrated and calibrated clusters)
 272         fHistos->CreateTHnSparse(Form("hClusterCalibHist%s", name.c_str()), Form("Calib. cluster-based histogram for %s events", title.c_str()), 4, clusteraxes, "s");
 273         fHistos->CreateTHnSparse(Form("hClusterUncalibHist%s", name.c_str()), Form("Uncalib. cluster-based histogram for %s events", title.c_str()), 4, clusteraxes, "s");
 274         if(fJetCollArray.GetEntries()){
 275           // Create histograms for particles in different jetfinder radii, only track-based
 276           for(int irad = 0; irad < kNJetRadii; irad++){
 277             fHistos->CreateTHnSparse(Form("hTrackHist%sRad%02d", name.c_str(), int(kJetRadii[irad]*10)),
 278                 Form("Track-based data for %s events for tracks in jets with Radius %.1f", title.c_str(), kJetRadii[irad]), 7, trackaxes, "s");
 279             fHistos->CreateTHnSparse(Form("hTrackInAcceptanceHist%sRad%02d", name.c_str(), int(kJetRadii[irad]*10)),
 280                 Form("Track-based data for %s events for tracks matched to EMCal clusters in jets with Radius %.1f", title.c_str(), kJetRadii[irad]),
 281                 7, trackaxes, "s");
 282             fHistos->CreateTHnSparse(Form("hMCTrackHist%sRad%02d", name.c_str(), int(kJetRadii[irad]*10)),
 283                 Form("Track-based data for %s events with MC kinematics for tracks in jets with Radius %.1f", title.c_str(), kJetRadii[irad]),
 284                 7, trackaxes, "s");
 285             fHistos->CreateTHnSparse(Form("hMCTrackInAcceptanceHist%sRad%02d", name.c_str(), int(kJetRadii[irad]*10)),
 286                 Form("Track-based data for %s events with MC kinematics for tracks matched to EMCal clusters in jets with Radius %.1f", title.c_str(), kJetRadii[irad]),
 287                 7, trackaxes, "s");
 288             fHistos->CreateTHnSparse(Form("hClusterCalibHist%sRad%02d", name.c_str(), int(kJetRadii[irad]*10)),
 289                 Form("Calib. cluster-based histogram for %s events for clusters in Jets with radius %.1f", title.c_str(), kJetRadii[irad]), 4, clusteraxes, "s");
 290           }
 291         }
 292         // Create also histograms for isolated particles
 293         fHistos->CreateTHnSparse(Form("hTrackHist%sIsolated", name.c_str()), Form("Track-based data for %s events for isolated tracks", title.c_str()), 7, trackaxes, "s");
 294         fHistos->CreateTHnSparse(Form("hTrackInAcceptanceHist%sIsolated", name.c_str()), Form("Track-based data for %s events for isolated tracks matched to EMCal clusters", title.c_str()), 7, trackaxes, "s");
 295         fHistos->CreateTHnSparse(Form("hMCTrackHist%sIsolated", name.c_str()), Form("Track-based data for %s events with MC kinematics for isolated tracks", title.c_str()), 7, trackaxes, "s");
 296         fHistos->CreateTHnSparse(Form("hMCTrackInAcceptanceHist%sIsolated", name.c_str()), Form("Track-based data for %s events with MC kinematics for isolated tracks matched to EMCal clusters", title.c_str()), 7, trackaxes, "s");
 297       }
 298       fHistos->CreateTHnSparse("hEventTriggers", "Trigger type per event", 5, triggeraxis);
 299       fHistos->CreateTHnSparse("hEventsTriggerbit", "Trigger bits for the different events", 4, bitaxes);
 300     } catch (HistoContainerContentException &e){
 301       std::stringstream errormessage;
 302       errormessage << "Creation of histogram failed: " << e.what();
 303       AliError(errormessage.str().c_str());
 304     }
 305     fOutput->Add(fHistos->GetListOfHistograms());
 306     if(fListTrackCuts && fListTrackCuts->GetEntries()){
 307       TIter cutIter(fListTrackCuts);
 308       AliEMCalPtTaskVTrackSelection *cutObject(NULL);
 309       while((cutObject = dynamic_cast<AliEMCalPtTaskVTrackSelection *>(cutIter()))){
 310         AliESDtrackCuts *cuts = dynamic_cast<AliESDtrackCuts *>(cutObject->GetTrackCuts());
 311         if(cuts){
 312           cuts->DefineHistograms();
 313           fOutput->Add(cuts);
 314         }
 315       }
 316     }
 317     PostData(1, fOutput);
 318   }
 319
 320   //______________________________________________________________________________
 321   Bool_t AliAnalysisTaskPtEMCalTrigger::Run(){
 322     /*
 323      * Runs the event loop
 324      *
 325      * @param option: Additional options
 326      */
 327     // Common checks: Have SPD vertex and primary vertex from tracks, and both need to have at least one contributor
 328     AliDebug(1,Form("Number of calibrated clusters: %d", fCaloClusters->GetEntries()));
 329     AliDebug(1,Form("Number of matched tracks: %d", fTracks->GetEntries()));
 330     if(fMCEvent){
 331       // Build always trigger strig from the trigger maker in case of MC
 332       fUseTriggersFromTriggerMaker = kTRUE;
 333     }
 334
 335     Bool_t emcalGood = fInputHandler->IsEventSelected() & AliEmcalPhysicsSelection::kEmcalOk;
 336
 337     // Loop over trigger patches, fill patch energy
 338     AliEmcalTriggerPatchInfo *triggerpatch(NULL);
 339     TIter patchIter(this->fTriggerPatchInfo);
 340     while((triggerpatch = dynamic_cast<AliEmcalTriggerPatchInfo *>(patchIter()))){
 341       double triggerpatchinfo[5] = {triggerpatch->GetPatchE(), triggerpatch->GetEtaGeo(), triggerpatch->GetPhiGeo(), triggerpatch->IsMainTrigger() ? 1. : 0., emcalGood ? 1. : 0.};
 342       double triggerpatchinfoamp[5] = {static_cast<double>(triggerpatch->GetADCAmp()), triggerpatch->GetEtaGeo(), triggerpatch->GetPhiGeo(), triggerpatch->IsMainTrigger() ? 1. : 0., emcalGood ? 1. : 0.};
 343       double triggerpatchinfoer[5] = {triggerpatch->GetADCAmpGeVRough(), triggerpatch->GetEtaGeo(), triggerpatch->GetPhiGeo(), triggerpatch->IsMainTrigger() ? 1. : 0., emcalGood ? 1. : 0.};
 344       if(triggerpatch->IsJetHigh()){
 345         fHistos->FillTHnSparse("EnergyJetHigh", triggerpatchinfo);
 346         fHistos->FillTHnSparse("AmplitudeJetHigh", triggerpatchinfoamp);
 347         fHistos->FillTHnSparse("EnergyRoughJetHigh", triggerpatchinfoer);
 348       }
 349       if(triggerpatch->IsJetLow()){
 350         fHistos->FillTHnSparse("EnergyJetLow", triggerpatchinfo);
 351         fHistos->FillTHnSparse("AmplitudeJetLow", triggerpatchinfoamp);
 352         fHistos->FillTHnSparse("EnergyRoughJetLow", triggerpatchinfoer);
 353       }
 354       if(triggerpatch->IsGammaHigh()){
 355         fHistos->FillTHnSparse("EnergyGammaHigh", triggerpatchinfo);
 356         fHistos->FillTHnSparse("AmplitudeGammaHigh", triggerpatchinfoamp);
 357         fHistos->FillTHnSparse("EnergyRoughGammaHigh", triggerpatchinfoer);
 358       }
 359       if(triggerpatch->IsGammaLow()){
 360         fHistos->FillTHnSparse("EnergyGammaLow", triggerpatchinfo);
 361         fHistos->FillTHnSparse("AmplitudeGammaLow", triggerpatchinfoamp);
 362         fHistos->FillTHnSparse("EnergyRoughGammaLow", triggerpatchinfoer);
 363       }
 364       if(triggerpatch->IsLevel0()){
 365         fHistos->FillTHnSparse("EnergyLevel0", triggerpatchinfo);
 366         fHistos->FillTHnSparse("AmplitudeLevel0", triggerpatchinfoamp);
 367         fHistos->FillTHnSparse("EnergyRoughLevel0", triggerpatchinfoer);
 368       }
 369     }
 370
 371     const AliVVertex *vtxTracks = fInputEvent->GetPrimaryVertex(),
 372         *vtxSPD = GetSPDVertex();
 373     if(!(vtxTracks && vtxSPD)) return false;
 374     if(!fVertexRange.IsInRange(vtxTracks->GetZ())) return false;
 375     if(vtxTracks->GetNContributors() < 1 || vtxSPD->GetNContributors() < 1) return false;
 376
 377     double triggers[5]; memset(triggers, 0, sizeof(double) * 5);
 378     double triggerbits[4]; memset(triggerbits, 0, sizeof(double) * 4);
 379     if(fInputHandler->IsEventSelected() & AliVEvent::kINT7){
 380       triggers[0] = 1.;
 381       triggerbits[0] = 1.;
 382     }
 383
 384     // check triggerbits
 385     if(fInputHandler->IsEventSelected() & AliVEvent::kEMC7){
 386       triggerbits[1] = 1.;
 387     }
 388     if(fInputHandler->IsEventSelected() & AliVEvent::kEMCEGA){
 389       triggerbits[2] = 1.;
 390     }
 391     if(fInputHandler->IsEventSelected() & AliVEvent::kEMCEJE){
 392       triggerbits[3] = 1.;
 393     }
 394     try{
 395       fHistos->FillTHnSparse("hEventsTriggerbit", triggerbits);
 396     } catch(HistoContainerContentException &e) {
 397       std::stringstream errormessage;
 398       errormessage << "Filling of histogram failed: " << e.what();
 399       AliError(errormessage.str().c_str());
 400     }
 401
 402     std::vector<std::string> triggerstrings;
 403     // EMCal-triggered event, distinguish types
 404     TString trgstr(fUseTriggersFromTriggerMaker ? BuildTriggerString() : fInputEvent->GetFiredTriggerClasses());
 405     AliDebug(1, Form("Triggerstring: %s\n", trgstr.Data()));
 406     if(trgstr.Contains("EJ1")){
 407       triggerstrings.push_back("EMCJHigh");
 408       triggers[1] = 1;
 409       if(trgstr.Contains("EG1"))
 410         triggerstrings.push_back("EMCHighBoth");
 411       else
 412         triggerstrings.push_back("EMCHighJetOnly");
 413     }
 414     if(trgstr.Contains("EJ2")){
 415       triggerstrings.push_back("EMCJLow");
 416       triggers[2] = 1;
 417       if(trgstr.Contains("EG2"))
 418         triggerstrings.push_back("EMCLowBoth");
 419       else
 420         triggerstrings.push_back("EMCLowJetOnly");
 421     }
 422     if(trgstr.Contains("EG1")){
 423       triggerstrings.push_back("EMCGHigh");
 424       triggers[3] = 1;
 425       if(!trgstr.Contains("EJ1"))
 426         triggerstrings.push_back("EMCHighGammaOnly");
 427     }
 428     if(trgstr.Contains("EG2")){
 429       triggerstrings.push_back("EMCGLow");
 430       triggers[4] = 1;
 431       if(!trgstr.Contains("EJ2"))
 432         triggerstrings.push_back("EMCLowGammaOnly");
 433     }
 434
 435     try{
 436       fHistos->FillTHnSparse("hEventTriggers", triggers);
 437     } catch (HistoContainerContentException &e){
 438       std::stringstream errormessage;
 439       errormessage << "Filling of histogram failed: " << e.what();
 440       AliError(errormessage.str().c_str());
 441     }
 442
 443     // apply event selection: Combine the Pileup cut from SPD with the other pA Vertex selection cuts.
 444     bool isPileupEvent = fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5.);
 445     isPileupEvent = isPileupEvent || (TMath::Abs(vtxTracks->GetZ() - vtxSPD->GetZ()) > 0.5);
 446     double covSPD[6]; vtxSPD->GetCovarianceMatrix(covSPD);
 447     isPileupEvent = isPileupEvent || (TString(vtxSPD->GetTitle()).Contains("vertexer:Z") && TMath::Sqrt(covSPD[5]) > 0.25);
 448
 449     // Fill event-based histogram
 450     const double &zv = vtxTracks->GetZ();
 451     if(triggers[0]) FillEventHist("MinBias", zv, isPileupEvent);
 452     if(!triggerstrings.size()) // Non-EMCal-triggered
 453       FillEventHist("NoEMCal", zv, isPileupEvent);
 454     else{
 455       // EMCal-triggered events
 456       for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
 457         FillEventHist(it->c_str(), zv, isPileupEvent);
 458     }
 459
 460     const AliEmcalJet *foundJet(NULL);
 461     char histname[1024];
 462     std::vector<double> coneradii;
 463     AliJetContainer *jetContMC(NULL), *jetContData(NULL);
 464
 465     // Fill MC truth
 466     AliVParticle *part(NULL);
 467     if(fMCEvent){
 468       if(fJetCollArray.GetEntries() &&
 469                   (jetContMC = dynamic_cast<AliJetContainer *>(fJetCollArray.FindObject((static_cast<TObjString *>(fJetContainersMC.At(0)))->String().Data())))){
 470         // In case we have a jet array we loop over all MC selected particles in the particle container of the jet array
 471         TIter particles(jetContMC->GetParticleContainer()->GetArray());
 472         while((part = dynamic_cast<AliVParticle *>(particles()))){
 473           if(part->Charge() == 0) continue;
 474           if(!fEtaRange.IsInRange(part->Eta())) continue;
 475           if(!fPtRange.IsInRange(part->Pt())) continue;
 476           FillMCParticleHist("hMCtrueParticles", part, zv, isPileupEvent);
 477
 478           /*
 479            * Jet part: Find track in jet container,
 480            * check according to number of particles in jet, and
 481            * check for different cone radii
 482            */
 483           foundJet = FoundTrackInJet(part, jetContMC);
 484           if(foundJet && foundJet->GetNumberOfConstituents() > 1){
 485             for(int irad = 0; irad < kNJetRadii; irad++){
 486               if(IsInRadius(part, foundJet, kJetRadii[irad])){
 487                 sprintf(histname, "hMCtrueParticlesRad%02d", int(kJetRadii[irad]*10));
 488                 FillMCParticleHist(histname,  part, zv, isPileupEvent);
 489               }
 490             }
 491           } else {
 492             // isolated track
 493             FillMCParticleHist("hMCtrueParticlesIsolated", part, zv, isPileupEvent);
 494           }
 495         }
 496       } else {
 497         // Use MC Event
 498         for(int ipart = 0; ipart < fMCEvent->GetNumberOfTracks(); ipart++){
 499           // Select only physical primary particles
 500           part = fMCEvent->GetTrack(ipart);
 501           if(part->Charge() == 0) continue;
 502           if(!fEtaRange.IsInRange(part->Eta())) continue;
 503           if(!fPtRange.IsInRange(part->Pt())) continue;
 504           if(!fMCEvent->IsPhysicalPrimary(ipart)) continue;
 505           FillMCParticleHist("hMCtrueParticles", part, zv, isPileupEvent);
 506         }
 507       }
 508     }
 509
 510     AliVTrack *track(NULL);
 511     AliPicoTrack *picoTrack(NULL);
 512     TObject *containerObject(NULL);
 513     // Loop over all tracks (No cuts applied)
 514     if(fSelectAllTracks){
 515       // loop over all tracks only if requested
 516       TIter allTrackIter(fTracks);
 517       while((containerObject = dynamic_cast<TObject *>(allTrackIter()))){
 518         if((picoTrack = dynamic_cast<AliPicoTrack *>(containerObject))){
 519           track = picoTrack->GetTrack();
 520         } else
 521           track = dynamic_cast<AliVTrack *>(containerObject);
 522         if(!IsTrueTrack(track)) continue;
 523         if(!fEtaRange.IsInRange(track->Eta())) continue;
 524         if(!fPtRange.IsInRange(track->Pt())) continue;
 525         if(triggers[0]) FillTrackHist("MinBias", track, zv, isPileupEvent, 0, triggers[0]);
 526         if(!triggerstrings.size()) // Non-EMCal-triggered
 527           FillTrackHist("NoEMCal", track, zv, isPileupEvent, 0, triggers[0]);
 528         else {
 529           // EMCal-triggered events
 530           for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
 531             FillTrackHist(it->c_str(), track, zv, isPileupEvent, 0, triggers[0]);
 532         }
 533       }
 534     }
 535
 536     // Now apply track selection cuts
 537     // allow for several track selections to be tested at the same time
 538     // each track selection gets a different cut ID starting from 1
 539     // cut ID 0 is reserved for the case of no cuts
 540     // In case we have a jet container attached we check whether the track is
 541     // found in a jet, and check for different cone radii around a jet
 542     if(fListTrackCuts && fListTrackCuts->GetEntries()){
 543       for(int icut = 0; icut < fListTrackCuts->GetEntries(); icut++){
 544         AliEMCalPtTaskVTrackSelection *trackSelection = static_cast<AliEMCalPtTaskVTrackSelection *>(fListTrackCuts->At(icut));
 545         TIter trackIter(trackSelection->GetAcceptedTracks(fTracks));
 546         while((track = dynamic_cast<AliVTrack *>(trackIter()))){
 547           if(fMCEvent && !IsTrueTrack(track)) continue;   // Reject fake tracks in case of MC
 548           if(!fEtaRange.IsInRange(track->Eta())) continue;
 549           if(!fPtRange.IsInRange(track->Pt())) continue;
 550           coneradii.clear();
 551           if(this->fJetCollArray.GetEntries() &&
 552               (jetContData = dynamic_cast<AliJetContainer *>(this->fJetCollArray.FindObject((static_cast<TObjString *>(this->fJetContainersData.At(0)))->String().Data())))){
 553             foundJet = FoundTrackInJet(track, jetContData);
 554             if(foundJet){
 555               for(int irad = 0; irad < kNJetRadii; irad++){
 556                 if(IsInRadius(track, foundJet, kJetRadii[irad])) coneradii.push_back(kJetRadii[irad]);
 557               }
 558             }
 559           }
 560           if(triggers[0]){
 561             FillTrackHist("MinBias", track, zv, isPileupEvent, icut + 1, triggers[0]);
 562             if(coneradii.size()){
 563               for(std::vector<double>::iterator radIter = coneradii.begin(); radIter != coneradii.end(); radIter++)
 564                 FillTrackHist("MinBias", track, zv, isPileupEvent, icut + 1, triggers[0], *radIter);
 565             }
 566           }
 567           if(!triggerstrings.size()){ // Non-EMCal-triggered
 568             FillTrackHist("NoEMCal", track, zv, isPileupEvent, icut + 1, triggers[0]);
 569             for(std::vector<double>::iterator radIter = coneradii.begin(); radIter != coneradii.end(); radIter++)
 570               FillTrackHist("NoEMCal", track, zv, isPileupEvent, icut + 1, triggers[0], *radIter);
 571           } else {
 572             // EMCal-triggered events
 573             for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it){
 574               FillTrackHist(it->c_str(), track, zv, isPileupEvent, icut + 1, triggers[0]);
 575               for(std::vector<double>::iterator radIter = coneradii.begin(); radIter != coneradii.end(); radIter++)
 576                 FillTrackHist(it->c_str(), track, zv, isPileupEvent, icut + 1, triggers[0], *radIter);
 577             }
 578           }
 579         }
 580       }
 581     }
 582
 583     // Next step we loop over the (uncalibrated) emcal clusters and fill histograms with the cluster energy
 584     const AliVCluster *clust(NULL);
 585     for(int icl = 0; icl < fInputEvent->GetNumberOfCaloClusters(); icl++){
 586       clust = fInputEvent->GetCaloCluster(icl);
 587       if(!clust->IsEMCAL()) continue;
 588       if(!fEnergyRange.IsInRange(clust->E())) continue;
 589       if(triggers[0]) FillClusterHist("hClusterUncalibHistMinBias", clust, zv, isPileupEvent, triggers[0]);
 590       if(!triggerstrings.size()){       // Non-EMCal-triggered
 591         FillClusterHist("hClusterUncalibHistNoEMCal", clust, zv, isPileupEvent, triggers[0]);
 592       } else{
 593         for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it){
 594           sprintf(histname, "hClusterUncalibHist%s", it->c_str());
 595           FillClusterHist(histname, clust, zv, isPileupEvent, triggers[0]);
 596         }
 597       }
 598     }
 599
 600     if(fCaloClusters){
 601       TIter clustIter(fCaloClusters);
 602       while((clust = dynamic_cast<const AliVCluster *>(clustIter()))){
 603         if(!clust->IsEMCAL()) continue;
 604         if(!fEnergyRange.IsInRange(clust->E())) continue;
 605         if(triggers[0]) FillClusterHist("hClusterCalibHistMinBias", clust, zv, isPileupEvent, triggers[0]);
 606         if(!triggerstrings.size())      // Non-EMCal-triggered
 607           FillClusterHist("hClusterCalibHistNoEMCal", clust, zv, isPileupEvent, triggers[0]);
 608         else{
 609           for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it){
 610             sprintf(histname, "hClusterCalibHist%s", it->c_str());
 611             FillClusterHist(histname, clust, zv, isPileupEvent, triggers[0]);
 612           }
 613         }
 614       }
 615     }
 616
 617     PostData(1, fOutput);
 618     return true;
 619   }
 620
 621   //______________________________________________________________________________
 622   void AliAnalysisTaskPtEMCalTrigger::CreateDefaultPtBinning(TArrayD &binning) const{
 623     /*
 624      * Creating the default pt binning.
 625      *
 626      * @param binning: Array where to store the results.
 627      */
 628     std::vector<double> mybinning;
 629     std::map<double,double> definitions;
 630     definitions.insert(std::pair<double,double>(2.5, 0.1));
 631     definitions.insert(std::pair<double,double>(7., 0.25));
 632     definitions.insert(std::pair<double,double>(15., 0.5));
 633     definitions.insert(std::pair<double,double>(25., 1.));
 634     definitions.insert(std::pair<double,double>(40., 2.5));
 635     definitions.insert(std::pair<double,double>(50., 5.));
 636     definitions.insert(std::pair<double,double>(100., 10.));
 637     double currentval = 0;
 638     for(std::map<double,double>::iterator id = definitions.begin(); id != definitions.end(); ++id){
 639       double limit = id->first, binwidth = id->second;
 640       while(currentval < limit){
 641         currentval += binwidth;
 642         mybinning.push_back(currentval);
 643       }
 644     }
 645     binning.Set(mybinning.size());
 646     int ib = 0;
 647     for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
 648       binning[ib++] = *it;
 649   }
 650
 651   //______________________________________________________________________________
 652   void AliAnalysisTaskPtEMCalTrigger::CreateDefaultZVertexBinning(TArrayD &binning) const {
 653     /*
 654      * Creating default z-Vertex binning.
 655      *
 656      * @param binning: Array where to store the results.
 657      */
 658     std::vector<double> mybinning;
 659     double currentval = -10;
 660     mybinning.push_back(currentval);
 661     while(currentval <= 10.){
 662       currentval += 5.;
 663       mybinning.push_back(currentval);
 664     }
 665     binning.Set(mybinning.size());
 666     int ib = 0;
 667     for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
 668       binning[ib++] = *it;
 669   }
 670
 671   //______________________________________________________________________________
 672   void AliAnalysisTaskPtEMCalTrigger::CreateDefaultEtaBinning(TArrayD& binning) const {
 673     /*
 674      * Creating default z-Vertex binning.
 675      *
 676      * @param binning: Array where to store the results.
 677      */
 678     std::vector<double> mybinning;
 679     double currentval = -0.8;
 680     mybinning.push_back(currentval);
 681     while(currentval <= 0.8){
 682       currentval += 0.1;
 683       mybinning.push_back(currentval);
 684     }
 685     binning.Set(mybinning.size());
 686     int ib = 0;
 687     for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
 688       binning[ib++] = *it;
 689   }
 690
 691   //______________________________________________________________________________
 692   void AliAnalysisTaskPtEMCalTrigger::DefineAxis(TAxis& axis, const char* name,
 693       const char* title, const TArrayD& binning, const char** labels) {
 694     /*
 695      * Define an axis with a given binning
 696      *
 697      * @param axis: Axis to be defined
 698      * @param name: Name of the axis
 699      * @param title: Title of the axis
 700      * @param binning: axis binning
 701      * @param labels (@optional): array of bin labels
 702      */
 703     axis.Set(binning.GetSize()-1, binning.GetArray());
 704     axis.SetName(name);
 705     axis.SetTitle(title);
 706     if(labels){
 707       for(int ib = 1; ib <= axis.GetNbins(); ++ib)
 708         axis.SetBinLabel(ib, labels[ib-1]);
 709     }
 710   }
 711
 712   //______________________________________________________________________________
 713   void AliAnalysisTaskPtEMCalTrigger::DefineAxis(TAxis& axis, const char* name,
 714       const char* title, int nbins, double min, double max,
 715       const char** labels) {
 716     /*
 717      * Define an axis with number of bins from min to max
 718      *
 719      * @param axis: Axis to be defined
 720      * @param name: Name of the axis
 721      * @param title: Title of the axis
 722      * @param nbins: Number of bins
 723      * @param min: lower limit of the axis
 724      * @param max: upper limit of the axis
 725      * @param labels (@optional): array of bin labels
 726      */
 727     axis.Set(nbins, min, max);
 728     axis.SetName(name);
 729     axis.SetTitle(title);
 730     if(labels){
 731       for(int ib = 1; ib <= axis.GetNbins(); ++ib)
 732         axis.SetBinLabel(ib, labels[ib-1]);
 733     }
 734   }
 735
 736
 737   //______________________________________________________________________________
 738   void AliAnalysisTaskPtEMCalTrigger::FillEventHist(const char* trigger,
 739       double vz, bool isPileup) {
 740     /*
 741      * Fill event-based histogram
 742      *
 743      * @param trigger: name of the trigger configuration to be processed
 744      * @param vz: z-position of the vertex
 745      * @param isPileup: signalises if the event is flagged as pileup event
 746      */
 747     char histname[1024];
 748     sprintf(histname, "hEventHist%s", trigger);
 749     try{
 750       fHistos->FillTH2(histname, 0., vz);
 751     } catch (HistoContainerContentException &e){
 752       std::stringstream errormessage;
 753       errormessage << "Filling of histogram failed: " << e.what();
 754       AliError(errormessage.str().c_str());
 755     }
 756     if(!isPileup){
 757       try{
 758         fHistos->FillTH2(histname, 1., vz);
 759       } catch(HistoContainerContentException &e){
 760         std::stringstream errormessage;
 761         errormessage << "Filling of histogram failed: " << e.what();
 762         AliError(errormessage.str().c_str());
 763       }
 764     }
 765   }
 766
 767   //______________________________________________________________________________
 768   void AliAnalysisTaskPtEMCalTrigger::FillTrackHist(const char* trigger,
 769       const AliVTrack* track, double vz, bool isPileup, int cut, bool isMinBias, double jetradius) {
 770     /*
 771      * Fill track-based histogram with corresponding information
 772      *
 773      * @param trigger: name of the trigger
 774      * @param track: ESD track selected
 775      * @param vz: z-position of the vertex
 776      * @param isPileup: flag event as pileup event
 777      * @param cut: id of the cut (0 = no cut)
 778      */
 779     double etasign = fSwapEta ? -1. : 1.;
 780     double data[7] = {TMath::Abs(track->Pt()), etasign * track->Eta(), track->Phi(), vz, 0, static_cast<double>(cut), isMinBias ? 1. : 0.};
 781     double dataMC[7] = {0., 0., 0., vz, 0, static_cast<double>(cut), isMinBias ? 1. : 0.};
 782     AliVParticle *assocMC(NULL);
 783     if(fMCEvent && (assocMC = fMCEvent->GetTrack(TMath::Abs(track->GetLabel())))){
 784       // Select track onl
 785       dataMC[0] = TMath::Abs(assocMC->Pt());
 786       dataMC[1] = etasign * assocMC->Eta();
 787       dataMC[2] = assocMC->Phi();
 788     }
 789     char histname[1024], histnameAcc[1024], histnameMC[1024], histnameMCAcc[1024];
 790     sprintf(histname, "hTrackHist%s", trigger);
 791     sprintf(histnameAcc, "hTrackInAcceptanceHist%s", trigger);
 792     sprintf(histnameMC, "hMCTrackHist%s", trigger);
 793     sprintf(histnameMCAcc, "hMCTrackInAcceptanceHist%s", trigger);
 794     if(jetradius > 0.){
 795       char *hnames[] = {histname, histnameAcc, histnameMC, histnameMCAcc};
 796       for(unsigned int iname = 0; iname < sizeof(hnames)/sizeof(char *); iname++){
 797         char *myhname = hnames[iname];
 798         sprintf(myhname, "%sRad%02d", myhname, int(jetradius * 10.));
 799       }
 800     }
 801     Bool_t isEMCAL = kFALSE;
 802     if(track->IsEMCAL()){
 803       // Check if the cluster is matched to only one track
 804       AliVCluster *emcclust(NULL);
 805       AliDebug(2, Form("cluster id: %d\n", track->GetEMCALcluster()));
 806       if(fCaloClusters) {
 807         AliDebug(2, "Using calibrated clusters");
 808         emcclust = dynamic_cast<AliVCluster *>(fCaloClusters->At(track->GetEMCALcluster()));
 809       } else {
 810         AliDebug(2, "Using uncalibrated clusters");
 811         emcclust = fInputEvent->GetCaloCluster(track->GetEMCALcluster());
 812       }
 813       if(!emcclust) AliError("Null pointer to EMCal cluster");
 814       if(emcclust && emcclust->GetNTracksMatched() <= 1){
 815         isEMCAL = kTRUE;
 816       }
 817     }
 818     try{
 819       fHistos->FillTHnSparse(histname, data);
 820       if(fMCEvent) fHistos->FillTHnSparse(histnameMC, dataMC);
 821       if(isEMCAL){
 822         fHistos->FillTHnSparse(histnameAcc, data);
 823         if(fMCEvent) fHistos->FillTHnSparse(histnameMCAcc, dataMC);
 824       }
 825     } catch (HistoContainerContentException &e){
 826       std::stringstream errormessage;
 827       errormessage << "Filling of histogram failed: " << e.what();
 828       AliError(errormessage.str().c_str());
 829     }
 830     if(!isPileup){
 831       data[4] = 1;
 832       dataMC[4] = 1;
 833       try{
 834         fHistos->FillTHnSparse(histname, data);
 835         if(fMCEvent) fHistos->FillTHnSparse(histnameMC, dataMC);
 836         if(isEMCAL){
 837           fHistos->FillTHnSparse(histnameAcc, data);
 838           if(fMCEvent) fHistos->FillTHnSparse(histnameMCAcc, dataMC);
 839         }
 840       } catch (HistoContainerContentException &e){
 841         std::stringstream errormessage;
 842         errormessage << "Filling of histogram failed: " << e.what();
 843         AliError(errormessage.str().c_str());
 844       }
 845     }
 846   }
 847
 848   //______________________________________________________________________________
 849   void AliAnalysisTaskPtEMCalTrigger::FillClusterHist(const char* histname,
 850       const AliVCluster* clust, double vz, bool isPileup, bool isMinBias) {
 851     /*
 852      * Fill cluster-based histogram with corresponding information
 853      *
 854      * @param trigger: name of the trigger
 855      * @param cluster: the EMCal cluster information
 856      * @param vz: z-position of the vertex
 857      * @param isPileup: flag event as pileup event
 858      */
 859     double data[4] =  {clust->E(), vz, 0, isMinBias ? 1. : 0.};
 860     try{
 861       fHistos->FillTHnSparse(histname, data);
 862     } catch (HistoContainerContentException &e){
 863       std::stringstream errormessage;
 864       errormessage << "Filling of histogram failed: " << e.what();
 865       AliError(errormessage.str().c_str());
 866     }
 867     if(!isPileup){
 868       data[2] = 1.;
 869       try{
 870         fHistos->FillTHnSparse(histname, data);
 871       } catch (HistoContainerContentException &e){
 872         std::stringstream errormessage;
 873         errormessage << "Filling of histogram failed: " << e.what();
 874         AliError(errormessage.str().c_str());
 875       }
 876     }
 877   }
 878
 879   //______________________________________________________________________________
 880   void AliAnalysisTaskPtEMCalTrigger::FillMCParticleHist(const char *histname, const AliVParticle * const track, double vz, bool isPileup){
 881     /*
 882      * Fill histogram for MC-true particles with the information pt, eta and phi
 883      *
 884      * @param track: the Monte-Carlo track
 885      */
 886     double data[5] = {TMath::Abs(track->Pt()), track->Eta(), track->Phi(), vz, 0.};
 887     fHistos->FillTHnSparse(histname, data);
 888     if(!isPileup){
 889       data[4] = 1.;
 890       fHistos->FillTHnSparse(histname, data);
 891     }
 892   }
 893
 894   //______________________________________________________________________________
 895   bool AliAnalysisTaskPtEMCalTrigger::IsTrueTrack(const AliVTrack *const track) const{
 896     /*
 897      * Check if the track has an associated MC particle, and that the particle is a physical primary
 898      * In case of data we do not do the selection at that step (always return true)
 899      *
 900      * @param track: Track to check
 901      * @result: true primary track (true or false)
 902      */
 903     if(!fMCEvent) return true;
 904     AliVParticle *mcassociate = fMCEvent->GetTrack(TMath::Abs(track->GetLabel()));
 905     if(!mcassociate) return false;
 906     return fMCEvent->IsPhysicalPrimary(TMath::Abs(track->GetLabel()));
 907   }
 908
 909   //______________________________________________________________________________
 910   void AliAnalysisTaskPtEMCalTrigger::AddESDTrackCuts(AliESDtrackCuts* trackCuts) {
 911     /*
 912      * Add new track cuts to the task
 913      *
 914      * @param trackCuts: Object of type AliESDtrackCuts
 915      */
 916     fListTrackCuts->AddLast(new AliEMCalPtTaskTrackSelectionESD(trackCuts));
 917   }
 918
 919   //______________________________________________________________________________
 920   void AliAnalysisTaskPtEMCalTrigger::AddCutsForAOD(AliESDtrackCuts* trackCuts, UInt_t filterbits) {
 921     /*
 922      * Add new track cuts to the task
 923      *
 924      * @param trackCuts: Object of type AliESDtrackCuts
 925      */
 926     fListTrackCuts->AddLast(new AliEMCalPtTaskTrackSelectionAOD(trackCuts, filterbits));
 927   }
 928
 929
 930   //______________________________________________________________________________
 931   TString AliAnalysisTaskPtEMCalTrigger::BuildTriggerString() {
 932     /*
 933      * Build trigger string from the trigger maker
 934      *
 935      * @return: blank-separated string of fired trigger classes
 936      */
 937     AliDebug(1, "trigger checking");
 938     TString result = "";
 939     if(HasTriggerType(kJ1)) result += "EJ1 ";
 940     if(HasTriggerType(kJ2)) result += "EJ2 ";
 941     if(HasTriggerType(kG1)) result += "EG1 ";
 942     if(HasTriggerType(kG2)) result += "EG2 ";
 943     return result;
 944   }
 945
 946   //______________________________________________________________________________
 947   const AliVVertex* AliAnalysisTaskPtEMCalTrigger::GetSPDVertex() const {
 948     /*
 949      * Accessor for the SPD vertex, creating transparency for ESDs and AODs
 950      *
 951      * @return: the spd vertex
 952      */
 953     AliESDEvent *esd = dynamic_cast<AliESDEvent *>(fInputEvent);
 954     if(esd){
 955       return esd->GetPrimaryVertexSPD();
 956     } else {
 957       AliAODEvent *aod = dynamic_cast<AliAODEvent *>(fInputEvent);
 958       if(aod){
 959         return aod->GetPrimaryVertexSPD();
 960       }
 961     }
 962     return NULL;
 963   }
 964
 965   //______________________________________________________________________________
 966   const AliEmcalJet * AliAnalysisTaskPtEMCalTrigger::FoundTrackInJet(
 967       const AliVParticle * const track, AliJetContainer *const jets) const
 968   {
 969     /*
 970      * Correlate track to reconstructed jet
 971      *
 972      * @param track: particle to be checked
 973      * @param jets: container of recontructed jets
 974      * @return: The matched jet (NULL if not found)
 975      */
 976     const AliEmcalJet *result = NULL;
 977     const AliEmcalJet *tmp = jets->GetNextAcceptJet(0);
 978     while(tmp){
 979       if(TrackInJet(track, tmp, jets->GetParticleContainer())){
 980         result = tmp;
 981         break;
 982       }
 983       tmp = jets->GetNextAcceptJet();
 984     }
 985     return result;
 986   }
 987
 988   //______________________________________________________________________________
 989   bool AliAnalysisTaskPtEMCalTrigger::IsInRadius(const AliVParticle *const track, const AliEmcalJet *reconstructedJet, Double_t radius) const {
 990     /*
 991      * Check if track is in radius around a given jet
 992      *
 993      * @param track: Track to check
 994      * @param reconstructed jet: jet to probe
 995      * @param radius: cone radius
 996      * @return: result of the test (true if track is inside the cone radius, false otherwise)
 997      */
 998     return reconstructedJet->DeltaR(track) < radius;
 999   }
1000
1001   //______________________________________________________________________________
1002   bool AliAnalysisTaskPtEMCalTrigger::IsInRadius(const AliVCluster *const clust, const AliEmcalJet *reconstructedJet, Double_t radius) const {
1003     /*
1004      * Check if track is in radius around a given jet
1005      *
1006      * @param track: Track to check
1007      * @param reconstructed jet: jet to probe
1008      * @param radius: cone radius
1009      * @return: result of the test (true if track is inside the cone radius, false otherwise)
1010      */
1011     double vertexPos[3];
1012     fInputEvent->GetPrimaryVertex()->GetXYZ(vertexPos);
1013     TLorentzVector clustVect;
1014     clust->GetMomentum(clustVect, vertexPos);
1015
1016     Double_t dPhi = reconstructedJet->Phi() - clustVect.Phi();
1017     Double_t dEta = reconstructedJet->Eta() - clustVect.Eta();
1018     dPhi = TVector2::Phi_mpi_pi ( dPhi );
1019
1020     double deltaR = TMath::Sqrt ( dPhi * dPhi + dEta * dEta );
1021     return deltaR < radius;
1022   }
1023
1024   //______________________________________________________________________________
1025   bool AliAnalysisTaskPtEMCalTrigger::TrackInJet(const AliVParticle* const track,
1026       const AliEmcalJet* reconstructedJet, const AliParticleContainer * const particles) const {
1027     /*
1028      * Check whether track is among the jet constituents
1029      *
1030      * @param track: track to check
1031      * @param reconstructedJet: reconstructed jet to check
1032      * @param tracks: container with tracks used for jetfinding
1033      * @return: true if found, false otherwise
1034      */
1035     bool found = false;
1036     const AliPicoTrack *picotmp(NULL);
1037     const AliVParticle *tmp(NULL);
1038     for(int ipart = 0; ipart < reconstructedJet->GetNumberOfTracks(); ipart++){
1039       tmp = dynamic_cast<const AliVParticle *>(reconstructedJet->TrackAt(ipart, particles->GetArray()));
1040       if((picotmp = dynamic_cast<const AliPicoTrack *>(tmp)))   // handle pico tracks
1041         tmp = picotmp->GetTrack();
1042       if(!tmp->Compare(track)){
1043         found = true;
1044         break;
1045       }
1046     }
1047     return found;
1048   }
1049
1050   //______________________________________________________________________________
1051   const AliEmcalJet* AliAnalysisTaskPtEMCalTrigger::FoundClusterInJet(const AliVCluster* const clust, AliJetContainer* const jets) const {
1052     /*
1053      * Check whether a cluster is in a radius around a given jet
1054      *
1055      * @param clust: The cluster to check
1056      * @param reconstructedJet: reconstructed jet to check
1057      * @return: the jet containing the cluster (null otherwise)
1058      */
1059     const AliEmcalJet *result = NULL;
1060     const AliEmcalJet *tmp = jets->GetNextAcceptJet(0);
1061     while(tmp){
1062       if(ClusterInJet(clust, tmp, jets->GetClusterContainer())){
1063         result = tmp;
1064         break;
1065       }
1066       tmp =jets->GetNextAcceptJet();
1067     }
1068     return result;
1069   }
1070
1071   //______________________________________________________________________________
1072   bool AliAnalysisTaskPtEMCalTrigger::ClusterInJet(const AliVCluster* const clust,
1073       const AliEmcalJet* reconstructedJet, const AliClusterContainer* const clusters) const {
1074     /*
1075      * Check whether cluster is among the jet constituents
1076      *
1077      * @param track: track to check
1078      * @param reconstructedJet: reconstructed jet to check
1079      * @param clusters: the cluster container
1080      * @return: true if found, false otherwise
1081      */
1082     bool found = false;
1083     const AliVCluster *tmp(NULL);
1084     for(int ipart = 0; ipart < reconstructedJet->GetNumberOfTracks(); ipart++){
1085       tmp = dynamic_cast<const AliVCluster *>(reconstructedJet->ClusterAt(ipart, clusters->GetArray()));
1086       if(!tmp->Compare(clust)){
1087         found = true;
1088         break;
1089       }
1090     }
1091     return found;
1092   }
1093
1094   //______________________________________________________________________________
1095   void EMCalTriggerPtAnalysis::AliAnalysisTaskPtEMCalTrigger::AddJetContainerName(const Char_t* contname, Bool_t isMC) {
1096     /*
1097      * Add new Jet input container to the analysis task
1098      *
1099      * @param contname: Name of the container
1100      * @param isMC: Defines whether the container is for MC truth or not
1101      */
1102     TList &mycontainer = isMC ? fJetContainersMC : fJetContainersData;
1103     mycontainer.Add(new TObjString(contname));
1104   }
1105 }