AliPhysics  2c6b7ad (2c6b7ad)
AliAnalysisTaskClusterQA.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Authors: Nicolas Schmidt *
5  * Version 1.0 *
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 
17 //---------------------------------------------
18 // Class used to create a tree for QA of clusters
19 //---------------------------------------------
21 
23 #include "TChain.h"
24 #include "TRandom.h"
25 #include "AliAnalysisManager.h"
26 #include "TParticle.h"
27 #include "TVectorF.h"
28 #include "AliPIDResponse.h"
29 #include "TFile.h"
30 #include "AliESDtrackCuts.h"
31 #include "AliAODMCParticle.h"
32 #include "AliAODMCHeader.h"
33 #include "AliAODEvent.h"
34 #include "AliMultSelection.h"
35 
36 class iostream;
37 
38 using namespace std;
39 
41 
42 //________________________________________________________________________
44  fV0Reader(NULL),
45  fV0ReaderName("V0ReaderV1"),
46  fCorrTaskSetting(""),
47  fConversionCuts(NULL),
48  fEventCuts(NULL),
49  fClusterCutsEMC(NULL),
50  fClusterCutsDMC(NULL),
51  fMesonCuts(NULL),
52  fMinNLMCut(1),
53  fMaxNLMCut(1),
54  fInputEvent(NULL),
55  fMCEvent(NULL),
56  fWeightJetJetMC(1),
57  fGeomEMCAL(NULL),
58  fClusterTree(NULL),
59  fIsHeavyIon(kFALSE),
60  ffillTree(-100),
61  ffillHistograms(kFALSE),
62  fOutputList(NULL),
63  fIsMC(kFALSE),
64  fCorrectForNonlinearity(kFALSE),
65  fSaveEventProperties(0),
66  fSaveCells(0),
67  fSaveSurroundingCells(0),
68  fSaveTracks(0),
69  fConeRadius(0),
70  fMinTrackPt(0),
71  fMinClusterEnergy(0),
72  fSaveMCInformation(0),
73  fSaveAdditionalHistos(0),
74  fExtractionPercentages(),
75  fExtractionPercentagePtBins(),
76  fBuffer_EventID(0),
77  fBuffer_ClusterE(0),
78  fBuffer_ClusterPhi(0),
79  fBuffer_ClusterEta(0),
80  fBuffer_ClusterIsEMCAL(kFALSE),
81  fBuffer_MC_Cluster_Flag(0),
82  fBuffer_ClusterNumCells(0),
83  fBuffer_LeadingCell_ID(0),
84  fBuffer_LeadingCell_E(0),
85  fBuffer_LeadingCell_Eta(0),
86  fBuffer_LeadingCell_Phi(0),
87  fBuffer_ClusterM02(0),
88  fBuffer_ClusterM20(0),
89  fBuffer_Event_Vertex_X(0),
90  fBuffer_Event_Vertex_Y(0),
91  fBuffer_Event_Vertex_Z(0),
92  fBuffer_Event_Multiplicity(0),
93  fBuffer_Event_NumActiveCells(0),
94  fBuffer_Cells_ID(0),
95  fBuffer_Cells_E(0),
96  fBuffer_Cells_RelativeEta(0),
97  fBuffer_Cells_RelativePhi(0),
98  fBuffer_Surrounding_NCells(0),
99  fBuffer_Surrounding_Cells_ID(0),
100  fBuffer_Surrounding_Cells_R(0),
101  fBuffer_Surrounding_Cells_E(0),
102  fBuffer_Surrounding_Cells_RelativeEta(0),
103  fBuffer_Surrounding_Cells_RelativePhi(0),
104  fBuffer_Surrounding_NTracks(10),
105  fBuffer_Surrounding_Tracks_R(0),
106  fBuffer_Surrounding_Tracks_Pt(0),
107  fBuffer_Surrounding_Tracks_RelativeEta(0),
108  fBuffer_Surrounding_Tracks_RelativePhi(0),
109  fBuffer_Cluster_MC_Label(-10),
110  fBuffer_Mother_MC_Label(-10),
111  hNCellsInClustersVsCentrality(NULL),
112  hNActiveCellsVsCentrality(NULL),
113  hNActiveCellsAbove50MeVVsCentrality(NULL),
114  hNActiveCellsAbove80MeVVsCentrality(NULL),
115  hNActiveCellsAbove100MeVVsCentrality(NULL),
116  hNActiveCellsAbove150MeVVsCentrality(NULL),
117  hECellsInClustersVsCentrality(NULL),
118  hEActiveCellsVsCentrality(NULL),
119  hEActiveCells50MeVVsCentrality(NULL),
120  hEActiveCells80MeVVsCentrality(NULL),
121  hEActiveCells100MeVVsCentrality(NULL),
122  hEActiveCells150MeVVsCentrality(NULL)
123 {
124  fBuffer_Cells_ID = new Int_t[kMaxActiveCells];
125  fBuffer_Cells_E = new Float_t[kMaxActiveCells];
126  fBuffer_Cells_RelativeEta = new Float_t[kMaxActiveCells];
127  fBuffer_Cells_RelativePhi = new Float_t[kMaxActiveCells];
128  fBuffer_Surrounding_Cells_ID = new Int_t[kMaxActiveCells];
129  fBuffer_Surrounding_Cells_R = new Float_t[kMaxActiveCells];
130  fBuffer_Surrounding_Cells_E = new Float_t[kMaxActiveCells];
131  fBuffer_Surrounding_Cells_RelativeEta = new Float_t[kMaxActiveCells];
132  fBuffer_Surrounding_Cells_RelativePhi = new Float_t[kMaxActiveCells];
133  fBuffer_Surrounding_Tracks_R = new Float_t[kMaxNTracks];
134  fBuffer_Surrounding_Tracks_Pt = new Float_t[kMaxNTracks];
135  fBuffer_Surrounding_Tracks_RelativeEta= new Float_t[kMaxNTracks];
136  fBuffer_Surrounding_Tracks_RelativePhi= new Float_t[kMaxNTracks];
137 }
138 
140  fV0Reader(NULL),
141  fV0ReaderName("V0ReaderV1"),
142  fCorrTaskSetting(""),
143  fConversionCuts(NULL),
144  fEventCuts(NULL),
145  fClusterCutsEMC(NULL),
146  fClusterCutsDMC(NULL),
147  fMesonCuts(NULL),
148  fMinNLMCut(1),
149  fMaxNLMCut(1),
150  fInputEvent(NULL),
151  fMCEvent(NULL),
152  fWeightJetJetMC(1),
153  fGeomEMCAL(NULL),
154  fClusterTree(NULL),
155  fIsHeavyIon(kFALSE),
156  ffillTree(-100),
157  ffillHistograms(kFALSE),
158  fOutputList(NULL),
159  fIsMC(kFALSE),
160  fCorrectForNonlinearity(kFALSE),
161  fSaveEventProperties(0),
162  fSaveCells(0),
163  fSaveSurroundingCells(0),
164  fSaveTracks(0),
165  fConeRadius(0),
166  fMinTrackPt(0),
167  fMinClusterEnergy(0),
168  fSaveMCInformation(0),
169  fSaveAdditionalHistos(0),
170  fExtractionPercentages(),
171  fExtractionPercentagePtBins(),
172  fBuffer_EventID(0),
173  fBuffer_ClusterE(0),
174  fBuffer_ClusterPhi(0),
175  fBuffer_ClusterEta(0),
176  fBuffer_ClusterIsEMCAL(kFALSE),
177  fBuffer_MC_Cluster_Flag(0),
178  fBuffer_ClusterNumCells(0),
179  fBuffer_LeadingCell_ID(0),
180  fBuffer_LeadingCell_E(0),
181  fBuffer_LeadingCell_Eta(0),
182  fBuffer_LeadingCell_Phi(0),
183  fBuffer_ClusterM02(0),
184  fBuffer_ClusterM20(0),
185  fBuffer_Event_Vertex_X(0),
186  fBuffer_Event_Vertex_Y(0),
187  fBuffer_Event_Vertex_Z(0),
188  fBuffer_Event_Multiplicity(0),
189  fBuffer_Event_NumActiveCells(0),
190  fBuffer_Cells_ID(0),
191  fBuffer_Cells_E(0),
192  fBuffer_Cells_RelativeEta(0),
193  fBuffer_Cells_RelativePhi(0),
194  fBuffer_Surrounding_NCells(0),
195  fBuffer_Surrounding_Cells_ID(0),
196  fBuffer_Surrounding_Cells_R(0),
197  fBuffer_Surrounding_Cells_E(0),
198  fBuffer_Surrounding_Cells_RelativeEta(0),
199  fBuffer_Surrounding_Cells_RelativePhi(0),
200  fBuffer_Surrounding_NTracks(10),
201  fBuffer_Surrounding_Tracks_R(0),
202  fBuffer_Surrounding_Tracks_Pt(0),
203  fBuffer_Surrounding_Tracks_RelativeEta(0),
204  fBuffer_Surrounding_Tracks_RelativePhi(0),
205  fBuffer_Cluster_MC_Label(-10),
206  fBuffer_Mother_MC_Label(-10),
207  hNCellsInClustersVsCentrality(NULL),
208  hNActiveCellsVsCentrality(NULL),
209  hNActiveCellsAbove50MeVVsCentrality(NULL),
210  hNActiveCellsAbove80MeVVsCentrality(NULL),
211  hNActiveCellsAbove100MeVVsCentrality(NULL),
212  hNActiveCellsAbove150MeVVsCentrality(NULL),
213  hECellsInClustersVsCentrality(NULL),
214  hEActiveCellsVsCentrality(NULL),
215  hEActiveCells50MeVVsCentrality(NULL),
216  hEActiveCells80MeVVsCentrality(NULL),
217  hEActiveCells100MeVVsCentrality(NULL),
218  hEActiveCells150MeVVsCentrality(NULL)
219 {
233  // Default constructor
234 
235  DefineInput(0, TChain::Class());
236  DefineOutput(1, TList::Class());
237  DefineOutput(2, TTree::Class());
238 }
239 
240 //________________________________________________________________________
242 {
243  // default deconstructor
244 
245 }
246 //________________________________________________________________________
248 {
249  // Create User Output Objects
250 
251  if(fOutputList != NULL){
252  delete fOutputList;
253  fOutputList = NULL;
254  }
255  if(fOutputList == NULL){
256  fOutputList = new TList();
257  fOutputList->SetOwner(kTRUE);
258  }
259 
260  if(ffillHistograms){
261 
262  if(((AliConvEventCuts*)fEventCuts)->GetCutHistograms()){
263  fOutputList->Add(((AliConvEventCuts*)fEventCuts)->GetCutHistograms());
264  }
265 
266  if(((AliCaloPhotonCuts*)fClusterCutsEMC)->GetCutHistograms()){
267  fOutputList->Add(((AliCaloPhotonCuts*)fClusterCutsEMC)->GetCutHistograms());
268  }
269  if(((AliCaloPhotonCuts*)fClusterCutsDMC)->GetCutHistograms()){
270  fOutputList->Add(((AliCaloPhotonCuts*)fClusterCutsDMC)->GetCutHistograms());
271  }
272  if(((AliConversionMesonCuts*)fMesonCuts)->GetCutHistograms()){
273  fOutputList->Add(((AliConversionMesonCuts*)fMesonCuts)->GetCutHistograms());
274  }
276  hNCellsInClustersVsCentrality = new TH2F("NCellsInClusters_Centrality","NCellsInClusters_Centrality",18000, 0, 18000, 100, 0, 100);
278  hNActiveCellsVsCentrality = new TH2F("NActiveCells_Centrality","NActiveCells_Centrality",18000, 0, 18000, 100, 0, 100);
280  hNActiveCellsAbove50MeVVsCentrality = new TH2F("NActiveCells50MeV_Centrality","NActiveCells50MeV_Centrality",18000, 0, 18000, 100, 0, 100);
282  hNActiveCellsAbove80MeVVsCentrality = new TH2F("NActiveCells80MeV_Centrality","NActiveCells80MeV_Centrality",18000, 0, 18000, 100, 0, 100);
284  hNActiveCellsAbove100MeVVsCentrality = new TH2F("NActiveCells100MeV_Centrality","NActiveCells100MeV_Centrality",18000, 0, 18000, 100, 0, 100);
286  hNActiveCellsAbove150MeVVsCentrality = new TH2F("NActiveCells150MeV_Centrality","NActiveCells150MeV_Centrality",18000, 0, 18000, 100, 0, 100);
288 
289  hECellsInClustersVsCentrality = new TH2F("EofCellsInClusters_Centrality","EofCellsInClusters_Centrality",2000, 0, 20, 100, 0, 100);
291  hEActiveCellsVsCentrality = new TH2F("EofActiveCells_Centrality","EofActiveCells_Centrality",2000, 0, 20, 100, 0, 100);
293  hEActiveCells50MeVVsCentrality = new TH2F("EofActiveCells50MeV_Centrality","EofActiveCells50MeV_Centrality",2000, 0, 20, 100, 0, 100);
295  hEActiveCells80MeVVsCentrality = new TH2F("EofActiveCells80MeV_Centrality","EofActiveCells80MeV_Centrality",2000, 0, 20, 100, 0, 100);
297  hEActiveCells100MeVVsCentrality = new TH2F("EofActiveCells100MeV_Centrality","EofActiveCells100MeV_Centrality",2000, 0, 20, 100, 0, 100);
299  hEActiveCells150MeVVsCentrality = new TH2F("EofActiveCells150MeV_Centrality","EofActiveCells150MeV_Centrality",2000, 0, 20, 100, 0, 100);
301  }
302  PostData(1, fOutputList);
303  }
304 
305  fClusterTree = new TTree(Form("ClusterQA_%s_%s_%s",(fEventCuts->GetCutNumber()).Data(),(fClusterCutsEMC->GetCutNumber()).Data(),(fClusterCutsDMC->GetCutNumber()).Data()),Form("ClusterQA_%s_%s_%s",(fEventCuts->GetCutNumber()).Data(),(fClusterCutsEMC->GetCutNumber()).Data(),(fClusterCutsDMC->GetCutNumber()).Data()));
306 
307  fClusterTree->Branch("Cluster_E", &fBuffer_ClusterE, "Cluster_E/F");
308  fClusterTree->Branch("Cluster_Eta", &fBuffer_ClusterEta, "Cluster_Eta/F");
309  fClusterTree->Branch("Cluster_Phi", &fBuffer_ClusterPhi, "Cluster_Phi/F");
310  fClusterTree->Branch("Cluster_IsEMCAL", &fBuffer_ClusterIsEMCAL, "Cluster_IsEMCAL/O");
311  fClusterTree->Branch("Cluster_NumCells", &fBuffer_ClusterNumCells, "Cluster_NumCells/I");
312  fClusterTree->Branch("Cluster_LeadingCell_ID", &fBuffer_LeadingCell_ID, "Cluster_LeadingCell_ID/I");
313  fClusterTree->Branch("Cluster_LeadingCell_E", &fBuffer_LeadingCell_E, "Cluster_LeadingCell_E/F");
314  fClusterTree->Branch("Cluster_LeadingCell_Eta", &fBuffer_LeadingCell_Eta, "Cluster_LeadingCell_Eta/F");
315  fClusterTree->Branch("Cluster_LeadingCell_Phi", &fBuffer_LeadingCell_Phi, "Cluster_LeadingCell_Phi/F");
316  fClusterTree->Branch("Cluster_M02", &fBuffer_ClusterM02, "Cluster_M02/F");
317  fClusterTree->Branch("Cluster_M20", &fBuffer_ClusterM20, "Cluster_M20/F");
319  {
320  fClusterTree->Branch("Event_ID", &fBuffer_EventID, "Event_ID/l");
321  fClusterTree->Branch("Event_Vertex_X", &fBuffer_Event_Vertex_X, "Event_Vertex_X/F");
322  fClusterTree->Branch("Event_Vertex_Y", &fBuffer_Event_Vertex_Y, "Event_Vertex_Y/F");
323  fClusterTree->Branch("Event_Vertex_Z", &fBuffer_Event_Vertex_Z, "Event_Vertex_Z/F");
324  fClusterTree->Branch("Event_Multiplicity", &fBuffer_Event_Multiplicity, "Event_Multiplicity/F");
325  fClusterTree->Branch("Event_NumActiveCells", &fBuffer_Event_NumActiveCells, "Event_NumActiveCells/I");
326  }
327 
328  if(fSaveCells)
329  {
330  fClusterTree->Branch("Cluster_Cells_ID", fBuffer_Cells_ID, "Cluster_Cells_ID[Cluster_NumCells]/I");
331  fClusterTree->Branch("Cluster_Cells_E", fBuffer_Cells_E, "Cluster_Cells_E[Cluster_NumCells]/F");
332  fClusterTree->Branch("Cluster_Cells_RelativeEta", fBuffer_Cells_RelativeEta, "Cluster_Cells_RelativeEta[Cluster_NumCells]/F");
333  fClusterTree->Branch("Cluster_Cells_RelativePhi", fBuffer_Cells_RelativePhi, "Cluster_Cells_RelativePhi[Cluster_NumCells]/F");
334  }
335 
337  {
338  fClusterTree->Branch("Surrounding_NCells", &fBuffer_Surrounding_NCells, "Surrounding_NCells/I");
339  fClusterTree->Branch("Surrounding_Cells_ID", fBuffer_Surrounding_Cells_ID, "Surrounding_Cells_ID[Surrounding_NCells]/I");
340  fClusterTree->Branch("Surrounding_Cells_R", fBuffer_Surrounding_Cells_R, "Surrounding_Cells_R[Surrounding_NCells]/F");
341  fClusterTree->Branch("Surrounding_Cells_E", fBuffer_Surrounding_Cells_E, "Surrounding_Cells_E[Surrounding_NCells]/F");
342  fClusterTree->Branch("Surrounding_Cells_RelativeEta", fBuffer_Surrounding_Cells_RelativeEta, "Surrounding_Cells_RelativeEta[Surrounding_NCells]/F");
343  fClusterTree->Branch("Surrounding_Cells_RelativePhi", fBuffer_Surrounding_Cells_RelativePhi, "Surrounding_Cells_RelativePhi[Surrounding_NCells]/F");
344  }
345  if(fSaveTracks)
346  {
347  fClusterTree->Branch("Surrounding_NTracks", &fBuffer_Surrounding_NTracks, "Surrounding_NTracks/I");
348  fClusterTree->Branch("Surrounding_Tracks_R", fBuffer_Surrounding_Tracks_R, "Surrounding_Tracks_R[Surrounding_NTracks]/F");
349  fClusterTree->Branch("Surrounding_Tracks_Pt", fBuffer_Surrounding_Tracks_Pt, "Surrounding_Tracks_Pt[Surrounding_NTracks]/F");
350  fClusterTree->Branch("Surrounding_Tracks_RelativeEta", fBuffer_Surrounding_Tracks_RelativeEta, "Surrounding_Tracks_RelativeEta[Surrounding_NTracks]/F");
351  fClusterTree->Branch("Surrounding_Tracks_RelativePhi", fBuffer_Surrounding_Tracks_RelativePhi, "Surrounding_Tracks_RelativePhi[Surrounding_NTracks]/F");
352  }
353 
355  {
356  fClusterTree->Branch("Cluster_MC_Label", &fBuffer_Cluster_MC_Label, "Cluster_MC_Label/I");
357  fClusterTree->Branch("Mother_MC_Label", &fBuffer_Mother_MC_Label, "Mother_MC_Label/I");
358  }
359 
360  fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data());
361  OpenFile(2);
362  PostData(2, fClusterTree);
363 }
364 //_____________________________________________________________________________
366 {
371  }
372 
373  return kTRUE;
374 }
375 //________________________________________________________________________
377 
378  Int_t eventQuality = ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEventQuality();
379  if(eventQuality != 0){// Event Not Accepted
380  return;
381  }
382  fInputEvent = InputEvent();
383  if(fIsMC) fMCEvent = MCEvent();
384 
386  if(eventNotAccepted) return; // Check Centrality, PileUp, SDD and V0AND --> Not Accepted => eventQuality = 1
387 
388  fGeomEMCAL = AliEMCALGeometry::GetInstance();
389  if(!fGeomEMCAL){ AliFatal("EMCal geometry not initialized!");}
390 
391  Int_t nclus = 0;
392  TClonesArray * arrClustersProcess = NULL;
393  // fNCurrentClusterBasic = 0;
394  if(!fCorrTaskSetting.CompareTo("")){
395  nclus = fInputEvent->GetNumberOfCaloClusters();
396  } else {
397  arrClustersProcess = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(Form("%sClustersBranch",fCorrTaskSetting.Data())));
398  if(!arrClustersProcess)
399  AliFatal(Form("%sClustersBranch was not found in AliAnalysisTaskGammaCalo! Check the correction framework settings!",fCorrTaskSetting.Data()));
400  nclus = arrClustersProcess->GetEntries();
401  }
402 
403  if(nclus == 0) return;
404 
405  ((AliCaloPhotonCuts*)fClusterCutsEMC)->FillHistogramsExtendedQA(fInputEvent,fIsMC);
406  // ((AliCaloPhotonCuts*)fClusterCutsDMC)->FillHistogramsExtendedQA(fInputEvent,fIsMC);
407 
408  // match tracks to clusters
409  ((AliCaloPhotonCuts*)fClusterCutsEMC)->MatchTracksToClusters(fInputEvent,fWeightJetJetMC,kTRUE, fMCEvent);
410  // ((AliCaloPhotonCuts*)fClusterCutsDMC)->MatchTracksToClusters(fInputEvent,fWeightJetJetMC,kTRUE, fMCEvent);
411 
412  AliVCaloCells* cells = fInputEvent->GetEMCALCells();
413  Double_t centrality = 0;
415  centrality = GetCentrality(fInputEvent);
416 
417  Int_t numberOfActiveCells50MeV =0;
418  Int_t numberOfActiveCells80MeV =0;
419  Int_t numberOfActiveCells100MeV =0;
420  Int_t numberOfActiveCells150MeV =0;
421  for(Int_t aCell=0;aCell<cells->GetNumberOfCells();aCell++){
422  // Define necessary variables
423  Short_t cellNumber = 0;
424  Double_t cellAmplitude = 0, cellTime = 0, cellEFrac = 0;
425  Int_t cellMCLabel = 0;
426  // Get Cell ID
427  cells->GetCell(aCell,cellNumber,cellAmplitude,cellTime,cellMCLabel,cellEFrac);
428  Double_t cellEnergy = cells->GetCellAmplitude(cellNumber);
429  hEActiveCellsVsCentrality->Fill(cellEnergy,centrality);
430  if(cellEnergy >0.050 ){
431  numberOfActiveCells50MeV++;
432  hEActiveCells50MeVVsCentrality->Fill(cellEnergy,centrality);
433  }
434  if(cellEnergy >0.080 ){
435  numberOfActiveCells80MeV++;
436  hEActiveCells80MeVVsCentrality->Fill(cellEnergy,centrality);
437  }
438  if(cellEnergy >0.100 ){
439  numberOfActiveCells100MeV++;
440  hEActiveCells100MeVVsCentrality->Fill(cellEnergy,centrality);
441  }
442  if(cellEnergy >0.150 ){
443  numberOfActiveCells150MeV++;
444  hEActiveCells150MeVVsCentrality->Fill(cellEnergy,centrality);
445  }
446  }
447  hNActiveCellsVsCentrality->Fill(cells->GetNumberOfCells(),centrality);
448  hNActiveCellsAbove50MeVVsCentrality->Fill(numberOfActiveCells50MeV,centrality);
449  hNActiveCellsAbove80MeVVsCentrality->Fill(numberOfActiveCells80MeV,centrality);
450  hNActiveCellsAbove100MeVVsCentrality->Fill(numberOfActiveCells100MeV,centrality);
451  hNActiveCellsAbove150MeVVsCentrality->Fill(numberOfActiveCells150MeV,centrality);
452  }
453 
454  // vertex
455  Double_t vertex[3] = {0};
456  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
457 
458  map<Long_t,Int_t> mapIsClusterAccepted;
459  map<Long_t,Int_t> mapIsClusterAcceptedWithoutTrackMatch;
460  Int_t nCellInCluster = 0;
461  // Loop over EMCal clusters
462  for(Long_t i = 0; i < nclus; i++){
463  Double_t tempClusterWeight = fWeightJetJetMC;
464  AliVCluster* clus = NULL;
465  if(fInputEvent->IsA()==AliESDEvent::Class()){
466  if(arrClustersProcess)
467  clus = new AliESDCaloCluster(*(AliESDCaloCluster*)arrClustersProcess->At(i));
468  else
469  clus = new AliESDCaloCluster(*(AliESDCaloCluster*)fInputEvent->GetCaloCluster(i));
470  } else if(fInputEvent->IsA()==AliAODEvent::Class()){
471  if(arrClustersProcess)
472  clus = new AliAODCaloCluster(*(AliAODCaloCluster*)arrClustersProcess->At(i));
473  else
474  clus = new AliAODCaloCluster(*(AliAODCaloCluster*)fInputEvent->GetCaloCluster(i));
475  }
476  if(!clus) continue;
477 
478 
479  if ( !clus->IsEMCAL()){ // for PHOS: cluster->GetType() == AliVCluster::kPHOSNeutral
480  delete clus;
481  continue;
482  }
484  nCellInCluster+=clus->GetNCells();
485  for(Int_t iCell=0;iCell<clus->GetNCells();iCell++){
486  hECellsInClustersVsCentrality->Fill(cells->GetCellAmplitude(clus->GetCellAbsId(iCell)),centrality);
487  }
488  }
489 
490  if ( clus->E() < fMinClusterEnergy){
491  delete clus;
492  continue;
493  }
494 
495  if(!((AliCaloPhotonCuts*)fClusterCutsEMC)->ClusterIsSelected(clus,fInputEvent,fMCEvent,fIsMC, tempClusterWeight,i) && !((AliCaloPhotonCuts*)fClusterCutsDMC)->ClusterIsSelected(clus,fInputEvent,fMCEvent,fIsMC, tempClusterWeight,i)){
496  delete clus;
497  continue;
498  }
500  }
502  hNCellsInClustersVsCentrality->Fill(nCellInCluster,centrality);
503 
504  PostData(1, fOutputList);
505 }
506 
508 void AliAnalysisTaskClusterQA::ProcessQATreeCluster(AliVEvent *event, AliVCluster* cluster, Long_t indexCluster){
509 
510  Float_t clusPos[3] = { 0,0,0 };
511  cluster->GetPosition(clusPos);
512  TVector3 clusterVector(clusPos[0],clusPos[1],clusPos[2]);
513  Double_t etaCluster = clusterVector.Eta();
514  Double_t phiCluster = clusterVector.Phi();
515  if (phiCluster < 0) phiCluster += 2*TMath::Pi();
516 
518  // Vertex position x, y, z
519  fBuffer_Event_Vertex_X = fInputEvent->GetPrimaryVertex()->GetX();
520  fBuffer_Event_Vertex_Y = fInputEvent->GetPrimaryVertex()->GetY();
521  fBuffer_Event_Vertex_Z = fInputEvent->GetPrimaryVertex()->GetZ();
522 
523  // Unique event ID of the current cluster
524  // AliESDEvent* esdEvent = dynamic_cast<AliESDEvent*>(event);
525  fBuffer_EventID = 0;
526  // if(esdEvent->GetHeader()){
527  // fBuffer_EventID = GetUniqueEventID(esdEvent->GetHeader());
528  // printf("event id: %lld\n",fBuffer_EventID);
529  // }
530 
531  // V0-based multiplicity of the event
532  fBuffer_Event_Multiplicity = fInputEvent->GetVZEROData()->GetMTotV0A()+fInputEvent->GetVZEROData()->GetMTotV0C();
533  }
534 
535  // Properties of the current cluster
536  fBuffer_ClusterE = cluster->E();
537 
538  if(etaCluster < 0.67 && etaCluster > -0.67){
539  if(phiCluster < 3.2 && phiCluster > 1.4){
540  fBuffer_ClusterIsEMCAL = kTRUE;
541  }
542  }
543  if(etaCluster < 0.67 && etaCluster > -0.67){
544  if(etaCluster > 0.23 || etaCluster < -0.23){
545  if(phiCluster < 5.6 && phiCluster > 4.6){
546  fBuffer_ClusterIsEMCAL = kFALSE;
547  }
548  }
549  }
550 
551  fBuffer_ClusterPhi = phiCluster;
552  fBuffer_ClusterEta = etaCluster;
553  fBuffer_ClusterM02 = cluster->GetM02();
554  fBuffer_ClusterM20 = cluster->GetM20();
555 
556  // Get all cells from the event
557  AliVCaloCells* cells = NULL;
558  if(cluster->IsEMCAL())
559  cells = event->GetEMCALCells();
560  else
561  return;
562 
563  // Determine all active cells in the event
565  fBuffer_Event_NumActiveCells = cells->GetNumberOfCells();
566 
567  // Get the number of cells from the current cluster
568  Int_t nCellCluster = cluster->GetNCells();
569  fBuffer_ClusterNumCells = nCellCluster;
570 
571  // Find the leading cell in the cluster and its position
573  fBuffer_LeadingCell_E = cells->GetCellAmplitude(fBuffer_LeadingCell_ID);
574  Float_t leadcelleta = 0.;
575  Float_t leadcellphi = 0.;
576  fGeomEMCAL->EtaPhiFromIndex(fBuffer_LeadingCell_ID, leadcelleta, leadcellphi);
577  if ( leadcellphi < 0 ) leadcellphi+=TMath::TwoPi();
578  fBuffer_LeadingCell_Eta = leadcelleta;
579  fBuffer_LeadingCell_Phi = leadcellphi;
580  // Treat the remaining cells of the cluster and get their relative position compared to the leading cell
581  if(fSaveCells){
582  for(Int_t iCell=0;iCell<nCellCluster;iCell++){
583  if(iCell<100){ // maximum number of cells for a cluster is set to 100
584  fBuffer_Cells_ID[iCell] = cluster->GetCellAbsId(iCell);
585  fBuffer_Cells_E[iCell] = cells->GetCellAmplitude(cluster->GetCellAbsId(iCell));
586  Float_t celleta = 0.;
587  Float_t cellphi = 0.;
588  fGeomEMCAL->EtaPhiFromIndex(fBuffer_Cells_ID[iCell], celleta, cellphi);
589  if ( cellphi < 0 ) cellphi+=TMath::TwoPi();
590  fBuffer_Cells_RelativeEta[iCell] = leadcelleta-celleta;
591  fBuffer_Cells_RelativePhi[iCell] = leadcellphi-cellphi;
592  }
593  }
594  }
596  Int_t nActiveCellsSurroundingInR = 0;
597  // fill surrounding cell buffer
598  for(Int_t aCell=0;aCell<cells->GetNumberOfCells();aCell++){
599  // Define necessary variables
600  Short_t cellNumber = 0;
601  Double_t cellAmplitude = 0, cellTime = 0, cellEFrac = 0;
602  Int_t cellMCLabel = 0;
603  Float_t surrcelleta = 0.;
604  Float_t surrcellphi = 0.;
605  // Get Cell ID
606  cells->GetCell(aCell,cellNumber,cellAmplitude,cellTime,cellMCLabel,cellEFrac);
607 
608  // Get eta and phi for the surounding cells
609  fGeomEMCAL->EtaPhiFromIndex(cellNumber, surrcelleta, surrcellphi);
610  if ( surrcellphi < 0 ) surrcellphi+=TMath::TwoPi();
611  Double_t dR2 = pow(leadcelleta-surrcelleta,2) + pow(leadcellphi-surrcellphi,2);
612  // Select those cells that are within fConeRadius of the leading cluster cell
613  if( dR2 < fConeRadius){
614  fBuffer_Surrounding_Cells_E[aCell] = cells->GetCellAmplitude(cellNumber);
615  fBuffer_Surrounding_Cells_ID[aCell] = cellNumber;
616  fBuffer_Surrounding_Cells_R[aCell] = dR2;
617  fBuffer_Surrounding_Cells_RelativeEta[aCell] = leadcelleta-surrcelleta;
618  fBuffer_Surrounding_Cells_RelativePhi[aCell] = leadcellphi-surrcellphi;
619  nActiveCellsSurroundingInR+=1;
620  }
621  }
622  fBuffer_Surrounding_NCells = nActiveCellsSurroundingInR;
623  }
624 
625  // write PDG code of mother of particle that mainly contributed to cluster to tree
626  if(fIsMC> 0){
627  if (cluster->GetNLabels()>0){
628  if((cluster->GetLabelAt(0)!=-1)){
629  fBuffer_Mother_MC_Label = (fMCEvent->Particle(cluster->GetLabelAt(0)))->GetPdgCode(); // mother of leading contribution
630  } else{
631  // mother is initial collision
633  }
634  } else{
635  // no mothers found (should not happen)
637  }
638  }
639  if(fIsMC) fBuffer_Cluster_MC_Label = MakePhotonCandidates(cluster, cells,indexCluster);
640  if(fSaveTracks) ProcessTracksAndMatching(cluster,indexCluster);
641  // Add everything to the tree
642  if (fClusterTree) fClusterTree->Fill();
643 }
644 
645 //________________________________________________________________________
647 
648 }
649 
650 //________________________________________________________________________
651 Int_t AliAnalysisTaskClusterQA::FindLargestCellInCluster(AliVCluster* cluster, AliVEvent* event){
652 
653  const Int_t nCells = cluster->GetNCells();
654  AliVCaloCells* cells = NULL;
655 
656  if(cluster->IsEMCAL())
657  cells = event->GetEMCALCells();
658  else if(cluster->IsPHOS())
659  cells = event->GetPHOSCells();
660 
661  Float_t eMax = 0.;
662  Int_t idMax = -1;
663 
664  if (nCells < 1) return idMax;
665  for (Int_t iCell = 0;iCell < nCells;iCell++){
666  Int_t cellAbsID = cluster->GetCellsAbsId()[iCell];
667  if (cells->GetCellAmplitude(cellAbsID)> eMax){
668  eMax = cells->GetCellAmplitude(cellAbsID);
669  idMax = cellAbsID;
670  }
671  }
672  return idMax;
673 }
674 
675 //________________________________________________________________________
677  Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0;
678  row=0;
679  column=0;
680  // Get SM number and relative row/column for SM
681  fGeomEMCAL->GetCellIndex(cellIndex, nSupMod,nModule,nIphi,nIeta);
682  fGeomEMCAL->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, row,column);
683  row += nSupMod/2 * 24;
684  column += nSupMod%2 * 48;
685 }
686 
687 //________________________________________________________________________
688 Int_t AliAnalysisTaskClusterQA::MakePhotonCandidates(AliVCluster* clus, AliVCaloCells* cells, Long_t indexCluster){
689 
690  Double_t vertex[3] = {0};
691  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
692 
693  TLorentzVector clusterVector;
694  clus->GetMomentum(clusterVector,vertex);
695 
696  TLorentzVector* tmpvec = new TLorentzVector();
697  tmpvec->SetPxPyPzE(clusterVector.Px(),clusterVector.Py(),clusterVector.Pz(),clusterVector.E());
698 
699  // convert to AODConversionPhoton
700  AliAODConversionPhoton *PhotonCandidate=new AliAODConversionPhoton(tmpvec);
701  if(!PhotonCandidate){
702  delete clus;
703  delete tmpvec;
704  if (PhotonCandidate) delete PhotonCandidate;
705  return -9;
706  }
707  // Flag Photon as CaloPhoton
708  PhotonCandidate->SetIsCaloPhoton();
709  PhotonCandidate->SetCaloClusterRef(indexCluster);
710 
711  // get MC label
712  if(fIsMC> 0){
713  Int_t* mclabelsCluster = clus->GetLabels();
714  Int_t nValidClusters = 0;
715  if (clus->GetNLabels()>0){
716  for (Int_t k =0; k< (Int_t)clus->GetNLabels(); k++){
717  if (mclabelsCluster[k]>0){
718  if (nValidClusters< 50)PhotonCandidate->SetCaloPhotonMCLabel(nValidClusters,mclabelsCluster[k]);
719  nValidClusters++;
720  }
721  }
722  }
723  PhotonCandidate->SetNCaloPhotonMCLabels(nValidClusters);
724  }
725 
726  AliAODCaloCluster* clusSub1 = new AliAODCaloCluster();
727  AliAODCaloCluster* clusSub2 = new AliAODCaloCluster();
728 
729 
730  // split clusters according to their shares in the cluster (NLM == 1) needs to be treated differently
731  if (fMinNLMCut == 1 && fMaxNLMCut == 1 ){
732  Int_t absCellIdFirst = ((AliCaloPhotonCuts*)fClusterCutsEMC)->FindLargestCellInCluster(clus, fInputEvent);
733  Int_t absCellIdSecond = ((AliCaloPhotonCuts*)fClusterCutsEMC)->FindSecondLargestCellInCluster(clus, fInputEvent);
734 
735  ((AliCaloPhotonCuts*)fClusterCutsEMC)->SplitEnergy(absCellIdFirst, absCellIdSecond, clus, fInputEvent, fIsMC, clusSub1, clusSub2);
736  } else if (fMinNLMCut > 1 ){
737  const Int_t nc = clus->GetNCells();
738  Int_t absCellIdList[nc];
739 
740  ((AliCaloPhotonCuts*)fClusterCutsEMC)->SplitEnergy(absCellIdList[0], absCellIdList[1], clus, fInputEvent, fIsMC, clusSub1, clusSub2);
741  }
742 
743  // TLorentzvector with sub cluster 1
744  TLorentzVector clusterVector1;
745  clusSub1->GetMomentum(clusterVector1,vertex);
746  TLorentzVector* tmpvec1 = new TLorentzVector();
747  tmpvec1->SetPxPyPzE(clusterVector1.Px(),clusterVector1.Py(),clusterVector1.Pz(),clusterVector1.E());
748  // convert to AODConversionPhoton
749  AliAODConversionPhoton *PhotonCandidate1=new AliAODConversionPhoton(tmpvec1);
750  if(!PhotonCandidate1){
751  delete clusSub1;
752  delete tmpvec1;
753  return -9;
754  }
755  // Flag Photon as CaloPhoton
756  PhotonCandidate1->SetIsCaloPhoton();
757  // TLorentzvector with sub cluster 2
758  TLorentzVector clusterVector2;
759  clusSub2->GetMomentum(clusterVector2,vertex);
760  TLorentzVector* tmpvec2 = new TLorentzVector();
761  tmpvec2->SetPxPyPzE(clusterVector2.Px(),clusterVector2.Py(),clusterVector2.Pz(),clusterVector2.E());
762  // convert to AODConversionPhoton
763  AliAODConversionPhoton *PhotonCandidate2=new AliAODConversionPhoton(tmpvec2);
764  if(!PhotonCandidate2){
765  delete clusSub2;
766  delete tmpvec2;
767  return -9;
768  }
769  // Flag Photon as CaloPhoton
770  PhotonCandidate2->SetIsCaloPhoton();
771  Int_t mclabel = -3;
772  // create pi0 candidate
773  AliAODConversionMother *pi0cand = new AliAODConversionMother(PhotonCandidate1,PhotonCandidate2);
774 // if((((AliConversionMesonCuts*)fMesonCuts)->MesonIsSelected(pi0cand,kTRUE,fEventCuts->GetEtaShift()))){
775  if(fIsMC> 0 && PhotonCandidate && PhotonCandidate1 && PhotonCandidate2 && fSaveMCInformation){
776  // if(fInputEvent->IsA()==AliESDEvent::Class())
777  mclabel = ProcessTrueClusterCandidates(PhotonCandidate, clus->GetM02(), PhotonCandidate1, PhotonCandidate2);
778  // if(fInputEvent->IsA()==AliAODEvent::Class())
779  // ProcessTrueClusterCandidatesAOD(PhotonCandidate, clus->GetM02(), PhotonCandidate1, PhotonCandidate2);
780  return mclabel;
781  }
782 // } else {
783 // return -7;
784 // }
785  return -1;
786 }
787 
788 //________________________________________________________________________
789 void AliAnalysisTaskClusterQA::ProcessTracksAndMatching(AliVCluster* clus, Long_t indexCluster){
790 
791 
792  Int_t nTracksInR = 0;
793  Int_t nModules = fGeomEMCAL->GetNumberOfSuperModules();
794 
795  AliESDEvent *esdev = dynamic_cast<AliESDEvent*>(fInputEvent);
796  AliESDtrackCuts *EsdTrackCuts = 0x0;
797  // Using standard function for setting Cuts
798  static int prevRun = -1;
799  // Using standard function for setting Cuts
800  if (esdev){
801  Int_t runNumber = fInputEvent->GetRunNumber();
802  if (prevRun!=runNumber) {
803  delete EsdTrackCuts;
804  EsdTrackCuts = 0;
805  prevRun = runNumber;
806  }
807  // if LHC11a or earlier or if LHC13g or if LHC12a-i -> use 2010 cuts
808  if( (runNumber<=146860) || (runNumber>=197470 && runNumber<=197692) || (runNumber>=172440 && runNumber<=193766) ){
809  EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
810  // else if run2 data use 2015 PbPb cuts
811  }else if (runNumber>=209122){
812  // hard coded track cuts for the moment, because AliESDtrackCuts::GetStandardITSTPCTrackCuts2015PbPb() gives spams warnings
813  EsdTrackCuts = new AliESDtrackCuts();
814  EsdTrackCuts->AliESDtrackCuts::SetMinNCrossedRowsTPC(70);
815  EsdTrackCuts->AliESDtrackCuts::SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
816  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2PerClusterTPC(4);
817  EsdTrackCuts->AliESDtrackCuts::SetAcceptKinkDaughters(kFALSE);
818  EsdTrackCuts->AliESDtrackCuts::SetRequireTPCRefit(kTRUE);
819  // ITS
820  EsdTrackCuts->AliESDtrackCuts::SetRequireITSRefit(kTRUE);
821  EsdTrackCuts->AliESDtrackCuts::SetClusterRequirementITS(AliESDtrackCuts::kSPD,
822  AliESDtrackCuts::kAny);
823  EsdTrackCuts->AliESDtrackCuts::SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.1");
824  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2TPCConstrainedGlobal(36);
825  EsdTrackCuts->AliESDtrackCuts::SetMaxDCAToVertexZ(2);
826  EsdTrackCuts->AliESDtrackCuts::SetDCAToVertex2D(kFALSE);
827  EsdTrackCuts->AliESDtrackCuts::SetRequireSigmaToVertex(kFALSE);
828  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2PerClusterITS(36);
829  // else use 2011 version of track cuts
830  }else{
831  EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
832  }
833  EsdTrackCuts->SetMaxDCAToVertexZ(2);
834  EsdTrackCuts->SetEtaRange(-0.8, 0.8);
835  EsdTrackCuts->SetPtRange(0.15);
836  }
837 
838  for (Int_t itr=0;itr<fInputEvent->GetNumberOfTracks();itr++){
839  AliVTrack *inTrack = esdev->GetTrack(itr);
840  if(!inTrack) continue;
841  if(inTrack->Pt()<fMinTrackPt) continue;
842  AliESDtrack *esdt = dynamic_cast<AliESDtrack*>(inTrack);
843  // check track cuts
844  if(!EsdTrackCuts->AcceptTrack(esdt)) continue;
845  const AliExternalTrackParam *in = esdt->GetInnerParam();
846  if (!in){AliError("Could not get InnerParam of Track, continue");continue;}
847  AliExternalTrackParam *trackParam =new AliExternalTrackParam(*in);
848  AliExternalTrackParam emcParam(*trackParam);
849  Float_t eta, phi, pt;
850 
851  //propagate tracks to emc surfaces
852  if (!AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&emcParam, 440., 0.139, 20., eta, phi, pt)) {
853  delete trackParam;
854  continue;
855  }
856  if( TMath::Abs(eta) > 0.75 ) {
857  delete trackParam;
858  continue;
859  }
860  // Save some time and memory in case of no DCal present
861  if( nModules < 13 && ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad())){
862  delete trackParam;
863  continue;
864  }
865  // Save some time and memory in case of run2
866  if( nModules > 12 ){
867  if (( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()) && ( phi < 250*TMath::DegToRad() || phi > 340*TMath::DegToRad())){
868  delete trackParam;
869  continue;
870  }
871  }
872  Float_t dEta=-999, dPhi=-999;
873  Double_t trkPos[3] = {0.,0.,0.};
874  if (!emcParam.GetXYZ(trkPos)){
875  delete trackParam;
876  continue;
877  }
878 
879  AliExternalTrackParam trackParamTmp(emcParam);//Retrieve the starting point every time before the extrapolation
880  if(!AliEMCALRecoUtils::ExtrapolateTrackToCluster(&trackParamTmp, clus, 0.139, 5., dEta, dPhi)) continue;
881 
882  Float_t dR2 = dPhi*dPhi + dEta*dEta;
883  if(dR2 < fConeRadius){
884  fBuffer_Surrounding_Tracks_R[nTracksInR]=dR2;
885  fBuffer_Surrounding_Tracks_Pt[nTracksInR]=inTrack->Pt();
888  nTracksInR+=1;
889  }
890  }
891  fBuffer_Surrounding_NTracks = nTracksInR;
892  if(nTracksInR==0){
893  fBuffer_Surrounding_Tracks_R[nTracksInR]=-1;
894  fBuffer_Surrounding_Tracks_Pt[nTracksInR]=-1;
897  }
898 
899  if(EsdTrackCuts){
900  delete EsdTrackCuts;
901  EsdTrackCuts=0x0;
902  }
903 }
904 
905 //________________________________________________________________________
907  AliAODConversionPhoton *TrueSubClusterCandidate1,
908  AliAODConversionPhoton *TrueSubClusterCandidate2)
909 {
910 
911  Int_t mcLabelCluster = -5;
912  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
913  Double_t mcProdVtxX = primVtxMC->GetX();
914  Double_t mcProdVtxY = primVtxMC->GetY();
915  Double_t mcProdVtxZ = primVtxMC->GetZ();
916 
917  TParticle *Photon = NULL;
918  if (!TrueClusterCandidate->GetIsCaloPhoton()) AliFatal("CaloPhotonFlag has not been set task will abort");
919  if (TrueClusterCandidate->GetCaloPhotonMCLabel(0) < 0){
920  mcLabelCluster = -10;
921  return mcLabelCluster;
922  }
923  if (TrueClusterCandidate->GetNCaloPhotonMCLabels()>0){
924  Photon = fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMCLabel(0));
925  } else{
926  mcLabelCluster = -11;
927  return mcLabelCluster;
928  }
929  if(Photon == NULL){
930  mcLabelCluster = -12;
931  return mcLabelCluster;
932  }
933  AliAODConversionMother *mesoncand = new AliAODConversionMother(TrueSubClusterCandidate1,TrueSubClusterCandidate2);
934 // Bool_t mesonIsSelected = (((AliConversionMesonCuts*)fMesonCuts)->MesonIsSelected(mesoncand,kTRUE,fEventCuts->GetEtaShift()));
935 // if (!mesonIsSelected){
936 // delete mesoncand;
937 // mcLabelCluster = -13;
938 // return mcLabelCluster;
939 // }
940 
941  TrueClusterCandidate->SetCaloPhotonMCFlags(fMCEvent, kFALSE);
942 
943  // cluster classification:
944  // 1 - nice merged cluster (2 gamma | contributions from 2 gamma) from pi0/eta
945  // 2 - contribution from only 1 partner (1 gamma, 1 fully coverted gamma) from pi0/eta
946  // 3 - contribution from part of 1 partner (1 electron) from pi0/eta
947 
948 
949  Int_t clusterClass = 0;
950  Bool_t isPrimary = fEventCuts->IsConversionPrimaryESD( fMCEvent, TrueClusterCandidate->GetCaloPhotonMCLabel(0), mcProdVtxX, mcProdVtxY, mcProdVtxZ);
951 
952 
953  Long_t motherLab = -1;
954  if (TrueClusterCandidate->IsMerged() || TrueClusterCandidate->IsMergedPartConv()){
955  clusterClass = 1;
956  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
957  } else if (TrueClusterCandidate->GetNCaloPhotonMotherMCLabels()> 0){
958  if (TrueClusterCandidate->IsLargestComponentElectron() || TrueClusterCandidate->IsLargestComponentPhoton()){
959  if (TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0) > -1 && (fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0))->GetPdgCode() == 111 || fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0))->GetPdgCode() == 221) ){
960  if ( TrueClusterCandidate->IsConversion() && !TrueClusterCandidate->IsConversionFullyContained() ){
961  clusterClass = 3;
962  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
963  } else {
964  clusterClass = 2;
965  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
966  }
967  }
968  } else if (TrueClusterCandidate->IsSubLeadingEM()){
969  if (TrueClusterCandidate->GetNCaloPhotonMotherMCLabels()> 1){
970  if ( TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1) > -1){
971  if (fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1))->GetPdgCode() == 111 || fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1))->GetPdgCode() == 221 ){
972  clusterClass = 2;
973  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1);
974  }
975  }
976  }
977  } else {
978  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
979  }
980  }
981 
982  // Get Mother particle
983  TParticle *mother = NULL;
984  Int_t motherPDG = -1;
985  if (motherLab > -1)
986  mother = fMCEvent->Particle(motherLab);
987  if (mother)
988  motherPDG = TMath::Abs(mother->GetPdgCode());
989 
990  //
991  if (clusterClass == 1 || clusterClass == 2 || clusterClass == 3 ){
992  // separate different components
993  if (clusterClass == 1 && TrueClusterCandidate->IsMerged()){
994  if (motherPDG == 111){
995  mcLabelCluster = 10; // NOTE merged pi0
996  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
997  mcLabelCluster = 11; // NOTE secondary merged pi0
998  }
999  if (motherPDG == 221)
1000  mcLabelCluster = 12; // NOTE merged eta
1001  } else if (clusterClass == 1 && TrueClusterCandidate->IsMergedPartConv()){
1002  if (motherPDG == 111){
1003  mcLabelCluster = 13; // NOTE merged pi0 with one converted gamma
1004  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
1005  mcLabelCluster = 14; // NOTE merged secondary pi0 with one converted gamma
1006  }
1007  if (motherPDG == 221)
1008  mcLabelCluster = 15; // NOTE merged eta with one converted gamma
1009  } else if (clusterClass == 2){
1010  if (motherPDG == 111){
1011  mcLabelCluster = 20; // NOTE decay photon from pi0
1012  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
1013  mcLabelCluster = 21; // NOTE decay photon from secondary pi0
1014  }
1015  if (motherPDG == 221)
1016  mcLabelCluster = 22; // NOTE decay photon from eta
1017  } else if (clusterClass == 3){
1018  if (motherPDG == 111) {
1019  mcLabelCluster = 30; // NOTE electron from decayed pi0
1020  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
1021  mcLabelCluster = 31; // NOTE electron from decayed secondary pi0
1022  }
1023  if (motherPDG == 221)
1024  mcLabelCluster = 32; // NOTE electron from decayed eta
1025  }
1026 
1027  // leading particle is a photon or the conversion is fully contained and its not from pi0 || eta
1028  } else if (TrueClusterCandidate->IsLargestComponentPhoton() || TrueClusterCandidate->IsConversionFullyContained()
1029  || TrueClusterCandidate->IsElectronFromFragPhoton()){
1030 
1031  if(TrueClusterCandidate->IsLargestComponentPhoton()){
1032  // cluster is produced by a photon that either has no mother or the mother is neither from a pi^0 nor eta, e.g. inital collision
1033 
1034  if (motherLab == -1) mcLabelCluster = 40; // direct photon from initial collision
1035  else if ((motherLab >0) && (motherLab<9)) mcLabelCluster = 41; // photon from quark
1036  else if (motherLab == 11) mcLabelCluster = 42; // photon from electron
1037  else if (motherLab == 22){ // check for frag photon
1038 
1039  TParticle *dummyMother = fMCEvent->Particle(motherLab);
1040  Bool_t originReached = kFALSE;
1041  Bool_t isFragPhoton = kFALSE;
1042 
1043  while (dummyMother->GetPdgCode() == 22 && !originReached){ // follow photon's history, as long as the mother is a photon
1044  if (dummyMother->GetMother(0) > -1){
1045  dummyMother = fMCEvent->Particle(dummyMother->GetMother(0));
1046  if (TMath::Abs(dummyMother->GetPdgCode()) == 22){ // if mother of mother is again a photon, continue
1047  if (dummyMother->GetMother(0) > -1){
1048  dummyMother = fMCEvent->Particle(dummyMother->GetMother(0));
1049  } else {
1050  originReached = kTRUE;
1051  }
1052  } else {
1053  originReached = kTRUE;
1054  }
1055  isFragPhoton = (TMath::Abs(dummyMother->GetPdgCode()) < 6);// photon stems from quark = fragmentation photon
1056  } else {
1057  originReached = kTRUE;
1058  }
1059  }
1060 
1061  if(isFragPhoton) mcLabelCluster = 43; // Fragmentation photon
1062  else{
1063  mcLabelCluster = 47; // something like cluster <- photon <- photon <- X (not photon)
1064  AliInfo(Form("Mother of photon is photon etc. but origin is not quark. ID: %i", motherLab));
1065  }
1066  }
1067  else{
1068  mcLabelCluster = 44; // other (e.g. from meson decays that are not pi0 or eta0)
1069  AliInfo(Form("Single cluster is mainly produced by a photon and mother is %i", motherLab));
1070  }
1071  } else if(TrueClusterCandidate->IsConversionFullyContained()){
1072  // cluster is from a fully contained conversion
1073  mcLabelCluster = 45;
1074  } else if(TrueClusterCandidate->IsElectronFromFragPhoton()){
1075  mcLabelCluster = 46; // electron from frac
1076  }
1077 
1078  // leading particle is an electron and its not from pi0 || eta and no electron from fragmentation photon conversion
1079  } else if (TrueClusterCandidate->IsLargestComponentElectron() && !TrueClusterCandidate->IsElectronFromFragPhoton()){
1080  mcLabelCluster = 50; // NOTE single electron
1081  } else {
1082  // leading particle from hadron
1083  mcLabelCluster = 60; // NOTE hadron cluster
1084  AliInfo(Form("Single cluster is mainly produced by hadron with id: %i", motherLab));
1085  }
1086 
1087  delete mesoncand;
1088  return mcLabelCluster;
1089 }
1090 
1091 //_____________________________________________________________________________
1092 ULong64_t AliAnalysisTaskClusterQA::GetUniqueEventID(AliVHeader* header)
1093 {
1094  // To have a unique id for each event in a run!
1095  // Modified from AliRawReader.h
1096  return ((ULong64_t)header->GetBunchCrossNumber()+
1097  (ULong64_t)header->GetOrbitNumber()*3564+
1098  (ULong64_t)header->GetPeriodNumber()*16777215*3564);
1099 }
1100 
1101 //-------------------------------------------------------------
1103 { // Get Event Centrality
1104 
1105  AliESDEvent *esdEvent=dynamic_cast<AliESDEvent*>(event);
1106  if(esdEvent){
1107  AliMultSelection *MultSelection = (AliMultSelection*)event->FindListObject("MultSelection");
1108  if(!MultSelection){
1109  AliWarning ("AliMultSelection object not found !");
1110  return -1;
1111  }else{
1112  return MultSelection->GetMultiplicityPercentile("V0M");// default
1113  }
1114  }
1115 
1116 
1117  return -1;
1118 }
void SetPeriodEnum(TString periodName)
Int_t fBuffer_Surrounding_NTracks
! array buffer
void SetCaloClusterRef(Long_t ref)
double Double_t
Definition: External.C:58
Float_t fBuffer_ClusterPhi
! array buffer
virtual void UserExec(Option_t *option)
Definition: External.C:236
void SetCaloPhotonMCFlags(AliMCEvent *mcEvent, Bool_t enableSort)
ULong64_t GetUniqueEventID(AliVHeader *header)
Float_t fBuffer_Event_Vertex_Y
! array buffer
Int_t MakePhotonCandidates(AliVCluster *clus, AliVCaloCells *cells, Long_t indexCluster)
PeriodVar GetPeriodEnum()
Bool_t fSaveMCInformation
save MC information
void SetPeriodEnumExplicit(PeriodVar periodEnum)
const Int_t kMaxNTracks
Float_t fBuffer_ClusterE
! array buffer
centrality
Float_t * fBuffer_Cells_RelativeEta
! array buffer
Bool_t fSaveEventProperties
save general event properties (centrality etc.)
Float_t fConeRadius
save arrays of all cells in event
TString GetPeriodName()
void SetCaloPhotonMCLabel(Int_t i, Int_t labelCaloPhoton)
Int_t fBuffer_Mother_MC_Label
! array buffer
Float_t * fBuffer_Surrounding_Tracks_R
! array buffer
Float_t fBuffer_Event_Vertex_Z
! array buffer
Int_t GetCaloPhotonMotherMCLabel(Int_t i)
Bool_t fBuffer_ClusterIsEMCAL
! array buffer
Float_t fBuffer_Event_Vertex_X
! array buffer
Int_t fBuffer_LeadingCell_ID
! array buffer
void GetRowAndColumnFromAbsCellID(Int_t cellIndex, Int_t &row, Int_t &column)
void ProcessTracksAndMatching(AliVCluster *clus, Long_t indexCluster)
Float_t fBuffer_Event_Multiplicity
! array buffer
Int_t * fBuffer_Surrounding_Cells_ID
! array buffer
Int_t FindLargestCellInCluster(AliVCluster *cluster, AliVEvent *event)
Float_t * fBuffer_Surrounding_Cells_E
! array buffer
Float_t fBuffer_LeadingCell_Phi
! array buffer
Int_t IsEventAcceptedByCut(AliConvEventCuts *ReaderCuts, AliVEvent *event, AliMCEvent *mcEvent, Int_t isHeavyIon, Bool_t isEMCALAnalysis)
int Int_t
Definition: External.C:63
Float_t * fBuffer_Surrounding_Cells_R
! array buffer
Class handling all kinds of selection cuts for Gamma Calo analysis.
void ProcessQATreeCluster(AliVEvent *event, AliVCluster *cluster, Long_t indexCluster)
Float_t * fBuffer_Cells_RelativePhi
! array buffer
Float_t * fBuffer_Surrounding_Tracks_RelativeEta
! array buffer
float Float_t
Definition: External.C:68
Float_t * fBuffer_Surrounding_Tracks_Pt
! array buffer
Bool_t IsConversionPrimaryESD(AliMCEvent *mcEvent, Long_t eventpos, Double_t prodVtxX, Double_t prodVtxY, Double_t prodVtxZ)
Int_t fBuffer_Cluster_MC_Label
! array buffer
Bool_t fSaveSurroundingCells
save arrays of all cells in event
Float_t fBuffer_ClusterM20
! array buffer
Bool_t fSaveAdditionalHistos
save MC information
Float_t fMinClusterEnergy
save arrays of all cells in event
Int_t * fBuffer_Cells_ID
! array buffer
Int_t fMaxNLMCut
save MC information
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
short Short_t
Definition: External.C:23
Float_t fBuffer_ClusterM02
! array buffer
virtual void Terminate(Option_t *)
Int_t fBuffer_Event_NumActiveCells
! array buffer
Float_t fBuffer_ClusterEta
! array buffer
Float_t GetCentrality(AliVEvent *event)
Bool_t fSaveCells
save arrays of cluster cells
static Bool_t ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, const AliVCluster *cluster, Double_t mass, Double_t step, Float_t &tmpEta, Float_t &tmpPhi)
ULong64_t fBuffer_EventID
! array buffer
Float_t * fBuffer_Surrounding_Cells_RelativePhi
! array buffer
Class handling all kinds of selection cuts for Gamma Conversion analysis.
Int_t fBuffer_ClusterNumCells
! array buffer
void SetNCaloPhotonMCLabels(Int_t nLabels)
Float_t fBuffer_LeadingCell_Eta
! array buffer
const Int_t kMaxActiveCells
Class handling all kinds of selection cuts for Gamma Conversion analysis.
Float_t * fBuffer_Surrounding_Cells_RelativeEta
! array buffer
AliConversionMesonCuts * fMesonCuts
Float_t fBuffer_LeadingCell_E
! array buffer
AliConvEventCuts * GetEventCuts()
Definition: AliV0ReaderV1.h:90
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
Float_t * fBuffer_Surrounding_Tracks_RelativePhi
! array buffer
Float_t * fBuffer_Cells_E
! array buffer
Int_t fMinNLMCut
save MC information
Int_t fBuffer_Surrounding_NCells
! array buffer
Bool_t fSaveTracks
save arrays of all cells in event
Int_t ProcessTrueClusterCandidates(AliAODConversionPhoton *TrueClusterCandidate, Float_t m02, AliAODConversionPhoton *TrueSubClusterCandidate1, AliAODConversionPhoton *TrueSubClusterCandidate2)
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
Float_t fMinTrackPt
save arrays of all cells in event