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