AliPhysics  59e0e03 (59e0e03)
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 
35 class iostream;
36 
37 using namespace std;
38 
40 
41 //________________________________________________________________________
43  fV0Reader(NULL),
44  fV0ReaderName("V0ReaderV1"),
45  fCorrTaskSetting(""),
46  fConversionCuts(NULL),
47  fEventCuts(NULL),
48  fClusterCutsEMC(NULL),
49  fClusterCutsDMC(NULL),
50  fMesonCuts(NULL),
51  fMinNLMCut(1),
52  fMaxNLMCut(1),
53  fInputEvent(NULL),
54  fMCEvent(NULL),
55  fWeightJetJetMC(1),
56  fGeomEMCAL(NULL),
57  fClusterTree(NULL),
58  fIsHeavyIon(kFALSE),
59  ffillTree(-100),
60  ffillHistograms(kFALSE),
61  fOutputList(NULL),
62  fIsMC(kFALSE),
63  fCorrectForNonlinearity(kFALSE),
64  fSaveEventProperties(0),
65  fSaveCells(0),
66  fSaveSurroundingCells(0),
67  fSaveMCInformation(0),
68  fNSurroundingCells(10),
69  fExtractionPercentages(),
70  fExtractionPercentagePtBins(),
71  fBuffer_ClusterE(0),
72  fBuffer_ClusterPhi(0),
73  fBuffer_ClusterEta(0),
74  fBuffer_ClusterIsEMCAL(kFALSE),
75  fBuffer_MC_Cluster_Flag(0),
76  fBuffer_ClusterNumCells(0),
77  fBuffer_LeadingCell_ID(0),
78  fBuffer_LeadingCell_Row(0),
79  fBuffer_LeadingCell_Column(0),
80  fBuffer_ClusterM02(0),
81  fBuffer_ClusterM20(0),
82  fBuffer_Event_Vertex_X(0),
83  fBuffer_Event_Vertex_Y(0),
84  fBuffer_Event_Vertex_Z(0),
85  fBuffer_Event_Multiplicity(0),
86  fBuffer_Event_NumActiveCells(0),
87  fBuffer_Cells_ID(0),
88  fBuffer_Cells_E(0),
89  fBuffer_Cells_RelativeRow(0),
90  fBuffer_Cells_RelativeColumn(0),
91  fBuffer_Surrounding_Cells_ID(0),
92  fBuffer_Surrounding_Cells_E(0),
93  fBuffer_Surrounding_Cells_RelativeRow(0),
94  fBuffer_Surrounding_Cells_RelativeColumn(0),
95  fBuffer_Cluster_MC_Label(-10)
96 {
97  fBuffer_Cells_ID = new Int_t[kMaxActiveCells];
98  fBuffer_Cells_E = new Float_t[kMaxActiveCells];
99  fBuffer_Cells_RelativeRow = new Int_t[kMaxActiveCells];
100  fBuffer_Cells_RelativeColumn = new Int_t[kMaxActiveCells];
101  fBuffer_Surrounding_Cells_ID = new Int_t[kMaxActiveCells];
102  fBuffer_Surrounding_Cells_E = new Float_t[kMaxActiveCells];
103  fBuffer_Surrounding_Cells_RelativeRow = new Int_t[kMaxActiveCells];
104  fBuffer_Surrounding_Cells_RelativeColumn = new Int_t[kMaxActiveCells];
105 }
106 
108  fV0Reader(NULL),
109  fV0ReaderName("V0ReaderV1"),
110  fCorrTaskSetting(""),
111  fConversionCuts(NULL),
112  fEventCuts(NULL),
113  fClusterCutsEMC(NULL),
114  fClusterCutsDMC(NULL),
115  fMesonCuts(NULL),
116  fMinNLMCut(1),
117  fMaxNLMCut(1),
118  fInputEvent(NULL),
119  fMCEvent(NULL),
120  fWeightJetJetMC(1),
121  fGeomEMCAL(NULL),
122  fClusterTree(NULL),
123  fIsHeavyIon(kFALSE),
124  ffillTree(-100),
125  ffillHistograms(kFALSE),
126  fOutputList(NULL),
127  fIsMC(kFALSE),
128  fCorrectForNonlinearity(kFALSE),
129  fSaveEventProperties(0),
130  fSaveCells(0),
131  fSaveSurroundingCells(0),
132  fSaveMCInformation(0),
133  fNSurroundingCells(10),
134  fExtractionPercentages(),
135  fExtractionPercentagePtBins(),
136  fBuffer_ClusterE(0),
137  fBuffer_ClusterPhi(0),
138  fBuffer_ClusterEta(0),
139  fBuffer_ClusterIsEMCAL(kFALSE),
140  fBuffer_MC_Cluster_Flag(0),
141  fBuffer_ClusterNumCells(0),
142  fBuffer_LeadingCell_ID(0),
143  fBuffer_LeadingCell_Row(0),
144  fBuffer_LeadingCell_Column(0),
145  fBuffer_ClusterM02(0),
146  fBuffer_ClusterM20(0),
147  fBuffer_Event_Vertex_X(0),
148  fBuffer_Event_Vertex_Y(0),
149  fBuffer_Event_Vertex_Z(0),
150  fBuffer_Event_Multiplicity(0),
151  fBuffer_Event_NumActiveCells(0),
152  fBuffer_Cells_ID(0),
153  fBuffer_Cells_E(0),
154  fBuffer_Cells_RelativeRow(0),
155  fBuffer_Cells_RelativeColumn(0),
156  fBuffer_Surrounding_Cells_ID(0),
157  fBuffer_Surrounding_Cells_E(0),
158  fBuffer_Surrounding_Cells_RelativeRow(0),
159  fBuffer_Surrounding_Cells_RelativeColumn(0),
160  fBuffer_Cluster_MC_Label(-10)
161 {
170  // Default constructor
171 
172  DefineInput(0, TChain::Class());
173  DefineOutput(1, TList::Class());
174  DefineOutput(2, TTree::Class());
175 }
176 
177 //________________________________________________________________________
179 {
180  // default deconstructor
181 
182 }
183 //________________________________________________________________________
185 {
186  // Create User Output Objects
187 
188  if(fOutputList != NULL){
189  delete fOutputList;
190  fOutputList = NULL;
191  }
192  if(fOutputList == NULL){
193  fOutputList = new TList();
194  fOutputList->SetOwner(kTRUE);
195  }
196 
197  if(ffillHistograms){
198 
199  if(((AliConvEventCuts*)fEventCuts)->GetCutHistograms()){
200  fOutputList->Add(((AliConvEventCuts*)fEventCuts)->GetCutHistograms());
201  }
202 
203  if(((AliCaloPhotonCuts*)fClusterCutsEMC)->GetCutHistograms()){
204  fOutputList->Add(((AliCaloPhotonCuts*)fClusterCutsEMC)->GetCutHistograms());
205  }
206  if(((AliCaloPhotonCuts*)fClusterCutsDMC)->GetCutHistograms()){
207  fOutputList->Add(((AliCaloPhotonCuts*)fClusterCutsDMC)->GetCutHistograms());
208  }
209  if(((AliConversionMesonCuts*)fMesonCuts)->GetCutHistograms()){
210  fOutputList->Add(((AliConversionMesonCuts*)fMesonCuts)->GetCutHistograms());
211  }
212  PostData(1, fOutputList);
213  }
214 
215  fClusterTree = new TTree("ClusterQA","ClusterQA");
216 
217  fClusterTree->Branch("Cluster_E", &fBuffer_ClusterE, "Cluster_E/F");
218  fClusterTree->Branch("Cluster_Phi", &fBuffer_ClusterPhi, "Cluster_Phi/F");
219  fClusterTree->Branch("Cluster_Eta", &fBuffer_ClusterEta, "Cluster_Eta/F");
220  fClusterTree->Branch("Cluster_IsEMCAL", &fBuffer_ClusterIsEMCAL, "Cluster_IsEMCAL/O");
221  fClusterTree->Branch("Cluster_NumCells", &fBuffer_ClusterNumCells, "Cluster_NumCells/I");
222  fClusterTree->Branch("Cluster_LeadingCell_ID", &fBuffer_LeadingCell_ID, "Cluster_LeadingCell_ID/I");
223  fClusterTree->Branch("Cluster_LeadingCell_Row", &fBuffer_LeadingCell_Row, "Cluster_LeadingCell_Row/I");
224  fClusterTree->Branch("Cluster_LeadingCell_Column", &fBuffer_LeadingCell_Column, "Cluster_LeadingCell_Column/I");
225  fClusterTree->Branch("Cluster_M02", &fBuffer_ClusterM02, "Cluster_M02/F");
226  fClusterTree->Branch("Cluster_M20", &fBuffer_ClusterM20, "Cluster_M20/F");
228  {
229  fClusterTree->Branch("Event_Vertex_X", &fBuffer_Event_Vertex_X, "Event_Vertex_X/F");
230  fClusterTree->Branch("Event_Vertex_Y", &fBuffer_Event_Vertex_Y, "Event_Vertex_Y/F");
231  fClusterTree->Branch("Event_Vertex_Z", &fBuffer_Event_Vertex_Z, "Event_Vertex_Z/F");
232  fClusterTree->Branch("Event_Multiplicity", &fBuffer_Event_Multiplicity, "Event_Multiplicity/F");
233  fClusterTree->Branch("Event_NumActiveCells", &fBuffer_Event_NumActiveCells, "Event_NumActiveCells/I");
234  }
235 
236  if(fSaveCells)
237  {
238  fClusterTree->Branch("Cluster_Cells_ID", fBuffer_Cells_ID, "Cluster_Cells_ID[Cluster_NumCells]/I");
239  fClusterTree->Branch("Cluster_Cells_E", fBuffer_Cells_E, "Cluster_Cells_E[Cluster_NumCells]/F");
240  fClusterTree->Branch("Cluster_Cells_RelativeRow", fBuffer_Cells_RelativeRow, "Cluster_Cells_RelativeRow[Cluster_NumCells]/I");
241  fClusterTree->Branch("Cluster_Cells_RelativeColumn", fBuffer_Cells_RelativeColumn, "Cluster_Cells_RelativeColumn[Cluster_NumCells]/I");
242  }
243 
245  {
246  fClusterTree->Branch("Surrounding_Cells_ID", fBuffer_Surrounding_Cells_ID, "Surrounding_Cells_ID[Event_NumActiveCells]/I");
247  fClusterTree->Branch("Surrounding_Cells_E", fBuffer_Surrounding_Cells_E, "Surrounding_Cells_E[Event_NumActiveCells]/F");
248  fClusterTree->Branch("Surrounding_Cells_RelativeRow", fBuffer_Surrounding_Cells_RelativeRow, "Surrounding_Cells_RelativeRow[Event_NumActiveCells]/I");
249  fClusterTree->Branch("Surrounding_Cells_RelativeColumn",fBuffer_Surrounding_Cells_RelativeColumn, "Surrounding_Cells_RelativeColumn[Event_NumActiveCells]/I");
250  }
251 
253  {
254  fClusterTree->Branch("Cluster_MC_Label", &fBuffer_Cluster_MC_Label, "Cluster_MC_Label/I");
255  }
256 
257  fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask(fV0ReaderName.Data());
258  OpenFile(2);
259  PostData(2, fClusterTree);
260 }
261 //_____________________________________________________________________________
263 {
268  }
269 
270  return kTRUE;
271 }
272 //________________________________________________________________________
274 
275  Int_t eventQuality = ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEventQuality();
276  if(eventQuality != 0){// Event Not Accepted
277  return;
278  }
279  fInputEvent = InputEvent();
280  if(fIsMC) fMCEvent = MCEvent();
281 
283  if(eventNotAccepted) return; // Check Centrality, PileUp, SDD and V0AND --> Not Accepted => eventQuality = 1
284 
285  fGeomEMCAL = AliEMCALGeometry::GetInstance();
286  if(!fGeomEMCAL){ AliFatal("EMCal geometry not initialized!");}
287 
288  Int_t nclus = 0;
289  TClonesArray * arrClustersProcess = NULL;
290  // fNCurrentClusterBasic = 0;
291  if(!fCorrTaskSetting.CompareTo("")){
292  nclus = fInputEvent->GetNumberOfCaloClusters();
293  } else {
294  arrClustersProcess = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(Form("%sClustersBranch",fCorrTaskSetting.Data())));
295  if(!arrClustersProcess)
296  AliFatal(Form("%sClustersBranch was not found in AliAnalysisTaskGammaCalo! Check the correction framework settings!",fCorrTaskSetting.Data()));
297  nclus = arrClustersProcess->GetEntries();
298  }
299 
300  if(nclus == 0) return;
301 
302  ((AliCaloPhotonCuts*)fClusterCutsEMC)->FillHistogramsExtendedQA(fInputEvent,fIsMC);
303  // ((AliCaloPhotonCuts*)fClusterCutsDMC)->FillHistogramsExtendedQA(fInputEvent,fIsMC);
304 
305  // match tracks to clusters
306  ((AliCaloPhotonCuts*)fClusterCutsEMC)->MatchTracksToClusters(fInputEvent,fWeightJetJetMC,kTRUE, fMCEvent);
307  // ((AliCaloPhotonCuts*)fClusterCutsDMC)->MatchTracksToClusters(fInputEvent,fWeightJetJetMC,kTRUE, fMCEvent);
308 
309  // vertex
310  Double_t vertex[3] = {0};
311  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
312 
313  map<Long_t,Int_t> mapIsClusterAccepted;
314  map<Long_t,Int_t> mapIsClusterAcceptedWithoutTrackMatch;
315  // Loop over EMCal clusters
316  for(Long_t i = 0; i < nclus; i++){
317  Double_t tempClusterWeight = fWeightJetJetMC;
318  AliVCluster* clus = NULL;
319  if(fInputEvent->IsA()==AliESDEvent::Class()){
320  if(arrClustersProcess)
321  clus = new AliESDCaloCluster(*(AliESDCaloCluster*)arrClustersProcess->At(i));
322  else
323  clus = new AliESDCaloCluster(*(AliESDCaloCluster*)fInputEvent->GetCaloCluster(i));
324  } else if(fInputEvent->IsA()==AliAODEvent::Class()){
325  if(arrClustersProcess)
326  clus = new AliAODCaloCluster(*(AliAODCaloCluster*)arrClustersProcess->At(i));
327  else
328  clus = new AliAODCaloCluster(*(AliAODCaloCluster*)fInputEvent->GetCaloCluster(i));
329  }
330  if(!clus) continue;
331 
332 
333  if ( !clus->IsEMCAL()){ // for PHOS: cluster->GetType() == AliVCluster::kPHOSNeutral
334  delete clus;
335  continue;
336  }
337 
338 
339  if(!((AliCaloPhotonCuts*)fClusterCutsEMC)->ClusterIsSelected(clus,fInputEvent,fMCEvent,fIsMC, tempClusterWeight,i) && !((AliCaloPhotonCuts*)fClusterCutsDMC)->ClusterIsSelected(clus,fInputEvent,fMCEvent,fIsMC, tempClusterWeight,i)){
340  delete clus;
341  continue;
342  }
344  }
345  PostData(1, fOutputList);
346 }
347 
349 void AliAnalysisTaskClusterQA::ProcessQATreeCluster(AliVEvent *event, AliVCluster* cluster, Long_t indexCluster){
350 
351  Float_t clusPos[3] = { 0,0,0 };
352  cluster->GetPosition(clusPos);
353  TVector3 clusterVector(clusPos[0],clusPos[1],clusPos[2]);
354  Double_t etaCluster = clusterVector.Eta();
355  Double_t phiCluster = clusterVector.Phi();
356  if (phiCluster < 0) phiCluster += 2*TMath::Pi();
357 
358  // Vertex position x, y, z
359  fBuffer_Event_Vertex_X = fInputEvent->GetPrimaryVertex()->GetX();
360  fBuffer_Event_Vertex_Y = fInputEvent->GetPrimaryVertex()->GetY();
361  fBuffer_Event_Vertex_Z = fInputEvent->GetPrimaryVertex()->GetZ();
362 
363  // V0-based multiplicity of the event
364  fBuffer_Event_Multiplicity = fInputEvent->GetVZEROData()->GetMTotV0A()+fInputEvent->GetVZEROData()->GetMTotV0C();
365 
366  // Properties of the current cluster
367  fBuffer_ClusterE = cluster->E();
368 
369  if(etaCluster < 0.67 && etaCluster > -0.67){
370  if(phiCluster < 3.2 && phiCluster > 1.4){
371  fBuffer_ClusterIsEMCAL = kTRUE;
372  }
373  } else if(etaCluster < 0.67 && etaCluster > -0.67){
374  if(etaCluster > 0.23 || etaCluster < -0.23){
375  if(phiCluster < 5.6 && phiCluster > 4.6){
376  fBuffer_ClusterIsEMCAL = kFALSE;
377  }
378  }
379  }
380 
381  fBuffer_ClusterPhi = phiCluster;
382  fBuffer_ClusterEta = etaCluster;
383  fBuffer_ClusterM02 = cluster->GetM02();
384  fBuffer_ClusterM20 = cluster->GetM20();
385 
386  // Get all cells from the event
387  AliVCaloCells* cells = NULL;
388  if(cluster->IsEMCAL())
389  cells = event->GetEMCALCells();
390  else
391  return;
392  // else if(cluster->IsPHOS())
393  // cells = event->GetPHOSCells();
394 
395  // Determine all cells in the event
396  fBuffer_Event_NumActiveCells = cells->GetNumberOfCells();
397 
398  // Get the number of cells from the current cluster
399  Int_t nCellCluster = cluster->GetNCells();
400  fBuffer_ClusterNumCells = nCellCluster;
401 
402  // Find the leading cell in the cluster and its position
403  Int_t rowLeading=0, columnLeading=0;
405  GetRowAndColumnFromAbsCellID(fBuffer_LeadingCell_ID, rowLeading, columnLeading);
406  fBuffer_LeadingCell_Row = rowLeading;
407  fBuffer_LeadingCell_Column = columnLeading;
408 
409  // Treat the remaining cells of the cluster and get their relative position compared to the leading cell
410  for(Int_t iCell=0;iCell<nCellCluster;iCell++){
411  if(iCell<100){ // maximum number of cells for a cluster is set to 100
412  fBuffer_Cells_ID[iCell] = cluster->GetCellAbsId(iCell);
413  fBuffer_Cells_E[iCell] = cells->GetCellAmplitude(cluster->GetCellAbsId(iCell));
414  Int_t row=0, column=0;
415  GetRowAndColumnFromAbsCellID(fBuffer_Cells_ID[iCell], row, column);
416  fBuffer_Cells_RelativeRow[iCell] = row-rowLeading;
417  fBuffer_Cells_RelativeColumn[iCell] = column-columnLeading;
418  }
419  }
420 
421  // fill surrounding cell buffer
422  // fBuffer_Event_NumActiveCellsAboveThreshold = 0;
423  for(Int_t aCell=0;aCell<cells->GetNumberOfCells();aCell++){
424  // Define necessary variables
425  Short_t cellNumber = 0;
426  Double_t cellAmplitude = 0, cellTime = 0, cellEFrac = 0;
427  Int_t cellMCLabel = 0, row = 0, column = 0;
428 
429  // Get Cell ID
430  cells->GetCell(aCell,cellNumber,cellAmplitude,cellTime,cellMCLabel,cellEFrac);
431 
432  // Get row and column for each cell
433  GetRowAndColumnFromAbsCellID(cellNumber, row, column);
434 
435  // Select those cells that are within fNSurroundingCells of the leading cluster cell
436  if( (TMath::Abs(row-rowLeading)<fNSurroundingCells) && (TMath::Abs(column-columnLeading))<fNSurroundingCells){
437  fBuffer_Surrounding_Cells_RelativeRow[aCell] = row-rowLeading;
438  fBuffer_Surrounding_Cells_RelativeColumn[aCell] = column-columnLeading;
439  fBuffer_Surrounding_Cells_E[aCell] = cells->GetCellAmplitude(cellNumber);
440  fBuffer_Surrounding_Cells_ID[aCell] = cellNumber;
441  }
442  }
443  if(fIsMC) fBuffer_Cluster_MC_Label = MakePhotonCandidates(cluster, cells,indexCluster);
444 
445  // Add everything to the tree
446  if (fClusterTree) fClusterTree->Fill();
447 }
448 
449 //________________________________________________________________________
451 
452 }
453 
454 //________________________________________________________________________
455 Int_t AliAnalysisTaskClusterQA::FindLargestCellInCluster(AliVCluster* cluster, AliVEvent* event){
456 
457  const Int_t nCells = cluster->GetNCells();
458  AliVCaloCells* cells = NULL;
459 
460  if(cluster->IsEMCAL())
461  cells = event->GetEMCALCells();
462  else if(cluster->IsPHOS())
463  cells = event->GetPHOSCells();
464 
465  Float_t eMax = 0.;
466  Int_t idMax = -1;
467 
468  if (nCells < 1) return idMax;
469  for (Int_t iCell = 0;iCell < nCells;iCell++){
470  Int_t cellAbsID = cluster->GetCellsAbsId()[iCell];
471  if (cells->GetCellAmplitude(cellAbsID)> eMax){
472  eMax = cells->GetCellAmplitude(cellAbsID);
473  idMax = cellAbsID;
474  }
475  }
476  return idMax;
477 }
478 
479 //________________________________________________________________________
481  Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0;
482  row=0;
483  column=0;
484  // Get SM number and relative row/column for SM
485  fGeomEMCAL->GetCellIndex(cellIndex, nSupMod,nModule,nIphi,nIeta);
486  fGeomEMCAL->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, row,column);
487  row += nSupMod/2 * 24;
488  column += nSupMod%2 * 48;
489 }
490 
491 //________________________________________________________________________
492 Int_t AliAnalysisTaskClusterQA::MakePhotonCandidates(AliVCluster* clus, AliVCaloCells* cells, Long_t indexCluster){
493 
494  Double_t vertex[3] = {0};
495  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
496 
497  TLorentzVector clusterVector;
498  clus->GetMomentum(clusterVector,vertex);
499 
500  TLorentzVector* tmpvec = new TLorentzVector();
501  tmpvec->SetPxPyPzE(clusterVector.Px(),clusterVector.Py(),clusterVector.Pz(),clusterVector.E());
502 
503  // convert to AODConversionPhoton
504  AliAODConversionPhoton *PhotonCandidate=new AliAODConversionPhoton(tmpvec);
505  if(!PhotonCandidate){
506  delete clus;
507  delete tmpvec;
508  if (PhotonCandidate) delete PhotonCandidate;
509  return -9;
510  }
511  // Flag Photon as CaloPhoton
512  PhotonCandidate->SetIsCaloPhoton();
513  PhotonCandidate->SetCaloClusterRef(indexCluster);
514 
515  // get MC label
516  if(fIsMC> 0){
517  Int_t* mclabelsCluster = clus->GetLabels();
518  Int_t nValidClusters = 0;
519  if (clus->GetNLabels()>0){
520  for (Int_t k =0; k< (Int_t)clus->GetNLabels(); k++){
521  if (mclabelsCluster[k]>0){
522  if (nValidClusters< 50)PhotonCandidate->SetCaloPhotonMCLabel(nValidClusters,mclabelsCluster[k]);
523  nValidClusters++;
524  }
525  }
526  }
527  PhotonCandidate->SetNCaloPhotonMCLabels(nValidClusters);
528  }
529 
530  AliAODCaloCluster* clusSub1 = new AliAODCaloCluster();
531  AliAODCaloCluster* clusSub2 = new AliAODCaloCluster();
532 
533 
534  // split clusters according to their shares in the cluster (NLM == 1) needs to be treated differently
535  if (fMinNLMCut == 1 && fMaxNLMCut == 1 ){
536  Int_t absCellIdFirst = ((AliCaloPhotonCuts*)fClusterCutsEMC)->FindLargestCellInCluster(clus, fInputEvent);
537  Int_t absCellIdSecond = ((AliCaloPhotonCuts*)fClusterCutsEMC)->FindSecondLargestCellInCluster(clus, fInputEvent);
538 
539  ((AliCaloPhotonCuts*)fClusterCutsEMC)->SplitEnergy(absCellIdFirst, absCellIdSecond, clus, fInputEvent, fIsMC, clusSub1, clusSub2);
540  } else if (fMinNLMCut > 1 ){
541  const Int_t nc = clus->GetNCells();
542  Int_t absCellIdList[nc];
543 
544  ((AliCaloPhotonCuts*)fClusterCutsEMC)->SplitEnergy(absCellIdList[0], absCellIdList[1], clus, fInputEvent, fIsMC, clusSub1, clusSub2);
545  }
546 
547  // TLorentzvector with sub cluster 1
548  TLorentzVector clusterVector1;
549  clusSub1->GetMomentum(clusterVector1,vertex);
550  TLorentzVector* tmpvec1 = new TLorentzVector();
551  tmpvec1->SetPxPyPzE(clusterVector1.Px(),clusterVector1.Py(),clusterVector1.Pz(),clusterVector1.E());
552  // convert to AODConversionPhoton
553  AliAODConversionPhoton *PhotonCandidate1=new AliAODConversionPhoton(tmpvec1);
554  if(!PhotonCandidate1){
555  delete clusSub1;
556  delete tmpvec1;
557  return -9;
558  }
559  // Flag Photon as CaloPhoton
560  PhotonCandidate1->SetIsCaloPhoton();
561  // TLorentzvector with sub cluster 2
562  TLorentzVector clusterVector2;
563  clusSub2->GetMomentum(clusterVector2,vertex);
564  TLorentzVector* tmpvec2 = new TLorentzVector();
565  tmpvec2->SetPxPyPzE(clusterVector2.Px(),clusterVector2.Py(),clusterVector2.Pz(),clusterVector2.E());
566  // convert to AODConversionPhoton
567  AliAODConversionPhoton *PhotonCandidate2=new AliAODConversionPhoton(tmpvec2);
568  if(!PhotonCandidate2){
569  delete clusSub2;
570  delete tmpvec2;
571  return -9;
572  }
573  // Flag Photon as CaloPhoton
574  PhotonCandidate2->SetIsCaloPhoton();
575  Int_t mclabel = -3;
576  // create pi0 candidate
577  AliAODConversionMother *pi0cand = new AliAODConversionMother(PhotonCandidate1,PhotonCandidate2);
578  if((((AliConversionMesonCuts*)fMesonCuts)->MesonIsSelected(pi0cand,kTRUE,fEventCuts->GetEtaShift()))){
579  if(fIsMC> 0 && PhotonCandidate && PhotonCandidate1 && PhotonCandidate2 && fSaveMCInformation){
580  // if(fInputEvent->IsA()==AliESDEvent::Class())
581  mclabel = ProcessTrueClusterCandidates(PhotonCandidate, clus->GetM02(), PhotonCandidate1, PhotonCandidate2);
582  // if(fInputEvent->IsA()==AliAODEvent::Class())
583  // ProcessTrueClusterCandidatesAOD(PhotonCandidate, clus->GetM02(), PhotonCandidate1, PhotonCandidate2);
584  return mclabel;
585  }
586  } else {
587  return -7;
588  }
589  return -1;
590 }
591 
592 
593 
594 //________________________________________________________________________
596  AliAODConversionPhoton *TrueSubClusterCandidate1,
597  AliAODConversionPhoton *TrueSubClusterCandidate2)
598 {
599  Int_t mcLabelCluster = -5;
600  const AliVVertex* primVtxMC = fMCEvent->GetPrimaryVertex();
601  Double_t mcProdVtxX = primVtxMC->GetX();
602  Double_t mcProdVtxY = primVtxMC->GetY();
603  Double_t mcProdVtxZ = primVtxMC->GetZ();
604 
605  TParticle *Photon = NULL;
606  if (!TrueClusterCandidate->GetIsCaloPhoton()) AliFatal("CaloPhotonFlag has not been set task will abort");
607  if (TrueClusterCandidate->GetCaloPhotonMCLabel(0) < 0) return mcLabelCluster;
608  if (TrueClusterCandidate->GetNCaloPhotonMCLabels()>0) Photon = fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMCLabel(0));
609  else return mcLabelCluster;
610 
611  if(Photon == NULL){
612  return mcLabelCluster;
613  }
614  AliAODConversionMother *mesoncand = new AliAODConversionMother(TrueSubClusterCandidate1,TrueSubClusterCandidate2);
615  Bool_t mesonIsSelected = (((AliConversionMesonCuts*)fMesonCuts)->MesonIsSelected(mesoncand,kTRUE,fEventCuts->GetEtaShift()));
616  if (!mesonIsSelected){
617  delete mesoncand;
618  return mcLabelCluster;
619  }
620 
621  TrueClusterCandidate->SetCaloPhotonMCFlags(fMCEvent, kFALSE);
622 
623  Int_t clusterClass = 0;
624  Bool_t isPrimary = fEventCuts->IsConversionPrimaryESD( fMCEvent, TrueClusterCandidate->GetCaloPhotonMCLabel(0), mcProdVtxX, mcProdVtxY, mcProdVtxZ);
625 
626  // cluster classification:
627  // 1 - nice merged cluster (2 gamma | contributions from 2 gamma) from pi0/eta
628  // 2 - contribution from only 1 partner (1 gamma, 1 fully coverted gamma) from pi0/eta
629  // 3 - contribution from part of 1 partner (1 electron) from pi0/eta
630  Long_t motherLab = -1;
631  if (TrueClusterCandidate->IsMerged() || TrueClusterCandidate->IsMergedPartConv()){
632  clusterClass = 1;
633  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
634  } else if (TrueClusterCandidate->GetNCaloPhotonMotherMCLabels()> 0){
635  if (TrueClusterCandidate->IsLargestComponentElectron() || TrueClusterCandidate->IsLargestComponentPhoton()){
636  if (TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0) > -1 && (fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0))->GetPdgCode() == 111 || fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0))->GetPdgCode() == 221) ){
637  if ( TrueClusterCandidate->IsConversion() && !TrueClusterCandidate->IsConversionFullyContained() ){
638  clusterClass = 3;
639  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
640  } else {
641  clusterClass = 2;
642  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
643  }
644  }
645  } else if (TrueClusterCandidate->IsSubLeadingEM()){
646  if (TrueClusterCandidate->GetNCaloPhotonMotherMCLabels()> 1){
647  if ( TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1) > -1){
648  if (fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1))->GetPdgCode() == 111 || fMCEvent->Particle(TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1))->GetPdgCode() == 221 ){
649  clusterClass = 2;
650  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(1);
651  }
652  }
653  }
654  } else {
655  motherLab = TrueClusterCandidate->GetCaloPhotonMotherMCLabel(0);
656  }
657  }
658 
659  // Get Mother particle
660  TParticle *mother = NULL;
661  Int_t motherPDG = -1;
662  if (motherLab > -1)
663  mother = fMCEvent->Particle(motherLab);
664  if (mother)
665  motherPDG = TMath::Abs(mother->GetPdgCode());
666 
667  //
668  if (clusterClass == 1 || clusterClass == 2 || clusterClass == 3 ){
669  // separate different components
670  if (clusterClass == 1 && TrueClusterCandidate->IsMerged()){
671  if (motherPDG == 111){
672  mcLabelCluster = 10; // NOTE merged pi0
673  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
674  mcLabelCluster = 11; // NOTE secondary merged pi0
675  }
676  if (motherPDG == 221)
677  mcLabelCluster = 12; // NOTE merged eta
678  } else if (clusterClass == 1 && TrueClusterCandidate->IsMergedPartConv()){
679  if (motherPDG == 111){
680  mcLabelCluster = 13; // NOTE merged pi0 with one converted gamma
681  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
682  mcLabelCluster = 14; // NOTE merged secondary pi0 with one converted gamma
683  }
684  if (motherPDG == 221)
685  mcLabelCluster = 15; // NOTE merged eta with one converted gamma
686  } else if (clusterClass == 2){
687  if (motherPDG == 111){
688  mcLabelCluster = 20; // NOTE decay photon from pi0
689  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
690  mcLabelCluster = 21; // NOTE decay photon from secondary pi0
691  }
692  if (motherPDG == 221)
693  mcLabelCluster = 22; // NOTE decay photon from eta
694  } else if (clusterClass == 3){
695  if (motherPDG == 111) {
696  mcLabelCluster = 30; // NOTE electron from decayed pi0
697  if (!isPrimary && m02 >= 0 && m02 <= 4.8 )
698  mcLabelCluster = 31; // NOTE electron from decayed secondary pi0
699  }
700  if (motherPDG == 221)
701  mcLabelCluster = 32; // NOTE electron from decayed eta
702  }
703  // leading particle is a photon or the conversion is fully contained and its not from pi0 || eta
704  } else if (TrueClusterCandidate->IsLargestComponentPhoton() || TrueClusterCandidate->IsConversionFullyContained()
705  || TrueClusterCandidate->IsElectronFromFragPhoton()){
706  if (motherLab == -1)
707  mcLabelCluster = 40; // NOTE direct photon
708  else
709  mcLabelCluster = 41;
710  // leading particle is an electron and its not from pi0 || eta and no electron from fragmentation photon conversion
711  } else if (TrueClusterCandidate->IsLargestComponentElectron() && !TrueClusterCandidate->IsElectronFromFragPhoton()){
712  mcLabelCluster = 50; // NOTE cluster from hadron
713  } else {
714  mcLabelCluster = 60; // NOTE BG cluster
715  }
716 
717  delete mesoncand;
718  return mcLabelCluster;
719 }
void SetPeriodEnum(TString periodName)
void SetCaloClusterRef(Long_t ref)
double Double_t
Definition: External.C:58
Float_t fBuffer_ClusterPhi
! array buffer
Int_t * fBuffer_Cells_RelativeColumn
! array buffer
virtual void UserExec(Option_t *option)
void SetCaloPhotonMCFlags(AliMCEvent *mcEvent, Bool_t enableSort)
Float_t fBuffer_Event_Vertex_Y
! array buffer
Int_t * fBuffer_Cells_RelativeRow
! array buffer
Int_t MakePhotonCandidates(AliVCluster *clus, AliVCaloCells *cells, Long_t indexCluster)
PeriodVar GetPeriodEnum()
Bool_t fSaveMCInformation
save MC information
void SetPeriodEnumExplicit(PeriodVar periodEnum)
Float_t fBuffer_ClusterE
! array buffer
Double_t GetEtaShift()
Bool_t fSaveEventProperties
save general event properties (centrality etc.)
TString GetPeriodName()
void SetCaloPhotonMCLabel(Int_t i, Int_t labelCaloPhoton)
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)
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
Int_t fNSurroundingCells
save MC information
Int_t IsEventAcceptedByCut(AliConvEventCuts *ReaderCuts, AliVEvent *event, AliMCEvent *mcEvent, Int_t isHeavyIon, Bool_t isEMCALAnalysis)
int Int_t
Definition: External.C:63
Class handling all kinds of selection cuts for Gamma Calo analysis.
void ProcessQATreeCluster(AliVEvent *event, AliVCluster *cluster, Long_t indexCluster)
float Float_t
Definition: External.C:68
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
Int_t * fBuffer_Surrounding_Cells_RelativeRow
! array buffer
Int_t fBuffer_LeadingCell_Row
! array buffer
Bool_t fSaveSurroundingCells
save arrays of all cells in event
Int_t fBuffer_LeadingCell_Column
! array buffer
Float_t fBuffer_ClusterM20
! array buffer
Int_t * fBuffer_Cells_ID
! array buffer
Int_t fMaxNLMCut
save MC information
short Short_t
Definition: External.C:23
Float_t fBuffer_ClusterM02
! array buffer
virtual void Terminate(Option_t *)
Int_t * fBuffer_Surrounding_Cells_RelativeColumn
! array buffer
Int_t fBuffer_Event_NumActiveCells
! array buffer
Float_t fBuffer_ClusterEta
! array buffer
Bool_t fSaveCells
save arrays of cluster cells
Class handling all kinds of selection cuts for Gamma Conversion analysis.
Int_t fBuffer_ClusterNumCells
! array buffer
void SetNCaloPhotonMCLabels(Int_t nLabels)
const Int_t kMaxActiveCells
Class handling all kinds of selection cuts for Gamma Conversion analysis.
AliConversionMesonCuts * fMesonCuts
AliConvEventCuts * GetEventCuts()
Definition: AliV0ReaderV1.h:90
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
Float_t * fBuffer_Cells_E
! array buffer
Int_t fMinNLMCut
save MC information
Int_t ProcessTrueClusterCandidates(AliAODConversionPhoton *TrueClusterCandidate, Float_t m02, AliAODConversionPhoton *TrueSubClusterCandidate1, AliAODConversionPhoton *TrueSubClusterCandidate2)
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65