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 <TDirectory.h>
  31 #include <TH1.h>
  32 #include <THashList.h>
  33 #include <TList.h>
  34 #include <TKey.h>
  35 #include <TMath.h>
  36 #include <TObjArray.h>
  37 #include <TString.h>
  38
  39 #include "AliESDEvent.h"
  40 #include "AliESDInputHandler.h"
  41 #include "AliESDtrack.h"
  42 #include "AliESDVertex.h"
  43
  44 #include "AliEMCalHistoContainer.h"
  45 #include "AliAnalysisTaskPtEMCalTrigger.h"
  46
  47 ClassImp(EMCalTriggerPtAnalysis::AliAnalysisTaskPtEMCalTrigger)
  48
  49 namespace EMCalTriggerPtAnalysis {
  50
  51         //______________________________________________________________________________
  52         AliAnalysisTaskPtEMCalTrigger::AliAnalysisTaskPtEMCalTrigger():
  53                                 AliAnalysisTaskSE(),
  54                                 fResults(NULL),
  55                                 fHistos(NULL),
  56                                 fListTrackCuts(NULL)
  57         {
  58                 /*
  59                  * Dummy constructor, initialising the values with default (NULL) values
  60                  */
  61         }
  62
  63         //______________________________________________________________________________
  64         AliAnalysisTaskPtEMCalTrigger::AliAnalysisTaskPtEMCalTrigger(const char *name):
  65                                 AliAnalysisTaskSE(name),
  66                                 fResults(NULL),
  67                                 fHistos(NULL),
  68                                 fListTrackCuts(NULL)
  69         {
  70                 /*
  71                  * Main constructor, setting default values for eta and zvertex cut
  72                  */
  73                 DefineOutput(1, TList::Class());
  74
  75                 fListTrackCuts = new TList;
  76                 fListTrackCuts->SetOwner(false);
  77
  78                 // Set default cuts
  79                 fEtaRange.SetLimits(-0.8, 0.8);
  80
  81         }
  82
  83         //______________________________________________________________________________
  84         AliAnalysisTaskPtEMCalTrigger::~AliAnalysisTaskPtEMCalTrigger(){
  85                 /*
  86                  * Destructor, deleting output
  87                  */
  88                 //if(fTrackSelection) delete fTrackSelection;
  89                 if(fHistos) delete fHistos;
  90                 if(fListTrackCuts) delete fListTrackCuts;
  91         }
  92
  93         //______________________________________________________________________________
  94         void AliAnalysisTaskPtEMCalTrigger::UserCreateOutputObjects(){
  95                 /*
  96                  * Create the list of output objects and define the histograms.
  97                  * Also adding the track cuts to the list of histograms.
  98                  */
  99                 fResults = new TList;
 100                 fResults->SetOwner();
 101
 102                 fHistos = new AliEMCalHistoContainer("PtEMCalTriggerHistograms");
 103                 fHistos->ReleaseOwner();
 104
 105                 std::map<std::string, std::string> triggerCombinations;
 106                 const char *triggernames[6] = {"MinBias", "EMCJHigh", "EMCJLow", "EMCGHigh", "EMCGLow", "NoEMCal"};
 107                 // Define axes for the trigger correlation histogram
 108                 const TAxis *triggeraxis[5]; memset(triggeraxis, 0, sizeof(const TAxis *) * 5);
 109                 const char *binlabels[2] = {"OFF", "ON"};
 110                 TAxis mytrgaxis[5];
 111                 for(int itrg = 0; itrg < 5; ++itrg){
 112                         DefineAxis(mytrgaxis[itrg], triggernames[itrg], triggernames[itrg], 2, -0.5, 1.5, binlabels);
 113                         triggeraxis[itrg] = mytrgaxis+itrg;
 114                 }
 115                 // Define names and titles for different triggers in the histogram container
 116                 triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[0], "min. bias events"));
 117                 triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[1], "jet-triggered events (high threshold)"));
 118                 triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[2], "jet-triggered events (low threshold)"));
 119                 triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[3], "jet-triggered events (high threshold)"));
 120                 triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[4], "jet-triggered events (low threshold)"));
 121                 triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[5], "non-EMCal-triggered events (low threshold)"));
 122                 // Define axes for the pt histogram
 123                 // Dimensions:
 124                 // 1. pt
 125                 // 2. eta
 126                 // 3. phi
 127                 // 4. vertex
 128                 // 5. pileup (0 = all events, 1 = after pileup rejection)
 129                 // 6. track cuts (0 = no cuts; 1 = after std cuts)
 130                 TArrayD ptbinning, zvertexBinning, etabinning, pileupaxis(3);
 131                 pileupaxis[0] = -0.5; pileupaxis[1] = 0.5; pileupaxis[2] = 1.5;
 132                 CreateDefaultPtBinning(ptbinning);
 133                 CreateDefaultZVertexBinning(zvertexBinning);
 134                 CreateDefaultEtaBinning(etabinning);
 135                 TAxis htrackaxes[6];
 136                 DefineAxis(htrackaxes[0], "pt", "p_{t} (GeV/c)", ptbinning);
 137                 DefineAxis(htrackaxes[1], "eta", "#eta", etabinning);
 138                 DefineAxis(htrackaxes[2], "phi", "#phi", 20, 0, 2 * TMath::Pi());
 139                 DefineAxis(htrackaxes[3], "zvertex", "z_{V} (cm)", zvertexBinning);
 140                 DefineAxis(htrackaxes[4], "pileup", "Pileup rejection", 2, -0.5, 1.5);
 141                 DefineAxis(htrackaxes[5], "trackcuts", "Track Cuts", (fListTrackCuts ? fListTrackCuts->GetEntries() : 0) + 1, -0.5, (fListTrackCuts ? fListTrackCuts->GetEntries() : 0) + 0.5);
 142                 const TAxis *trackaxes[6];
 143                 for(int iaxis = 0; iaxis < 6; ++iaxis) trackaxes[iaxis] = htrackaxes + iaxis;
 144                 try{
 145                         for(std::map<std::string,std::string>::iterator it = triggerCombinations.begin(); it != triggerCombinations.end(); ++it){
 146                                 const std::string name = it->first, &title = it->second;
 147                                 // Create event-based histogram
 148                                 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);
 149                                 // Create track-based histogram
 150                                 fHistos->CreateTHnSparse(Form("hTrackHist%s", name.c_str()), Form("Track-based data for %s events", title.c_str()), 6, trackaxes);
 151                         }
 152                         fHistos->CreateTHnSparse("hEventTriggers", "Trigger type per event", 5, triggeraxis);
 153                 } catch (HistoContainerContentException &e){
 154                         std::stringstream errormessage;
 155                         errormessage << "Creation of histogram failed: " << e.what();
 156                         AliError(errormessage.str().c_str());
 157                 }
 158                 fResults->Add(fHistos->GetListOfHistograms());
 159                 if(fListTrackCuts && fListTrackCuts->GetEntries()){
 160                         TIter cutIter(fListTrackCuts);
 161                         AliESDtrackCuts *cutObject(NULL);
 162                         while((cutObject = dynamic_cast<AliESDtrackCuts *>(cutIter()))){
 163                                 cutObject->DefineHistograms();
 164                                 fResults->Add(cutObject);
 165                         }
 166                 }
 167                 PostData(1, fResults);
 168         }
 169
 170         //______________________________________________________________________________
 171         void AliAnalysisTaskPtEMCalTrigger::UserExec(Option_t* /*option*/){
 172                 /*
 173                  * Runs the event loop
 174                  *
 175                  * @param option: Additional options
 176                  */
 177
 178                 // Common checks: Have SPD vertex and primary vertex from tracks, and both need to have at least one contributor
 179                 AliESDEvent *esd = static_cast<AliESDEvent *>(fInputEvent);
 180                 const AliESDVertex *vtxTracks = esd->GetPrimaryVertex(),
 181                                 *vtxSPD = esd->GetPrimaryVertexSPD();
 182                 if(!(vtxTracks && vtxSPD)) return;
 183                 if(vtxTracks->GetNContributors() < 1 || vtxSPD->GetNContributors() < 1) return;
 184
 185                 double triggers[5]; memset(triggers, 0, sizeof(double) *5);
 186                 if(fInputHandler->IsEventSelected() & AliVEvent::kINT7)
 187                         triggers[0] = 1.;
 188
 189                 std::vector<std::string> triggerstrings;
 190                 if(fInputHandler->IsEventSelected() & AliVEvent::kEMC7){
 191                         // EMCal-triggered event, distinguish types
 192                         TString trgstr(fInputEvent->GetFiredTriggerClasses());
 193                         if(trgstr.Contains("EJ1")){
 194                                 triggerstrings.push_back("EMCJHigh");
 195                                 triggers[1] = 1;
 196                         }
 197                         if(trgstr.Contains("EJ2")){
 198                                 triggerstrings.push_back("EMCJLow");
 199                                 triggers[2] = 1;
 200                         }
 201                         if(trgstr.Contains("EG1")){
 202                                 triggerstrings.push_back("EMCGHigh");
 203                                 triggers[3] = 1;
 204                         }
 205                         if(trgstr.Contains("EG2")){
 206                                 triggerstrings.push_back("EMCGLow");
 207                                 triggers[4] = 1;
 208                         }
 209                 }
 210                 try{
 211                         fHistos->FillTHnSparse("hEventTriggers", triggers);
 212                 } catch (HistoContainerContentException &e){
 213                         std::stringstream errormessage;
 214                         errormessage << "Filling of histogram failed: " << e.what();
 215                         AliError(errormessage.str().c_str());
 216                 }
 217
 218                 // apply event selection: Combine the Pileup cut from SPD with the other pA Vertex selection cuts.
 219                 bool isPileupEvent = esd->IsPileupFromSPD();
 220                 isPileupEvent = isPileupEvent || (TMath::Abs(vtxTracks->GetZ() - vtxSPD->GetZ()) > 0.5);
 221                 double covSPD[6]; vtxSPD->GetCovarianceMatrix(covSPD);
 222                 isPileupEvent = isPileupEvent || (TString(vtxSPD->GetTitle()).Contains("vertexer:Z") && TMath::Sqrt(covSPD[5]) > 0.25);
 223
 224                 // Fill event-based histogram
 225                 const double &zv = vtxTracks->GetZ();
 226                 if(triggers[0]) FillEventHist("MinBias", zv, isPileupEvent);
 227                 if(!triggerstrings.size()) // Non-EMCal-triggered
 228                         FillEventHist("NoEMCal", zv, isPileupEvent);
 229                 else{
 230                         // EMCal-triggered events
 231                         for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
 232                                 FillEventHist(it->c_str(), zv, isPileupEvent);
 233                 }
 234
 235                 AliESDtrack *track(NULL);
 236                 // Loop over all tracks (No cuts applied)
 237                 for(int itrk = 0; itrk < fInputEvent->GetNumberOfTracks(); ++itrk){
 238                         track = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(itrk));
 239                         // first fill without pielup cut
 240                         if(fEtaRange.IsInRange(track->Eta())) continue;
 241                         if(triggers[0]) FillTrackHist("MinBias", track, zv, isPileupEvent, 0);
 242                         if(!triggerstrings.size()) // Non-EMCal-triggered
 243                                 FillTrackHist("NoEMCal", track, zv, isPileupEvent, 0);
 244                         else {
 245                                 // EMCal-triggered events
 246                                 for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
 247                                         FillTrackHist(it->c_str(), track, zv, isPileupEvent, 0);
 248                         }
 249                 }
 250
 251                 // Now apply track selection cuts
 252                 // allow for several track selections to be tested at the same time
 253                 // each track selection gets a different cut ID starting from 1
 254                 // cut ID 0 is reserved for the case of no cuts
 255                 if(fListTrackCuts && fListTrackCuts->GetEntries()){
 256                         for(int icut = 0; icut < fListTrackCuts->GetEntries(); icut++){
 257                                 AliESDtrackCuts *trackSelection = static_cast<AliESDtrackCuts *>(fListTrackCuts->At(icut));
 258                                 std::auto_ptr<TObjArray> acceptedTracks(trackSelection->GetAcceptedTracks(esd));
 259                                 TIter trackIter(acceptedTracks.get());
 260                                 while((track = dynamic_cast<AliESDtrack *>(trackIter()))){
 261                                         if(!fEtaRange.IsInRange(track->Eta())) continue;
 262                                         if(triggers[0]) FillTrackHist("MinBias", track, zv, isPileupEvent, icut + 1);
 263                                         if(!triggerstrings.size()) // Non-EMCal-triggered
 264                                                 FillTrackHist("NoEMCal", track, zv, isPileupEvent, icut + 1);
 265                                         else {
 266                                                 // EMCal-triggered events
 267                                                 for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
 268                                                         FillTrackHist(it->c_str(), track, zv, isPileupEvent, icut + 1);
 269                                         }
 270                                 }
 271                         }
 272                 }
 273                 PostData(1, fResults);
 274         }
 275
 276         //______________________________________________________________________________
 277         void AliAnalysisTaskPtEMCalTrigger::CreateDefaultPtBinning(TArrayD &binning) const{
 278                 /*
 279                  * Creating the default pt binning.
 280                  *
 281                  * @param binning: Array where to store the results.
 282                  */
 283                 std::vector<double> mybinning;
 284                 std::map<double,double> definitions;
 285                 definitions.insert(std::pair<double,double>(2.5, 0.1));
 286                 definitions.insert(std::pair<double,double>(7., 0.25));
 287                 definitions.insert(std::pair<double,double>(15., 0.5));
 288                 definitions.insert(std::pair<double,double>(25., 1.));
 289                 definitions.insert(std::pair<double,double>(40., 2.5));
 290                 definitions.insert(std::pair<double,double>(60., 5.));
 291                 definitions.insert(std::pair<double,double>(100., 5.));
 292                 double currentval = 0;
 293                 for(std::map<double,double>::iterator id = definitions.begin(); id != definitions.end(); ++id){
 294                         double limit = id->first, binwidth = id->second;
 295                         while(currentval <= limit){
 296                                 currentval += binwidth;
 297                                 mybinning.push_back(currentval);
 298                         }
 299                 }
 300                 binning.Set(mybinning.size());
 301                 int ib = 0;
 302                 for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
 303                         binning[ib++] = *it;
 304         }
 305
 306         //______________________________________________________________________________
 307         void AliAnalysisTaskPtEMCalTrigger::CreateDefaultZVertexBinning(TArrayD &binning) const {
 308                 /*
 309                  * Creating default z-Vertex binning.
 310                  *
 311                  * @param binning: Array where to store the results.
 312                  */
 313                 std::vector<double> mybinning;
 314                 double currentval = -40;
 315                 mybinning.push_back(currentval);
 316                 while(currentval <= 40.){
 317                         currentval += 1.;
 318                         mybinning.push_back(currentval);
 319                 }
 320                 binning.Set(mybinning.size());
 321                 int ib = 0;
 322                 for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
 323                         binning[ib++] = *it;
 324         }
 325
 326         //______________________________________________________________________________
 327         void AliAnalysisTaskPtEMCalTrigger::CreateDefaultEtaBinning(TArrayD& binning) const {
 328                 /*
 329                  * Creating default z-Vertex binning.
 330                  *
 331                  * @param binning: Array where to store the results.
 332                  */
 333                 std::vector<double> mybinning;
 334                 double currentval = -0.8;
 335                 mybinning.push_back(currentval);
 336                 while(currentval <= 0.8){
 337                         currentval += 0.1;
 338                         mybinning.push_back(currentval);
 339                 }
 340                 binning.Set(mybinning.size());
 341                 int ib = 0;
 342                 for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
 343                         binning[ib++] = *it;
 344         }
 345
 346         //______________________________________________________________________________
 347         void AliAnalysisTaskPtEMCalTrigger::DefineAxis(TAxis& axis, const char* name,
 348                         const char* title, const TArrayD& binning, const char** labels) {
 349                 /*
 350                  * Define an axis with a given binning
 351                  *
 352                  * @param axis: Axis to be defined
 353                  * @param name: Name of the axis
 354                  * @param title: Title of the axis
 355                  * @param binning: axis binning
 356                  * @param labels (@optional): array of bin labels
 357                  */
 358                 axis.Set(binning.GetSize()-1, binning.GetArray());
 359                 axis.SetName(name);
 360                 axis.SetTitle(title);
 361                 if(labels){
 362                         for(int ib = 1; ib <= axis.GetNbins(); ++ib)
 363                                 axis.SetBinLabel(ib, labels[ib]);
 364                 }
 365         }
 366
 367         //______________________________________________________________________________
 368         void AliAnalysisTaskPtEMCalTrigger::DefineAxis(TAxis& axis, const char* name,
 369                         const char* title, int nbins, double min, double max,
 370                         const char** labels) {
 371                 /*
 372                  * Define an axis with number of bins from min to max
 373                  *
 374                  * @param axis: Axis to be defined
 375                  * @param name: Name of the axis
 376                  * @param title: Title of the axis
 377                  * @param nbins: Number of bins
 378                  * @param min: lower limit of the axis
 379                  * @param max: upper limit of the axis
 380                  * @param labels (@optional): array of bin labels
 381                  */
 382                 axis.Set(nbins, min, max);
 383                 axis.SetName(name);
 384                 axis.SetTitle(title);
 385                 if(labels){
 386                         for(int ib = 1; ib <= axis.GetNbins(); ++ib)
 387                                 axis.SetBinLabel(ib, labels[ib]);
 388                 }
 389         }
 390
 391         //______________________________________________________________________________
 392         void AliAnalysisTaskPtEMCalTrigger::FillEventHist(const char* trigger,
 393                         double vz, bool isPileup) {
 394                 /*
 395                  * Fill event-based histogram
 396                  *
 397                  * @param trigger: name of the trigger configuration to be processed
 398                  * @param vz: z-position of the vertex
 399                  * @param isPileup: signalises if the event is flagged as pileup event
 400                  */
 401                 char histname[1024];
 402                 sprintf(histname, "hEventHist%s", trigger);
 403                 try{
 404                         fHistos->FillTH2(histname, vz, 0);
 405                 } catch (HistoContainerContentException &e){
 406                         std::stringstream errormessage;
 407                         errormessage << "Filling of histogram failed: " << e.what();
 408                         AliError(errormessage.str().c_str());
 409                 }
 410                 if(!isPileup){
 411                         try{
 412                                 fHistos->FillTH2(histname, vz, 1);
 413                         } catch(HistoContainerContentException &e){
 414                                 std::stringstream errormessage;
 415                                 errormessage << "Filling of histogram failed: " << e.what();
 416                                 AliError(errormessage.str().c_str());
 417                         }
 418                 }
 419         }
 420
 421         //______________________________________________________________________________
 422         void AliAnalysisTaskPtEMCalTrigger::FillTrackHist(const char* trigger,
 423                         const AliESDtrack* track, double vz, bool isPileup, int cut) {
 424                 /*
 425                  * Fill track-based histogram with corresponding information
 426                  *
 427                  * @param trigger: name of the trigger
 428                  * @param track: ESD track selected
 429                  * @param vz: z-position of the vertex
 430                  * @param isPileup: flag event as pileup event
 431                  * @param cut: id of the cut (0 = no cut)
 432                  */
 433                 double data[6] = {track->Pt(), track->Eta(), track->Phi(), vz, 0, cut};
 434                 char histname[1024];
 435                 sprintf(histname, "hTrackHist%s", trigger);
 436                 try{
 437                         fHistos->FillTHnSparse(histname, data);
 438                 } catch (HistoContainerContentException &e){
 439                         std::stringstream errormessage;
 440                         errormessage << "Filling of histogram failed: " << e.what();
 441                         AliError(errormessage.str().c_str());
 442                 }
 443                 if(!isPileup){
 444                         data[4] = 1;
 445                         try{
 446                                 fHistos->FillTHnSparse(histname, data);
 447                         } catch (HistoContainerContentException &e){
 448                                 std::stringstream errormessage;
 449                                 errormessage << "Filling of histogram failed: " << e.what();
 450                                 AliError(errormessage.str().c_str());
 451                         }
 452                 }
 453         }
 454
 455 }
 456