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