PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskPtEMCalTrigger.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2007, 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 <TString.h>
  39
  40 #include "AliAODEvent.h"
  41 #include "AliESDEvent.h"
  42 #include "AliInputEventHandler.h"
  43 #include "AliMCEvent.h"
  44 #include "AliParticleContainer.h"
  45 #include "AliVCluster.h"
  46 #include "AliVParticle.h"
  47 #include "AliVTrack.h"
  48 #include "AliVVertex.h"
  49
  50 #include "AliEmcalTriggerPatchInfo.h"
  51 #include "AliEMCalHistoContainer.h"
  52 #include "AliEMCalPtTaskVTrackSelection.h"
  53 #include "AliEMCalPtTaskTrackSelectionAOD.h"
  54 #include "AliEMCalPtTaskTrackSelectionESD.h"
  55 #include "AliAnalysisTaskPtEMCalTrigger.h"
  56
  57 ClassImp(EMCalTriggerPtAnalysis::AliAnalysisTaskPtEMCalTrigger)
  58
  59 namespace EMCalTriggerPtAnalysis {
  60
  61   //______________________________________________________________________________
  62   AliAnalysisTaskPtEMCalTrigger::AliAnalysisTaskPtEMCalTrigger():
  63                                     AliAnalysisTaskEmcal(),
  64                                     fHistos(NULL),
  65                                     fListTrackCuts(NULL),
  66                                     fEtaRange(),
  67                                     fPtRange(),
  68                                     fSwapEta(kFALSE),
  69                                     fUseTriggersFromTriggerMaker(kFALSE)
  70   {
  71     /*
  72      * Dummy constructor, initialising the values with default (NULL) values
  73      */
  74   }
  75
  76   //______________________________________________________________________________
  77   AliAnalysisTaskPtEMCalTrigger::AliAnalysisTaskPtEMCalTrigger(const char *name):
  78                                     AliAnalysisTaskEmcal(name, kTRUE),
  79                                     fHistos(NULL),
  80                                     fListTrackCuts(NULL),
  81                                     fEtaRange(),
  82                                     fPtRange(),
  83                                     fSwapEta(kFALSE),
  84                                     fUseTriggersFromTriggerMaker(kFALSE)
  85   {
  86     /*
  87      * Main constructor, setting default values for eta and zvertex cut
  88      */
  89
  90     fListTrackCuts = new TList;
  91     fListTrackCuts->SetOwner(false);
  92
  93     // Set default cuts
  94     fEtaRange.SetLimits(-0.8, 0.8);
  95     fPtRange.SetLimits(0.15, 100.);
  96     SetMakeGeneralHistograms(kTRUE);
  97   }
  98
  99   //______________________________________________________________________________
 100   AliAnalysisTaskPtEMCalTrigger::~AliAnalysisTaskPtEMCalTrigger(){
 101     /*
 102      * Destructor, deleting output
 103      */
 104     //if(fTrackSelection) delete fTrackSelection;
 105     if(fHistos) delete fHistos;
 106     if(fListTrackCuts) delete fListTrackCuts;
 107   }
 108
 109   //______________________________________________________________________________
 110   void AliAnalysisTaskPtEMCalTrigger::UserCreateOutputObjects(){
 111     /*
 112      * Create the list of output objects and define the histograms.
 113      * Also adding the track cuts to the list of histograms.
 114      */
 115     AliAnalysisTaskEmcal::UserCreateOutputObjects();
 116     TString trackContainerName = "ESDFilterTracks", clusterContainerName = "EmcCaloClusters";
 117     if(!fIsEsd){
 118       trackContainerName = "AODFilterTracks";
 119       clusterContainerName = "EmcCaloClusters";
 120     }
 121     AliParticleContainer *trackContainer = this->AddParticleContainer(trackContainerName.Data());
 122     trackContainer->SetClassName("AliVTrack");
 123     this->AddClusterContainer(clusterContainerName.Data());
 124     this->SetCaloTriggerPatchInfoName("EmcalTriggers");
 125     fHistos = new AliEMCalHistoContainer("PtEMCalTriggerHistograms");
 126     fHistos->ReleaseOwner();
 127
 128     std::map<std::string, std::string> triggerCombinations;
 129     const char *triggernames[12] = {"MinBias", "EMCJHigh", "EMCJLow", "EMCGHigh", "EMCGLow", "NoEMCal", "EMCHighBoth", "EMCHighGammaOnly", "EMCHighJetOnly", "EMCLowBoth", "EMCLowGammaOnly", "EMCLowJetOnly"},
 130         *bitnames[4] = {"CINT7", "EMC7", "kEMCEGA", "kEMCEJE"};
 131     // Define axes for the trigger correlation histogram
 132     const TAxis *triggeraxis[5]; memset(triggeraxis, 0, sizeof(const TAxis *) * 5);
 133     const TAxis *bitaxes[4]; memset(bitaxes, 0, sizeof(TAxis *) * 4);
 134     const char *binlabels[2] = {"OFF", "ON"};
 135     TAxis mytrgaxis[5], mybitaxis[4];
 136     for(int itrg = 0; itrg < 5; ++itrg){
 137       DefineAxis(mytrgaxis[itrg], triggernames[itrg], triggernames[itrg], 2, -0.5, 1.5, binlabels);
 138       triggeraxis[itrg] = mytrgaxis+itrg;
 139       if(itrg < 4){
 140         DefineAxis(mybitaxis[itrg], bitnames[itrg], bitnames[itrg], 2, -0.5, 1.5, binlabels);
 141         bitaxes[itrg] = mybitaxis+itrg;
 142       }
 143     }
 144     // Define names and titles for different triggers in the histogram container
 145     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[0], "min. bias events"));
 146     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[1], "jet-triggered events (high threshold)"));
 147     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[2], "jet-triggered events (low threshold)"));
 148     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[3], "gamma-triggered events (high threshold)"));
 149     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[4], "gamma-triggered events (low threshold)"));
 150     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[5], "non-EMCal-triggered events"));
 151     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[6], "jet and gamma triggered events (high threshold)"));
 152     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[7], "exclusively gamma-triggered events (high threshold)"));
 153     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[8], "exclusively jet-triggered events (high threshold)"));
 154     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[9], "jet and gamma triggered events (low threshold)"));
 155     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[10], "exclusively gamma-triggered events (low threshold)"));
 156     triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[11], "exclusively-triggered events (low threshold)"));
 157     // Define axes for the pt histogram
 158     // Dimensions:
 159     // 1. pt
 160     // 2. eta
 161     // 3. phi
 162     // 4. vertex
 163     // 5. pileup (0 = all events, 1 = after pileup rejection)
 164     // 6. track cuts (0 = no cuts; 1 = after std cuts)
 165     TArrayD ptbinning, zvertexBinning, etabinning, pileupaxis(3);
 166     pileupaxis[0] = -0.5; pileupaxis[1] = 0.5; pileupaxis[2] = 1.5;
 167     CreateDefaultPtBinning(ptbinning);
 168     CreateDefaultZVertexBinning(zvertexBinning);
 169     CreateDefaultEtaBinning(etabinning);
 170     TAxis htrackaxes[7];
 171     DefineAxis(htrackaxes[0], "pt", "p_{t} (GeV/c)", ptbinning);
 172     DefineAxis(htrackaxes[1], "eta", "#eta", etabinning);
 173     DefineAxis(htrackaxes[2], "phi", "#phi", 20, 0, 2 * TMath::Pi());
 174     DefineAxis(htrackaxes[3], "zvertex", "z_{V} (cm)", zvertexBinning);
 175     DefineAxis(htrackaxes[4], "pileup", "Pileup rejection", 2, -0.5, 1.5);
 176     DefineAxis(htrackaxes[5], "trackcuts", "Track Cuts", (fListTrackCuts ? fListTrackCuts->GetEntries() : 0) + 1, -0.5, (fListTrackCuts ? fListTrackCuts->GetEntries() : 0) + 0.5);
 177     DefineAxis(htrackaxes[6], "mbtrigger", "Has MB trigger", 2, -0.5, 1.5);
 178     const TAxis *trackaxes[7];
 179     for(int iaxis = 0; iaxis < 7; ++iaxis) trackaxes[iaxis] = htrackaxes + iaxis;
 180     TAxis hclusteraxes[4];
 181     DefineAxis(hclusteraxes[0], "energy", "E (GeV)", ptbinning);
 182     DefineAxis(hclusteraxes[1], "zvertex", "z_{V} (cm)", zvertexBinning);
 183     DefineAxis(hclusteraxes[2], "pileup", "Pileup rejection", 2,1 -0.5, 1.5);
 184     DefineAxis(hclusteraxes[3], "mbtrigger", "Has MB trigger", 2, -0.5, 1.5);
 185     const TAxis *clusteraxes[4];
 186     for(int iaxis = 0; iaxis < 4; ++iaxis) clusteraxes[iaxis] = hclusteraxes + iaxis;
 187     TAxis hpatchenergyaxes[4];
 188     DefineAxis(hpatchenergyaxes[0], "energy", "Patch energy (GeV)", 100, 0., 100.);
 189     DefineAxis(hpatchenergyaxes[1], "eta", "#eta", etabinning);
 190     DefineAxis(hpatchenergyaxes[2], "phi", "#phi",  20, 0, 2 * TMath::Pi());
 191     DefineAxis(hpatchenergyaxes[3], "isMain", "Main trigger", 2, -0.5, 1.5);
 192     const TAxis *patchenergyaxes[4];
 193     for(int iaxis = 0; iaxis < 4; ++iaxis) patchenergyaxes[iaxis] = hpatchenergyaxes + iaxis;
 194     TAxis hpatchampaxes[4];
 195     DefineAxis(hpatchampaxes[0], "amplitude", "Patch energy (GeV)", 2000, 0., 2000.);
 196     DefineAxis(hpatchampaxes[1], "eta", "#eta", etabinning);
 197     DefineAxis(hpatchampaxes[2], "phi", "#phi",  20, 0, 2 * TMath::Pi());
 198     DefineAxis(hpatchampaxes[3], "isMain", "Main trigger", 2, -0.5, 1.5);
 199     const TAxis *patchampaxes[4];
 200     for(int iaxis = 0; iaxis < 4; ++iaxis) patchampaxes[iaxis] = hpatchampaxes + iaxis;
 201     try{
 202       std::string patchnames[] = {"Level0", "JetHigh", "JetLow", "GammaHigh", "GammaLow"};
 203       for(std::string * triggerpatch = patchnames; triggerpatch < patchnames + sizeof(patchnames)/sizeof(std::string); ++triggerpatch){
 204         fHistos->CreateTHnSparse(Form("Energy%s", triggerpatch->c_str()), Form("Patch energy for %s trigger patches", triggerpatch->c_str()), 4, patchenergyaxes);
 205         fHistos->CreateTHnSparse(Form("EnergyRough%s", triggerpatch->c_str()), Form("Rough patch energy for %s trigger patches", triggerpatch->c_str()), 4, patchenergyaxes);
 206         fHistos->CreateTHnSparse(Form("Amplitude%s", triggerpatch->c_str()), Form("Patch amplitude for %s trigger patches", triggerpatch->c_str()), 4, patchampaxes);
 207       }
 208
 209       // Create histogram for MC-truth
 210       fHistos->CreateTHnSparse("hMCtrueParticles", "Particle-based histogram for MC-true particles", 5, trackaxes);
 211       for(std::map<std::string,std::string>::iterator it = triggerCombinations.begin(); it != triggerCombinations.end(); ++it){
 212         const std::string name = it->first, &title = it->second;
 213         // Create event-based histogram
 214         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);
 215         // Create track-based histogram
 216         fHistos->CreateTHnSparse(Form("hTrackHist%s", name.c_str()), Form("Track-based data for %s events", title.c_str()), 7, trackaxes);
 217         fHistos->CreateTHnSparse(Form("hTrackInAcceptanceHist%s", name.c_str()), Form("Track-based data for %s events", title.c_str()), 7, trackaxes);
 218         fHistos->CreateTHnSparse(Form("hMCTrackHist%s", name.c_str()), Form("Track-based data for %s events with MC kinematics", title.c_str()), 7, trackaxes);
 219         fHistos->CreateTHnSparse(Form("hMCTrackInAcceptanceHist%s", name.c_str()), Form("Track-based data for %s events with MC kinematics", title.c_str()), 7, trackaxes);
 220         // Create cluster-based histogram (Uncalibrated and calibrated clusters)
 221         fHistos->CreateTHnSparse(Form("hClusterCalibHist%s", name.c_str()), Form("Calib. cluster-based histogram for %s events", title.c_str()), 4, clusteraxes);
 222         fHistos->CreateTHnSparse(Form("hClusterUncalibHist%s", name.c_str()), Form("Uncalib. cluster-based histogram for %s events", title.c_str()), 4, clusteraxes);
 223       }
 224       fHistos->CreateTHnSparse("hEventTriggers", "Trigger type per event", 5, triggeraxis);
 225       fHistos->CreateTHnSparse("hEventsTriggerbit", "Trigger bits for the different events", 4, bitaxes);
 226     } catch (HistoContainerContentException &e){
 227       std::stringstream errormessage;
 228       errormessage << "Creation of histogram failed: " << e.what();
 229       AliError(errormessage.str().c_str());
 230     }
 231     fOutput->Add(fHistos->GetListOfHistograms());
 232     if(fListTrackCuts && fListTrackCuts->GetEntries()){
 233       TIter cutIter(fListTrackCuts);
 234       AliEMCalPtTaskVTrackSelection *cutObject(NULL);
 235       while((cutObject = dynamic_cast<AliEMCalPtTaskVTrackSelection *>(cutIter()))){
 236         AliESDtrackCuts *cuts = dynamic_cast<AliESDtrackCuts *>(cutObject->GetTrackCuts());
 237         if(cuts){
 238           cuts->DefineHistograms();
 239           fOutput->Add(cuts);
 240         }
 241       }
 242     }
 243     PostData(1, fOutput);
 244   }
 245
 246   //______________________________________________________________________________
 247   Bool_t AliAnalysisTaskPtEMCalTrigger::Run(){
 248     /*
 249      * Runs the event loop
 250      *
 251      * @param option: Additional options
 252      */
 253     // Common checks: Have SPD vertex and primary vertex from tracks, and both need to have at least one contributor
 254     AliDebug(1,Form("Number of calibrated clusters: %d", fCaloClusters->GetEntries()));
 255     AliDebug(1,Form("Number of matched tracks: %d", fTracks->GetEntries()));
 256     if(fMCEvent){
 257       // Build always trigger strig from the trigger maker in case of MC
 258       fUseTriggersFromTriggerMaker = kTRUE;
 259     }
 260
 261     // Loop over trigger patches, fill patch energy
 262     AliEmcalTriggerPatchInfo *triggerpatch(NULL);
 263     TIter patchIter(this->fTriggerPatchInfo);
 264     while((triggerpatch = dynamic_cast<AliEmcalTriggerPatchInfo *>(patchIter()))){
 265       double triggerpatchinfo[4] = {triggerpatch->GetPatchE(), triggerpatch->GetEtaGeo(), triggerpatch->GetPhiGeo(), triggerpatch->IsMainTrigger() ? 1 : 0};
 266       double triggerpatchinfoamp[4] = {triggerpatch->GetADCAmp(), triggerpatch->GetEtaGeo(), triggerpatch->GetPhiGeo(), triggerpatch->IsMainTrigger() ? 1 : 0};
 267       double triggerpatchinfoer[4] = {triggerpatch->GetADCAmpGeVRough(), triggerpatch->GetEtaGeo(), triggerpatch->GetPhiGeo(), triggerpatch->IsMainTrigger() ? 1 : 0};
 268       if(triggerpatch->IsJetHigh()){
 269         fHistos->FillTHnSparse("EnergyJetHigh", triggerpatchinfo);
 270         fHistos->FillTHnSparse("AmplitudeJetHigh", triggerpatchinfoamp);
 271         fHistos->FillTHnSparse("EnergyRoughJetHigh", triggerpatchinfoer);
 272       }
 273       if(triggerpatch->IsJetLow()){
 274         fHistos->FillTHnSparse("EnergyJetLow", triggerpatchinfo);
 275         fHistos->FillTHnSparse("AmplitudeJetLow", triggerpatchinfoamp);
 276         fHistos->FillTHnSparse("EnergyRoughJetLow", triggerpatchinfoer);
 277       }
 278       if(triggerpatch->IsGammaHigh()){
 279         fHistos->FillTHnSparse("EnergyGammaHigh", triggerpatchinfo);
 280         fHistos->FillTHnSparse("AmplitudeGammaHigh", triggerpatchinfoamp);
 281         fHistos->FillTHnSparse("EnergyRoughGammaHigh", triggerpatchinfoer);
 282       }
 283       if(triggerpatch->IsGammaLow()){
 284         fHistos->FillTHnSparse("EnergyGammaLow", triggerpatchinfo);
 285         fHistos->FillTHnSparse("AmplitudeGammaLow", triggerpatchinfoamp);
 286         fHistos->FillTHnSparse("EnergyRoughGammaLow", triggerpatchinfoer);
 287       }
 288       if(triggerpatch->IsLevel0()){
 289         fHistos->FillTHnSparse("EnergyLevel0", triggerpatchinfo);
 290         fHistos->FillTHnSparse("AmplitudeLevel0", triggerpatchinfoamp);
 291         fHistos->FillTHnSparse("EnergyRoughLevel0", triggerpatchinfoer);
 292       }
 293     }
 294
 295     const AliVVertex *vtxTracks = fInputEvent->GetPrimaryVertex(),
 296         *vtxSPD = GetSPDVertex();
 297     if(!(vtxTracks && vtxSPD)) return false;
 298     if(vtxTracks->GetNContributors() < 1 || vtxSPD->GetNContributors() < 1) return false;
 299
 300     double triggers[5]; memset(triggers, 0, sizeof(double) * 5);
 301     double triggerbits[4]; memset(triggerbits, 0, sizeof(double) * 4);
 302     if(fInputHandler->IsEventSelected() & AliVEvent::kINT7){
 303       triggers[0] = 1.;
 304       triggerbits[0] = 1.;
 305     }
 306
 307     // check triggerbits
 308     if(fInputHandler->IsEventSelected() & AliVEvent::kEMC7){
 309       triggerbits[1] = 1.;
 310     }
 311     if(fInputHandler->IsEventSelected() & AliVEvent::kEMCEGA){
 312       triggerbits[2] = 1.;
 313     }
 314     if(fInputHandler->IsEventSelected() & AliVEvent::kEMCEJE){
 315       triggerbits[3] = 1.;
 316     }
 317     try{
 318       fHistos->FillTHnSparse("hEventsTriggerbit", triggerbits);
 319     } catch(HistoContainerContentException &e) {
 320       std::stringstream errormessage;
 321       errormessage << "Filling of histogram failed: " << e.what();
 322       AliError(errormessage.str().c_str());
 323     }
 324
 325     std::vector<std::string> triggerstrings;
 326     // EMCal-triggered event, distinguish types
 327     TString trgstr(fUseTriggersFromTriggerMaker ? BuildTriggerString() : fInputEvent->GetFiredTriggerClasses());
 328     AliDebug(1, Form("Triggerstring: %s\n", trgstr.Data()));
 329     if(trgstr.Contains("EJ1")){
 330       triggerstrings.push_back("EMCJHigh");
 331       triggers[1] = 1;
 332       if(trgstr.Contains("EG1"))
 333         triggerstrings.push_back("EMCHighBoth");
 334       else
 335         triggerstrings.push_back("EMCHighJetOnly");
 336     }
 337     if(trgstr.Contains("EJ2")){
 338       triggerstrings.push_back("EMCJLow");
 339       triggers[2] = 1;
 340       if(trgstr.Contains("EG2"))
 341         triggerstrings.push_back("EMCLowBoth");
 342       else
 343         triggerstrings.push_back("EMCLowJetOnly");
 344     }
 345     if(trgstr.Contains("EG1")){
 346       triggerstrings.push_back("EMCGHigh");
 347       triggers[3] = 1;
 348       if(!trgstr.Contains("EJ1"))
 349         triggerstrings.push_back("EMCHighGammaOnly");
 350     }
 351     if(trgstr.Contains("EG2")){
 352       triggerstrings.push_back("EMCGLow");
 353       triggers[4] = 1;
 354       if(!trgstr.Contains("EJ2"))
 355         triggerstrings.push_back("EMCLowGammaOnly");
 356     }
 357
 358     try{
 359       fHistos->FillTHnSparse("hEventTriggers", triggers);
 360     } catch (HistoContainerContentException &e){
 361       std::stringstream errormessage;
 362       errormessage << "Filling of histogram failed: " << e.what();
 363       AliError(errormessage.str().c_str());
 364     }
 365
 366     // apply event selection: Combine the Pileup cut from SPD with the other pA Vertex selection cuts.
 367     bool isPileupEvent = fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5.);
 368     isPileupEvent = isPileupEvent || (TMath::Abs(vtxTracks->GetZ() - vtxSPD->GetZ()) > 0.5);
 369     double covSPD[6]; vtxSPD->GetCovarianceMatrix(covSPD);
 370     isPileupEvent = isPileupEvent || (TString(vtxSPD->GetTitle()).Contains("vertexer:Z") && TMath::Sqrt(covSPD[5]) > 0.25);
 371
 372     // Fill event-based histogram
 373     const double &zv = vtxTracks->GetZ();
 374     if(triggers[0]) FillEventHist("MinBias", zv, isPileupEvent);
 375     if(!triggerstrings.size()) // Non-EMCal-triggered
 376       FillEventHist("NoEMCal", zv, isPileupEvent);
 377     else{
 378       // EMCal-triggered events
 379       for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
 380         FillEventHist(it->c_str(), zv, isPileupEvent);
 381     }
 382
 383     // Fill MC truth
 384     if(fMCEvent){
 385       for(int ipart = 0; ipart < fMCEvent->GetNumberOfTracks(); ipart++){
 386         // Select only physical primary particles
 387         AliVParticle *part = fMCEvent->GetTrack(ipart);
 388         if(!fEtaRange.IsInRange(part->Eta())) continue;
 389         if(!fPtRange.IsInRange(part->Pt())) continue;
 390         if(!fMCEvent->IsPhysicalPrimary(ipart)) continue;
 391         FillMCParticleHist(part, zv, isPileupEvent);
 392       }
 393     }
 394
 395     AliVTrack *track(NULL);
 396     // Loop over all tracks (No cuts applied)
 397     TIter allTrackIter(fTracks);
 398     while((track = dynamic_cast<AliVTrack *>(allTrackIter()))){
 399       if(!IsTrueTrack(track)) continue;
 400       if(!fEtaRange.IsInRange(track->Eta())) continue;
 401       if(!fPtRange.IsInRange(track->Pt())) continue;
 402       if(triggers[0]) FillTrackHist("MinBias", track, zv, isPileupEvent, 0, triggers[0]);
 403       if(!triggerstrings.size()) // Non-EMCal-triggered
 404         FillTrackHist("NoEMCal", track, zv, isPileupEvent, 0, triggers[0]);
 405       else {
 406         // EMCal-triggered events
 407         for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
 408           FillTrackHist(it->c_str(), track, zv, isPileupEvent, 0, triggers[0]);
 409       }
 410     }
 411
 412     // Now apply track selection cuts
 413     // allow for several track selections to be tested at the same time
 414     // each track selection gets a different cut ID starting from 1
 415     // cut ID 0 is reserved for the case of no cuts
 416     if(fListTrackCuts && fListTrackCuts->GetEntries()){
 417       for(int icut = 0; icut < fListTrackCuts->GetEntries(); icut++){
 418         AliEMCalPtTaskVTrackSelection *trackSelection = static_cast<AliEMCalPtTaskVTrackSelection *>(fListTrackCuts->At(icut));
 419         TIter trackIter(trackSelection->GetAcceptedTracks(fTracks));
 420         while((track = dynamic_cast<AliVTrack *>(trackIter()))){
 421           //if(!IsTrueTrack(track)) continue;
 422           if(!fEtaRange.IsInRange(track->Eta())) continue;
 423           if(!fPtRange.IsInRange(track->Pt())) continue;
 424           if(triggers[0]) FillTrackHist("MinBias", track, zv, isPileupEvent, icut + 1, triggers[0]);
 425           if(!triggerstrings.size()) // Non-EMCal-triggered
 426             FillTrackHist("NoEMCal", track, zv, isPileupEvent, icut + 1, triggers[0]);
 427           else {
 428             // EMCal-triggered events
 429             for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
 430               FillTrackHist(it->c_str(), track, zv, isPileupEvent, icut + 1, triggers[0]);
 431           }
 432         }
 433       }
 434     }
 435
 436     // Next step we loop over the (uncalibrated) emcal clusters and fill histograms with the cluster energy
 437     const AliVCluster *clust(NULL);
 438     for(int icl = 0; icl < fInputEvent->GetNumberOfCaloClusters(); icl++){
 439       clust = fInputEvent->GetCaloCluster(icl);
 440       if(!clust->IsEMCAL()) continue;
 441       if(triggers[0]) FillClusterHist("MinBias", clust, false, zv, isPileupEvent, triggers[0]);
 442       if(!triggerstrings.size())        // Non-EMCal-triggered
 443         FillClusterHist("NoEMCal", clust, false, zv, isPileupEvent, triggers[0]);
 444       else{
 445         for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it){
 446           FillClusterHist(it->c_str(), clust, false, zv, isPileupEvent, triggers[0]);
 447         }
 448       }
 449     }
 450
 451     if(fCaloClusters){
 452       TIter clustIter(fCaloClusters);
 453       while((clust = dynamic_cast<const AliVCluster *>(clustIter()))){
 454         if(!clust->IsEMCAL()) continue;
 455         if(triggers[0]) FillClusterHist("MinBias", clust, true, zv, isPileupEvent, triggers[0]);
 456         if(!triggerstrings.size())      // Non-EMCal-triggered
 457           FillClusterHist("NoEMCal", clust, true, zv, isPileupEvent, triggers[0]);
 458         else{
 459           for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it){
 460             FillClusterHist(it->c_str(), clust, true, zv, isPileupEvent, triggers[0]);
 461           }
 462         }
 463       }
 464     }
 465
 466     PostData(1, fOutput);
 467     return true;
 468   }
 469
 470   //______________________________________________________________________________
 471   void AliAnalysisTaskPtEMCalTrigger::CreateDefaultPtBinning(TArrayD &binning) const{
 472     /*
 473      * Creating the default pt binning.
 474      *
 475      * @param binning: Array where to store the results.
 476      */
 477     std::vector<double> mybinning;
 478     std::map<double,double> definitions;
 479     definitions.insert(std::pair<double,double>(2.5, 0.1));
 480     definitions.insert(std::pair<double,double>(7., 0.25));
 481     definitions.insert(std::pair<double,double>(15., 0.5));
 482     definitions.insert(std::pair<double,double>(25., 1.));
 483     definitions.insert(std::pair<double,double>(40., 2.5));
 484     definitions.insert(std::pair<double,double>(60., 5.));
 485     definitions.insert(std::pair<double,double>(100., 5.));
 486     double currentval = 0;
 487     for(std::map<double,double>::iterator id = definitions.begin(); id != definitions.end(); ++id){
 488       double limit = id->first, binwidth = id->second;
 489       while(currentval <= limit){
 490         currentval += binwidth;
 491         mybinning.push_back(currentval);
 492       }
 493     }
 494     binning.Set(mybinning.size());
 495     int ib = 0;
 496     for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
 497       binning[ib++] = *it;
 498   }
 499
 500   //______________________________________________________________________________
 501   void AliAnalysisTaskPtEMCalTrigger::CreateDefaultZVertexBinning(TArrayD &binning) const {
 502     /*
 503      * Creating default z-Vertex binning.
 504      *
 505      * @param binning: Array where to store the results.
 506      */
 507     std::vector<double> mybinning;
 508     double currentval = -40;
 509     mybinning.push_back(currentval);
 510     while(currentval <= 40.){
 511       currentval += 1.;
 512       mybinning.push_back(currentval);
 513     }
 514     binning.Set(mybinning.size());
 515     int ib = 0;
 516     for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
 517       binning[ib++] = *it;
 518   }
 519
 520   //______________________________________________________________________________
 521   void AliAnalysisTaskPtEMCalTrigger::CreateDefaultEtaBinning(TArrayD& binning) const {
 522     /*
 523      * Creating default z-Vertex binning.
 524      *
 525      * @param binning: Array where to store the results.
 526      */
 527     std::vector<double> mybinning;
 528     double currentval = -0.8;
 529     mybinning.push_back(currentval);
 530     while(currentval <= 0.8){
 531       currentval += 0.1;
 532       mybinning.push_back(currentval);
 533     }
 534     binning.Set(mybinning.size());
 535     int ib = 0;
 536     for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
 537       binning[ib++] = *it;
 538   }
 539
 540   //______________________________________________________________________________
 541   void AliAnalysisTaskPtEMCalTrigger::DefineAxis(TAxis& axis, const char* name,
 542       const char* title, const TArrayD& binning, const char** labels) {
 543     /*
 544      * Define an axis with a given binning
 545      *
 546      * @param axis: Axis to be defined
 547      * @param name: Name of the axis
 548      * @param title: Title of the axis
 549      * @param binning: axis binning
 550      * @param labels (@optional): array of bin labels
 551      */
 552     axis.Set(binning.GetSize()-1, binning.GetArray());
 553     axis.SetName(name);
 554     axis.SetTitle(title);
 555     if(labels){
 556       for(int ib = 1; ib <= axis.GetNbins(); ++ib)
 557         axis.SetBinLabel(ib, labels[ib-1]);
 558     }
 559   }
 560
 561   //______________________________________________________________________________
 562   void AliAnalysisTaskPtEMCalTrigger::DefineAxis(TAxis& axis, const char* name,
 563       const char* title, int nbins, double min, double max,
 564       const char** labels) {
 565     /*
 566      * Define an axis with number of bins from min to max
 567      *
 568      * @param axis: Axis to be defined
 569      * @param name: Name of the axis
 570      * @param title: Title of the axis
 571      * @param nbins: Number of bins
 572      * @param min: lower limit of the axis
 573      * @param max: upper limit of the axis
 574      * @param labels (@optional): array of bin labels
 575      */
 576     axis.Set(nbins, min, max);
 577     axis.SetName(name);
 578     axis.SetTitle(title);
 579     if(labels){
 580       for(int ib = 1; ib <= axis.GetNbins(); ++ib)
 581         axis.SetBinLabel(ib, labels[ib-1]);
 582     }
 583   }
 584
 585
 586   //______________________________________________________________________________
 587   void AliAnalysisTaskPtEMCalTrigger::FillEventHist(const char* trigger,
 588       double vz, bool isPileup) {
 589     /*
 590      * Fill event-based histogram
 591      *
 592      * @param trigger: name of the trigger configuration to be processed
 593      * @param vz: z-position of the vertex
 594      * @param isPileup: signalises if the event is flagged as pileup event
 595      */
 596     char histname[1024];
 597     sprintf(histname, "hEventHist%s", trigger);
 598     try{
 599       fHistos->FillTH2(histname, 0., vz);
 600     } catch (HistoContainerContentException &e){
 601       std::stringstream errormessage;
 602       errormessage << "Filling of histogram failed: " << e.what();
 603       AliError(errormessage.str().c_str());
 604     }
 605     if(!isPileup){
 606       try{
 607         fHistos->FillTH2(histname, 1., vz);
 608       } catch(HistoContainerContentException &e){
 609         std::stringstream errormessage;
 610         errormessage << "Filling of histogram failed: " << e.what();
 611         AliError(errormessage.str().c_str());
 612       }
 613     }
 614   }
 615
 616   //______________________________________________________________________________
 617   void AliAnalysisTaskPtEMCalTrigger::FillTrackHist(const char* trigger,
 618       const AliVTrack* track, double vz, bool isPileup, int cut, bool isMinBias) {
 619     /*
 620      * Fill track-based histogram with corresponding information
 621      *
 622      * @param trigger: name of the trigger
 623      * @param track: ESD track selected
 624      * @param vz: z-position of the vertex
 625      * @param isPileup: flag event as pileup event
 626      * @param cut: id of the cut (0 = no cut)
 627      */
 628     double etasign = fSwapEta ? -1. : 1.;
 629     double data[7] = {TMath::Abs(track->Pt()), etasign * track->Eta(), track->Phi(), vz, 0, static_cast<double>(cut), isMinBias ? 1. : 0.};
 630     double dataMC[7] = {0., 0., 0., vz, 0, static_cast<double>(cut), isMinBias ? 1. : 0.};
 631     AliVParticle *assocMC(NULL);
 632     if(fMCEvent && (assocMC = fMCEvent->GetTrack(TMath::Abs(track->GetLabel())))){
 633       dataMC[0] = TMath::Abs(assocMC->Pt());
 634       dataMC[1] = etasign * assocMC->Eta();
 635       dataMC[2] = assocMC->Phi();
 636     }
 637     char histname[1024], histnameAcc[1024], histnameMC[1024], histnameMCAcc[1024];
 638     sprintf(histname, "hTrackHist%s", trigger);
 639     sprintf(histnameAcc, "hTrackInAcceptanceHist%s", trigger);
 640     sprintf(histnameMC, "hMCTrackHist%s", trigger);
 641     sprintf(histnameMCAcc, "hMCTrackInAcceptanceHist%s", trigger);
 642     Bool_t isEMCAL = kFALSE;
 643     if(track->IsEMCAL()){
 644       // Check if the cluster is matched to only one track
 645       AliVCluster *emcclust(NULL);
 646       AliDebug(2, Form("cluster id: %d\n", track->GetEMCALcluster()));
 647       if(fCaloClusters) {
 648         AliDebug(2, "Using calibrated clusters");
 649         emcclust = dynamic_cast<AliVCluster *>(fCaloClusters->At(track->GetEMCALcluster()));
 650       } else {
 651         AliDebug(2, "Using uncalibrated clusters");
 652         emcclust = fInputEvent->GetCaloCluster(track->GetEMCALcluster());
 653       }
 654       if(!emcclust) AliError("Null pointer to EMCal cluster");
 655       if(emcclust && emcclust->GetNTracksMatched() <= 1){
 656         isEMCAL = kTRUE;
 657       }
 658     }
 659     try{
 660       fHistos->FillTHnSparse(histname, data);
 661       if(fMCEvent) fHistos->FillTHnSparse(histnameMC, dataMC);
 662       if(isEMCAL){
 663         fHistos->FillTHnSparse(histnameAcc, data);
 664         if(fMCEvent) fHistos->FillTHnSparse(histnameMCAcc, dataMC);
 665       }
 666     } catch (HistoContainerContentException &e){
 667       std::stringstream errormessage;
 668       errormessage << "Filling of histogram failed: " << e.what();
 669       AliError(errormessage.str().c_str());
 670     }
 671     if(!isPileup){
 672       data[4] = 1;
 673       dataMC[4] = 1;
 674       try{
 675         fHistos->FillTHnSparse(histname, data);
 676         if(fMCEvent) fHistos->FillTHnSparse(histnameMC, dataMC);
 677         if(isEMCAL){
 678           fHistos->FillTHnSparse(histnameAcc, data);
 679           if(fMCEvent) fHistos->FillTHnSparse(histnameMCAcc, dataMC);
 680         }
 681       } catch (HistoContainerContentException &e){
 682         std::stringstream errormessage;
 683         errormessage << "Filling of histogram failed: " << e.what();
 684         AliError(errormessage.str().c_str());
 685       }
 686     }
 687   }
 688
 689   //______________________________________________________________________________
 690   void AliAnalysisTaskPtEMCalTrigger::FillClusterHist(const char* trigger,
 691       const AliVCluster* clust, bool isCalibrated, double vz, bool isPileup, bool isMinBias) {
 692     /*
 693      * Fill cluster-based histogram with corresponding information
 694      *
 695      * @param trigger: name of the trigger
 696      * @param cluster: the EMCal cluster information
 697      * @param vz: z-position of the vertex
 698      * @param isPileup: flag event as pileup event
 699      */
 700     double data[4] =  {clust->E(), vz, 0, isMinBias ? 1. : 0.};
 701     char histname[1024];
 702     sprintf(histname, "hCluster%sHist%s", isCalibrated ? "Calib" : "Uncalib", trigger);
 703     try{
 704       fHistos->FillTHnSparse(histname, data);
 705     } catch (HistoContainerContentException &e){
 706       std::stringstream errormessage;
 707       errormessage << "Filling of histogram failed: " << e.what();
 708       AliError(errormessage.str().c_str());
 709     }
 710     if(!isPileup){
 711       data[2] = 1.;
 712       try{
 713         fHistos->FillTHnSparse(histname, data);
 714       } catch (HistoContainerContentException &e){
 715         std::stringstream errormessage;
 716         errormessage << "Filling of histogram failed: " << e.what();
 717         AliError(errormessage.str().c_str());
 718       }
 719     }
 720   }
 721
 722   //______________________________________________________________________________
 723   void AliAnalysisTaskPtEMCalTrigger::FillMCParticleHist(const AliVParticle * const track, double vz, bool isPileup){
 724     /*
 725      * Fill histogram for MC-true particles with the information pt, eta and phi
 726      *
 727      * @param track: the Monte-Carlo track
 728      */
 729     double data[5] = {TMath::Abs(track->Pt()), track->Eta(), track->Phi(), vz, 0.};
 730     fHistos->FillTHnSparse("hMCtrueParticles", data);
 731     if(!isPileup){
 732       data[4] = 1.;
 733       fHistos->FillTHnSparse("hMCtrueParticles", data);
 734     }
 735   }
 736
 737   //______________________________________________________________________________
 738   bool AliAnalysisTaskPtEMCalTrigger::IsTrueTrack(const AliVTrack *const track) const{
 739     /*
 740      * Check if the track has an associated MC particle, and that the particle is a physical primary
 741      * In case of data we do not do the selection at that step (always return true)
 742      *
 743      * @param track: Track to check
 744      * @result: true primary track (true or false)
 745      */
 746     if(!fMCEvent) return true;
 747     AliVParticle *mcassociate = fMCEvent->GetTrack(TMath::Abs(track->GetLabel()));
 748     if(!mcassociate) return false;
 749     return fMCEvent->IsPhysicalPrimary(TMath::Abs(track->GetLabel()));
 750   }
 751
 752   //______________________________________________________________________________
 753   void AliAnalysisTaskPtEMCalTrigger::AddESDTrackCuts(AliESDtrackCuts* trackCuts) {
 754     /*
 755      * Add new track cuts to the task
 756      *
 757      * @param trackCuts: Object of type AliESDtrackCuts
 758      */
 759     fListTrackCuts->AddLast(new AliEMCalPtTaskTrackSelectionESD(trackCuts));
 760   }
 761
 762   //______________________________________________________________________________
 763   void AliAnalysisTaskPtEMCalTrigger::AddCutsForAOD(AliESDtrackCuts* trackCuts, UInt_t filterbits) {
 764     /*
 765      * Add new track cuts to the task
 766      *
 767      * @param trackCuts: Object of type AliESDtrackCuts
 768      */
 769     fListTrackCuts->AddLast(new AliEMCalPtTaskTrackSelectionAOD(trackCuts, filterbits));
 770   }
 771
 772
 773   //______________________________________________________________________________
 774   TString AliAnalysisTaskPtEMCalTrigger::BuildTriggerString() {
 775     /*
 776      * Build trigger string from the trigger maker
 777      *
 778      * @return: blank-separated string of fired trigger classes
 779      */
 780     AliDebug(1, "trigger checking");
 781     TString result = "";
 782     if(HasTriggerType(kJ1)) result += "EJ1 ";
 783     if(HasTriggerType(kJ2)) result += "EJ2 ";
 784     if(HasTriggerType(kG1)) result += "EG1 ";
 785     if(HasTriggerType(kG2)) result += "EG2 ";
 786     return result;
 787   }
 788
 789   //______________________________________________________________________________
 790   const AliVVertex* AliAnalysisTaskPtEMCalTrigger::GetSPDVertex() const {
 791     /*
 792      * Accessor for the SPD vertex, creating transparency for ESDs and AODs
 793      *
 794      * @return: the spd vertex
 795      */
 796     AliESDEvent *esd = dynamic_cast<AliESDEvent *>(fInputEvent);
 797     if(esd){
 798       return esd->GetPrimaryVertexSPD();
 799     } else {
 800       AliAODEvent *aod = dynamic_cast<AliAODEvent *>(fInputEvent);
 801       if(aod){
 802         return aod->GetPrimaryVertexSPD();
 803       }
 804     }
 805     return NULL;
 806   }
 807
 808 }
 809