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 hpatchaxes[3];
 188                 DefineAxis(hpatchaxes[0], "energy", "Patch energy (GeV)", 100, 0., 100.);
 189                 DefineAxis(hpatchaxes[1], "eta", "#eta", etabinning);
 190                 DefineAxis(hpatchaxes[2], "phi", "#phi",  20, 0, 2 * TMath::Pi());
 191                 const TAxis *patchaxes[3];
 192                 for(int iaxis = 0; iaxis < 3; ++iaxis) patchaxes[iaxis] = hpatchaxes + iaxis;
 193                 try{
 194                         std::string patchnames[] = {"Level0", "JetHigh", "JetLow", "GammaHigh", "GammaLow"};
 195                         for(std::string * triggerpatch = patchnames; triggerpatch < patchnames + sizeof(patchnames)/sizeof(std::string); ++triggerpatch){
 196                                 fHistos->CreateTHnSparse(Form("Energy%s", triggerpatch->c_str()), Form("Patch energy for %s trigger patches", triggerpatch->c_str()), 3, patchaxes);
 197                                 fHistos->CreateTHnSparse(Form("EnergyMain%s", triggerpatch->c_str()), Form("Patch energy for main %s trigger patches", triggerpatch->c_str()), 3, patchaxes);
 198                         }
 199
 200                         // Create histogram for MC-truth
 201                         fHistos->CreateTHnSparse("hMCtrueParticles", "Particle-based histogram for MC-true particles", 5, trackaxes);
 202                         for(std::map<std::string,std::string>::iterator it = triggerCombinations.begin(); it != triggerCombinations.end(); ++it){
 203                                 const std::string name = it->first, &title = it->second;
 204                                 // Create event-based histogram
 205                                 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);
 206                                 // Create track-based histogram
 207                                 fHistos->CreateTHnSparse(Form("hTrackHist%s", name.c_str()), Form("Track-based data for %s events", title.c_str()), 7, trackaxes);
 208                                 fHistos->CreateTHnSparse(Form("hTrackInAcceptanceHist%s", name.c_str()), Form("Track-based data for %s events", title.c_str()), 7, trackaxes);
 209                                 fHistos->CreateTHnSparse(Form("hMCTrackHist%s", name.c_str()), Form("Track-based data for %s events with MC kinematics", title.c_str()), 7, trackaxes);
 210                                 fHistos->CreateTHnSparse(Form("hMCTrackInAcceptanceHist%s", name.c_str()), Form("Track-based data for %s events with MC kinematics", title.c_str()), 7, trackaxes);
 211                                 // Create cluster-based histogram (Uncalibrated and calibrated clusters)
 212                                 fHistos->CreateTHnSparse(Form("hClusterCalibHist%s", name.c_str()), Form("Calib. cluster-based histogram for %s events", title.c_str()), 4, clusteraxes);
 213                                 fHistos->CreateTHnSparse(Form("hClusterUncalibHist%s", name.c_str()), Form("Uncalib. cluster-based histogram for %s events", title.c_str()), 4, clusteraxes);
 214                         }
 215                         fHistos->CreateTHnSparse("hEventTriggers", "Trigger type per event", 5, triggeraxis);
 216                         fHistos->CreateTHnSparse("hEventsTriggerbit", "Trigger bits for the different events", 4, bitaxes);
 217                 } catch (HistoContainerContentException &e){
 218                         std::stringstream errormessage;
 219                         errormessage << "Creation of histogram failed: " << e.what();
 220                         AliError(errormessage.str().c_str());
 221                 }
 222                 fOutput->Add(fHistos->GetListOfHistograms());
 223                 if(fListTrackCuts && fListTrackCuts->GetEntries()){
 224                         TIter cutIter(fListTrackCuts);
 225                         AliEMCalPtTaskVTrackSelection *cutObject(NULL);
 226                         while((cutObject = dynamic_cast<AliEMCalPtTaskVTrackSelection *>(cutIter()))){
 227                                 AliESDtrackCuts *cuts = dynamic_cast<AliESDtrackCuts *>(cutObject->GetTrackCuts());
 228                                 if(cuts){
 229                                         cuts->DefineHistograms();
 230                                         fOutput->Add(cuts);
 231                                 }
 232                         }
 233                 }
 234                 PostData(1, fOutput);
 235         }
 236
 237         //______________________________________________________________________________
 238         Bool_t AliAnalysisTaskPtEMCalTrigger::Run(){
 239                 /*
 240                  * Runs the event loop
 241                  *
 242                  * @param option: Additional options
 243                  */
 244                 // Common checks: Have SPD vertex and primary vertex from tracks, and both need to have at least one contributor
 245                 AliDebug(1,Form("Number of calibrated clusters: %d", fCaloClusters->GetEntries()));
 246                 AliDebug(1,Form("Number of matched tracks: %d", fTracks->GetEntries()));
 247                 if(fMCEvent){
 248                         // Build always trigger strig from the trigger maker in case of MC
 249                         fUseTriggersFromTriggerMaker = kTRUE;
 250                 }
 251
 252                 // Loop over trigger patches, fill patch energy
 253                 AliEmcalTriggerPatchInfo *triggerpatch(NULL);
 254                 TIter patchIter(this->fTriggerPatchInfo);
 255                 while((triggerpatch = dynamic_cast<AliEmcalTriggerPatchInfo *>(patchIter()))){
 256                         double triggerpatchinfo[3] = {triggerpatch->GetPatchE(), triggerpatch->GetEtaCM(), triggerpatch->GetPhiCM()};
 257                         if(triggerpatch->IsJetHigh()){
 258                                 fHistos->FillTHnSparse("EnergyJetHigh", triggerpatchinfo);
 259                                 if(triggerpatch->IsMainTrigger())
 260                                         fHistos->FillTHnSparse("EnergyMainJetHigh", triggerpatchinfo);
 261                         }
 262                         if(triggerpatch->IsJetLow()){
 263                                 fHistos->FillTHnSparse("EnergyJetLow", triggerpatchinfo);
 264                                 if(triggerpatch->IsMainTrigger())
 265                                         fHistos->FillTHnSparse("EnergyMainJetLow", triggerpatchinfo);
 266                         }
 267                         if(triggerpatch->IsGammaHigh()){
 268                                 fHistos->FillTHnSparse("EnergyGammaHigh", triggerpatchinfo);
 269                                 if(triggerpatch->IsMainTrigger())
 270                                         fHistos->FillTHnSparse("EnergyMainGammaHigh", triggerpatchinfo);
 271                         }
 272                         if(triggerpatch->IsGammaLow()){
 273                                 fHistos->FillTHnSparse("EnergyGammaLow", triggerpatchinfo);
 274                                 if(triggerpatch->IsMainTrigger())
 275                                         fHistos->FillTHnSparse("EnergyMainGammaLow", triggerpatchinfo);
 276                         }
 277                         if(triggerpatch->IsLevel0()){
 278                                 fHistos->FillTHnSparse("EnergyLevel0", triggerpatchinfo);
 279                                 if(triggerpatch->IsMainTrigger())
 280                                         fHistos->FillTHnSparse("EnergyMainLevel0", triggerpatchinfo);
 281                         }
 282                 }
 283
 284                 const AliVVertex *vtxTracks = fInputEvent->GetPrimaryVertex(),
 285                                 *vtxSPD = GetSPDVertex();
 286                 if(!(vtxTracks && vtxSPD)) return false;
 287                 if(vtxTracks->GetNContributors() < 1 || vtxSPD->GetNContributors() < 1) return false;
 288
 289                 double triggers[5]; memset(triggers, 0, sizeof(double) * 5);
 290                 double triggerbits[4]; memset(triggerbits, 0, sizeof(double) * 4);
 291                 if(fInputHandler->IsEventSelected() & AliVEvent::kINT7){
 292                         triggers[0] = 1.;
 293                         triggerbits[0] = 1.;
 294                 }
 295
 296                 // check triggerbits
 297                 if(fInputHandler->IsEventSelected() & AliVEvent::kEMC7){
 298                         triggerbits[1] = 1.;
 299                 }
 300                 if(fInputHandler->IsEventSelected() & AliVEvent::kEMCEGA){
 301                         triggerbits[2] = 1.;
 302                 }
 303                 if(fInputHandler->IsEventSelected() & AliVEvent::kEMCEJE){
 304                         triggerbits[3] = 1.;
 305                 }
 306                 try{
 307                         fHistos->FillTHnSparse("hEventsTriggerbit", triggerbits);
 308                 } catch(HistoContainerContentException &e) {
 309                         std::stringstream errormessage;
 310                         errormessage << "Filling of histogram failed: " << e.what();
 311                         AliError(errormessage.str().c_str());
 312                 }
 313
 314                 std::vector<std::string> triggerstrings;
 315                 // EMCal-triggered event, distinguish types
 316                 TString trgstr(fUseTriggersFromTriggerMaker ? BuildTriggerString() : fInputEvent->GetFiredTriggerClasses());
 317                 AliDebug(1, Form("Triggerstring: %s\n", trgstr.Data()));
 318                 if(trgstr.Contains("EJ1")){
 319                         triggerstrings.push_back("EMCJHigh");
 320                         triggers[1] = 1;
 321                         if(trgstr.Contains("EG1"))
 322                                 triggerstrings.push_back("EMCHighBoth");
 323                         else
 324                                 triggerstrings.push_back("EMCHighJetOnly");
 325                 }
 326                 if(trgstr.Contains("EJ2")){
 327                         triggerstrings.push_back("EMCJLow");
 328                         triggers[2] = 1;
 329                         if(trgstr.Contains("EG2"))
 330                                 triggerstrings.push_back("EMCLowBoth");
 331                         else
 332                                 triggerstrings.push_back("EMCLowJetOnly");
 333                 }
 334                 if(trgstr.Contains("EG1")){
 335                         triggerstrings.push_back("EMCGHigh");
 336                         triggers[3] = 1;
 337                         if(!trgstr.Contains("EJ1"))
 338                                 triggerstrings.push_back("EMCHighGammaOnly");
 339                 }
 340                 if(trgstr.Contains("EG2")){
 341                         triggerstrings.push_back("EMCGLow");
 342                         triggers[4] = 1;
 343                         if(!trgstr.Contains("EJ2"))
 344                                 triggerstrings.push_back("EMCLowGammaOnly");
 345                 }
 346
 347                 try{
 348                         fHistos->FillTHnSparse("hEventTriggers", triggers);
 349                 } catch (HistoContainerContentException &e){
 350                         std::stringstream errormessage;
 351                         errormessage << "Filling of histogram failed: " << e.what();
 352                         AliError(errormessage.str().c_str());
 353                 }
 354
 355                 // apply event selection: Combine the Pileup cut from SPD with the other pA Vertex selection cuts.
 356                 bool isPileupEvent = fInputEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5.);
 357                 isPileupEvent = isPileupEvent || (TMath::Abs(vtxTracks->GetZ() - vtxSPD->GetZ()) > 0.5);
 358                 double covSPD[6]; vtxSPD->GetCovarianceMatrix(covSPD);
 359                 isPileupEvent = isPileupEvent || (TString(vtxSPD->GetTitle()).Contains("vertexer:Z") && TMath::Sqrt(covSPD[5]) > 0.25);
 360
 361                 // Fill event-based histogram
 362                 const double &zv = vtxTracks->GetZ();
 363                 if(triggers[0]) FillEventHist("MinBias", zv, isPileupEvent);
 364                 if(!triggerstrings.size()) // Non-EMCal-triggered
 365                         FillEventHist("NoEMCal", zv, isPileupEvent);
 366                 else{
 367                         // EMCal-triggered events
 368                         for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
 369                                 FillEventHist(it->c_str(), zv, isPileupEvent);
 370                 }
 371
 372                 // Fill MC truth
 373                 if(fMCEvent){
 374                         for(int ipart = 0; ipart < fMCEvent->GetNumberOfTracks(); ipart++){
 375                                 // Select only physical primary particles
 376                                 AliVParticle *part = fMCEvent->GetTrack(ipart);
 377                                 if(!fEtaRange.IsInRange(part->Eta())) continue;
 378                                 if(!fPtRange.IsInRange(part->Pt())) continue;
 379                                 if(!fMCEvent->IsPhysicalPrimary(ipart)) continue;
 380                                 FillMCParticleHist(part, zv, isPileupEvent);
 381                         }
 382                 }
 383
 384                 AliVTrack *track(NULL);
 385                 // Loop over all tracks (No cuts applied)
 386                 TIter allTrackIter(fTracks);
 387                 while((track = dynamic_cast<AliVTrack *>(allTrackIter()))){
 388                         if(!IsTrueTrack(track)) continue;
 389                         if(!fEtaRange.IsInRange(track->Eta())) continue;
 390                         if(!fPtRange.IsInRange(track->Pt())) continue;
 391                         if(triggers[0]) FillTrackHist("MinBias", track, zv, isPileupEvent, 0, triggers[0]);
 392                         if(!triggerstrings.size()) // Non-EMCal-triggered
 393                                 FillTrackHist("NoEMCal", track, zv, isPileupEvent, 0, triggers[0]);
 394                         else {
 395                                 // EMCal-triggered events
 396                                 for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
 397                                         FillTrackHist(it->c_str(), track, zv, isPileupEvent, 0, triggers[0]);
 398                         }
 399                 }
 400
 401                 // Now apply track selection cuts
 402                 // allow for several track selections to be tested at the same time
 403                 // each track selection gets a different cut ID starting from 1
 404                 // cut ID 0 is reserved for the case of no cuts
 405                 if(fListTrackCuts && fListTrackCuts->GetEntries()){
 406                         for(int icut = 0; icut < fListTrackCuts->GetEntries(); icut++){
 407                                 AliEMCalPtTaskVTrackSelection *trackSelection = static_cast<AliEMCalPtTaskVTrackSelection *>(fListTrackCuts->At(icut));
 408                                 TIter trackIter(trackSelection->GetAcceptedTracks(fTracks));
 409                                 while((track = dynamic_cast<AliVTrack *>(trackIter()))){
 410                                         //if(!IsTrueTrack(track)) continue;
 411                                         if(!fEtaRange.IsInRange(track->Eta())) continue;
 412                                         if(!fPtRange.IsInRange(track->Pt())) continue;
 413                                         if(triggers[0]) FillTrackHist("MinBias", track, zv, isPileupEvent, icut + 1, triggers[0]);
 414                                         if(!triggerstrings.size()) // Non-EMCal-triggered
 415                                                 FillTrackHist("NoEMCal", track, zv, isPileupEvent, icut + 1, triggers[0]);
 416                                         else {
 417                                                 // EMCal-triggered events
 418                                                 for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
 419                                                         FillTrackHist(it->c_str(), track, zv, isPileupEvent, icut + 1, triggers[0]);
 420                                         }
 421                                 }
 422                         }
 423                 }
 424
 425                 // Next step we loop over the (uncalibrated) emcal clusters and fill histograms with the cluster energy
 426                 const AliVCluster *clust(NULL);
 427                 for(int icl = 0; icl < fInputEvent->GetNumberOfCaloClusters(); icl++){
 428                         clust = fInputEvent->GetCaloCluster(icl);
 429                         if(!clust->IsEMCAL()) continue;
 430                         if(triggers[0]) FillClusterHist("MinBias", clust, false, zv, isPileupEvent, triggers[0]);
 431                         if(!triggerstrings.size())      // Non-EMCal-triggered
 432                                 FillClusterHist("NoEMCal", clust, false, zv, isPileupEvent, triggers[0]);
 433                         else{
 434                                 for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it){
 435                                         FillClusterHist(it->c_str(), clust, false, zv, isPileupEvent, triggers[0]);
 436                                 }
 437                         }
 438                 }
 439
 440                 if(fCaloClusters){
 441                         TIter clustIter(fCaloClusters);
 442                         while((clust = dynamic_cast<const AliVCluster *>(clustIter()))){
 443                                 if(!clust->IsEMCAL()) continue;
 444                                 if(triggers[0]) FillClusterHist("MinBias", clust, true, zv, isPileupEvent, triggers[0]);
 445                                 if(!triggerstrings.size())      // Non-EMCal-triggered
 446                                         FillClusterHist("NoEMCal", clust, true, zv, isPileupEvent, triggers[0]);
 447                                 else{
 448                                         for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it){
 449                                                 FillClusterHist(it->c_str(), clust, true, zv, isPileupEvent, triggers[0]);
 450                                         }
 451                                 }
 452                         }
 453                 }
 454
 455                 PostData(1, fOutput);
 456                 return true;
 457         }
 458
 459         //______________________________________________________________________________
 460         void AliAnalysisTaskPtEMCalTrigger::CreateDefaultPtBinning(TArrayD &binning) const{
 461                 /*
 462                  * Creating the default pt binning.
 463                  *
 464                  * @param binning: Array where to store the results.
 465                  */
 466                 std::vector<double> mybinning;
 467                 std::map<double,double> definitions;
 468                 definitions.insert(std::pair<double,double>(2.5, 0.1));
 469                 definitions.insert(std::pair<double,double>(7., 0.25));
 470                 definitions.insert(std::pair<double,double>(15., 0.5));
 471                 definitions.insert(std::pair<double,double>(25., 1.));
 472                 definitions.insert(std::pair<double,double>(40., 2.5));
 473                 definitions.insert(std::pair<double,double>(60., 5.));
 474                 definitions.insert(std::pair<double,double>(100., 5.));
 475                 double currentval = 0;
 476                 for(std::map<double,double>::iterator id = definitions.begin(); id != definitions.end(); ++id){
 477                         double limit = id->first, binwidth = id->second;
 478                         while(currentval <= limit){
 479                                 currentval += binwidth;
 480                                 mybinning.push_back(currentval);
 481                         }
 482                 }
 483                 binning.Set(mybinning.size());
 484                 int ib = 0;
 485                 for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
 486                         binning[ib++] = *it;
 487         }
 488
 489         //______________________________________________________________________________
 490         void AliAnalysisTaskPtEMCalTrigger::CreateDefaultZVertexBinning(TArrayD &binning) const {
 491                 /*
 492                  * Creating default z-Vertex binning.
 493                  *
 494                  * @param binning: Array where to store the results.
 495                  */
 496                 std::vector<double> mybinning;
 497                 double currentval = -40;
 498                 mybinning.push_back(currentval);
 499                 while(currentval <= 40.){
 500                         currentval += 1.;
 501                         mybinning.push_back(currentval);
 502                 }
 503                 binning.Set(mybinning.size());
 504                 int ib = 0;
 505                 for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
 506                         binning[ib++] = *it;
 507         }
 508
 509         //______________________________________________________________________________
 510         void AliAnalysisTaskPtEMCalTrigger::CreateDefaultEtaBinning(TArrayD& binning) const {
 511                 /*
 512                  * Creating default z-Vertex binning.
 513                  *
 514                  * @param binning: Array where to store the results.
 515                  */
 516                 std::vector<double> mybinning;
 517                 double currentval = -0.8;
 518                 mybinning.push_back(currentval);
 519                 while(currentval <= 0.8){
 520                         currentval += 0.1;
 521                         mybinning.push_back(currentval);
 522                 }
 523                 binning.Set(mybinning.size());
 524                 int ib = 0;
 525                 for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
 526                         binning[ib++] = *it;
 527         }
 528
 529         //______________________________________________________________________________
 530         void AliAnalysisTaskPtEMCalTrigger::DefineAxis(TAxis& axis, const char* name,
 531                         const char* title, const TArrayD& binning, const char** labels) {
 532                 /*
 533                  * Define an axis with a given binning
 534                  *
 535                  * @param axis: Axis to be defined
 536                  * @param name: Name of the axis
 537                  * @param title: Title of the axis
 538                  * @param binning: axis binning
 539                  * @param labels (@optional): array of bin labels
 540                  */
 541                 axis.Set(binning.GetSize()-1, binning.GetArray());
 542                 axis.SetName(name);
 543                 axis.SetTitle(title);
 544                 if(labels){
 545                         for(int ib = 1; ib <= axis.GetNbins(); ++ib)
 546                                 axis.SetBinLabel(ib, labels[ib-1]);
 547                 }
 548         }
 549
 550         //______________________________________________________________________________
 551         void AliAnalysisTaskPtEMCalTrigger::DefineAxis(TAxis& axis, const char* name,
 552                         const char* title, int nbins, double min, double max,
 553                         const char** labels) {
 554                 /*
 555                  * Define an axis with number of bins from min to max
 556                  *
 557                  * @param axis: Axis to be defined
 558                  * @param name: Name of the axis
 559                  * @param title: Title of the axis
 560                  * @param nbins: Number of bins
 561                  * @param min: lower limit of the axis
 562                  * @param max: upper limit of the axis
 563                  * @param labels (@optional): array of bin labels
 564                  */
 565                 axis.Set(nbins, min, max);
 566                 axis.SetName(name);
 567                 axis.SetTitle(title);
 568                 if(labels){
 569                         for(int ib = 1; ib <= axis.GetNbins(); ++ib)
 570                                 axis.SetBinLabel(ib, labels[ib-1]);
 571                 }
 572         }
 573
 574
 575         //______________________________________________________________________________
 576         void AliAnalysisTaskPtEMCalTrigger::FillEventHist(const char* trigger,
 577                         double vz, bool isPileup) {
 578                 /*
 579                  * Fill event-based histogram
 580                  *
 581                  * @param trigger: name of the trigger configuration to be processed
 582                  * @param vz: z-position of the vertex
 583                  * @param isPileup: signalises if the event is flagged as pileup event
 584                  */
 585                 char histname[1024];
 586                 sprintf(histname, "hEventHist%s", trigger);
 587                 try{
 588                         fHistos->FillTH2(histname, 0., vz);
 589                 } catch (HistoContainerContentException &e){
 590                         std::stringstream errormessage;
 591                         errormessage << "Filling of histogram failed: " << e.what();
 592                         AliError(errormessage.str().c_str());
 593                 }
 594                 if(!isPileup){
 595                         try{
 596                                 fHistos->FillTH2(histname, 1., vz);
 597                         } catch(HistoContainerContentException &e){
 598                                 std::stringstream errormessage;
 599                                 errormessage << "Filling of histogram failed: " << e.what();
 600                                 AliError(errormessage.str().c_str());
 601                         }
 602                 }
 603         }
 604
 605         //______________________________________________________________________________
 606         void AliAnalysisTaskPtEMCalTrigger::FillTrackHist(const char* trigger,
 607                         const AliVTrack* track, double vz, bool isPileup, int cut, bool isMinBias) {
 608                 /*
 609                  * Fill track-based histogram with corresponding information
 610                  *
 611                  * @param trigger: name of the trigger
 612                  * @param track: ESD track selected
 613                  * @param vz: z-position of the vertex
 614                  * @param isPileup: flag event as pileup event
 615                  * @param cut: id of the cut (0 = no cut)
 616                  */
 617                 double etasign = fSwapEta ? -1. : 1.;
 618         double data[7] = {TMath::Abs(track->Pt()), etasign * track->Eta(), track->Phi(), vz, 0, static_cast<double>(cut), isMinBias ? 1. : 0.};
 619         double dataMC[7] = {0., 0., 0., vz, 0, static_cast<double>(cut), isMinBias ? 1. : 0.};
 620         AliVParticle *assocMC(NULL);
 621         if(fMCEvent && (assocMC = fMCEvent->GetTrack(TMath::Abs(track->GetLabel())))){
 622                 dataMC[0] = TMath::Abs(assocMC->Pt());
 623                 dataMC[1] = etasign * assocMC->Eta();
 624                 dataMC[2] = assocMC->Phi();
 625         }
 626                 char histname[1024], histnameAcc[1024], histnameMC[1024], histnameMCAcc[1024];
 627                 sprintf(histname, "hTrackHist%s", trigger);
 628                 sprintf(histnameAcc, "hTrackInAcceptanceHist%s", trigger);
 629                 sprintf(histnameMC, "hMCTrackHist%s", trigger);
 630                 sprintf(histnameMCAcc, "hMCTrackInAcceptanceHist%s", trigger);
 631                 Bool_t isEMCAL = kFALSE;
 632                 if(track->IsEMCAL()){
 633                         // Check if the cluster is matched to only one track
 634                         AliVCluster *emcclust(NULL);
 635                         AliDebug(2, Form("cluster id: %d\n", track->GetEMCALcluster()));
 636                         if(fCaloClusters) {
 637                                 AliDebug(2, "Using calibrated clusters");
 638                                 emcclust = dynamic_cast<AliVCluster *>(fCaloClusters->At(track->GetEMCALcluster()));
 639                         } else {
 640                                 AliDebug(2, "Using uncalibrated clusters");
 641                                 emcclust = fInputEvent->GetCaloCluster(track->GetEMCALcluster());
 642                         }
 643                         if(!emcclust) AliError("Null pointer to EMCal cluster");
 644                         if(emcclust && emcclust->GetNTracksMatched() <= 1){
 645                                 isEMCAL = kTRUE;
 646                         }
 647                 }
 648                 try{
 649                         fHistos->FillTHnSparse(histname, data);
 650                         if(fMCEvent) fHistos->FillTHnSparse(histnameMC, dataMC);
 651                         if(isEMCAL){
 652                                 fHistos->FillTHnSparse(histnameAcc, data);
 653                                 if(fMCEvent) fHistos->FillTHnSparse(histnameMCAcc, dataMC);
 654                         }
 655                 } catch (HistoContainerContentException &e){
 656                         std::stringstream errormessage;
 657                         errormessage << "Filling of histogram failed: " << e.what();
 658                         AliError(errormessage.str().c_str());
 659                 }
 660                 if(!isPileup){
 661                         data[4] = 1;
 662                         dataMC[4] = 1;
 663                         try{
 664                                 fHistos->FillTHnSparse(histname, data);
 665                                 if(fMCEvent) fHistos->FillTHnSparse(histnameMC, dataMC);
 666                                 if(isEMCAL){
 667                                         fHistos->FillTHnSparse(histnameAcc, data);
 668                                         if(fMCEvent) fHistos->FillTHnSparse(histnameMCAcc, dataMC);
 669                                 }
 670                         } catch (HistoContainerContentException &e){
 671                                 std::stringstream errormessage;
 672                                 errormessage << "Filling of histogram failed: " << e.what();
 673                                 AliError(errormessage.str().c_str());
 674                         }
 675                 }
 676         }
 677
 678         //______________________________________________________________________________
 679         void AliAnalysisTaskPtEMCalTrigger::FillClusterHist(const char* trigger,
 680                         const AliVCluster* clust, bool isCalibrated, double vz, bool isPileup, bool isMinBias) {
 681                 /*
 682                  * Fill cluster-based histogram with corresponding information
 683                  *
 684                  * @param trigger: name of the trigger
 685                  * @param cluster: the EMCal cluster information
 686                  * @param vz: z-position of the vertex
 687                  * @param isPileup: flag event as pileup event
 688                  */
 689                 double data[4] =  {clust->E(), vz, 0, isMinBias ? 1. : 0.};
 690                 char histname[1024];
 691                 sprintf(histname, "hCluster%sHist%s", isCalibrated ? "Calib" : "Uncalib", trigger);
 692                 try{
 693                         fHistos->FillTHnSparse(histname, data);
 694                 } catch (HistoContainerContentException &e){
 695                         std::stringstream errormessage;
 696                         errormessage << "Filling of histogram failed: " << e.what();
 697                         AliError(errormessage.str().c_str());
 698                 }
 699                 if(!isPileup){
 700                         data[2] = 1.;
 701                         try{
 702                                 fHistos->FillTHnSparse(histname, data);
 703                         } catch (HistoContainerContentException &e){
 704                                 std::stringstream errormessage;
 705                                 errormessage << "Filling of histogram failed: " << e.what();
 706                                 AliError(errormessage.str().c_str());
 707                         }
 708                 }
 709         }
 710
 711         //______________________________________________________________________________
 712         void AliAnalysisTaskPtEMCalTrigger::FillMCParticleHist(const AliVParticle * const track, double vz, bool isPileup){
 713                 /*
 714                  * Fill histogram for MC-true particles with the information pt, eta and phi
 715                  *
 716                  * @param track: the Monte-Carlo track
 717                  */
 718                 double data[5] = {TMath::Abs(track->Pt()), track->Eta(), track->Phi(), vz, 0.};
 719                 fHistos->FillTHnSparse("hMCtrueParticles", data);
 720                 if(!isPileup){
 721                         data[4] = 1.;
 722                         fHistos->FillTHnSparse("hMCtrueParticles", data);
 723                 }
 724         }
 725
 726         //______________________________________________________________________________
 727         bool AliAnalysisTaskPtEMCalTrigger::IsTrueTrack(const AliVTrack *const track) const{
 728                 /*
 729                  * Check if the track has an associated MC particle, and that the particle is a physical primary
 730                  * In case of data we do not do the selection at that step (always return true)
 731                  *
 732                  * @param track: Track to check
 733                  * @result: true primary track (true or false)
 734                  */
 735                 if(!fMCEvent) return true;
 736                 AliVParticle *mcassociate = fMCEvent->GetTrack(TMath::Abs(track->GetLabel()));
 737                 if(!mcassociate) return false;
 738                 return fMCEvent->IsPhysicalPrimary(TMath::Abs(track->GetLabel()));
 739         }
 740
 741         //______________________________________________________________________________
 742         void AliAnalysisTaskPtEMCalTrigger::AddESDTrackCuts(AliESDtrackCuts* trackCuts) {
 743                 /*
 744                  * Add new track cuts to the task
 745                  *
 746                  * @param trackCuts: Object of type AliESDtrackCuts
 747                  */
 748                 fListTrackCuts->AddLast(new AliEMCalPtTaskTrackSelectionESD(trackCuts));
 749         }
 750
 751         //______________________________________________________________________________
 752         void AliAnalysisTaskPtEMCalTrigger::AddCutsForAOD(AliESDtrackCuts* trackCuts, UInt_t filterbits) {
 753                 /*
 754                  * Add new track cuts to the task
 755                  *
 756                  * @param trackCuts: Object of type AliESDtrackCuts
 757                  */
 758                 fListTrackCuts->AddLast(new AliEMCalPtTaskTrackSelectionAOD(trackCuts, filterbits));
 759         }
 760
 761
 762         //______________________________________________________________________________
 763         TString AliAnalysisTaskPtEMCalTrigger::BuildTriggerString() {
 764                 /*
 765                  * Build trigger string from the trigger maker
 766                  *
 767                  * @return: blank-separated string of fired trigger classes
 768                  */
 769                 AliDebug(1, "trigger checking");
 770                 TString result = "";
 771                 if(HasTriggerType(kJ1)) result += "EJ1 ";
 772                 if(HasTriggerType(kJ2)) result += "EJ2 ";
 773                 if(HasTriggerType(kG1)) result += "EG1 ";
 774                 if(HasTriggerType(kG2)) result += "EG2 ";
 775                 return result;
 776         }
 777
 778         //______________________________________________________________________________
 779         const AliVVertex* AliAnalysisTaskPtEMCalTrigger::GetSPDVertex() const {
 780                 /*
 781                  * Accessor for the SPD vertex, creating transparency for ESDs and AODs
 782                  *
 783                  * @return: the spd vertex
 784                  */
 785                 AliESDEvent *esd = dynamic_cast<AliESDEvent *>(fInputEvent);
 786                 if(esd){
 787                         return esd->GetPrimaryVertexSPD();
 788                 } else {
 789                         AliAODEvent *aod = dynamic_cast<AliAODEvent *>(fInputEvent);
 790                         if(aod){
 791                                 return aod->GetPrimaryVertexSPD();
 792                         }
 793                 }
 794                 return NULL;
 795         }
 796
 797 }
 798