AliPhysics  master (3d17d9d)
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  fReaderGammas(NULL),
47  fPIDResponse(NULL),
48  fCorrTaskSetting(""),
49  fConversionCuts(NULL),
50  fEventCuts(NULL),
51  fClusterCutsEMC(NULL),
52  // fMesonCuts(NULL),
53  fMinNLMCut(1),
54  fMaxNLMCut(1),
55  fInputEvent(NULL),
56  fMCEvent(NULL),
57  fWeightJetJetMC(1),
58  fGeomEMCAL(NULL),
59  fClusterTree(NULL),
60  fIsHeavyIon(kFALSE),
61  ffillTree(-100),
62  ffillHistograms(kFALSE),
63  fOutputList(NULL),
64  fIsMC(0),
65  fCorrectForNonlinearity(kFALSE),
66  fSaveEventProperties(0),
67  fSaveCells(0),
68  fSaveSurroundingCells(0),
69  fSaveTracks(0),
70  fConeRadius(0),
71  fMinTrackPt(0),
72  fMinClusterEnergy(0),
73  fSaveMCInformation(0),
74  fSaveAdditionalHistos(0),
75  fSaveEventsInVector(0),
76  fExtractionPercentages(),
77  fExtractionPercentagePtBins(),
78  fBuffer_EventWeight(1),
79  fBuffer_ClusterE(0),
80  fBuffer_ClusterPhi(0),
81  fBuffer_ClusterEta(0),
82  fBuffer_ClusterIsEMCAL(kFALSE),
83  fBuffer_ClusterSupMod(-1),
84  fBuffer_MC_Cluster_Flag(0),
85  fBuffer_ClusterNumCells(0),
86  // fBuffer_ClusterNLM(0),
87  fBuffer_LeadingCell_ID(0),
88  fBuffer_LeadingCell_E(0),
89  fBuffer_LeadingCell_Eta(0),
90  fBuffer_LeadingCell_Phi(0),
91  fBuffer_ClusterM02(0),
92  fBuffer_ClusterM20(0),
93  fBuffer_Event_Vertex_X(0),
94  fBuffer_Event_Vertex_Y(0),
95  fBuffer_Event_Vertex_Z(0),
96  fBuffer_Event_Multiplicity(0),
97  fBuffer_Event_NumActiveCells(0),
98  // fBuffer_ClusterNLM_ID(0),
99  // fBuffer_ClusterNLM_E(0),
100  fBuffer_Cells_ID(0),
101  fBuffer_Cells_E(0),
102  fBuffer_Cells_RelativeEta(0),
103  fBuffer_Cells_RelativePhi(0),
104  fBuffer_Surrounding_NCells(0),
105  fBuffer_Surrounding_Cells_ID(0),
106  fBuffer_Surrounding_Cells_R(0),
107  fBuffer_Surrounding_Cells_E(0),
108  fBuffer_Surrounding_Cells_RelativeEta(0),
109  fBuffer_Surrounding_Cells_RelativePhi(0),
110  fBuffer_Surrounding_NTracks(10),
111  fBuffer_Surrounding_Tracks_R(0),
112  fBuffer_Surrounding_Tracks_Pt(0),
113  fBuffer_Surrounding_Tracks_P(0),
114  fBuffer_Surrounding_Tracks_nSigdEdxE(0),
115  fBuffer_Surrounding_Tracks_RelativeEta(0),
116  fBuffer_Surrounding_Tracks_RelativePhi(0),
117  fBuffer_Surrounding_Tracks_V0Flag(0),
118  fBuffer_Cluster_MC_Label(-10),
119  fBuffer_Mother_MC_Label(-10),
120  fBuffer_Cluster_MC_EFracFirstLabel(-10),
121  fBuffer_Cluster_MC_EFracLeadingPi0(-10),
122  fBuffer_Cluster_MC_LeadingPi0_Pt(-10),
123  fBuffer_Cluster_MC_LeadingPi0_E(-10),
124  fVBuffer_Cluster_E(0),
125  fVBuffer_Cluster_Eta(0),
126  fVBuffer_Cluster_Phi(0),
127  fVBuffer_Cluster_isEMCal(0),
128  fVTrueNeutralPionDaughterIndex(0)
129 {
130  // fBuffer_ClusterNLM_ID = new Int_t[kMaxActiveCells];
131  // fBuffer_ClusterNLM_E = new Float_t[kMaxActiveCells];
132  fBuffer_Cells_ID = new Int_t[kMaxActiveCells];
133  fBuffer_Cells_E = new Float_t[kMaxActiveCells];
134  fBuffer_Cells_RelativeEta = new Float_t[kMaxActiveCells];
135  fBuffer_Cells_RelativePhi = new Float_t[kMaxActiveCells];
136  fBuffer_Surrounding_Cells_ID = new Int_t[kMaxActiveCells];
137  fBuffer_Surrounding_Cells_R = new Float_t[kMaxActiveCells];
138  fBuffer_Surrounding_Cells_E = new Float_t[kMaxActiveCells];
139  fBuffer_Surrounding_Cells_RelativeEta = new Float_t[kMaxActiveCells];
140  fBuffer_Surrounding_Cells_RelativePhi = new Float_t[kMaxActiveCells];
141  fBuffer_Surrounding_Tracks_R = new Float_t[kMaxNTracks];
142  fBuffer_Surrounding_Tracks_Pt = new Float_t[kMaxNTracks];
143  fBuffer_Surrounding_Tracks_P = new Float_t[kMaxNTracks];
144  fBuffer_Surrounding_Tracks_nSigdEdxE = new Float_t[kMaxNTracks];
145  fBuffer_Surrounding_Tracks_RelativeEta= new Float_t[kMaxNTracks];
146  fBuffer_Surrounding_Tracks_RelativePhi= new Float_t[kMaxNTracks];
147  fBuffer_Surrounding_Tracks_V0Flag= new Bool_t[kMaxNTracks];
148 }
149 
151  fV0Reader(NULL),
152  fV0ReaderName("V0ReaderV1"),
153  fReaderGammas(NULL),
154  fPIDResponse(NULL),
155  fCorrTaskSetting(""),
156  fConversionCuts(NULL),
157  fEventCuts(NULL),
158  fClusterCutsEMC(NULL),
159  // fMesonCuts(NULL),
160  fMinNLMCut(1),
161  fMaxNLMCut(1),
162  fInputEvent(NULL),
163  fMCEvent(NULL),
164  fWeightJetJetMC(1),
165  fGeomEMCAL(NULL),
166  fClusterTree(NULL),
167  fIsHeavyIon(kFALSE),
168  ffillTree(-100),
169  ffillHistograms(kFALSE),
170  fOutputList(NULL),
171  fIsMC(0),
172  fCorrectForNonlinearity(kFALSE),
173  fSaveEventProperties(0),
174  fSaveCells(0),
175  fSaveSurroundingCells(0),
176  fSaveTracks(0),
177  fConeRadius(0),
178  fMinTrackPt(0),
179  fMinClusterEnergy(0),
180  fSaveMCInformation(0),
181  fSaveAdditionalHistos(0),
182  fSaveEventsInVector(0),
183  fExtractionPercentages(),
184  fExtractionPercentagePtBins(),
185  fBuffer_EventWeight(1),
186  fBuffer_ClusterE(0),
187  fBuffer_ClusterPhi(0),
188  fBuffer_ClusterEta(0),
189  fBuffer_ClusterIsEMCAL(kFALSE),
190  fBuffer_ClusterSupMod(-1),
191  fBuffer_MC_Cluster_Flag(0),
192  fBuffer_ClusterNumCells(0),
193  // fBuffer_ClusterNLM(0),
194  fBuffer_LeadingCell_ID(0),
195  fBuffer_LeadingCell_E(0),
196  fBuffer_LeadingCell_Eta(0),
197  fBuffer_LeadingCell_Phi(0),
198  fBuffer_ClusterM02(0),
199  fBuffer_ClusterM20(0),
200  fBuffer_Event_Vertex_X(0),
201  fBuffer_Event_Vertex_Y(0),
202  fBuffer_Event_Vertex_Z(0),
203  fBuffer_Event_Multiplicity(0),
204  fBuffer_Event_NumActiveCells(0),
205  // fBuffer_ClusterNLM_ID(0),
206  // fBuffer_ClusterNLM_E(0),
207  fBuffer_Cells_ID(0),
208  fBuffer_Cells_E(0),
209  fBuffer_Cells_RelativeEta(0),
210  fBuffer_Cells_RelativePhi(0),
211  fBuffer_Surrounding_NCells(0),
212  fBuffer_Surrounding_Cells_ID(0),
213  fBuffer_Surrounding_Cells_R(0),
214  fBuffer_Surrounding_Cells_E(0),
215  fBuffer_Surrounding_Cells_RelativeEta(0),
216  fBuffer_Surrounding_Cells_RelativePhi(0),
217  fBuffer_Surrounding_NTracks(10),
218  fBuffer_Surrounding_Tracks_R(0),
219  fBuffer_Surrounding_Tracks_Pt(0),
220  fBuffer_Surrounding_Tracks_P(0),
221  fBuffer_Surrounding_Tracks_nSigdEdxE(0),
222  fBuffer_Surrounding_Tracks_RelativeEta(0),
223  fBuffer_Surrounding_Tracks_RelativePhi(0),
224  fBuffer_Surrounding_Tracks_V0Flag(0),
225  fBuffer_Cluster_MC_Label(-10),
226  fBuffer_Mother_MC_Label(-10),
227  fBuffer_Cluster_MC_EFracFirstLabel(-10),
228  fBuffer_Cluster_MC_EFracLeadingPi0(-10),
229  fBuffer_Cluster_MC_LeadingPi0_Pt(-10),
230  fBuffer_Cluster_MC_LeadingPi0_E(-10),
231  fVBuffer_Cluster_E(0),
232  fVBuffer_Cluster_Eta(0),
233  fVBuffer_Cluster_Phi(0),
234  fVBuffer_Cluster_isEMCal(0),
235  fVTrueNeutralPionDaughterIndex(0)
236 {
237  // fBuffer_ClusterNLM_ID = new Int_t[kMaxActiveCells];
238  // fBuffer_ClusterNLM_E = new Float_t[kMaxActiveCells];
255  // Default constructor
256 
257  DefineInput(0, TChain::Class());
258  DefineOutput(1, TList::Class());
259  DefineOutput(2, TTree::Class());
260 }
261 
262 //________________________________________________________________________
264 {
265  // default deconstructor
266 
267 }
268 //________________________________________________________________________
270 {
271  // Create User Output Objects
272 
273  fV0Reader = (AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data());
274  if(!fV0Reader){printf("Error: No V0 Reader");return;}// GetV0Reader
275 
276  if(fOutputList != NULL){
277  delete fOutputList;
278  fOutputList = NULL;
279  }
280  if(fOutputList == NULL){
281  fOutputList = new TList();
282  fOutputList->SetOwner(kTRUE);
283  }
284 
285  if(ffillHistograms){
286 
287  if(((AliConvEventCuts*)fEventCuts)->GetCutHistograms()){
288  fOutputList->Add(((AliConvEventCuts*)fEventCuts)->GetCutHistograms());
289  }
290 
291  if(((AliCaloPhotonCuts*)fClusterCutsEMC)->GetCutHistograms()){
292  fOutputList->Add(((AliCaloPhotonCuts*)fClusterCutsEMC)->GetCutHistograms());
293  }
294  // if(((AliConversionMesonCuts*)fMesonCuts)->GetCutHistograms()){
295  // fOutputList->Add(((AliConversionMesonCuts*)fMesonCuts)->GetCutHistograms());
296  // }
297  PostData(1, fOutputList);
298  }
299  if(!fCorrTaskSetting.CompareTo("")){
300  fClusterTree = new TTree(Form("ClusterQA_%s_%s",(fEventCuts->GetCutNumber()).Data(),(fClusterCutsEMC->GetCutNumber()).Data()),Form("ClusterQA_%s_%s",(fEventCuts->GetCutNumber()).Data(),(fClusterCutsEMC->GetCutNumber()).Data()));
301  } else {
302  fClusterTree = new TTree(Form("ClusterQA_%s_%s_%s",(fEventCuts->GetCutNumber()).Data(),(fClusterCutsEMC->GetCutNumber()).Data(),fCorrTaskSetting.Data()),Form("ClusterQA_%s_%s_%s",(fEventCuts->GetCutNumber()).Data(),(fClusterCutsEMC->GetCutNumber()).Data(),fCorrTaskSetting.Data()));
303  }
305  {
306  fClusterTree->Branch("Cluster_E", &fBuffer_ClusterE, "Cluster_E/F");
307  fClusterTree->Branch("Cluster_Eta", &fBuffer_ClusterEta, "Cluster_Eta/F");
308  fClusterTree->Branch("Cluster_Phi", &fBuffer_ClusterPhi, "Cluster_Phi/F");
309  fClusterTree->Branch("Cluster_IsEMCAL", &fBuffer_ClusterIsEMCAL, "Cluster_IsEMCAL/O");
310  fClusterTree->Branch("Cluster_SM", &fBuffer_ClusterSupMod, "Cluster_SM/I");
311  fClusterTree->Branch("Cluster_NumCells", &fBuffer_ClusterNumCells, "Cluster_NumCells/I");
312  // fClusterTree->Branch("Cluster_NLM", &fBuffer_ClusterNLM, "Cluster_NLM/I");
313  // fClusterTree->Branch("Cluster_NLM_ID", fBuffer_ClusterNLM_ID, "Cluster_NLM_ID[Cluster_NLM]/I");
314  // fClusterTree->Branch("Cluster_NLM_E", fBuffer_ClusterNLM_E, "Cluster_NLM_E[Cluster_NLM]/F");
315  fClusterTree->Branch("Cluster_LeadingCell_ID", &fBuffer_LeadingCell_ID, "Cluster_LeadingCell_ID/I");
316  fClusterTree->Branch("Cluster_LeadingCell_E", &fBuffer_LeadingCell_E, "Cluster_LeadingCell_E/F");
317  fClusterTree->Branch("Cluster_LeadingCell_Eta", &fBuffer_LeadingCell_Eta, "Cluster_LeadingCell_Eta/F");
318  fClusterTree->Branch("Cluster_LeadingCell_Phi", &fBuffer_LeadingCell_Phi, "Cluster_LeadingCell_Phi/F");
319  fClusterTree->Branch("Cluster_M02", &fBuffer_ClusterM02, "Cluster_M02/F");
320  fClusterTree->Branch("Cluster_M20", &fBuffer_ClusterM20, "Cluster_M20/F");
321 
322  fClusterTree->Branch("Event_Weight", &fBuffer_EventWeight, "Event_Weight/F");
323  }
324  else
325  {
326  fClusterTree->Branch("Cluster_E", &fVBuffer_Cluster_E);
327  fClusterTree->Branch("Cluster_Eta", &fVBuffer_Cluster_Eta);
328  fClusterTree->Branch("Cluster_Phi", &fVBuffer_Cluster_Phi);
329  fClusterTree->Branch("Cluster_isEMCal", &fVBuffer_Cluster_isEMCal);
330  if(fIsMC) fClusterTree->Branch("Cluster_MotherPi0", &fVTrueNeutralPionDaughterIndex);
331  }
332 
334  {
335  fClusterTree->Branch("Event_Vertex_X", &fBuffer_Event_Vertex_X, "Event_Vertex_X/F");
336  fClusterTree->Branch("Event_Vertex_Y", &fBuffer_Event_Vertex_Y, "Event_Vertex_Y/F");
337  fClusterTree->Branch("Event_Vertex_Z", &fBuffer_Event_Vertex_Z, "Event_Vertex_Z/F");
338  fClusterTree->Branch("Event_Multiplicity", &fBuffer_Event_Multiplicity, "Event_Multiplicity/F");
339  fClusterTree->Branch("Event_NumActiveCells", &fBuffer_Event_NumActiveCells, "Event_NumActiveCells/I");
340  }
341 
343  {
344  fClusterTree->Branch("Cluster_Cells_ID", fBuffer_Cells_ID, "Cluster_Cells_ID[Cluster_NumCells]/I");
345  fClusterTree->Branch("Cluster_Cells_E", fBuffer_Cells_E, "Cluster_Cells_E[Cluster_NumCells]/F");
346  fClusterTree->Branch("Cluster_Cells_RelativeEta", fBuffer_Cells_RelativeEta, "Cluster_Cells_RelativeEta[Cluster_NumCells]/F");
347  fClusterTree->Branch("Cluster_Cells_RelativePhi", fBuffer_Cells_RelativePhi, "Cluster_Cells_RelativePhi[Cluster_NumCells]/F");
348  }
349 
351  {
352  fClusterTree->Branch("Surrounding_NCells", &fBuffer_Surrounding_NCells, "Surrounding_NCells/I");
353  fClusterTree->Branch("Surrounding_Cells_ID", fBuffer_Surrounding_Cells_ID, "Surrounding_Cells_ID[Surrounding_NCells]/I");
354  fClusterTree->Branch("Surrounding_Cells_R", fBuffer_Surrounding_Cells_R, "Surrounding_Cells_R[Surrounding_NCells]/F");
355  fClusterTree->Branch("Surrounding_Cells_E", fBuffer_Surrounding_Cells_E, "Surrounding_Cells_E[Surrounding_NCells]/F");
356  fClusterTree->Branch("Surrounding_Cells_RelativeEta", fBuffer_Surrounding_Cells_RelativeEta, "Surrounding_Cells_RelativeEta[Surrounding_NCells]/F");
357  fClusterTree->Branch("Surrounding_Cells_RelativePhi", fBuffer_Surrounding_Cells_RelativePhi, "Surrounding_Cells_RelativePhi[Surrounding_NCells]/F");
358  }
360  {
361  fClusterTree->Branch("Surrounding_NTracks", &fBuffer_Surrounding_NTracks, "Surrounding_NTracks/I");
362  fClusterTree->Branch("Surrounding_Tracks_R", fBuffer_Surrounding_Tracks_R, "Surrounding_Tracks_R[Surrounding_NTracks]/F");
363  fClusterTree->Branch("Surrounding_Tracks_Pt", fBuffer_Surrounding_Tracks_Pt, "Surrounding_Tracks_Pt[Surrounding_NTracks]/F");
364  fClusterTree->Branch("Surrounding_Tracks_P", fBuffer_Surrounding_Tracks_P, "Surrounding_Tracks_P[Surrounding_NTracks]/F");
365  fClusterTree->Branch("Surrounding_Tracks_nSigdEdxE", fBuffer_Surrounding_Tracks_nSigdEdxE, "Surrounding_Tracks_nSigdEdxE[Surrounding_NTracks]/F");
366  fClusterTree->Branch("Surrounding_Tracks_RelativeEta", fBuffer_Surrounding_Tracks_RelativeEta, "Surrounding_Tracks_RelativeEta[Surrounding_NTracks]/F");
367  fClusterTree->Branch("Surrounding_Tracks_RelativePhi", fBuffer_Surrounding_Tracks_RelativePhi, "Surrounding_Tracks_RelativePhi[Surrounding_NTracks]/F");
368  fClusterTree->Branch("Surrounding_Tracks_V0Flag", fBuffer_Surrounding_Tracks_V0Flag, "Surrounding_Tracks_V0Flag[Surrounding_NTracks]/O");
369  }
370 
372  {
373  fClusterTree->Branch("Cluster_MC_Label", &fBuffer_Cluster_MC_Label, "Cluster_MC_Label/I");
374  fClusterTree->Branch("Mother_MC_Label", &fBuffer_Mother_MC_Label, "Mother_MC_Label/I");
375  fClusterTree->Branch("Cluster_MC_EFracFirstLabel", &fBuffer_Cluster_MC_EFracFirstLabel, "Cluster_MC_EFracFirstLabel/F");
376  fClusterTree->Branch("Cluster_MC_EFracLeadingPi0", &fBuffer_Cluster_MC_EFracLeadingPi0, "Cluster_MC_EFracLeadingPi0/F");
377  fClusterTree->Branch("Cluster_MC_LeadingPi0_Pt", &fBuffer_Cluster_MC_LeadingPi0_Pt, "Cluster_MC_LeadingPi0_Pt/F");
378  fClusterTree->Branch("Cluster_MC_LeadingPi0_E", &fBuffer_Cluster_MC_LeadingPi0_E, "Cluster_MC_LeadingPi0_E/F");
379  }
380 
381 
382 
383  fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data());
384  OpenFile(2);
385  PostData(2, fClusterTree);
386 }
387 //_____________________________________________________________________________
389 {
390  if (fEventCuts->GetPeriodEnum() == AliConvEventCuts::kNoPeriod && ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetPeriodEnum() != AliConvEventCuts::kNoPeriod){
391  fEventCuts->SetPeriodEnumExplicit(((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetPeriodEnum());
392  } else if (fEventCuts->GetPeriodEnum() == AliConvEventCuts::kNoPeriod ){
393  fEventCuts->SetPeriodEnum(fV0Reader->GetPeriodName());
394  }
395 
396  return kTRUE;
397 }
398 //________________________________________________________________________
400 
401  Int_t eventQuality = ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEventQuality();
402  if(eventQuality != 0){// Event Not Accepted
403  return;
404  }
405 
406  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
407  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
408  fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
409 
410  fInputEvent = InputEvent();
411  if(fIsMC>0) fMCEvent = MCEvent();
412 
413  Int_t eventNotAccepted = fEventCuts->IsEventAcceptedByCut(fV0Reader->GetEventCuts(),fInputEvent,fMCEvent,fIsHeavyIon,kFALSE);
414  if(eventNotAccepted) return; // Check Centrality, PileUp, SDD and V0AND --> Not Accepted => eventQuality = 1
415 
416  fGeomEMCAL = AliEMCALGeometry::GetInstance();
417  if(!fGeomEMCAL){ AliFatal("EMCal geometry not initialized!");}
418 
419  Int_t nclus = 0;
420  TClonesArray * arrClustersProcess = NULL;
421  // fNCurrentClusterBasic = 0;
422  if(!fCorrTaskSetting.CompareTo("")){
423  nclus = fInputEvent->GetNumberOfCaloClusters();
424  } else {
425  arrClustersProcess = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(Form("%sClustersBranch",fCorrTaskSetting.Data())));
426  if(!arrClustersProcess)
427  AliFatal(Form("%sClustersBranch was not found in AliAnalysisTaskGammaCalo! Check the correction framework settings!",fCorrTaskSetting.Data()));
428  nclus = arrClustersProcess->GetEntries();
429  }
430 
431  if(nclus == 0) return;
432 
433  ((AliCaloPhotonCuts*)fClusterCutsEMC)->FillHistogramsExtendedQA(fInputEvent,fIsMC);
434 
435  fReaderGammas = fV0Reader->GetReconstructedGammas(); // Gammas from default Cut
436 
437  // if(fIsMC==2){
438  // Float_t xsection = -1.; Float_t ntrials = -1.;
439  // ((AliConvEventCuts*)fEventCuts)->GetXSectionAndNTrials(fMCEvent,xsection,ntrials,fInputEvent);
440  // if((xsection==-1.) || (ntrials==-1.)) AliFatal(Form("ERROR: GetXSectionAndNTrials returned invalid xsection/ntrials, periodName from V0Reader: '%s'",fV0Reader->GetPeriodName().Data()));
441  // // fProfileJetJetXSection[iCut]->Fill(0.,xsection);
442  // // fHistoJetJetNTrials[iCut]->Fill("#sum{NTrials}",ntrials);
443  // }
444 
445  fWeightJetJetMC = 1;
446  Float_t pthard = -1;
447  Bool_t isMCJet = ((AliConvEventCuts*)fEventCuts)->IsJetJetMCEventAccepted( fMCEvent, fWeightJetJetMC,pthard, fInputEvent );
448  if (!isMCJet){
449  return;
450  }
452  // ((AliCaloPhotonCuts*)fClusterCutsDMC)->FillHistogramsExtendedQA(fInputEvent,fIsMC);
453 
454  // match tracks to clusters
455  // ((AliCaloPhotonCuts*)fClusterCutsEMC)->MatchTracksToClusters(fInputEvent,fWeightJetJetMC,kTRUE, fMCEvent);
456  // ((AliCaloPhotonCuts*)fClusterCutsDMC)->MatchTracksToClusters(fInputEvent,fWeightJetJetMC,kTRUE, fMCEvent);
457 
458  // AliVCaloCells* cells = fInputEvent->GetEMCALCells();
459 
460  // vertex
461  Double_t vertex[3] = {0};
462  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
463 
465 
466  // map<Long_t,Int_t> mapIsClusterAccepted;
467  // map<Long_t,Int_t> mapIsClusterAcceptedWithoutTrackMatch;
468  // Int_t nCellInCluster = 0;
469  // Loop over EMCal clusters
470  for(Long_t i = 0; i < nclus; i++){
471  Double_t tempClusterWeight = fWeightJetJetMC;
472  AliVCluster* clus = NULL;
473  if(fInputEvent->IsA()==AliESDEvent::Class()){
474  if(arrClustersProcess)
475  clus = new AliESDCaloCluster(*(AliESDCaloCluster*)arrClustersProcess->At(i));
476  else
477  clus = new AliESDCaloCluster(*(AliESDCaloCluster*)fInputEvent->GetCaloCluster(i));
478  } else if(fInputEvent->IsA()==AliAODEvent::Class()){
479  if(arrClustersProcess)
480  clus = new AliAODCaloCluster(*(AliAODCaloCluster*)arrClustersProcess->At(i));
481  else
482  clus = new AliAODCaloCluster(*(AliAODCaloCluster*)fInputEvent->GetCaloCluster(i));
483  }
484  if(!clus) continue;
485 
486 
487  if ( !clus->IsEMCAL()){ // for PHOS: cluster->GetType() == AliVCluster::kPHOSNeutral
488  delete clus;
489  continue;
490  }
491 
492  if ( clus->E() < fMinClusterEnergy){
493  delete clus;
494  continue;
495  }
496 
497  if(!((AliCaloPhotonCuts*)fClusterCutsEMC)->ClusterIsSelected(clus,fInputEvent,fMCEvent,fIsMC, tempClusterWeight,i)){
498  delete clus;
499  continue;
500  }
501  ResetBuffer();
503  }
504  if(fSaveEventsInVector) fClusterTree->Fill();
505 
506  PostData(1, fOutputList);
507 }
508 
510 void AliAnalysisTaskClusterQA::ProcessQATreeCluster(AliVEvent *event, AliVCluster* cluster, Long_t indexCluster){
511  Float_t clusPos[3] = { 0,0,0 };
512  cluster->GetPosition(clusPos);
513  TVector3 clusterVector(clusPos[0],clusPos[1],clusPos[2]);
514  Double_t etaCluster = clusterVector.Eta();
515  Double_t phiCluster = clusterVector.Phi();
516  if (phiCluster < 0) phiCluster += 2*TMath::Pi();
517 
519  // Vertex position x, y, z
520  fBuffer_Event_Vertex_X = fInputEvent->GetPrimaryVertex()->GetX();
521  fBuffer_Event_Vertex_Y = fInputEvent->GetPrimaryVertex()->GetY();
522  fBuffer_Event_Vertex_Z = fInputEvent->GetPrimaryVertex()->GetZ();
523 
524  // V0-based multiplicity of the event
525  fBuffer_Event_Multiplicity = fInputEvent->GetVZEROData()->GetMTotV0A()+fInputEvent->GetVZEROData()->GetMTotV0C();
526  }
528 
529  // Properties of the current cluster
530  fBuffer_ClusterE = cluster->E();
531 
532  if(etaCluster < 0.67 && etaCluster > -0.67){
533  if(phiCluster < 3.2 && phiCluster > 1.4){
534  fBuffer_ClusterIsEMCAL = kTRUE;
535  }
536  }
537  if(etaCluster < 0.67 && etaCluster > -0.67){
538  if(etaCluster > 0.23 || etaCluster < -0.23){
539  if(phiCluster < 5.6 && phiCluster > 4.6){
540  fBuffer_ClusterIsEMCAL = kFALSE;
541  }
542  }
543  }
544 
545  fBuffer_ClusterPhi = phiCluster;
546  fBuffer_ClusterEta = etaCluster;
547  fBuffer_ClusterM02 = cluster->GetM02();
548  fBuffer_ClusterM20 = cluster->GetM20();
549 
550  // Get all cells from the event
551  AliVCaloCells* cells = NULL;
552  if(cluster->IsEMCAL())
553  cells = event->GetEMCALCells();
554  else
555  return;
556 
557  // Determine all active cells in the event
559  fBuffer_Event_NumActiveCells = cells->GetNumberOfCells();
560 
561  // Get the number of cells from the current cluster
562  Int_t nCellCluster = cluster->GetNCells();
563  fBuffer_ClusterNumCells = nCellCluster;
564 
565  // Int_t nLM = GetNumberOfLocalMaxima(cluster, event);
566  // fBuffer_ClusterNLM = nLM;
567  // Find the leading cell in the cluster and its position
569 
570  Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0;
571  // Get SM number for the leading cell in the cluster
572  fGeomEMCAL->GetCellIndex(fBuffer_LeadingCell_ID, nSupMod,nModule,nIphi,nIeta);
573  fBuffer_ClusterSupMod = nSupMod;
574 
575  fBuffer_LeadingCell_E = cells->GetCellAmplitude(fBuffer_LeadingCell_ID);
576  Float_t leadcelleta = 0.;
577  Float_t leadcellphi = 0.;
578  fGeomEMCAL->EtaPhiFromIndex(fBuffer_LeadingCell_ID, leadcelleta, leadcellphi);
579  if ( leadcellphi < 0 ) leadcellphi+=TMath::TwoPi();
580  fBuffer_LeadingCell_Eta = leadcelleta;
581  fBuffer_LeadingCell_Phi = leadcellphi;
582  // Treat the remaining cells of the cluster and get their relative position compared to the leading cell
583  if(fSaveCells){
584  for(Int_t iCell=0;iCell<nCellCluster;iCell++){
585  if(iCell<100){ // maximum number of cells for a cluster is set to 100
586  fBuffer_Cells_ID[iCell] = cluster->GetCellAbsId(iCell);
587  fBuffer_Cells_E[iCell] = cells->GetCellAmplitude(cluster->GetCellAbsId(iCell));
588  Float_t celleta = 0.;
589  Float_t cellphi = 0.;
590  fGeomEMCAL->EtaPhiFromIndex(fBuffer_Cells_ID[iCell], celleta, cellphi);
591  if ( cellphi < 0 ) cellphi+=TMath::TwoPi();
592  fBuffer_Cells_RelativeEta[iCell] = leadcelleta-celleta;
593  fBuffer_Cells_RelativePhi[iCell] = leadcellphi-cellphi;
594  }
595  }
596  }
598  Int_t nActiveCellsSurroundingInR = 0;
599  // fill surrounding cell buffer
600  for(Int_t aCell=0;aCell<cells->GetNumberOfCells();aCell++){
601  // Define necessary variables
602  Short_t cellNumber = 0;
603  Double_t cellAmplitude = 0, cellTime = 0, cellEFrac = 0;
604  Int_t cellMCLabel = 0;
605  Float_t surrcelleta = 0.;
606  Float_t surrcellphi = 0.;
607  // Get Cell ID
608  cells->GetCell(aCell,cellNumber,cellAmplitude,cellTime,cellMCLabel,cellEFrac);
609 
610  // Get eta and phi for the surounding cells
611  fGeomEMCAL->EtaPhiFromIndex(cellNumber, surrcelleta, surrcellphi);
612  if ( surrcellphi < 0 ) surrcellphi+=TMath::TwoPi();
613  Double_t dR2 = pow(leadcelleta-surrcelleta,2) + pow(leadcellphi-surrcellphi,2);
614  // Select those cells that are within fConeRadius of the leading cluster cell
615  if( dR2 < fConeRadius){
616  fBuffer_Surrounding_Cells_E[nActiveCellsSurroundingInR] = cells->GetCellAmplitude(cellNumber);
617  fBuffer_Surrounding_Cells_ID[nActiveCellsSurroundingInR] = cellNumber;
618  fBuffer_Surrounding_Cells_R[nActiveCellsSurroundingInR] = dR2;
619  fBuffer_Surrounding_Cells_RelativeEta[nActiveCellsSurroundingInR] = leadcelleta-surrcelleta;
620  fBuffer_Surrounding_Cells_RelativePhi[nActiveCellsSurroundingInR] = leadcellphi-surrcellphi;
621  nActiveCellsSurroundingInR+=1;
622  }
623  }
624  fBuffer_Surrounding_NCells = nActiveCellsSurroundingInR;
625  }
626 
627  // write PDG code of mother of particle that mainly contributed to cluster to tree
628  if(fIsMC> 0){
629  if (cluster->GetNLabels()>0){
630  if((cluster->GetLabelAt(0)!=-1)){
631  if(fInputEvent->IsA()==AliESDEvent::Class()){
632  fBuffer_Mother_MC_Label = (fMCEvent->Particle(cluster->GetLabelAt(0)))->GetPdgCode(); // mother of leading contribution
633  } else if(fInputEvent->IsA()==AliAODEvent::Class()){
634  TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
635  if (AODMCTrackArray == NULL) return;
636 
637  AliAODMCParticle* particle = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(cluster->GetLabelAt(0)));
638  fBuffer_Mother_MC_Label = particle->GetPdgCode();
639  }
640 
641  } else{
642  // mother is initial collision
644  }
645  } else{
646  // no mothers found (should not happen)
648  }
649  }
650  if(fIsMC) fBuffer_Cluster_MC_Label = MakePhotonCandidates(cluster, cells,indexCluster);
651  if(fSaveTracks) ProcessTracksAndMatching(cluster,indexCluster);
652  // Add everything to the tree
655  fVBuffer_Cluster_E.push_back(cluster->E());
656  fVBuffer_Cluster_Eta.push_back(etaCluster);
657  fVBuffer_Cluster_Phi.push_back(phiCluster);
659  }
660 
661 }
662 
663 //________________________________________________________________________
665 
666 }
667 
668 //________________________________________________________________________
669 Int_t AliAnalysisTaskClusterQA::FindLargestCellInCluster(AliVCluster* cluster, AliVEvent* event){
670 
671  const Int_t nCells = cluster->GetNCells();
672  AliVCaloCells* cells = NULL;
673 
674  if(cluster->IsEMCAL())
675  cells = event->GetEMCALCells();
676  else if(cluster->IsPHOS())
677  cells = event->GetPHOSCells();
678 
679  Float_t eMax = 0.;
680  Int_t idMax = -1;
681 
682  if (nCells < 1) return idMax;
683  for (Int_t iCell = 0;iCell < nCells;iCell++){
684  Int_t cellAbsID = cluster->GetCellsAbsId()[iCell];
685  if (cells->GetCellAmplitude(cellAbsID)> eMax){
686  eMax = cells->GetCellAmplitude(cellAbsID);
687  idMax = cellAbsID;
688  }
689  }
690  return idMax;
691 }
692 
693 //________________________________________________________________________
695  Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0;
696  row=0;
697  column=0;
698  // Get SM number and relative row/column for SM
699  fGeomEMCAL->GetCellIndex(cellIndex, nSupMod,nModule,nIphi,nIeta);
700  fGeomEMCAL->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, row,column);
701  row += nSupMod/2 * 24;
702  column += nSupMod%2 * 48;
703 }
704 
705 //________________________________________________________________________
706 Int_t AliAnalysisTaskClusterQA::MakePhotonCandidates(AliVCluster* clus, AliVCaloCells* cells, Long_t indexCluster){
707 
708  Double_t vertex[3] = {0};
709  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
710 
711  TLorentzVector clusterVector;
712  clus->GetMomentum(clusterVector,vertex);
713 
714  TLorentzVector* tmpvec = new TLorentzVector();
715  tmpvec->SetPxPyPzE(clusterVector.Px(),clusterVector.Py(),clusterVector.Pz(),clusterVector.E());
716 
717  // convert to AODConversionPhoton
718  AliAODConversionPhoton *PhotonCandidate=new AliAODConversionPhoton(tmpvec);
719  if(!PhotonCandidate){
720  delete clus;
721  delete tmpvec;
722  if (PhotonCandidate) delete PhotonCandidate;
723  return -9;
724  }
725  // Flag Photon as CaloPhoton
726  PhotonCandidate->SetIsCaloPhoton(((AliCaloPhotonCuts*)fClusterCutsEMC)->GetClusterType());
727  PhotonCandidate->SetCaloClusterRef(indexCluster);
728 
729  // get MC label
730  if(fIsMC> 0){
731  Int_t* mclabelsCluster = clus->GetLabels();
732  Int_t nValidClusters = 0;
733  if (clus->GetNLabels()>0){
734  for (Int_t k =0; k< (Int_t)clus->GetNLabels(); k++){
735  if (mclabelsCluster[k]>0){
736  if (nValidClusters< 50)PhotonCandidate->SetCaloPhotonMCLabel(nValidClusters,mclabelsCluster[k]);
737  nValidClusters++;
738  }
739  }
740  }
741  PhotonCandidate->SetNCaloPhotonMCLabels(nValidClusters);
742  }
743 
744  AliAODCaloCluster* clusSub1 = new AliAODCaloCluster();
745  AliAODCaloCluster* clusSub2 = new AliAODCaloCluster();
746 
747 
748  // split clusters according to their shares in the cluster (NLM == 1) needs to be treated differently
749  if (fMinNLMCut == 1 && fMaxNLMCut == 1 ){
750  Int_t absCellIdFirst = ((AliCaloPhotonCuts*)fClusterCutsEMC)->FindLargestCellInCluster(clus, fInputEvent);
751  Int_t absCellIdSecond = ((AliCaloPhotonCuts*)fClusterCutsEMC)->FindSecondLargestCellInCluster(clus, fInputEvent);
752 
753  ((AliCaloPhotonCuts*)fClusterCutsEMC)->SplitEnergy(absCellIdFirst, absCellIdSecond, clus, fInputEvent, fIsMC, clusSub1, clusSub2);
754  } else if (fMinNLMCut > 1 ){
755  const Int_t nc = clus->GetNCells();
756  Int_t absCellIdList[nc];
757 
758  ((AliCaloPhotonCuts*)fClusterCutsEMC)->SplitEnergy(absCellIdList[0], absCellIdList[1], clus, fInputEvent, fIsMC, clusSub1, clusSub2);
759  }
760 
761  // TLorentzvector with sub cluster 1
762  TLorentzVector clusterVector1;
763  clusSub1->GetMomentum(clusterVector1,vertex);
764  TLorentzVector* tmpvec1 = new TLorentzVector();
765  tmpvec1->SetPxPyPzE(clusterVector1.Px(),clusterVector1.Py(),clusterVector1.Pz(),clusterVector1.E());
766  // convert to AODConversionPhoton
767  AliAODConversionPhoton *PhotonCandidate1=new AliAODConversionPhoton(tmpvec1);
768  if(!PhotonCandidate1){
769  delete clusSub1;
770  delete tmpvec1;
771  return -9;
772  }
773  // Flag Photon as CaloPhoton
774  PhotonCandidate1->SetIsCaloPhoton(((AliCaloPhotonCuts*)fClusterCutsEMC)->GetClusterType());
775  // TLorentzvector with sub cluster 2
776  TLorentzVector clusterVector2;
777  clusSub2->GetMomentum(clusterVector2,vertex);
778  TLorentzVector* tmpvec2 = new TLorentzVector();
779  tmpvec2->SetPxPyPzE(clusterVector2.Px(),clusterVector2.Py(),clusterVector2.Pz(),clusterVector2.E());
780  // convert to AODConversionPhoton
781  AliAODConversionPhoton *PhotonCandidate2=new AliAODConversionPhoton(tmpvec2);
782  if(!PhotonCandidate2){
783  delete clusSub2;
784  delete tmpvec2;
785  return -9;
786  }
787  // Flag Photon as CaloPhoton
788  PhotonCandidate2->SetIsCaloPhoton(((AliCaloPhotonCuts*)fClusterCutsEMC)->GetClusterType());
789  Int_t mclabel = -3;
790  if(fIsMC> 0 && PhotonCandidate && PhotonCandidate1 && PhotonCandidate2 && fSaveMCInformation){
791  if(fInputEvent->IsA()==AliESDEvent::Class())
792  mclabel = ProcessTrueClusterCandidates(PhotonCandidate, clus, PhotonCandidate1, PhotonCandidate2);
793  if(fInputEvent->IsA()==AliAODEvent::Class())
794  mclabel = ProcessTrueClusterCandidatesAOD(PhotonCandidate, clus, PhotonCandidate1, PhotonCandidate2);
795  return mclabel;
796  } else {
798  return -7;
799  }
800  return -1;
801 }
802 
803 //________________________________________________________________________
804 void AliAnalysisTaskClusterQA::ProcessTracksAndMatching(AliVCluster* clus, Long_t indexCluster){
805 
806 
807  Int_t nTracksInR = 0;
808  Int_t nModules = fGeomEMCAL->GetNumberOfSuperModules();
809 
810  AliESDEvent *esdev = dynamic_cast<AliESDEvent*>(fInputEvent);
811  AliAODEvent *aodev = 0;
812  if (!esdev) {
813  aodev = dynamic_cast<AliAODEvent*>(fInputEvent);
814  if (!aodev) {
815  AliError("Task needs AOD or ESD event, returning");
816  return;
817  }
818  }
819  AliESDtrackCuts *EsdTrackCuts = 0x0;
820  // Using standard function for setting Cuts
821  static int prevRun = -1;
822  // Using standard function for setting Cuts
823  if (esdev){
824  Int_t runNumber = fInputEvent->GetRunNumber();
825  if (prevRun!=runNumber) {
826  delete EsdTrackCuts;
827  EsdTrackCuts = 0;
828  prevRun = runNumber;
829  }
830  // if LHC11a or earlier or if LHC13g or if LHC12a-i -> use 2010 cuts
831  if( (runNumber<=146860) || (runNumber>=197470 && runNumber<=197692) || (runNumber>=172440 && runNumber<=193766) ){
832  EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
833  // else if run2 data use 2015 PbPb cuts
834  }else if (runNumber>=209122){
835  // hard coded track cuts for the moment, because AliESDtrackCuts::GetStandardITSTPCTrackCuts2015PbPb() gives spams warnings
836  EsdTrackCuts = new AliESDtrackCuts();
837  EsdTrackCuts->AliESDtrackCuts::SetMinNCrossedRowsTPC(70);
838  EsdTrackCuts->AliESDtrackCuts::SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
839  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2PerClusterTPC(4);
840  EsdTrackCuts->AliESDtrackCuts::SetAcceptKinkDaughters(kFALSE);
841  EsdTrackCuts->AliESDtrackCuts::SetRequireTPCRefit(kTRUE);
842  // ITS
843  EsdTrackCuts->AliESDtrackCuts::SetRequireITSRefit(kTRUE);
844  EsdTrackCuts->AliESDtrackCuts::SetClusterRequirementITS(AliESDtrackCuts::kSPD,
845  AliESDtrackCuts::kAny);
846  EsdTrackCuts->AliESDtrackCuts::SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.1");
847  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2TPCConstrainedGlobal(36);
848  EsdTrackCuts->AliESDtrackCuts::SetMaxDCAToVertexZ(2);
849  EsdTrackCuts->AliESDtrackCuts::SetDCAToVertex2D(kFALSE);
850  EsdTrackCuts->AliESDtrackCuts::SetRequireSigmaToVertex(kFALSE);
851  EsdTrackCuts->AliESDtrackCuts::SetMaxChi2PerClusterITS(36);
852  // else use 2011 version of track cuts
853  }else{
854  EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
855  }
856  EsdTrackCuts->SetMaxDCAToVertexZ(2);
857  EsdTrackCuts->SetEtaRange(-0.8, 0.8);
858  EsdTrackCuts->SetPtRange(0.15);
859  }
860 
861  for (Int_t itr=0;itr<fInputEvent->GetNumberOfTracks();itr++){
862  AliExternalTrackParam *trackParam = 0;
863  AliVTrack *inTrack = 0x0;
864  if(esdev){
865  inTrack = esdev->GetTrack(itr);
866  if(!inTrack) continue;
867  if(inTrack->Pt()<fMinTrackPt) continue;
868  AliESDtrack *esdt = dynamic_cast<AliESDtrack*>(inTrack);
869  // check track cuts
870  if(!EsdTrackCuts->AcceptTrack(esdt)) continue;
871  const AliExternalTrackParam *in = esdt->GetInnerParam();
872  if (!in){AliError("Could not get InnerParam of Track, continue");continue;}
873  trackParam =new AliExternalTrackParam(*in);
874  } else if(aodev){
875  inTrack = dynamic_cast<AliVTrack*>(aodev->GetTrack(itr));
876  if(!inTrack) continue;
877  if(inTrack->Pt()<fMinTrackPt) continue;
878  AliAODTrack *aodt = dynamic_cast<AliAODTrack*>(inTrack);
879  // check track cuts
880  if(!aodt->IsHybridGlobalConstrainedGlobal()) continue;
881 
882  Double_t xyz[3] = {0}, pxpypz[3] = {0}, cv[21] = {0};
883  aodt->GetPxPyPz(pxpypz);
884  aodt->GetXYZ(xyz);
885  aodt->GetCovarianceXYZPxPyPz(cv);
886 
887  trackParam = new AliExternalTrackParam(xyz,pxpypz,cv,aodt->Charge());
888  }
889  AliExternalTrackParam emcParam(*trackParam);
890  Float_t eta, phi, pt;
891 
892  //propagate tracks to emc surfaces
893  if (!AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&emcParam, 440., 0.139, 20., eta, phi, pt)) {
894  delete trackParam;
895  continue;
896  }
897  if( TMath::Abs(eta) > 0.75 ) {
898  delete trackParam;
899  continue;
900  }
901  // Save some time and memory in case of no DCal present
902  if( nModules < 13 && ( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad())){
903  delete trackParam;
904  continue;
905  }
906  // Save some time and memory in case of run2
907  if( nModules > 12 ){
908  if (( phi < 70*TMath::DegToRad() || phi > 190*TMath::DegToRad()) && ( phi < 250*TMath::DegToRad() || phi > 340*TMath::DegToRad())){
909  delete trackParam;
910  continue;
911  }
912  }
913  Float_t dEta=-999, dPhi=-999;
914  Double_t trkPos[3] = {0.,0.,0.};
915  if (!emcParam.GetXYZ(trkPos)){
916  delete trackParam;
917  continue;
918  }
919 
920  Bool_t trackIsFromV0 = kFALSE;
921  for(Int_t i = 0; i < fReaderGammas->GetEntriesFast(); i++){
922  AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(i);
923  if(!PhotonCandidate) continue;
924  //apply cuts to maximize electron purity
925  // if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelected(PhotonCandidate,fInputEvent))continue;
926 
927  for (Int_t iElec = 0;iElec < 2;iElec++){
928  Int_t tracklabel = PhotonCandidate->GetLabel(iElec);
929  if(tracklabel==itr){
930  trackIsFromV0 = kTRUE;
931  } else {
932  trackIsFromV0 = kFALSE;
933  }
934  }
935 
936  }
937 
938 
939  AliExternalTrackParam trackParamTmp(emcParam);//Retrieve the starting point every time before the extrapolation
940  if(!AliEMCALRecoUtils::ExtrapolateTrackToCluster(&trackParamTmp, clus, 0.139, 5., dEta, dPhi)) continue;
941 
942  Float_t dR2 = dPhi*dPhi + dEta*dEta;
943  if(dR2 < fConeRadius){
944  fBuffer_Surrounding_Tracks_R[nTracksInR]=dR2;
945  fBuffer_Surrounding_Tracks_Pt[nTracksInR]=inTrack->Pt();
946  fBuffer_Surrounding_Tracks_P[nTracksInR]=inTrack->P();
947  fBuffer_Surrounding_Tracks_nSigdEdxE[nTracksInR]= fPIDResponse ? fPIDResponse->NumberOfSigmasTPC(inTrack,AliPID::kElectron) : -100;
950  fBuffer_Surrounding_Tracks_V0Flag[nTracksInR]=trackIsFromV0;
951  nTracksInR+=1;
952  }
953  }
954  fBuffer_Surrounding_NTracks = nTracksInR;
955  if(nTracksInR==0){
956  fBuffer_Surrounding_Tracks_R[nTracksInR]=-100;
957  fBuffer_Surrounding_Tracks_Pt[nTracksInR]=-1;
958  fBuffer_Surrounding_Tracks_P[nTracksInR]=-1;
959  fBuffer_Surrounding_Tracks_nSigdEdxE[nTracksInR]=-100;
962  fBuffer_Surrounding_Tracks_V0Flag[nTracksInR]=kFALSE;
963  }
964 
965  if(EsdTrackCuts){
966  delete EsdTrackCuts;
967  EsdTrackCuts=0x0;
968  }
969 }
970 
971 
972 
973 //________________________________________________________________________
974 Int_t AliAnalysisTaskClusterQA::ProcessTrueClusterCandidates(AliAODConversionPhoton *TrueClusterCandidate, AliVCluster* cluster,
975  AliAODConversionPhoton *TrueSubClusterCandidate1,
976  AliAODConversionPhoton *TrueSubClusterCandidate2)
977 {
978 
979  Int_t mcLabelCluster = -5;
980  Float_t m02 = cluster->GetM02();
981  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
982  Double_t mcProdVtxX = primVtxMC->GetX();
983  Double_t mcProdVtxY = primVtxMC->GetY();
984  Double_t mcProdVtxZ = primVtxMC->GetZ();
985 
986  TParticle *Photon = NULL;
987  TParticle *Pi0Dummy = NULL;
988  if (TrueClusterCandidate->GetIsCaloPhoton() == 0) AliFatal("CaloPhotonFlag has not been set task will abort");
989  if (TrueClusterCandidate->GetCaloPhotonMCLabel(0) < 0){
990  mcLabelCluster = -10;
991  return mcLabelCluster;
992  }
993 
994  // Setting all MC Flags (IsMerged, etc)
995  TrueClusterCandidate->SetCaloPhotonMCFlags(fMCEvent, kFALSE, kTRUE,cluster);
996 
997  if (TrueClusterCandidate->GetNCaloPhotonMCLabels()>0){
998  fBuffer_Cluster_MC_EFracFirstLabel = cluster->GetClusterMCEdepFraction(0);
999  if(TrueClusterCandidate->GetNNeutralPionMCLabels()>0){
1000  if(TrueClusterCandidate->GetLeadingNeutralPionIndex()>-1){
1001  fBuffer_Cluster_MC_EFracLeadingPi0 = TrueClusterCandidate->GetNeutralPionEnergyFraction(TrueClusterCandidate->GetLeadingNeutralPionIndex());
1002  Pi0Dummy = fMCEvent->Particle(TrueClusterCandidate->GetNeutralPionMCLabel(TrueClusterCandidate->GetLeadingNeutralPionIndex()));
1003  if(Pi0Dummy){
1004  fBuffer_Cluster_MC_LeadingPi0_Pt = Pi0Dummy->Pt();
1005  fBuffer_Cluster_MC_LeadingPi0_E = Pi0Dummy->Energy();
1006  if(fSaveEventsInVector) fVTrueNeutralPionDaughterIndex.push_back(TrueClusterCandidate->GetNeutralPionMCLabel(TrueClusterCandidate->GetLeadingNeutralPionIndex()));
1007  } else {
1009  }
1010  } else {
1012  }
1013  } else {
1015  }
1016  // check if leading pi0 comes not from label 0 in cluster
1017  // for this do:
1018  // -> check if neutral pions were found in cluster
1019  // -> if the leading daughter index is not 0
1020  // -> the leading neutral pion has a larger cluster energy fraction than the cluster label 0
1021  if( TrueClusterCandidate->GetNNeutralPionMCLabels()>0 && TrueClusterCandidate->GetLeadingNeutralPionDaughterIndex()!=0 && TrueClusterCandidate->GetNeutralPionEnergyFraction(TrueClusterCandidate->GetLeadingNeutralPionIndex())>cluster->GetClusterMCEdepFraction(0)) {
1022  // load particle corresponding to largest daughter of leading pi0
1023  Photon = fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMCLabel(TrueClusterCandidate->GetLeadingNeutralPionDaughterIndex()));
1024  } else {
1025  // load particle corresponding to MC label 0 in cluster
1026  Photon = fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMCLabel(0));
1027  }
1028  } else {
1029  // return if there are no MC labels in the cluster
1030  mcLabelCluster = -11;
1032  return mcLabelCluster;
1033  }
1034  if(Photon == NULL){
1035  mcLabelCluster = -12;
1037  return mcLabelCluster;
1038  }
1039 
1040  // AliAODConversionMother *mesoncand = new AliAODConversionMother(TrueSubClusterCandidate1,TrueSubClusterCandidate2);
1041  // //apply rapidity cut on clusters
1042  // Bool_t mesonIsSelected = (((AliConversionMesonCuts*)fMesonCuts)->MesonIsSelected(mesoncand,kTRUE,fEventCuts->GetEtaShift()));
1043  // if (!mesonIsSelected){
1044  // delete mesoncand;
1045  // mcLabelCluster = -13;
1046  // return mcLabelCluster;
1047  // }
1048 
1049 
1050  // cluster classification:
1051  // 1 - nice merged cluster (2 gamma | contributions from 2 gamma) from pi0/eta
1052  // 2 - contribution from only 1 partner (1 gamma, 1 fully coverted gamma) from pi0/eta
1053  // 3 - contribution from part of 1 partner (1 electron) from pi0/eta
1054 
1055 
1056  Int_t clusterClass = 0;
1057  Bool_t isPrimary = fEventCuts->IsConversionPrimaryESD( fMCEvent, TrueClusterCandidate->GetCaloPhotonMCLabel(0), mcProdVtxX, mcProdVtxY, mcProdVtxZ);
1058 
1059 
1060  Long_t motherLab = -1;
1061  if (TrueClusterCandidate->IsMerged() || TrueClusterCandidate->IsMergedPartConv()){
1062  clusterClass = 1;
1063  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
1064  } else if (TrueClusterCandidate->GetNCaloPhotonMotherMCLabels()> 0){
1065  if (TrueClusterCandidate->IsLargestComponentElectron() || TrueClusterCandidate->IsLargestComponentPhoton()){
1066  if (TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0) > -1 && (fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0))->GetPdgCode() == 111 || fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0))->GetPdgCode() == 221) ){
1067  if ( TrueClusterCandidate->IsConversion() && !TrueClusterCandidate->IsConversionFullyContained() ){
1068  clusterClass = 3;
1069  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
1070  } else {
1071  clusterClass = 2;
1072  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
1073  }
1074  }
1075  } else if (TrueClusterCandidate->IsSubLeadingEM()){
1076  if (TrueClusterCandidate->GetNCaloPhotonMotherMCLabels()> 1){
1077  if ( TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1) > -1){
1078  if (fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1))->GetPdgCode() == 111 || fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1))->GetPdgCode() == 221 ){
1079  clusterClass = 2;
1080  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1);
1081  }
1082  }
1083  }
1084  } else {
1085  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
1086  }
1087  }
1088 
1089  // Get Mother particle
1090  TParticle *mother = NULL;
1091  Int_t motherPDG = -1;
1092  if (motherLab > -1)
1093  mother = fMCEvent->Particle(motherLab);
1094  if (mother)
1095  motherPDG = TMath::Abs(mother->GetPdgCode());
1096 
1097  //
1098  if (clusterClass == 1 || clusterClass == 2 || clusterClass == 3 ){
1099  // separate different components
1100  if (clusterClass == 1 && TrueClusterCandidate->IsMerged()){
1101  if (motherPDG == 111){
1102  mcLabelCluster = 10; // NOTE merged pi0
1103  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
1104  mcLabelCluster = 11; // NOTE secondary merged pi0
1105  }
1106  if (motherPDG == 221)
1107  mcLabelCluster = 12; // NOTE merged eta
1108  } else if (clusterClass == 1 && TrueClusterCandidate->IsMergedPartConv()){
1109  if (motherPDG == 111){
1110  mcLabelCluster = 13; // NOTE merged pi0 with one converted gamma
1111  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
1112  mcLabelCluster = 14; // NOTE merged secondary pi0 with one converted gamma
1113  }
1114  if (motherPDG == 221)
1115  mcLabelCluster = 15; // NOTE merged eta with one converted gamma
1116  } else if (clusterClass == 2){
1117  if (motherPDG == 111){
1118  mcLabelCluster = 20; // NOTE decay photon from pi0
1119  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
1120  mcLabelCluster = 21; // NOTE decay photon from secondary pi0
1121  }
1122  if (motherPDG == 221)
1123  mcLabelCluster = 22; // NOTE decay photon from eta
1124  } else if (clusterClass == 3){
1125  if (motherPDG == 111) {
1126  mcLabelCluster = 30; // NOTE electron from decayed pi0
1127  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
1128  mcLabelCluster = 31; // NOTE electron from decayed secondary pi0
1129  }
1130  if (motherPDG == 221)
1131  mcLabelCluster = 32; // NOTE electron from decayed eta
1132  }
1133 
1134  // leading particle is a photon or the conversion is fully contained and its not from pi0 || eta
1135  } else if (TrueClusterCandidate->IsLargestComponentPhoton() || TrueClusterCandidate->IsConversionFullyContained()
1136  || TrueClusterCandidate->IsElectronFromFragPhoton()){
1137 
1138  if(TrueClusterCandidate->IsLargestComponentPhoton()){
1139  // 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
1140 
1141  if (motherLab == -1) mcLabelCluster = 40; // direct photon from initial collision
1142  else if ((motherLab >0) && (motherLab<9)) mcLabelCluster = 41; // photon from quark
1143  else if (motherLab == 11) mcLabelCluster = 42; // photon from electron
1144  else if (motherLab == 22){ // check for frag photon
1145 
1146  TParticle *dummyMother = fMCEvent->Particle(motherLab);
1147  Bool_t originReached = kFALSE;
1148  Bool_t isFragPhoton = kFALSE;
1149 
1150  while (dummyMother->GetPdgCode() == 22 && !originReached){ // follow photon's history, as long as the mother is a photon
1151  if (dummyMother->GetMother(0) > -1){
1152  dummyMother = fMCEvent->Particle(dummyMother->GetMother(0));
1153  if (TMath::Abs(dummyMother->GetPdgCode()) == 22){ // if mother of mother is again a photon, continue
1154  if (dummyMother->GetMother(0) > -1){
1155  dummyMother = fMCEvent->Particle(dummyMother->GetMother(0));
1156  } else {
1157  originReached = kTRUE;
1158  }
1159  } else {
1160  originReached = kTRUE;
1161  }
1162  isFragPhoton = (TMath::Abs(dummyMother->GetPdgCode()) < 6);// photon stems from quark = fragmentation photon
1163  } else {
1164  originReached = kTRUE;
1165  }
1166  }
1167 
1168  if(isFragPhoton) mcLabelCluster = 43; // Fragmentation photon
1169  else{
1170  mcLabelCluster = 47; // something like cluster <- photon <- photon <- X (not photon)
1171  // AliInfo(Form("Mother of photon is photon etc. but origin is not quark. ID: %li", motherLab));
1172  }
1173  }
1174  else{
1175  mcLabelCluster = 44; // other (e.g. from meson decays that are not pi0 or eta0)
1176  // AliInfo(Form("Single cluster is mainly produced by a photon and mother is %li", motherLab));
1177  }
1178  } else if(TrueClusterCandidate->IsConversionFullyContained()){
1179  // cluster is from a fully contained conversion
1180  mcLabelCluster = 45;
1181  } else if(TrueClusterCandidate->IsElectronFromFragPhoton()){
1182  mcLabelCluster = 46; // electron from frac
1183  }
1184 
1185  // leading particle is an electron and its not from pi0 || eta and no electron from fragmentation photon conversion
1186  } else if (TrueClusterCandidate->IsLargestComponentElectron() && !TrueClusterCandidate->IsElectronFromFragPhoton()){
1187  mcLabelCluster = 50; // NOTE single electron
1188  } else {
1189  // leading particle from hadron
1190  mcLabelCluster = 60; // NOTE hadron cluster
1191  // AliInfo(Form("Single cluster is mainly produced by hadron with id: %li", motherLab));
1192  }
1193 
1194  // delete mesoncand;
1195  return mcLabelCluster;
1196 }
1197 
1198 //________________________________________________________________________
1199 Int_t AliAnalysisTaskClusterQA::ProcessTrueClusterCandidatesAOD(AliAODConversionPhoton *TrueClusterCandidate, AliVCluster* cluster,
1200  AliAODConversionPhoton *TrueSubClusterCandidate1,
1201  AliAODConversionPhoton *TrueSubClusterCandidate2)
1202 {
1203  Int_t mcLabelCluster = -5;
1204  Float_t m02 = cluster->GetM02();
1205  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
1206  Double_t mcProdVtxX = primVtxMC->GetX();
1207  Double_t mcProdVtxY = primVtxMC->GetY();
1208  Double_t mcProdVtxZ = primVtxMC->GetZ();
1209 
1210  AliAODMCParticle *Photon = NULL;
1211  AliAODMCParticle *Pi0Dummy = NULL;
1212  TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1213 
1214  if (AODMCTrackArray){
1215  if (TrueClusterCandidate->GetIsCaloPhoton() == 0) AliFatal("CaloPhotonFlag has not been set task will abort");
1216  if (TrueClusterCandidate->GetCaloPhotonMCLabel(0) < 0) {
1217  mcLabelCluster = -10;
1218  return mcLabelCluster;
1219  }
1220  } else {
1221  AliInfo("AODMCTrackArray could not be loaded");
1222  mcLabelCluster = -90;
1223  return mcLabelCluster;
1224  }
1225 
1226  TrueClusterCandidate->SetCaloPhotonMCFlagsAOD(fInputEvent, kFALSE, kTRUE,cluster);
1227 
1228  if (TrueClusterCandidate->GetNCaloPhotonMCLabels()>0){
1229  fBuffer_Cluster_MC_EFracFirstLabel = cluster->GetClusterMCEdepFraction(0);
1230  if(TrueClusterCandidate->GetNNeutralPionMCLabels()>0){
1231  if(TrueClusterCandidate->GetLeadingNeutralPionIndex()>-1){
1232  fBuffer_Cluster_MC_EFracLeadingPi0 = TrueClusterCandidate->GetNeutralPionEnergyFraction(TrueClusterCandidate->GetLeadingNeutralPionIndex());
1233  Pi0Dummy = (AliAODMCParticle*) AODMCTrackArray->At(TrueClusterCandidate->GetNeutralPionMCLabel(TrueClusterCandidate->GetLeadingNeutralPionIndex()));
1234  if(Pi0Dummy){
1235  fBuffer_Cluster_MC_LeadingPi0_Pt = Pi0Dummy->Pt();
1236  fBuffer_Cluster_MC_LeadingPi0_E = Pi0Dummy->E();
1237  if(fSaveEventsInVector) fVTrueNeutralPionDaughterIndex.push_back(TrueClusterCandidate->GetNeutralPionMCLabel(TrueClusterCandidate->GetLeadingNeutralPionIndex()));
1238  } else {
1240  }
1241  } else {
1243  }
1244  } else {
1246  }
1247  // check if leading pi0 comes not from label 0 in cluster
1248  // for this do:
1249  // -> check if neutral pions were found in cluster
1250  // -> if the leading daughter index is not 0
1251  // -> the leading neutral pion has a larger cluster energy fraction than the cluster label 0
1252  if( TrueClusterCandidate->GetNNeutralPionMCLabels()>0 && TrueClusterCandidate->GetLeadingNeutralPionDaughterIndex()!=0 && TrueClusterCandidate->GetNeutralPionEnergyFraction(TrueClusterCandidate->GetLeadingNeutralPionIndex())>cluster->GetClusterMCEdepFraction(0)) {
1253  // load particle corresponding to largest daughter of leading pi0
1254  Photon = (AliAODMCParticle*) AODMCTrackArray->At(TrueClusterCandidate->GetCaloPhotonMCLabel(TrueClusterCandidate->GetLeadingNeutralPionDaughterIndex()));
1255  } else {
1256  // load particle corresponding to MC label 0 in cluster
1257  Photon = (AliAODMCParticle*) AODMCTrackArray->At(TrueClusterCandidate->GetCaloPhotonMCLabel(0));
1258  }
1259  } else {
1260  // return if there are no MC labels in the cluster
1261  mcLabelCluster = -11;
1263  return mcLabelCluster;
1264  }
1265 
1266  if(Photon == NULL){
1267  mcLabelCluster = -12;
1269  return mcLabelCluster;
1270  }
1271 
1272  // AliAODConversionMother *mesoncand = new AliAODConversionMother(TrueSubClusterCandidate1,TrueSubClusterCandidate2);
1273  // //apply rapidity cut on clusters
1274  // Bool_t mesonIsSelected = (((AliConversionMesonCuts*)fMesonCuts)->MesonIsSelected(mesoncand,kTRUE,fEventCuts->GetEtaShift()));
1275  // if (!mesonIsSelected){
1276  // delete mesoncand;
1277  // mcLabelCluster = -13;
1278  // return mcLabelCluster;
1279  // }
1280 
1281  Int_t clusterClass = 0;
1282  Bool_t isPrimary = ((AliConvEventCuts*)fEventCuts)->IsConversionPrimaryAOD( fInputEvent, Photon, mcProdVtxX, mcProdVtxY, mcProdVtxZ);
1283 
1284  Long_t motherLab = -1;
1285  if (TrueClusterCandidate->IsMerged() || TrueClusterCandidate->IsMergedPartConv()){
1286  clusterClass = 1;
1287  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
1288  } else if (TrueClusterCandidate->GetNCaloPhotonMotherMCLabels()> 0){
1289  if (TrueClusterCandidate->IsLargestComponentElectron() || TrueClusterCandidate->IsLargestComponentPhoton()){
1290  if (TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0) > -1 && (((AliAODMCParticle*) AODMCTrackArray->At(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0)))->GetPdgCode() == 111 || ((AliAODMCParticle*) AODMCTrackArray->At(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0)))->GetPdgCode() == 221) ){
1291  if ( TrueClusterCandidate->IsConversion() && !TrueClusterCandidate->IsConversionFullyContained() ){
1292  clusterClass = 3;
1293  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
1294  } else {
1295  clusterClass = 2;
1296  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
1297  }
1298  }
1299  } else if (TrueClusterCandidate->IsSubLeadingEM()){
1300  if (TrueClusterCandidate->GetNCaloPhotonMotherMCLabels()> 1){
1301  if ( TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1) > -1){
1302  if (TMath::Abs(((AliAODMCParticle*) AODMCTrackArray->At(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1)))->GetPdgCode()) == 111 || TMath::Abs(((AliAODMCParticle*) AODMCTrackArray->At(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1)))->GetPdgCode()) == 221 ){
1303  clusterClass = 2;
1304  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1);
1305  }
1306  }
1307  }
1308  } else {
1309  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
1310  }
1311  }
1312 
1313  // Get Mother particle
1314  AliAODMCParticle *mother = NULL;
1315  Int_t motherPDG = -1;
1316  if (motherLab > -1)
1317  mother = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(motherLab));
1318  if (mother)
1319  motherPDG = TMath::Abs(mother->GetPdgCode());
1320 
1321  if (clusterClass == 1 || clusterClass == 2 || clusterClass == 3 ){
1322  // separate different components
1323  if (clusterClass == 1 && TrueClusterCandidate->IsMerged()){
1324  if (motherPDG == 111){
1325  mcLabelCluster = 10; // NOTE merged pi0
1326  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
1327  mcLabelCluster = 11; // NOTE secondary merged pi0
1328  }
1329  if (motherPDG == 221)
1330  mcLabelCluster = 12; // NOTE merged eta
1331  } else if (clusterClass == 1 && TrueClusterCandidate->IsMergedPartConv()){
1332  if (motherPDG == 111){
1333  mcLabelCluster = 13; // NOTE merged pi0 with one converted gamma
1334  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
1335  mcLabelCluster = 14; // NOTE merged secondary pi0 with one converted gamma
1336  }
1337  if (motherPDG == 221)
1338  mcLabelCluster = 15; // NOTE merged eta with one converted gamma
1339  } else if (clusterClass == 2){
1340  if (motherPDG == 111){
1341  mcLabelCluster = 20; // NOTE decay photon from pi0
1342  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
1343  mcLabelCluster = 21; // NOTE decay photon from secondary pi0
1344  }
1345  if (motherPDG == 221)
1346  mcLabelCluster = 22; // NOTE decay photon from eta
1347  } else if (clusterClass == 3){
1348  if (motherPDG == 111) {
1349  mcLabelCluster = 30; // NOTE electron from decayed pi0
1350  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
1351  mcLabelCluster = 31; // NOTE electron from decayed secondary pi0
1352  }
1353  if (motherPDG == 221)
1354  mcLabelCluster = 32; // NOTE electron from decayed eta
1355  }
1356  // leading particle is a photon or the conversion is fully contained and its not from pi0 || eta
1357  } else if (TrueClusterCandidate->IsLargestComponentPhoton() || TrueClusterCandidate->IsConversionFullyContained()
1358  || TrueClusterCandidate->IsElectronFromFragPhoton()){
1359 
1360  if(TrueClusterCandidate->IsLargestComponentPhoton()){
1361  // 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
1362 
1363  if (motherLab == -1) mcLabelCluster = 40; // direct photon from initial collision
1364  else if ((motherLab >0) && (motherLab<9)) mcLabelCluster = 41; // photon from quark
1365  else if (motherLab == 11) mcLabelCluster = 42; // photon from electron
1366  else if (motherLab == 22){ // check for frag photon
1367 
1368  AliAODMCParticle *dummyMother = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(motherLab));
1369  Bool_t originReached = kFALSE;
1370  Bool_t isFragPhoton = kFALSE;
1371 
1372  while (dummyMother->GetPdgCode() == 22 && !originReached){ // follow photon's history, as long as the mother is a photon
1373  if (dummyMother->GetMother() > -1){
1374  dummyMother = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(dummyMother->GetMother()));
1375 
1376  if (TMath::Abs(dummyMother->GetPdgCode()) == 22){ // if mother of mother is again a photon, continue
1377  if (dummyMother->GetMother() > -1){
1378  dummyMother = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(dummyMother->GetMother()));
1379  } else {
1380  originReached = kTRUE;
1381  }
1382  } else {
1383  originReached = kTRUE;
1384  }
1385  isFragPhoton = (TMath::Abs(dummyMother->GetPdgCode()) < 6);// photon stems from quark = fragmentation photon
1386  } else {
1387  originReached = kTRUE;
1388  }
1389  }
1390 
1391  if(isFragPhoton) mcLabelCluster = 43; // Fragmentation photon
1392  else{
1393  mcLabelCluster = 47; // something like cluster <- photon <- photon <- X (not photon)
1394  // AliInfo(Form("Mother of photon is photon etc. but origin is not quark. ID: %li", motherLab));
1395  }
1396  }
1397  else{
1398  mcLabelCluster = 44; // other (e.g. from meson decays that are not pi0 or eta0)
1399  // AliInfo(Form("Single cluster is mainly produced by a photon and mother is %li", motherLab));
1400  }
1401  } else if(TrueClusterCandidate->IsConversionFullyContained()){
1402  // cluster is from a fully contained conversion
1403  mcLabelCluster = 45;
1404  } else if(TrueClusterCandidate->IsElectronFromFragPhoton()){
1405  mcLabelCluster = 46; // electron from frac
1406  }
1407 
1408  // leading particle is an electron and its not from pi0 || eta and no electron from fragmentation photon conversion
1409  } else if (TrueClusterCandidate->IsLargestComponentElectron() && !TrueClusterCandidate->IsElectronFromFragPhoton()){
1410  mcLabelCluster = 50; // NOTE single electron
1411  } else {
1412  // leading particle from hadron
1413  mcLabelCluster = 60; // NOTE hadron cluster
1414  // AliInfo(Form("Single cluster is mainly produced by hadron with id: %li", motherLab));
1415  }
1416 
1417  // delete mesoncand;
1418  return mcLabelCluster;
1419 }
1420 //_____________________________________________________________________________
1421 ULong64_t AliAnalysisTaskClusterQA::GetUniqueEventID(AliVHeader* header)
1422 {
1423  // To have a unique id for each event in a run!
1424  // Modified from AliRawReader.h
1425  return ((ULong64_t)header->GetBunchCrossNumber()+
1426  (ULong64_t)header->GetOrbitNumber()*3564+
1427  (ULong64_t)header->GetPeriodNumber()*16777215*3564);
1428 }
1429 
1430 //________________________________________________________________________
1431 //************** Find number of local maxima in cluster ******************
1432 //* derived from G. Conesa Balbastre's AliCalorimeterUtils *******************
1433 //************************************************************************
1434 // Int_t AliAnalysisTaskClusterQA::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVEvent * event){
1435 
1436 // const Int_t nc = cluster->GetNCells();
1437 
1438 // Int_t absCellIdList[nc];
1439 // Float_t maxEList[nc];
1440 
1441 // Int_t nMax = GetNumberOfLocalMaxima(cluster, event, absCellIdList, maxEList);
1442 
1443 // return nMax;
1444 // }
1445 //________________________________________________________________________
1446 //************** Find number of local maxima in cluster ******************
1447 //* derived from G. Conesa Balbastre's AliCalorimeterUtils ***************
1448 //************************************************************************
1449 // Int_t AliAnalysisTaskClusterQA::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVEvent * event, Int_t *absCellIdList, Float_t* maxEList){
1450 
1451 // Int_t absCellId1 = -1;
1452 // Int_t absCellId2 = -1;
1453 // const Int_t nCells = cluster->GetNCells();
1454 // AliVCaloCells* cells = NULL;
1455 
1456 // // if (fClusterType == 1 || fClusterType == 3 || fClusterType == 4)
1457 // cells = event->GetEMCALCells();
1458 // // else if (fClusterType ==2 )
1459 // // cells = event->GetPHOSCells();
1460 
1461 // Float_t eMax = 0.;
1462 // Int_t idMax = -1;
1463 
1464 // for (Int_t iCell = 0;iCell < nCells;iCell++){
1465 // absCellIdList[iCell] = cluster->GetCellsAbsId()[iCell];
1466 // Int_t imod = -1, icol = -1, irow = -1;
1467 // imod = ((AliCaloPhotonCuts*)fClusterCutsEMC)->GetModuleNumberAndCellPosition(absCellIdList[iCell], icol, irow);
1468 // if (cells->GetCellAmplitude(absCellIdList[iCell])> eMax){
1469 // eMax = cells->GetCellAmplitude(absCellIdList[iCell]);
1470 // idMax = absCellIdList[iCell];
1471 // }
1472 // }
1473 
1474 // // find the largest separated cells in cluster
1475 // for (Int_t iCell = 0;iCell < nCells;iCell++){
1476 // // check whether valid cell number is selected
1477 // if (absCellIdList[iCell] >= 0){
1478 // // store current energy and cell id
1479 // absCellId1 = cluster->GetCellsAbsId()[iCell];
1480 // Float_t en1 = cells->GetCellAmplitude(absCellId1);
1481 
1482 // // loop over other cells in cluster
1483 // for (Int_t iCellN = 0;iCellN < nCells;iCellN++){
1484 // // jump out if array has changed in the meantime
1485 // if (absCellIdList[iCell] == -1) continue;
1486 // // get cell id & check whether its valid
1487 // absCellId2 = cluster->GetCellsAbsId()[iCellN];
1488 
1489 // // don't compare to yourself
1490 // if (absCellId2 == -1) continue;
1491 // if (absCellId1 == absCellId2) continue;
1492 
1493 // // get cell energy of second cell
1494 // Float_t en2 = cells->GetCellAmplitude(absCellId2);
1495 
1496 // // check if cells are Neighbours
1497 // if (((AliCaloPhotonCuts*)fClusterCutsEMC)->AreNeighbours(absCellId1, absCellId2)){
1498 // // determine which cell has larger energy, mask the other
1499 // if (en1 > en2 ){
1500 // absCellIdList[iCellN] = -1;
1501 // if (en1 < en2 )
1502 // absCellIdList[iCell] = -1;
1503 // } else {
1504 // absCellIdList[iCell] = -1;
1505 // if (en1 > en2 )
1506 // absCellIdList[iCellN] = -1;
1507 // }
1508 // }
1509 // }
1510 // }
1511 // }
1512 
1513 // // shrink list of cells to only maxima
1514 // Int_t nMaximaNew = 0;
1515 // for (Int_t iCell = 0;iCell < nCells;iCell++){
1516 // if (absCellIdList[iCell] > -1){
1517 // Float_t en = cells->GetCellAmplitude(absCellIdList[iCell]);
1518 // // check whether cell energy is larger than required seed
1519 // if (en < 0.1) continue;
1520 // absCellIdList[nMaximaNew] = absCellIdList[iCell];
1521 // maxEList[nMaximaNew] = en;
1522 // fBuffer_ClusterNLM_ID[nMaximaNew] = absCellIdList[nMaximaNew];
1523 // fBuffer_ClusterNLM_E[nMaximaNew] = maxEList[nMaximaNew];
1524 // nMaximaNew++;
1525 // }
1526 // }
1527 
1528 // // check whether a local maximum was found
1529 // // if no maximum was found use highest cell as maximum
1530 // if (nMaximaNew == 0){
1531 // nMaximaNew = 1;
1532 // maxEList[0] = eMax;
1533 // absCellIdList[0] = idMax;
1534 // }
1535 
1536 // return nMaximaNew;
1537 // }
1538 
1539 //-------------------------------------------------------------
1541 { // Get Event Centrality
1542 
1543  AliESDEvent *esdEvent=dynamic_cast<AliESDEvent*>(event);
1544  if(esdEvent){
1545  AliMultSelection *MultSelection = (AliMultSelection*)event->FindListObject("MultSelection");
1546  if(!MultSelection){
1547  AliWarning ("AliMultSelection object not found !");
1548  return -1;
1549  }else{
1550  return MultSelection->GetMultiplicityPercentile("V0M");// default
1551  }
1552  }
1553 
1554 
1555  return -1;
1556 }
1558  fBuffer_EventWeight = 1;
1559  fBuffer_ClusterE = 0;
1560  fBuffer_ClusterPhi = 0;
1561  fBuffer_ClusterEta = 0;
1562  fBuffer_ClusterIsEMCAL = kFALSE;
1563  fBuffer_ClusterSupMod = -1;
1570  fBuffer_ClusterM02 = 0;
1571  fBuffer_ClusterM20 = 0;
1583  for(Int_t cell = 0; cell < kMaxActiveCells; cell++){
1584  fBuffer_Cells_E[cell] = 0;
1585  fBuffer_Cells_RelativeEta[cell] = 0;
1586  fBuffer_Cells_RelativePhi[cell] = 0;
1587  fBuffer_Surrounding_Cells_ID[cell] = 0;
1588  fBuffer_Surrounding_Cells_R[cell] = 0;
1589  fBuffer_Surrounding_Cells_E[cell] = 0;
1592  }
1593  for(Int_t track = 0; track < kMaxNTracks; track++){
1594  fBuffer_Surrounding_Tracks_R[track] = 100;
1595  fBuffer_Surrounding_Tracks_Pt[track] = 0;
1596  fBuffer_Surrounding_Tracks_P[track] = 0;
1600  fBuffer_Surrounding_Tracks_V0Flag[track] = kFALSE;
1601  }
1602 }
1603 
1604 //________________________________________________________________________
1606 
1607  fVBuffer_Cluster_E.clear();
1608  fVBuffer_Cluster_Eta.clear();
1609  fVBuffer_Cluster_Phi.clear();
1610  fVBuffer_Cluster_isEMCal.clear();
1612 
1613  fVBuffer_Cluster_E.resize(0);
1614  fVBuffer_Cluster_Eta.resize(0);
1615  fVBuffer_Cluster_Phi.resize(0);
1616  fVBuffer_Cluster_isEMCal.resize(0);
1618 
1619 }
Int_t fBuffer_Surrounding_NTracks
! array buffer
Bool_t * fBuffer_Surrounding_Tracks_V0Flag
! array buffer
double Double_t
Definition: External.C:58
Float_t fBuffer_ClusterPhi
! array buffer
virtual void UserExec(Option_t *option)
Int_t ProcessTrueClusterCandidatesAOD(AliAODConversionPhoton *TrueClusterCandidate, AliVCluster *cluster, AliAODConversionPhoton *TrueSubClusterCandidate1, AliAODConversionPhoton *TrueSubClusterCandidate2)
ULong64_t GetUniqueEventID(AliVHeader *header)
Float_t fBuffer_Event_Vertex_Y
! array buffer
Int_t MakePhotonCandidates(AliVCluster *clus, AliVCaloCells *cells, Long_t indexCluster)
Bool_t fSaveMCInformation
save MC information
const Int_t kMaxNTracks
Float_t fBuffer_ClusterE
! array buffer
Float_t * fBuffer_Cells_RelativeEta
! array buffer
AliPIDResponse * fPIDResponse
PID response.
Bool_t fSaveEventProperties
save general event properties (centrality etc.)
std::vector< Float_t > fVBuffer_Cluster_Eta
! vector buffer
Float_t fConeRadius
save arrays of all cells in event
Int_t fBuffer_ClusterSupMod
! array buffer
Int_t fBuffer_Mother_MC_Label
! array buffer
Float_t * fBuffer_Surrounding_Tracks_R
! array buffer
Float_t fBuffer_Event_Vertex_Z
! array buffer
Bool_t fBuffer_ClusterIsEMCAL
! array buffer
Float_t fBuffer_EventWeight
! array buffer
Float_t * fBuffer_Surrounding_Tracks_nSigdEdxE
! array buffer
Float_t fBuffer_Event_Vertex_X
! array buffer
Int_t fBuffer_LeadingCell_ID
! array buffer
std::vector< Float_t > fVBuffer_Cluster_E
! vector 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 Int_t
Definition: External.C:63
Float_t * fBuffer_Surrounding_Cells_R
! array buffer
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
Float_t fBuffer_Cluster_MC_EFracFirstLabel
! array buffer
Int_t fBuffer_Cluster_MC_Label
! array buffer
Float_t fBuffer_Cluster_MC_EFracLeadingPi0
! array buffer
Bool_t fSaveSurroundingCells
save arrays of all cells in event
Float_t fBuffer_ClusterM20
! array buffer
Float_t fMinClusterEnergy
save arrays of all cells in event
Int_t * fBuffer_Cells_ID
! array buffer
Int_t fMaxNLMCut
save MC information
std::vector< Bool_t > fVBuffer_Cluster_isEMCal
! vector buffer
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)
Float_t fBuffer_Cluster_MC_LeadingPi0_E
! array buffer
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)
std::vector< Int_t > fVTrueNeutralPionDaughterIndex
! vector buffer store the MC stack ID of mother pi0 for true information
Float_t fBuffer_Cluster_MC_LeadingPi0_Pt
! array buffer
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)
Float_t * fBuffer_Surrounding_Tracks_P
! array buffer
Float_t * fBuffer_Surrounding_Cells_RelativePhi
! array buffer
Int_t fBuffer_ClusterNumCells
! array buffer
Int_t ProcessTrueClusterCandidates(AliAODConversionPhoton *TrueClusterCandidate, AliVCluster *cluster, AliAODConversionPhoton *TrueSubClusterCandidate1, AliAODConversionPhoton *TrueSubClusterCandidate2)
Float_t fBuffer_LeadingCell_Eta
! array buffer
const Int_t kMaxActiveCells
Int_t fBuffer_MC_Cluster_Flag
! array buffer
Bool_t fSaveEventsInVector
save cluster information in event vectors information
Float_t * fBuffer_Surrounding_Cells_RelativeEta
! array buffer
Float_t fBuffer_LeadingCell_E
! array buffer
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
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
std::vector< Float_t > fVBuffer_Cluster_Phi
! vector buffer
Float_t fMinTrackPt
save arrays of all cells in event