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