AliPhysics  master (3d17d9d)
AliCaloTrackReader.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // --- ROOT system ---
17 #include <TFile.h>
18 #include <TGeoManager.h>
19 #include <TStreamerInfo.h>
20 
21 // ---- ANALYSIS system ----
22 #include "AliMCEvent.h"
23 #include "AliAODMCHeader.h"
24 #include "AliGenPythiaEventHeader.h"
25 #include "AliGenEventHeader.h"
26 #include "AliESDEvent.h"
27 #include "AliAODEvent.h"
28 #include "AliVTrack.h"
29 #include "AliVParticle.h"
30 #include "AliMixedEvent.h"
31 //#include "AliTriggerAnalysis.h"
32 #include "AliESDVZERO.h"
33 #include "AliVCaloCells.h"
34 #include "AliAnalysisManager.h"
35 #include "AliInputEventHandler.h"
36 #include "AliAODMCParticle.h"
37 #include "AliMultSelection.h"
38 
39 // ---- Detectors ----
40 #include "AliPHOSGeoUtils.h"
41 #include "AliEMCALGeometry.h"
42 #include "AliEMCALRecoUtils.h"
43 
44 // ---- CaloTrackCorr ---
45 #include "AliCalorimeterUtils.h"
46 #include "AliCaloTrackReader.h"
47 #include "AliMCAnalysisUtils.h"
48 
49 // ---- Jets ----
50 #include "AliAODJet.h"
51 //#include "AliAODJetEventBackground.h"
52 
54 ClassImp(AliCaloTrackReader) ;
56 
57 //________________________________________
59 //________________________________________
61 TObject(), fEventNumber(-1), //fCurrentFileName(""),
62 fDataType(0), fDebug(0),
63 fFiducialCut(0x0), fCheckFidCut(kFALSE),
64 fComparePtHardAndJetPt(0), fPtHardAndJetPtFactor(0),
65 fComparePtHardAndClusterPt(0),fPtHardAndClusterPtFactor(0),
66 fCTSPtMin(0), fEMCALPtMin(0), fPHOSPtMin(0),
67 fCTSPtMax(0), fEMCALPtMax(0), fPHOSPtMax(0),
68 fEMCALBadChMinDist(0), fPHOSBadChMinDist (0),
69 fEMCALNCellsCut(0), fPHOSNCellsCut(0),
70 fEMCALHighEnergyNdiffCut(0), fEMCALMinCellEnNdiffCut(0),
71 fUseEMCALTimeCut(1), fUseParamTimeCut(0),
72 fUseTrackTimeCut(0), fAccessTrackTOF(0),
73 fEMCALTimeCutMin(-10000), fEMCALTimeCutMax(10000),
74 fEMCALParamTimeCutMin(), fEMCALParamTimeCutMax(),
75 fTrackTimeCutMin(-10000), fTrackTimeCutMax(10000),
76 fUseTrackDCACut(0),
77 fAODBranchList(0x0),
78 fCTSTracks(0x0), fEMCALClusters(0x0),
79 fDCALClusters(0x0), fPHOSClusters(0x0),
80 fEMCALCells(0x0), fPHOSCells(0x0),
81 fInputEvent(0x0), fOutputEvent(0x0), fMC(0x0),
82 fFillCTS(0), fFillEMCAL(0),
83 fFillDCAL(0), fFillPHOS(0),
84 fFillEMCALCells(0), fFillPHOSCells(0),
85 fRecalculateClusters(kFALSE),fCorrectELinearity(kTRUE),
86 fSelectEmbeddedClusters(kFALSE),
87 fScaleEPerSM(kFALSE),
88 fSmearShowerShape(0), fSmearShowerShapeWidth(0), fRandom(),
89 fSmearingFunction(0), fSmearNLMMin(0), fSmearNLMMax(0),
90 fTrackStatus(0), fSelectSPDHitTracks(0),
91 fTrackMultNPtCut(0), fTrackMultEtaCut(0.9),
92 fDeltaAODFileName(""), fFiredTriggerClassName(""),
93 
94 fEventTriggerMask(0), fMixEventTriggerMask(0), fEventTriggerAtSE(0),
95 fEventTrigMinBias(0), fEventTrigCentral(0),
96 fEventTrigSemiCentral(0), fEventTrigEMCALL0(0),
97 fEventTrigEMCALL1Gamma1(0), fEventTrigEMCALL1Gamma2(0),
98 fEventTrigEMCALL1Jet1(0), fEventTrigEMCALL1Jet2(0),
99 fBitEGA(0), fBitEJE(0),
100 
101 fEventType(-1),
102 fTaskName(""), fCaloUtils(0x0), fMCUtils(0x0),
103 fWeightUtils(0x0), fEventWeight(1),
104 fMixedEvent(NULL), fNMixedEvent(0), fVertex(NULL),
105 fListMixedTracksEvents(), fListMixedCaloEvents(),
106 fLastMixedTracksEvent(-1), fLastMixedCaloEvent(-1),
107 fWriteOutputDeltaAOD(kFALSE),
108 fEMCALClustersListName(""), fEMCALCellsListName(""),
109 fZvtxCut(0.),
110 fAcceptFastCluster(kFALSE), fRemoveLEDEvents(0),
111 fLEDHighEnergyCutSM(0), fLEDHighNCellsCutSM(0),
112 fLEDLowEnergyCutSM3(0), fLEDLowNCellsCutSM3(0),
113 fLEDMinCellEnergy(0), fLEDMaxCellEnergy(0),
114 fRemoveLEDStripEvents(0), fLEDEventMaxNumberOfStrips(0),
115 fLEDLowEnergyCutSM3Strip(0), fLEDLowNCellsCutSM3Strip(0),
116 
117 //Trigger rejection
118 fRemoveBadTriggerEvents(0), fTriggerPatchClusterMatch(0),
119 fTriggerPatchTimeWindow(), fTriggerL0EventThreshold(0),
120 fTriggerL1EventThreshold(0), fTriggerL1EventThresholdFix(0),
121 fTriggerClusterBC(0), fTriggerClusterIndex(0), fTriggerClusterId(0),
122 fIsExoticEvent(0), fIsBadCellEvent(0), fIsBadMaxCellEvent(0),
123 fIsTriggerMatch(0), fIsTriggerMatchOpenCut(),
124 fTriggerClusterTimeRecal(kTRUE), fRemoveUnMatchedTriggers(kTRUE),
125 fDoPileUpEventRejection(kFALSE), fDoV0ANDEventSelection(kFALSE),
126 fDoVertexBCEventSelection(kFALSE),
127 fDoRejectNoTrackEvents(kFALSE),
128 fUseEventsWithPrimaryVertex(kFALSE),
129 //fTriggerAnalysis (0x0),
130 fTimeStampEventSelect(0),
131 fTimeStampEventFracMin(0), fTimeStampEventFracMax(0),
132 fTimeStampRunMin(0), fTimeStampRunMax(0),
133 fTimeStampEventCTPBCCorrExclude(0),
134 fTimeStampEventCTPBCCorrMin(0), fTimeStampEventCTPBCCorrMax(0),
135 fNPileUpClusters(-1), fNNonPileUpClusters(-1), fNPileUpClustersCut(3),
136 fVertexBC(-200), fRecalculateVertexBC(0),
137 fUseAliCentrality(0), fMultWithEventSel(0),
138 fCentralityClass(""), fCentralityOpt(0),
139 fEventPlaneMethod(""),
140 fFillInputNonStandardJetBranch(kFALSE),
141 fNonStandardJets(new TClonesArray("AliAODJet",100)), fInputNonStandardJetBranchName("jets"),
142 fFillInputBackgroundJetBranch(kFALSE),
143 //fBackgroundJets(0x0),
144 fBackgroundJets(new TClonesArray("AliAODJet",100)),
145 fInputBackgroundJetBranchName("jets"),
146 fAcceptEventsWithBit(0), fRejectEventsWithBit(0), fRejectEMCalTriggerEventsWith2Tresholds(0),
147 fMomentum(), fParRun(kFALSE), fCurrentParIndex(0),
148 fOutputContainer(0x0), fhEMCALClusterEtaPhi(0), fhEMCALClusterEtaPhiFidCut(0),
149 fhEMCALClusterTimeE(0),
150 fhEMCALNSumEnCellsPerSM(0), fhEMCALNSumEnCellsPerSMAfter(0), fhEMCALNSumEnCellsPerSMAfterStripCut(0),
151 fhEMCALNSumEnCellsPerStrip(0), fhEMCALNSumEnCellsPerStripAfter(0),
152 fEnergyHistogramNbins(0),
153 fhNEventsAfterCut(0), fNMCGenerToAccept(0), fMCGenerEventHeaderToAccept(""),
154 fGenEventHeader(0), fGenPythiaEventHeader(0)
155 {
156  for(Int_t i = 0; i < 8; i++) fhEMCALClusterCutsE [i]= 0x0 ;
157  for(Int_t i = 0; i < 7; i++) fhPHOSClusterCutsE [i]= 0x0 ;
158  for(Int_t i = 0; i < 6; i++) fhCTSTrackCutsPt [i]= 0x0 ;
159  for(Int_t j = 0; j < 5; j++) { fMCGenerToAccept [j] = ""; fMCGenerIndexToAccept[j] = -1; }
160 
161  InitParameters();
162 }
163 
164 //_______________________________________
166 //_______________________________________
168 {
169  DeletePointers();
170 }
171 
172 //_______________________________________
175 //_______________________________________
177 {
178  delete fFiducialCut ;
179 
180  if(fAODBranchList)
181  {
182  fAODBranchList->Delete();
183  delete fAODBranchList ;
184  }
185 
186  if(fCTSTracks)
187  {
188  if(fDataType!=kMC)fCTSTracks->Clear() ;
189  else fCTSTracks->Delete() ;
190  delete fCTSTracks ;
191  }
192 
193  if(fEMCALClusters)
194  {
195  if(fDataType!=kMC)fEMCALClusters->Clear("C") ;
196  else fEMCALClusters->Delete() ;
197  delete fEMCALClusters ;
198  }
199 
200  if(fDCALClusters)
201  {
202  if(fDataType!=kMC)fDCALClusters->Clear("C") ;
203  else fDCALClusters->Delete() ;
204  delete fDCALClusters ;
205  }
206 
207  if(fPHOSClusters)
208  {
209  if(fDataType!=kMC)fPHOSClusters->Clear("C") ;
210  else fPHOSClusters->Delete() ;
211  delete fPHOSClusters ;
212  }
213 
214  if(fVertex)
215  {
216  for (Int_t i = 0; i < fNMixedEvent; i++)
217  {
218  delete [] fVertex[i] ;
219 
220  }
221  delete [] fVertex ;
222  }
223 
224  //delete fTriggerAnalysis;
225 
226  if(fNonStandardJets)
227  {
228  if(fDataType!=kMC) fNonStandardJets->Clear("C") ;
229  else fNonStandardJets->Delete() ;
230  delete fNonStandardJets ;
231  }
232 
233  if(fBackgroundJets)
234  {
235  if(fDataType!=kMC) fBackgroundJets->Clear("C") ;
236  else fBackgroundJets->Delete() ;
237  delete fBackgroundJets;
238  }
239  // delete fBackgroundJets ;
240 
241  fRejectEventsWithBit.Reset();
242  fAcceptEventsWithBit.Reset();
243 
244  if ( fWeightUtils ) delete fWeightUtils ;
245 
246  if ( fMCUtils ) delete fMCUtils ;
247 
248  // Pointers not owned, done by the analysis frame
249  // if(fInputEvent) delete fInputEvent ;
250  // if(fOutputEvent) delete fOutputEvent ;
251  // if(fMC) delete fMC ;
252  // Pointer not owned, deleted by maker
253  // if (fCaloUtils) delete fCaloUtils ;
254 }
255 
256 //____________________________________________________________
260 //____________________________________________________________
262 {
263  Float_t cut = fTrackDCACut[0]+fTrackDCACut[1]/TMath::Power(pt,fTrackDCACut[2]);
264 
265  if(TMath::Abs(dca) < cut)
266  return kTRUE;
267  else
268  return kFALSE;
269 }
270 
271 //_____________________________________________________
274 //_____________________________________________________
276 {
277  Int_t nAccept = fAcceptEventsWithBit.GetSize();
278 
279  //printf("N accept %d\n", nAccept);
280 
281  if( nAccept <= 0 )
282  return kTRUE ; // accept the event
283 
284  UInt_t trigFired = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
285 
286  for(Int_t ibit = 0; ibit < nAccept; ibit++)
287  {
288  Bool_t accept = (trigFired & fAcceptEventsWithBit.At(ibit));
289 
290  //printf("accept %d, ibit %d, bit %d \n",accept, ibit,fAcceptEventsWithBit.At(ibit));
291  if(accept) return kTRUE ; // accept the event
292  }
293 
294  return kFALSE ; // reject the event
295 }
296 
297 //_____________________________________________________
302 //_____________________________________________________
304 {
305  if( !fMC || fNMCGenerToAccept <= 0 ) return kTRUE;
306 
307  TString genName;
308  Int_t genIndex;
309  genIndex = GetCocktailGeneratorAndIndex(mcLabel, genName);
310  //fMC->GetCocktailGenerator(mcLabel,genName);
311 
312  Bool_t generOK = kFALSE;
313  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
314  {
315  if ( fMCGenerToAccept[ig].Contains(genName) ) generOK = kTRUE;
316 
317  if ( generOK && fMCGenerIndexToAccept[ig] >= 0 && fMCGenerToAccept[ig] != genIndex) generOK = kFALSE;
318  }
319 
320  if ( !generOK ) AliDebug(1, Form("skip label %d, gen %s",mcLabel,genName.Data()) );
321 
322  return generOK;
323 }
324 
325 //_____________________________________________________________________
333 //_____________________________________________________________________
335 {
336  //method that gives the generator for a given particle with label index (or that of the corresponding primary)
337  AliVParticle* mcpart0 = (AliVParticle*) GetMC()->GetTrack(index);
338  Int_t genIndex = -1;
339 
340  if(!mcpart0)
341  {
342  printf("AliMCEvent-BREAK: No valid AliMCParticle at label %i\n",index);
343  return -1;
344  }
345 
346  nameGen = GetGeneratorNameAndIndex(index,genIndex);
347 
348  if(nameGen.Contains("nococktailheader") ) return -1;
349 
350  Int_t lab=index;
351 
352  while(nameGen.IsWhitespace())
353  {
354  AliVParticle* mcpart = (AliVParticle*) GetMC()->GetTrack(lab);
355 
356  if(!mcpart)
357  {
358  printf("AliMCEvent-BREAK: No valid AliMCParticle at label %i\n",lab);
359  break;
360  }
361 
362  Int_t mother=0;
363  mother = mcpart->GetMother();
364 
365  if(mother<0)
366  {
367  printf("AliMCEvent - BREAK: Reached primary particle without valid mother\n");
368  break;
369  }
370 
371  AliVParticle* mcmom = (AliVParticle*) GetMC()->GetTrack(mother);
372  if(!mcmom)
373  {
374  printf("AliMCEvent-BREAK: No valid AliMCParticle mother at label %i\n",mother);
375  break;
376  }
377 
378  lab=mother;
379 
380  nameGen = GetGeneratorNameAndIndex(mother,genIndex);
381  }
382 
383  return genIndex;
384 }
385 //_____________________________________________________________________
393 //_____________________________________________________________________
395 {
396  Int_t nsumpart = GetMC()->GetNumberOfPrimaries();
397 
398  genIndex = -1;
399 
400  TList* lh = GetMC()->GetCocktailList();
401  if(!lh)
402  {
403  TString noheader="nococktailheader";
404  return noheader;
405  }
406 
407  Int_t nh = lh->GetEntries();
408 
409  for (Int_t i = nh-1; i >= 0; i--)
410  {
411  AliGenEventHeader* gh = (AliGenEventHeader*)lh->At(i);
412 
413  TString genname = gh->GetName();
414 
415  Int_t npart=gh->NProduced();
416 
417  if (i == 0) npart = nsumpart;
418 
419  if(index < nsumpart && index >= (nsumpart-npart))
420  {
421  genIndex = i ;
422  return genname;
423  }
424 
425  nsumpart-=npart;
426  }
427 
428  TString empty="";
429  return empty;
430 }
431 
432 
433 //_____________________________________________________
436 //_____________________________________________________
438 {
439  Int_t nReject = fRejectEventsWithBit.GetSize();
440 
441  //printf("N reject %d\n", nReject);
442 
443  if( nReject <= 0 )
444  return kTRUE ; // accept the event
445 
446  UInt_t trigFired = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
447 
448  for(Int_t ibit = 0; ibit < nReject; ibit++)
449  {
450  Bool_t reject = (trigFired & fRejectEventsWithBit.At(ibit));
451 
452  //printf("reject %d, ibit %d, bit %d \n",reject, ibit,fRejectEventsWithBit.At(ibit));
453  if(reject) return kFALSE ; // reject the event
454  }
455 
456  return kTRUE ; // accept the event
457 }
458 
459 //_____________________________________________
463 //_____________________________________________
465 {
466  //-----------------------------------------------------------
467  // Reject events depending on the event species type
468  //-----------------------------------------------------------
469 
470  // Event types:
471  // kStartOfRun = 1, // START_OF_RUN
472  // kEndOfRun = 2, // END_OF_RUN
473  // kStartOfRunFiles = 3, // START_OF_RUN_FILES
474  // kEndOfRunFiles = 4, // END_OF_RUN_FILES
475  // kStartOfBurst = 5, // START_OF_BURST
476  // kEndOfBurst = 6, // END_OF_BURST
477  // kPhysicsEvent = 7, // PHYSICS_EVENT
478  // kCalibrationEvent = 8, // CALIBRATION_EVENT
479  // kFormatError = 9, // EVENT_FORMAT_ERROR
480  // kStartOfData = 10, // START_OF_DATA
481  // kEndOfData = 11, // END_OF_DATA
482  // kSystemSoftwareTriggerEvent = 12, // SYSTEM_SOFTWARE_TRIGGER_EVENT
483  // kDetectorSoftwareTriggerEvent = 13 // DETECTOR_SOFTWARE_TRIGGER_EVENT
484 
485  Int_t eventType = 0;
486  if(fInputEvent->GetHeader())
487  eventType = ((AliVHeader*)fInputEvent->GetHeader())->GetEventType();
488 
489  AliDebug(3,Form("Event type %d",eventType));
490 
491  // Select only Physics events in data, eventType = 7,
492  // usually done by physics selection
493  // MC not set, eventType =0, I think, so do not apply a value to fEventType by default
494  // LED events have eventType = 8, implemented a selection cut in the past
495  // but not really useful for EMCal data since LED are not reconstructed (unless wrongly assumed as physics)
496 
497  if ( fEventType >= 0 && eventType != fEventType ) return kFALSE ;
498 
499  AliDebug(1,"Pass event species selection");
500 
501  fhNEventsAfterCut->Fill(1.5);
502 
503  //-----------------------------------------------------------------
504  // In case of mixing analysis, select here the trigger of the event
505  //-----------------------------------------------------------------
506 
507  UInt_t isTrigger = kFALSE;
508  UInt_t isMB = kFALSE;
509 
510  if(!fEventTriggerAtSE)
511  {
512  // In case of mixing analysis, accept MB events, not only Trigger
513  // Track and cluster arrays filled for MB in order to create the pool in the corresponding analysis
514  // via de method in the base class FillMixedEventPool()
515 
516  AliAnalysisManager *manager = AliAnalysisManager::GetAnalysisManager();
517  AliInputEventHandler *inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
518 
519  if(!inputHandler) return kFALSE ; // to content coverity
520 
521  isTrigger = inputHandler->IsEventSelected() & fEventTriggerMask;
522  isMB = inputHandler->IsEventSelected() & fMixEventTriggerMask;
523 
524  if(!isTrigger && !isMB) return kFALSE;
525 
526  //printf("Selected triggered event : %s\n",GetFiredTriggerClasses().Data());
527  AliDebug(1,"Pass uninteresting triggered events rejection in case of mixing analysis");
528 
529  fhNEventsAfterCut->Fill(2.5);
530  }
531 
532  //-----------------------------------------------------------
533  // Reject events depending on the trigger name
534  // Careful!, if a special MB event string is selected but the option
535  // to select events via the mask in the reader is done, it will not
536  // be taken into account.
537  //-----------------------------------------------------------
538 
539  AliDebug(1,Form("FiredTriggerClass <%s>, selected class <%s>, compare name %d",
542 
543  if ( fFiredTriggerClassName != "" && !isMB )
544  {
546  return kFALSE;
547 
548  AliDebug(1,"Accepted triggered event");
549 
550  fhNEventsAfterCut->Fill(3.5);
551  }
552 
553  //-------------------------------------------------------------------------------------
554  // Reject or accept events depending on the trigger bit
555  //-------------------------------------------------------------------------------------
556 
559 
560  //printf("AliCaloTrackReader::FillInputEvent() - Accept event? %d, Reject event %d? \n",okA,okR);
561 
562  if(!okA || !okR) return kFALSE;
563 
564  AliDebug(1,"Pass event bit rejection");
565 
566  fhNEventsAfterCut->Fill(4.5);
567 
568  //----------------------------------------------------------------------
569  // Do not count events that were likely triggered by an exotic cluster
570  // or out BC cluster
571  //----------------------------------------------------------------------
572 
573  // Set a bit with the event kind, MB, L0, L1 ...
575 
576  // In case of Mixing, avoid checking the triggers in the min bias events
577  if(!fEventTriggerAtSE && (isMB && !isTrigger)) return kTRUE;
578 
580  {
582  {
583  // Reject triggered events when there is coincidence on both EMCal trigger thresholds,
584  // but the requested trigger is the low trigger threshold
585  if(IsEventEMCALL1Jet1 () && IsEventEMCALL1Jet2 () && fFiredTriggerClassName.Contains("EJ2")) return kFALSE;
586  if(IsEventEMCALL1Gamma1() && IsEventEMCALL1Gamma2() && fFiredTriggerClassName.Contains("EG2")) return kFALSE;
587  }
588 
589  //Get Patches that triggered
591 
592  MatchTriggerCluster(patches);
593 
594  patches.Reset();
595 
596  // If requested, remove badly triggeed events, but only when the EMCal trigger bit is set
598  {
599  AliDebug(1,Form("ACCEPT triggered event? \n exotic? %d - bad cell %d - bad Max cell %d - BC %d - Matched %d\n",
601 
602  if (fIsExoticEvent) return kFALSE;
603  else if(fIsBadCellEvent) return kFALSE;
604  else if(fRemoveUnMatchedTriggers && !fIsTriggerMatch) return kFALSE ;
605  else if(fTriggerClusterBC != 0) return kFALSE;
606 
607  AliDebug(1,Form("\t *** YES for %s",GetFiredTriggerClasses().Data()));
608  }
609 
610  AliDebug(1,"Pass EMCal triggered event rejection \n");
611 
612  fhNEventsAfterCut->Fill(5.5);
613  }
614 
615  //-------------------------------------------------------------------------------------
616  // Select events only fired by a certain trigger configuration if it is provided
617  //-------------------------------------------------------------------------------------
618 
619  if (GetFiredTriggerClasses().Contains("FAST") && !GetFiredTriggerClasses().Contains("ALL") && !fAcceptFastCluster)
620  {
621  AliDebug(1,Form("Do not count events from fast cluster, trigger name %s\n",fFiredTriggerClassName.Data()));
622  return kFALSE;
623 
624  fhNEventsAfterCut->Fill(6.5);
625  }
626 
627  //-------------------------------------------------------------------------------------
628  // Reject event if large clusters with large energy
629  // Use only for LHC11a data for the moment, and if input is clusterizer V1 or V1+unfolding
630  // If clusterzer NxN or V2 it does not help
631  //-------------------------------------------------------------------------------------
632 
633  //Int_t run = fInputEvent->GetRunNumber();
634 
635  if ( fRemoveLEDEvents )
636  {
638 
639  if(reject) return kFALSE;
640 
641  AliDebug(1,"Pass LED event rejection");
642 
643  fhNEventsAfterCut->Fill(7.5);
644  } // Remove LED events
645 
646  // All selection criteria passed, accept the event
647  return kTRUE ;
648 }
649 
650 //________________________________________________
659 //________________________________________________
661 {
662  if ( !fGenEventHeader )
663  {
664  AliError("Skip event, event header is not available!");
665  return kFALSE;
666  }
667 
668  if ( fGenPythiaEventHeader )
669  {
670  // Do this check only for jet-jet and gamma-jet productions
671 
672  if ( !processName.Contains("Jet") ) return kTRUE;
673 
674  Int_t nTriggerJets = fGenPythiaEventHeader->NTriggerJets();
675  Float_t ptHard = fGenPythiaEventHeader->GetPtHard();
676  TParticle * jet = 0;
677 
678  AliDebug(1,Form("Njets: %d, pT Hard %f",nTriggerJets, ptHard));
679 
680  Float_t tmpjet[]={0,0,0,0};
681  for(Int_t ijet = 0; ijet< nTriggerJets; ijet++)
682  {
683  fGenPythiaEventHeader->TriggerJet(ijet, tmpjet);
684  jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
685 
686  AliDebug(1,Form("jet %d; pycell jet pT %f",ijet, jet->Pt()));
687 
688  //Compare jet pT and pt Hard
689  if(jet->Pt() > fPtHardAndJetPtFactor * ptHard)
690  {
691  AliInfo(Form("Reject jet event with : process %d <%s>, pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f",
692  process, processName.Data(), ptHard, jet->Pt(), fPtHardAndJetPtFactor));
693  return kFALSE;
694  }
695  } // jet loop
696 
697  if(jet) delete jet;
698  } // pythia header
699 
700  return kTRUE ;
701 }
702 
703 //____________________________________________________
712 //____________________________________________________
714 {
715  if ( !fGenEventHeader )
716  {
717  AliError("Skip event, event header is not available!");
718  return kFALSE;
719  }
720 
721  if ( fGenPythiaEventHeader )
722  {
723  // Do this check only for gamma-jet productions
724  if(processName !="Gamma-Jet") return kTRUE;
725 
726  Float_t ptHard = fGenPythiaEventHeader->GetPtHard();
727 
728  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
729  for (Int_t iclus = 0; iclus < nclusters; iclus++)
730  {
731  AliVCluster * clus = fInputEvent->GetCaloCluster(iclus) ;
732  Float_t ecluster = clus->E();
733 
734  if(ecluster > fPtHardAndClusterPtFactor*ptHard)
735  {
736  AliInfo(Form("Reject : process %d <%s>, ecluster %2.2f, calo %d, factor %2.2f, ptHard %f",
737  process, processName.Data(), ecluster ,clus->GetType(), fPtHardAndClusterPtFactor,ptHard));
738 
739  return kFALSE;
740  }
741  } // cluster loop
742 
743  } // pythia header
744 
745  return kTRUE ;
746 }
747 
748 //___________________________________________________
751 //___________________________________________________
753 {
754  fhNEventsAfterCut = new TH1I("hNEventsAfterCut", "Number of analyzed events", 20, 0, 20) ;
755  //fhNEventsAfterCut->SetXTitle("Selection");
756  fhNEventsAfterCut->SetYTitle("# events");
757  fhNEventsAfterCut->GetXaxis()->SetBinLabel(1 ,"1=Input");
758  fhNEventsAfterCut->GetXaxis()->SetBinLabel(2 ,"2=Event Type");
759  fhNEventsAfterCut->GetXaxis()->SetBinLabel(3 ,"3=Mixing Event");
760  fhNEventsAfterCut->GetXaxis()->SetBinLabel(4 ,"4=Trigger string");
761  fhNEventsAfterCut->GetXaxis()->SetBinLabel(5 ,"5=Trigger Bit");
762  fhNEventsAfterCut->GetXaxis()->SetBinLabel(6 ,"6=Good EMC Trigger");
763  fhNEventsAfterCut->GetXaxis()->SetBinLabel(7 ,"7=!Fast Cluster");
764  fhNEventsAfterCut->GetXaxis()->SetBinLabel(8 ,"8=!LED");
765  fhNEventsAfterCut->GetXaxis()->SetBinLabel(9 ,"9=Time stamp");
766  fhNEventsAfterCut->GetXaxis()->SetBinLabel(10,"10=Primary vertex");
767  fhNEventsAfterCut->GetXaxis()->SetBinLabel(11,"11=Null 3 vertex");
768  fhNEventsAfterCut->GetXaxis()->SetBinLabel(12,"12=Z vertex window");
769  fhNEventsAfterCut->GetXaxis()->SetBinLabel(13,"13=Pile-up");
770  fhNEventsAfterCut->GetXaxis()->SetBinLabel(14,"14=V0AND");
771  fhNEventsAfterCut->GetXaxis()->SetBinLabel(15,"15=Centrality");
772  fhNEventsAfterCut->GetXaxis()->SetBinLabel(16,"16=GenHeader");
773  fhNEventsAfterCut->GetXaxis()->SetBinLabel(17,"17=PtHard-Jet");
774  fhNEventsAfterCut->GetXaxis()->SetBinLabel(18,"18=PtHard-Cluster");
775  fhNEventsAfterCut->GetXaxis()->SetBinLabel(19,"19=N Track>0");
776  fhNEventsAfterCut->GetXaxis()->SetBinLabel(20,"20=TOF BC");
778 
779  if ( fFillEMCAL )
780  {
781  for(Int_t i = 0; i < 9; i++)
782  {
783  TString names[] =
784  { "NoCut", "Corrected", "GoodCluster", "NonLinearity",
785  "EnergyAndFidutial", "NCells", "BadDist", "Time","NcellsDiff" } ;
786 
787  fhEMCALClusterCutsE[i] = new TH1F(Form("hEMCALReaderClusterCuts_%d_%s",i,names[i].Data()),
788  Form("EMCal %d, %s",i,names[i].Data()),
790  fhEMCALClusterCutsE[i]->SetYTitle("# clusters");
791  fhEMCALClusterCutsE[i]->SetXTitle("#it{E} (GeV)");
793  }
794 
796  ("hEMCALReaderTimeE","#it{time}_{cluster} vs #it{E}_{cluster} after cuts", 250,0,250,1201,-1201,1201);
797  fhEMCALClusterTimeE->SetXTitle("#it{E}_{cluster} (GeV)");
798  fhEMCALClusterTimeE->SetYTitle("#it{time}_{cluster} (ns)");
800 
802  ("hEMCALReaderEtaPhi","#eta vs #varphi",80,-2, 2,100, 0,10);
803  // Very open limits to check problems
804  fhEMCALClusterEtaPhi->SetXTitle("#eta");
805  fhEMCALClusterEtaPhi->SetYTitle("#varphi (rad)");
807 
809  ("hEMCALReaderEtaPhiFidCut","#eta vs #varphi after fidutial cut",80,-2, 2,100, 0,10);
810  fhEMCALClusterEtaPhiFidCut->SetXTitle("#eta");
811  fhEMCALClusterEtaPhiFidCut->SetYTitle("#varphi (rad)");
813 
814  if ( fRemoveLEDEvents > 1 )
815  {
817  ("hEMCALNSumEnCellsPerSM","Total number of cells and energy in any SM",144,0,1152,250,0,5000);
818  fhEMCALNSumEnCellsPerSM->SetXTitle("#it{n}_{cells}^{SM}");
819  fhEMCALNSumEnCellsPerSM->SetYTitle("#Sigma #it{E}_{cells}^{SM} (GeV)");
821 
823  ("hEMCALNSumEnCellsPerSMAfter","Total number of cells and energy in any SM, after LED SM event rejection",144,0,1152,250,0,5000);
824  fhEMCALNSumEnCellsPerSMAfter->SetXTitle("#it{n}_{cells}^{SM}");
825  fhEMCALNSumEnCellsPerSMAfter->SetYTitle("#Sigma #it{E}_{cells}^{SM} (GeV)");
827 
828  if ( fRemoveLEDStripEvents )
829  {
831  ("hEMCALNSumEnCellsPerSMAfterStripCut","Total number of cells and energy in any SM, after LED SM and strip event rejection ",144,0,1152,250,0,5000);
832  fhEMCALNSumEnCellsPerSMAfterStripCut->SetXTitle("#it{n}_{cells}^{SM}");
833  fhEMCALNSumEnCellsPerSMAfterStripCut->SetYTitle("#Sigma #it{E}_{cells}^{SM} (GeV)");
835  }
836  }
837 
839  {
841  ("hEMCALNSumEnCellsPerStrip","Total number of cells and energy in any strip, after LED SM event rejection",48,0,48,100,0,500);
842  fhEMCALNSumEnCellsPerStrip->SetXTitle("#it{n}_{cells}^{strip}");
843  fhEMCALNSumEnCellsPerStrip->SetYTitle("#Sigma #it{E}_{cells}^{strip} (GeV)");
845 
847  ("hEMCALNSumEnCellsPerStripAfter","Total number of cells and energy in any strip, after LED SM event rejection",48,0,48,100,0,500);
848  fhEMCALNSumEnCellsPerStripAfter->SetXTitle("#it{n}_{cells}^{strip}");
849  fhEMCALNSumEnCellsPerStripAfter->SetYTitle("#Sigma #it{E}_{cells}^{strip} (GeV)");
851  }
852  }
853 
854  if(fFillPHOS)
855  {
856  for(Int_t i = 0; i < 7; i++)
857  {
858  TString names[] = {"NoCut", "ExcludeCPV", "BorderCut", "FiducialCut", "EnergyCut", "NCells", "BadDist"};
859 
860  fhPHOSClusterCutsE[i] = new TH1F(Form("hPHOSReaderClusterCuts_%d_%s",i,names[i].Data()),
861  Form("PHOS Cut %d, %s",i,names[i].Data()),
863  fhPHOSClusterCutsE[i]->SetYTitle("# clusters");
864  fhPHOSClusterCutsE[i]->SetXTitle("#it{E} (GeV)");
866  }
867  }
868 
869  if(fFillCTS)
870  {
871  for(Int_t i = 0; i < 6; i++)
872  {
873  TString names[] = {"NoCut", "Status", "ESD_AOD", "TOF", "DCA","PtAcceptanceMult"};
874 
875  fhCTSTrackCutsPt[i] = new TH1F(Form("hCTSReaderClusterCuts_%d_%s",i,names[i].Data()),
876  Form("CTS Cut %d, %s",i,names[i].Data()),
878  fhCTSTrackCutsPt[i]->SetYTitle("# tracks");
879  fhCTSTrackCutsPt[i]->SetXTitle("#it{p}_{T} (GeV)");
881  }
882  }
883 
884  return fOutputContainer ;
885 }
886 
887 //_____________________________________________________
889 //_____________________________________________________
891 {
892  TString parList ; //this will be list of parameters used for this analysis.
893 
894  const Int_t buffersize = 255;
895  char onePar[buffersize] ;
896 
897  snprintf(onePar,buffersize,"--- Reader ---:") ;
898  parList+=onePar ;
899  snprintf(onePar,buffersize,"Data type: %d; zvertex cut %2.2f; EMC cluster name: <%s> ",fDataType, fZvtxCut, fEMCALClustersListName.Data()) ;
900  parList+=onePar ;
901  snprintf(onePar,buffersize,"Use detector: EMC %d, DCA %d, PHOS %d, CTS %d, EMCcells %d, PHOScells %d ; ",
903  snprintf(onePar,buffersize,"E-pT window: EMC (%2.1f,%2.1f), PHOS (%2.1f,%2.1f), CTS (%2.1f,%2.1f); ",
905  parList+=onePar ;
906  snprintf(onePar,buffersize,"Dist to bad channel: EMC > %2.1f, PHOS > %2.1f; ",fEMCALBadChMinDist,fPHOSBadChMinDist) ;
907  parList+=onePar ;
908  snprintf(onePar,buffersize,"N cells: EMC > %d, PHOS > %d; ",fEMCALNCellsCut,fPHOSNCellsCut) ;
909  parList+=onePar ;
910  snprintf(onePar,buffersize,"N cells diff TCard: E > %2.2f, Ecell > %1.2f; ",fEMCALHighEnergyNdiffCut,fEMCALMinCellEnNdiffCut) ;
911  parList+=onePar ;
912  snprintf(onePar,buffersize,"EMC time cut single window (%2.2f,%2.2f); ",fEMCALTimeCutMin,fEMCALTimeCutMax) ;
913  parList+=onePar ;
914  snprintf(onePar,buffersize,"Check: calo fid cut %d; ",fCheckFidCut) ;
915  parList+=onePar ;
916  snprintf(onePar,buffersize,"Track: status %d, SPD hit %d; ",(Int_t) fTrackStatus, fSelectSPDHitTracks) ;
917  parList+=onePar ;
918  snprintf(onePar,buffersize,"multip. eta cut %1.1f; npt cuts %d;",fTrackMultEtaCut, fTrackMultNPtCut) ;
919  parList+=onePar ;
920 
921  if(fUseTrackDCACut)
922  {
923  snprintf(onePar,buffersize,"DCA cut ON, param (%2.4f,%2.4f,%2.4f); ",fTrackDCACut[0],fTrackDCACut[1],fTrackDCACut[2]) ;
924  parList+=onePar ;
925  }
926 
927  snprintf(onePar,buffersize,"Recalculate Clusters = %d, E linearity = %d; ",fRecalculateClusters, fCorrectELinearity) ;
928  parList+=onePar ;
929 
930  snprintf(onePar,buffersize,"SE trigger sel. %d, not? trigger Mask? %d, MB Trigger Mask for mixed %d; ",
932  parList+=onePar ;
933 
934  snprintf(onePar,buffersize,"Select fired trigger %s; Remove Bad trigger event %d, unmatched %d; Accept fastcluster %d",
936  parList+=onePar ;
937 
938  if ( fRemoveLEDEvents > 0 )
939  {
940  snprintf(onePar,buffersize,"Remove LED %d, %2.1f < Ecell < %1.2f : SM - nCell > %d, Sum E > %2.0f; SM3 - nCell < %d, Sum E < %2.0f;",
943  parList+=onePar ;
944 
945  if ( fRemoveLEDStripEvents > 0 )
946  {
947  snprintf(onePar,buffersize,"Remove LED strip? %d, with n strip > %d: "
948  "Full SM, nCell > %d, Sum E > %2.0f; "
949  "1/3 SM, nCell > %d, Sum E > %2.0f; "
950  "SM3, nCell < %d, Sum E < %2.0f;",
955  parList+=onePar ;
956  }
957  }
958 
959 
961  {
962  snprintf(onePar,buffersize,"Accept only labels from: ");
963  parList+=onePar ;
964  for(Int_t i = 0; i< fNMCGenerToAccept; i++)
965  parList+=(fMCGenerToAccept[i]+" ") ;
966 
967  snprintf(onePar,buffersize,"; ");
968  parList+=onePar ;
969  }
970 
972  {
973  snprintf(onePar,buffersize,"EMC M02 smear ON, function %d, param %2.4f, %d<=NLM<=%d; ",
975  parList+=onePar ;
976  }
977 
979  {
980  snprintf(onePar,buffersize,"jet pt / pt hard < %2.1f; ",fPtHardAndJetPtFactor);
981  parList+=onePar ;
982  }
983 
985  {
986  snprintf(onePar,buffersize,"cluster pt / pt hard < %2.2f",fPtHardAndClusterPtFactor);
987  parList+=onePar ;
988  }
989 
990  snprintf(onePar,buffersize,"Centrality: Class %s, Option %d, Bin [%d,%d]; New centrality %d; Mult. PS %d; Event plane method %s; ",
991  fCentralityClass.Data(),fCentralityOpt,fCentralityBin[0], fCentralityBin[1],
993  parList+=onePar ;
994 
995  return new TObjString(parList) ;
996 }
997 
998 
999 //______________________________________________
1001 //______________________________________________
1003 {
1004  if(fMC)
1005  {
1006  return fMC->Header();
1007  }
1008  else
1009  {
1010  AliInfo("Header is not available");
1011  return 0x0 ;
1012  }
1013 }
1014 
1015 //_________________________________________________________
1018 //_________________________________________________________
1020 {
1021  AliInfo("Input are not AODs");
1022 
1023  return NULL ;
1024 }
1025 
1026 //________________________________________________________
1029 //________________________________________________________
1030 AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader() const
1031 {
1032  AliInfo("Input are not AODs");
1033 
1034  return NULL ;
1035 }
1036 
1037 //___________________________________________________________
1045 //___________________________________________________________
1047 {
1048  Int_t vertexBC=vtx->GetBC();
1049 
1050  if(!fRecalculateVertexBC) return vertexBC;
1051 
1052  // Value not available, recalculate it.
1053 
1054  Double_t bz = fInputEvent->GetMagneticField();
1055  Bool_t bc0 = kFALSE;
1056  Int_t ntr = GetCTSTracks()->GetEntriesFast();
1057  //printf("N Tracks %d\n",ntr);
1058 
1059  for(Int_t i = 0 ; i < ntr ; i++)
1060  {
1061  AliVTrack * track = (AliVTrack*) (GetCTSTracks()->At(i));
1062 
1063  //Check if has TOF info, if not skip
1064  ULong_t status = track->GetStatus();
1065  Bool_t okTOF = (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ;
1066  vertexBC = track->GetTOFBunchCrossing(bz);
1067  Float_t pt = track->Pt();
1068 
1069  if(!okTOF) continue;
1070 
1071  // Get DCA x, y
1072  Double_t dca[2] = {1e6,1e6};
1073  Double_t covar[3] = {1e6,1e6,1e6};
1074  track->PropagateToDCA(vtx,bz,100.,dca,covar);
1075 
1076  if(AcceptDCA(pt,dca[0]))
1077  {
1078  if (vertexBC !=0 && fVertexBC != AliVTrack::kTOFBCNA) return vertexBC;
1079  else if(vertexBC == 0) bc0 = kTRUE;
1080  }
1081  }
1082 
1083  if( bc0 ) vertexBC = 0 ;
1084  else vertexBC = AliVTrack::kTOFBCNA ;
1085 
1086  return vertexBC;
1087 }
1088 
1089 //_____________________________
1096 //_____________________________
1098 {
1099  return track->GetID();
1100 }
1101 
1102 //_____________________________
1105 //_____________________________
1107 {
1108  // Activate debug level in reader
1109  if( fDebug >= 0 )
1110  (AliAnalysisManager::GetAnalysisManager())->AddClassDebug(this->ClassName(),fDebug);
1111 
1112  // Activate debug level in AliAnaWeights
1113  if( fWeightUtils->GetDebug() >= 0 )
1114  (AliAnalysisManager::GetAnalysisManager())->AddClassDebug(fWeightUtils->ClassName(), fWeightUtils->GetDebug());
1115 
1116  // Activate debug level in AliMCAnalysisUtils
1117  if( GetMCAnalysisUtils()->GetDebug() >= 0 )
1118  (AliAnalysisManager::GetAnalysisManager())->AddClassDebug(GetMCAnalysisUtils()->ClassName(),GetMCAnalysisUtils()->GetDebug());
1119 
1120  //printf("Debug levels: Reader %d, Neutral Sel %d, Iso %d\n",fDebug,GetMCAnalysisUtils()->GetDebug(),fWeightUtils->GetDebug());
1121 }
1122 
1123 //_______________________________________
1125 //_______________________________________
1127 {
1128  fDataType = kESD ;
1129  fCTSPtMin = 0.1 ;
1130  fEMCALPtMin = 0.1 ;
1131  fPHOSPtMin = 0.1 ;
1132  fCTSPtMax = 1000. ;
1133  fEMCALPtMax = 1000. ;
1134  fPHOSPtMax = 1000. ;
1135 
1136  fEMCALBadChMinDist = 0; // open, 2; // standard
1137  fPHOSBadChMinDist = 0; // open, 2; // standard
1138 
1139  fEMCALNCellsCut = 0; // open, 1; // standard
1140  fPHOSNCellsCut = 0; // open, 2; // standard
1141 
1142  // For clusters with energy above fEMCALHighEnergyNdiffCut, count cells
1143  // with E > fEMCALMinCellEnNdiffCut in different T-Card than highest energy cell.
1144  // Reject if 0.
1147 
1148  //Track DCA cuts
1149  // dca_xy cut = 0.0105+0.0350/TMath::Power(pt,1.1);
1150  fTrackDCACut[0] = 0.0105;
1151  fTrackDCACut[1] = 0.0350;
1152  fTrackDCACut[2] = 1.1;
1153 
1154  //Do not filter the detectors input by default.
1155  fFillEMCAL = kFALSE;
1156  fFillDCAL = kFALSE;
1157  fFillPHOS = kFALSE;
1158  fFillCTS = kFALSE;
1159  fFillEMCALCells = kFALSE;
1160  fFillPHOSCells = kFALSE;
1161 
1162  fDeltaAODFileName = "deltaAODPartCorr.root";
1164  fEventTriggerMask = AliVEvent::kAny;
1165  fMixEventTriggerMask = AliVEvent::kAnyINT;
1166  fEventTriggerAtSE = kTRUE; // Use only events that pass event selection at SE base class
1167 
1168  fAcceptFastCluster = kTRUE;
1169  fEventType = -1;
1170 
1171  fRemoveLEDEvents = 0;
1174  fLEDMinCellEnergy = 0.5;
1175  fLEDMaxCellEnergy = 15.;
1176 
1177  fRemoveLEDStripEvents = 0 ;
1179  fLEDHighEnergyCutStrip[0] = 80 ; fLEDHighEnergyCutStrip[1] = 55 ;
1181  fLEDLowEnergyCutSM3Strip = 100; // open
1182  fLEDLowNCellsCutSM3Strip = 100; // open
1183 
1184  //We want tracks fitted in the detectors:
1185  //fTrackStatus=AliESDtrack::kTPCrefit;
1186  //fTrackStatus|=AliESDtrack::kITSrefit;
1187  fTrackStatus = 0;
1188 
1189  fV0ADC[0] = 0; fV0ADC[1] = 0;
1190  fV0Mul[0] = 0; fV0Mul[1] = 0;
1191 
1192  fZvtxCut = 10.;
1193 
1194  fNMixedEvent = 1;
1195 
1196  fPtHardAndJetPtFactor = 4. ;
1198 
1199  //Centrality
1200  fUseAliCentrality = kFALSE;
1201  fMultWithEventSel = kTRUE;
1202  fCentralityClass = "V0M";
1203  fCentralityOpt = 100;
1204  fCentralityBin[0] = fCentralityBin[1]=-1;
1205 
1206  fEventPlaneMethod = "V0";
1207 
1208  // Allocate memory (not sure this is the right place)
1209  fCTSTracks = new TObjArray();
1210  fEMCALClusters = new TObjArray();
1211  fDCALClusters = new TObjArray();
1212  fPHOSClusters = new TObjArray();
1213  //fTriggerAnalysis = new AliTriggerAnalysis;
1214  fAODBranchList = new TList ;
1215  fOutputContainer = new TList ;
1216 
1217  fEnergyHistogramNbins = 500 ;
1218  fEnergyHistogramLimit[0] = 0. ;
1219  fEnergyHistogramLimit[1] = 250.;
1220 
1221  fPileUpParamSPD[0] = 3 ; fPileUpParamSPD[1] = 0.8 ;
1222  fPileUpParamSPD[2] = 3.0 ; fPileUpParamSPD[3] = 2.0 ; fPileUpParamSPD[4] = 5.0;
1223 
1224  // Parametrized time cut (LHC11d)
1227 
1228  // Parametrized time cut (LHC11c)
1229  //fEMCALParamTimeCutMin[0] =-5; fEMCALParamTimeCutMin[1] =-1 ; fEMCALParamTimeCutMin[2] = 1.87; fEMCALParamTimeCutMin[3] = 0.4;
1230  //fEMCALParamTimeCutMax[0] = 3.5; fEMCALParamTimeCutMax[1] = 50; fEMCALParamTimeCutMax[2] = 0.15; fEMCALParamTimeCutMax[3] = 1.6;
1231 
1232  fTimeStampRunMin = -1;
1233  fTimeStampRunMax = 1e12;
1236 
1239 
1240  for(Int_t i = 0; i < 19; i++)
1241  {
1242  fEMCalBCEvent [i] = 0;
1243  fEMCalBCEventCut[i] = 0;
1244  fTrackBCEvent [i] = 0;
1245  fTrackBCEventCut[i] = 0;
1246  }
1247 
1248  // Trigger match-rejection
1249  fTriggerPatchTimeWindow[0] = 8;
1250  fTriggerPatchTimeWindow[1] = 9;
1251 
1252  fTriggerClusterBC = -10000 ;
1255  fTriggerClusterIndex = -1;
1256  fTriggerClusterId = -1;
1257 
1258  //Jets
1261  //if(!fNonStandardJets) fNonStandardJets = new TClonesArray("AliAODJet",100);
1262  if(!fNonStandardJets) fNonStandardJets = new TClonesArray("AliEmcalJet",100);
1265  //if(!fBackgroundJets) fBackgroundJets = new AliAODJetEventBackground();
1266  if(!fBackgroundJets) fBackgroundJets = new TClonesArray("AliEmcalJet",100);
1267 
1268  fSmearShowerShapeWidth = 0.005;
1269  fSmearNLMMin = 1;
1270  fSmearNLMMax = 1;
1271 
1272  fMCUtils = new AliMCAnalysisUtils() ;
1273 
1274  fWeightUtils = new AliAnaWeights() ;
1275  fEventWeight = 1 ;
1276 
1277  fTrackMultNPtCut = 8;
1278  fTrackMultPtCut[0] = 0.15; fTrackMultPtCut[1] = 0.5; fTrackMultPtCut[2] = 1.0;
1279  fTrackMultPtCut[3] = 2.0 ; fTrackMultPtCut[4] = 4.0; fTrackMultPtCut[5] = 6.0;
1280  fTrackMultPtCut[6] = 8.0 ; fTrackMultPtCut[7] = 10.;
1281  fTrackMultPtCut[8] = 15.0; fTrackMultPtCut[9] = 20.;
1282 
1283  for(Int_t ism = 0; ism < 22; ism++) fScaleFactorPerSM[ism] = 1. ;
1284 }
1285 
1286 //__________________________________________________________________________
1289 //__________________________________________________________________________
1291 {
1292  // Parametrized cut depending on E
1293  if(fUseParamTimeCut)
1294  {
1297  //printf("tof %f, minCut %f, maxCut %f\n",tof,minCut,maxCut);
1298  if( tof < minCut || tof > maxCut ) return kFALSE ;
1299  }
1300 
1301  //In any case, the time should to be larger than the fixed window ...
1302  if( tof < fEMCALTimeCutMin || tof > fEMCALTimeCutMax ) return kFALSE ;
1303 
1304  return kTRUE ;
1305 }
1306 
1307 //________________________________________________
1310 //________________________________________________
1312 {
1313  return fInputEvent->IsPileupFromSPD((Int_t) fPileUpParamSPD[0] , fPileUpParamSPD[1] ,
1314  fPileUpParamSPD[2] , fPileUpParamSPD[3] , fPileUpParamSPD[4] );
1315  //printf("Param : %d, %2.2f, %2.2f, %2.2f, %2.2f\n",(Int_t) fPileUpParamSPD[0], fPileUpParamSPD[1], fPileUpParamSPD[2], fPileUpParamSPD[3], fPileUpParamSPD[4]);
1316 }
1317 
1318 //__________________________________________________
1320 //__________________________________________________
1322 {
1323  if(fNPileUpClusters > fNPileUpClustersCut) return kTRUE ;
1324  else return kFALSE;
1325 }
1326 
1327 //________________________________________________________
1329 //________________________________________________________
1331 {
1332  if( IsPileUpFromSPD() && IsPileUpFromEMCal()) return kTRUE ;
1333  else return kFALSE;
1334 }
1335 
1336 //_______________________________________________________
1338 //_______________________________________________________
1340 {
1341  if( IsPileUpFromSPD() || IsPileUpFromEMCal()) return kTRUE ;
1342  else return kFALSE;
1343 }
1344 
1345 //___________________________________________________________
1347 //___________________________________________________________
1349 {
1350  if( IsPileUpFromSPD() && !IsPileUpFromEMCal()) return kTRUE ;
1351  else return kFALSE;
1352 }
1353 
1354 //___________________________________________________________
1356 //___________________________________________________________
1358 {
1359  if( !IsPileUpFromSPD() && IsPileUpFromEMCal()) return kTRUE ;
1360  else return kFALSE;
1361 }
1362 
1363 //______________________________________________________________
1365 //______________________________________________________________
1367 {
1368  if( !IsPileUpFromSPD() && !IsPileUpFromEMCal()) return kTRUE ;
1369  else return kFALSE;
1370 }
1371 
1372 //___________________________________________________________________________________
1380 //___________________________________________________________________________________
1381 Bool_t AliCaloTrackReader::FillInputEvent(Int_t iEntry, const char * /*curFileName*/)
1382 {
1383  fEventNumber = iEntry;
1384  fTriggerClusterIndex = -1;
1385  fTriggerClusterId = -1;
1386  fIsTriggerMatch = kFALSE;
1387  fTriggerClusterBC = -10000;
1388  fIsExoticEvent = kFALSE;
1389  fIsBadCellEvent = kFALSE;
1390  fIsBadMaxCellEvent = kFALSE;
1391 
1392  fIsTriggerMatchOpenCut[0] = kFALSE ;
1393  fIsTriggerMatchOpenCut[1] = kFALSE ;
1394  fIsTriggerMatchOpenCut[2] = kFALSE ;
1395 
1396  fCurrentParIndex = 0;
1397  if(IsParRun())
1398  {
1399  ULong64_t globalEventID = (ULong64_t)fInputEvent->GetBunchCrossNumber() + (ULong64_t)fInputEvent->GetOrbitNumber() * (ULong64_t)3564 + (ULong64_t)fInputEvent->GetPeriodNumber() * (ULong64_t)59793994260;
1400  for(Short_t ipar=0;ipar<GetCaloUtils()->GetEMCALRecoUtils()->GetNPars();ipar++)
1401  {
1402  if(globalEventID >= GetCaloUtils()->GetEMCALRecoUtils()->GetGlobalIDPar(ipar))
1403  {
1404  fCurrentParIndex++;
1405  }
1406  }
1407  }
1409 
1410  //fCurrentFileName = TString(currentFileName);
1411  if(!fInputEvent)
1412  {
1413  AliInfo("Input event not available, skip event analysis");
1414  return kFALSE;
1415  }
1416 
1417  fhNEventsAfterCut->Fill(0.5);
1418 
1419  //-----------------------------------------------
1420  // Select the event depending on the trigger type
1421  // and other event characteristics
1422  // like the goodness of the EMCal trigger
1423  //-----------------------------------------------
1424 
1425  Bool_t accept = CheckEventTriggers();
1426  if(!accept) return kFALSE;
1427 
1428  AliDebug(1,"Pass Event trigger selection");
1429 
1430 
1431  //------------------------------------------------------
1432  // Event rejection depending on time stamp
1433  //------------------------------------------------------
1434 
1435  if ( fTimeStampEventSelect )
1436  {
1437  Int_t timeStamp = fInputEvent->GetTimeStamp();
1438  Float_t timeStampFrac = 1.*(timeStamp-fTimeStampRunMin) / (fTimeStampRunMax-fTimeStampRunMin);
1439 
1440  //printf("stamp0 %d, max0 %d, frac %f\n", timeStamp-fTimeStampRunMin,fTimeStampRunMax-fTimeStampRunMin, timeStampFrac);
1441 
1442  if(timeStampFrac < fTimeStampEventFracMin || timeStampFrac > fTimeStampEventFracMax) return kFALSE;
1443 
1444  AliDebug(1,"Pass Time Stamp rejection");
1445 
1446  fhNEventsAfterCut->Fill(8.5);
1447  }
1448 
1450  {
1451  AliESDEvent* esd = dynamic_cast<AliESDEvent*> (fInputEvent);
1452  if(esd)
1453  {
1454  Int_t timeStamp = esd->GetTimeStampCTPBCCorr();
1455 
1456  if(timeStamp > fTimeStampEventCTPBCCorrMin && timeStamp <= fTimeStampEventCTPBCCorrMax) return kFALSE;
1457  }
1458 
1459  AliDebug(1,"Pass Time Stamp CTPBCCorr rejection");
1460 
1461  fhNEventsAfterCut->Fill(8.5);
1462  }
1463 
1464 
1465  //------------------------------------------------------
1466  // Event rejection depending on vertex, pileup, v0and
1467  //------------------------------------------------------
1468 
1469  FillVertexArray();
1470 
1472  {
1473  if( !CheckForPrimaryVertex() ) return kFALSE; // algorithm in ESD/AOD Readers
1474 
1475  fhNEventsAfterCut->Fill(9.5);
1476 
1477  if( TMath::Abs(fVertex[0][0] ) < 1.e-6 &&
1478  TMath::Abs(fVertex[0][1] ) < 1.e-6 &&
1479  TMath::Abs(fVertex[0][2] ) < 1.e-6 ) return kFALSE;
1480 
1481  AliDebug(1,"Pass primary vertex/null rejection");
1482 
1483  fhNEventsAfterCut->Fill(10.5);
1484  }
1485 
1486  //Reject events with Z vertex too large, only for SE analysis, if not, cut on the analysis code
1487  if(!GetMixedEvent() && TMath::Abs(fVertex[0][2]) > fZvtxCut) return kFALSE;
1488 
1489  fhNEventsAfterCut->Fill(11.5);
1490 
1491  AliDebug(1,"Pass z vertex rejection");
1492 
1493 
1494  //printf("Reader : IsPileUp %d, Multi %d\n",IsPileUpFromSPD(),fInputEvent->IsPileupFromSPDInMultBins());
1495 
1497  {
1498  // Do not analyze events with pileup
1499  Bool_t bPileup = IsPileUpFromSPD();
1500  //IsPileupFromSPDInMultBins() // method to try
1501  //printf("pile-up %d, %d, %2.2f, %2.2f, %2.2f, %2.2f\n",bPileup, (Int_t) fPileUpParamSPD[0], fPileUpParamSPD[1], fPileUpParamSPD[2], fPileUpParamSPD[3], fPileUpParamSPD[4]);
1502  if(bPileup) return kFALSE;
1503 
1504  AliDebug(1,"Pass Pile-Up event rejection");
1505 
1506  fhNEventsAfterCut->Fill(12.5);
1507  }
1508 
1510  {
1511  AliVVZERO* v0 = fInputEvent->GetVZEROData();
1512 
1513  Bool_t bV0AND = ((v0->GetV0ADecision()==1) && (v0->GetV0CDecision()==1));
1514  //bV0AND = fTriggerAnalysis->IsOfflineTriggerFired((AliESDEvent*)fInputEvent, AliTriggerAnalysis::kV0AND);
1515  //printf("V0AND event? %d\n",bV0AND);
1516 
1517  if(!bV0AND)
1518  {
1519  AliDebug(1,"Reject event by V0AND");
1520  return kFALSE;
1521  }
1522 
1523  AliDebug(1,"Pass V0AND event rejection");
1524 
1525  fhNEventsAfterCut->Fill(13.5);
1526  }
1527 
1528  //------------------------------------------------------
1529  // Check if there is a centrality value, PbPb analysis,
1530  // and if a centrality bin selection is requested
1531  //------------------------------------------------------
1532 
1533  //If we need a centrality bin, we select only those events in the corresponding bin.
1534  Int_t cen = -1;
1535  if ( fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
1536  {
1537  cen = GetEventCentrality();
1538 
1539  AliDebug(1,Form("Centrality %d in [%d,%d]?", cen, fCentralityBin[0], fCentralityBin[1]));
1540 
1541  if ( cen >= fCentralityBin[1] || cen < fCentralityBin[0] ) return kFALSE; //reject events out of bin.
1542 
1543  AliDebug(1,"Pass centrality rejection");
1544 
1545  fhNEventsAfterCut->Fill(14.5);
1546  }
1547 
1548  //----------------------------------------------------------------
1549  // MC events selections
1550  //----------------------------------------------------------------
1551  if ( GetMC() )
1552  {
1553  //----------------------------------------------------------------
1554  // Get the event headers
1555  //----------------------------------------------------------------
1556 
1557  // Main header
1558  // Init it first to 0 to tell the method to recover it.
1559  fGenEventHeader = 0;
1561 
1562  if ( fGenEventHeader )
1563  {
1564  AliDebug(1,Form("Selected event header class <%s>, name <%s>; cocktail %p",
1565  fGenEventHeader->ClassName(),
1566  fGenEventHeader->GetName(),
1567  GetMC()->GetCocktailList()));
1568  }
1569 
1570  //----------------------------------------------------------------
1571  // Reject the event if the event header name is not
1572  // the one requested among the possible generators.
1573  // Needed in case of cocktail MC generation with multiple options.
1574  //----------------------------------------------------------------
1576  {
1577  if ( !fGenEventHeader ) return kFALSE;
1578 
1579  AliDebug(1,"Pass Event header selection");
1580 
1581  fhNEventsAfterCut->Fill(15.5);
1582  }
1583 
1584  // Pythia header
1585  TString pyGenName = "";
1586  TString pyProcessName = "";
1587  Int_t pyProcess = 0;
1588  Int_t pyFirstGenPart = 0;
1589  Int_t pythiaVersion = 0;
1590 
1591  // Init it first to 0 to tell the method to recover it.
1594  pyGenName,pyProcessName,pyProcess,pyFirstGenPart,pythiaVersion);
1595 
1597  {
1598  AliDebug(2,Form("Pythia v%d name <%s>, process %d <%s>, first generated particle %d",
1599  pythiaVersion, pyGenName.Data(), pyProcess, pyProcessName.Data(), pyFirstGenPart));
1600 
1601  //---------------------------------------------------------------------------
1602  // In case of analysis of events with jets, skip those with jet pt > 5 pt hard
1603  // To be used on for MC data in pT hard bins
1604  //---------------------------------------------------------------------------
1605 
1607  {
1608  if ( !ComparePtHardAndJetPt(pyProcess, pyProcessName) ) return kFALSE ;
1609 
1610  AliDebug(1,"Pass Pt Hard - Jet rejection");
1611 
1612  fhNEventsAfterCut->Fill(16.5);
1613  }
1614 
1616  {
1617  if ( !ComparePtHardAndClusterPt(pyProcess, pyProcessName) ) return kFALSE ;
1618 
1619  AliDebug(1,"Pass Pt Hard - Cluster rejection");
1620 
1621  fhNEventsAfterCut->Fill(17.5);
1622  }
1623  } // pythia header
1624  } // MC
1625 
1626  //------------------------------------------------------------------
1627  // Recover the weight assigned to the event, if provided
1628  // right now only for pT-hard bins and centrality dependent weights
1629  //------------------------------------------------------------------
1631  {
1633 
1635 
1637  }
1638 
1639  //-------------------------------------------------------
1640  // Get the main vertex BC, in case not available
1641  // it is calculated in FillCTS checking the BC of tracks
1642  //------------------------------------------------------
1643  fVertexBC = fInputEvent->GetPrimaryVertex()->GetBC();
1644 
1645  //-----------------------------------------------
1646  // Fill the arrays with cluster/tracks/cells data
1647  //-----------------------------------------------
1648 
1649  if(fFillCTS)
1650  {
1651  FillInputCTS();
1652  //Accept events with at least one track
1653  if(fTrackMult[0] == 0 && fDoRejectNoTrackEvents) return kFALSE ;
1654 
1655  AliDebug(1,"Pass rejection of null track events");
1656 
1657  fhNEventsAfterCut->Fill(18.5);
1658  }
1659 
1661  {
1662  if(fVertexBC != 0 && fVertexBC != AliVTrack::kTOFBCNA) return kFALSE ;
1663 
1664  AliDebug(1,"Pass rejection of events with vertex at BC!=0");
1665 
1666  fhNEventsAfterCut->Fill(19.5);
1667  }
1668 
1669  if(fFillEMCALCells)
1671 
1672  if(fFillPHOSCells)
1674 
1675  if(fFillEMCAL || fFillDCAL)
1676  FillInputEMCAL();
1677 
1678  if(fFillPHOS)
1679  FillInputPHOS();
1680 
1681  FillInputVZERO();
1682 
1683  //one specified jet branch
1688 
1689  AliDebug(1,"Event accepted for analysis");
1690 
1691  return kTRUE ;
1692 }
1693 
1694 //__________________________________________________
1697 //__________________________________________________
1699 {
1700  if(fUseAliCentrality)
1701  {
1702  if ( !GetCentrality() ) return -1;
1703 
1704  AliDebug(1,Form("Cent. Percentile: V0M %2.2f, CL0 %2.2f, CL1 %2.2f; selected class %s",
1705  GetCentrality()->GetCentralityPercentile("V0M"),
1706  GetCentrality()->GetCentralityPercentile("CL0"),
1707  GetCentrality()->GetCentralityPercentile("CL1"),
1708  fCentralityClass.Data()));
1709 
1710  if (fCentralityOpt == 100) return (Int_t) GetCentrality()->GetCentralityPercentile(fCentralityClass); // 100 bins max
1711  else if(fCentralityOpt == 10) return GetCentrality()->GetCentralityClass10(fCentralityClass);// 10 bins max
1712  else if(fCentralityOpt == 20) return GetCentrality()->GetCentralityClass5(fCentralityClass); // 20 bins max
1713  else
1714  {
1715  AliInfo(Form("Unknown centrality option %d, use 10, 20 or 100",fCentralityOpt));
1716  return -1;
1717  }
1718  }
1719  else
1720  {
1721  if ( !GetMultSelCen() ) return -1;
1722 
1723  AliDebug(1,Form("Mult. Percentile: V0M %2.2f, CL0 %2.2f, CL1 %2.2f; selected class %s",
1724  GetMultSelCen()->GetMultiplicityPercentile("V0M",1),
1725  GetMultSelCen()->GetMultiplicityPercentile("CL0",1),
1726  GetMultSelCen()->GetMultiplicityPercentile("CL1",1),
1727  fCentralityClass.Data()));
1728 
1729  return (Int_t) GetMultSelCen()->GetMultiplicityPercentile(fCentralityClass, fMultWithEventSel); // returns centrality only for events used in calibration
1730 
1731  // equivalent to
1732  //GetMultSelCen()->GetMultiplicityPercentile("V0M", kFALSE); // returns centrality for any event
1733  //Int_t qual = GetMultSelCen()->GetEvSelCode(); if (qual ! = 0) cent = qual;
1734  }
1735 }
1736 
1737 //_____________________________________________________
1740 //_____________________________________________________
1742 {
1743  if( !GetEventPlane() ) return -1000;
1744 
1745  Float_t ep = GetEventPlane()->GetEventplane(GetEventPlaneMethod(), GetInputEvent());
1746 
1747  if(GetEventPlaneMethod()=="Q" && (ep < 0 || ep > TMath::Pi()))
1748  {
1749  AliDebug(1,Form("Bad EP for <Q> method : %f",ep));
1750  return -1000;
1751  }
1752  else if(GetEventPlaneMethod().Contains("V0") )
1753  {
1754  if((ep > TMath::Pi()/2 || ep < -TMath::Pi()/2))
1755  {
1756  AliDebug(1,Form("Bad EP for <%s> method : %f",GetEventPlaneMethod().Data(), ep));
1757  return -1000;
1758  }
1759 
1760  ep+=TMath::Pi()/2; // put same range as for <Q> method
1761  }
1762 
1763  AliDebug(3,Form("Event plane angle %f",ep));
1764 
1765 // if(fDebug > 0 )
1766 // {
1767 // if (ep > TMath::Pi()) printf("AliCaloTrackReader::GetEventPlaneAngle() - Too large angle = %f\n",ep);
1768 // else if(ep < 0 ) printf("AliCaloTrackReader::GetEventPlaneAngle() - Negative angle = %f\n" ,ep);
1769 // }
1770 
1771  return ep;
1772 }
1773 
1774 //__________________________________________________________
1776 //__________________________________________________________
1778 {
1779  vertex[0] = fVertex[0][0];
1780  vertex[1] = fVertex[0][1];
1781  vertex[2] = fVertex[0][2];
1782 }
1783 
1784 //__________________________________________________________________________
1786 //__________________________________________________________________________
1787 void AliCaloTrackReader::GetVertex(Double_t vertex[3], Int_t evtIndex) const
1788 {
1789  vertex[0] = fVertex[evtIndex][0];
1790  vertex[1] = fVertex[evtIndex][1];
1791  vertex[2] = fVertex[evtIndex][2];
1792 }
1793 
1794 //________________________________________
1797 //________________________________________
1799 {
1800  // Delete previous vertex
1801  if(fVertex)
1802  {
1803  for (Int_t i = 0; i < fNMixedEvent; i++)
1804  {
1805  delete [] fVertex[i] ;
1806  }
1807  delete [] fVertex ;
1808  }
1809 
1810  fVertex = new Double_t*[fNMixedEvent] ;
1811  for (Int_t i = 0; i < fNMixedEvent; i++)
1812  {
1813  fVertex[i] = new Double_t[3] ;
1814  fVertex[i][0] = 0.0 ;
1815  fVertex[i][1] = 0.0 ;
1816  fVertex[i][2] = 0.0 ;
1817  }
1818 
1819  if (!fMixedEvent)
1820  { // Single event analysis
1821  if(fDataType!=kMC)
1822  {
1823 
1824  if(fInputEvent->GetPrimaryVertex())
1825  {
1826  fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]);
1827  }
1828  else
1829  {
1830  AliWarning("NULL primary vertex");
1831  fVertex[0][0]=0.; fVertex[0][1]=0.; fVertex[0][2]=0.;
1832  }//Primary vertex pointer do not exist
1833 
1834  } else
1835  {// MC read event
1836  fVertex[0][0]=0.; fVertex[0][1]=0.; fVertex[0][2]=0.;
1837  }
1838 
1839  AliDebug(1,Form("Single Event Vertex : %f,%f,%f",fVertex[0][0],fVertex[0][1],fVertex[0][2]));
1840 
1841  } else
1842  { // MultiEvent analysis
1843  for (Int_t iev = 0; iev < fNMixedEvent; iev++)
1844  {
1845  if (fMixedEvent->GetVertexOfEvent(iev))
1846  fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
1847  else
1848  AliWarning("No vertex found");
1849 
1850  AliDebug(1,Form("Multi Event %d Vertex : %f,%f,%f",iev,fVertex[iev][0],fVertex[iev][1],fVertex[iev][2]));
1851  }
1852  }
1853 }
1854 
1855 //_____________________________________
1861 //_____________________________________
1863 {
1864  AliDebug(1,"Begin");
1865 
1866  Double_t pTrack[3] = {0,0,0};
1867 
1868  Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
1869  Int_t nstatus = 0;
1870  Double_t bz = GetInputEvent()->GetMagneticField();
1871 
1872  for(Int_t i = 0; i < 19; i++)
1873  {
1874  fTrackBCEvent [i] = 0;
1875  fTrackBCEventCut[i] = 0;
1876  }
1877 
1878  for(Int_t iptCut = 0; iptCut < fTrackMultNPtCut; iptCut++ )
1879  {
1880  fTrackMult [iptCut] = 0;
1881  fTrackSumPt[iptCut] = 0;
1882  }
1883 
1884  Bool_t bc0 = kFALSE;
1885  if(fRecalculateVertexBC) fVertexBC = AliVTrack::kTOFBCNA;
1886 
1887  for (Int_t itrack = 0; itrack < nTracks; itrack++)
1888  {
1889  AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
1890 
1891  if ( !AcceptParticleMCLabel( TMath::Abs(track->GetLabel()) ) ) continue ;
1892 
1893  fhCTSTrackCutsPt[0]->Fill(track->Pt());
1894 
1895  //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
1896  ULong_t status = track->GetStatus();
1897 
1898  if (fTrackStatus && !((status & fTrackStatus) == fTrackStatus))
1899  continue ;
1900 
1901  fhCTSTrackCutsPt[1]->Fill(track->Pt());
1902 
1903  nstatus++;
1904 
1905  //-------------------------
1906  // Select the tracks depending on cuts of AOD or ESD
1907  if(!SelectTrack(track, pTrack)) continue ;
1908 
1909  fhCTSTrackCutsPt[2]->Fill(track->Pt());
1910 
1911  //-------------------------
1912  // TOF cuts
1913  Bool_t okTOF = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
1914  Double_t tof = -1000;
1915  Int_t trackBC = -1000 ;
1916 
1917  if(fAccessTrackTOF)
1918  {
1919  if(okTOF)
1920  {
1921  trackBC = track->GetTOFBunchCrossing(bz);
1922  SetTrackEventBC(trackBC+9);
1923 
1924  tof = track->GetTOFsignal()*1e-3;
1925 
1926  //After selecting tracks with small DCA, pointing to vertex, set vertex BC depeding on tracks BC
1928  {
1929  if (trackBC != 0 && trackBC != AliVTrack::kTOFBCNA) fVertexBC = trackBC;
1930  else if(trackBC == 0) bc0 = kTRUE;
1931  }
1932 
1933  //In any case, the time should to be larger than the fixed window ...
1934  if( fUseTrackTimeCut && (trackBC !=0 || tof < fTrackTimeCutMin || tof > fTrackTimeCutMax) )
1935  {
1936  //printf("Remove track time %f and bc = %d\n",tof,trackBC);
1937  continue ;
1938  }
1939  //else printf("Accept track time %f and bc = %d\n",tof,trackBC);
1940  }
1941  }
1942 
1943  fhCTSTrackCutsPt[3]->Fill(track->Pt());
1944 
1945  //---------------------
1946  // DCA cuts
1947  //
1948  fMomentum.SetPxPyPzE(pTrack[0],pTrack[1],pTrack[2],0);
1949 
1950  if(fUseTrackDCACut)
1951  {
1952  Float_t dcaTPC =-999;
1953  //In case of AODs, TPC tracks cannot be propagated back to primary vertex,
1954  if( fDataType == kAOD ) dcaTPC = ((AliAODTrack*) track)->DCA();
1955 
1956  //normal way to get the dca, cut on dca_xy
1957  if(dcaTPC==-999)
1958  {
1959  Double_t dca[2] = {1e6,1e6};
1960  Double_t covar[3] = {1e6,1e6,1e6};
1961  Bool_t okDCA = track->PropagateToDCA(fInputEvent->GetPrimaryVertex(),bz,100.,dca,covar);
1962  if( okDCA) okDCA = AcceptDCA(fMomentum.Pt(),dca[0]);
1963  if(!okDCA)
1964  {
1965  //printf("AliCaloTrackReader::FillInputCTS() - Reject track pt %2.2f, dca_xy %2.4f\n",fMomentum.Pt(),dca[0]);
1966  continue ;
1967  }
1968  }
1969  }// DCA cuts
1970 
1971  fhCTSTrackCutsPt[4]->Fill(track->Pt());
1972 
1973  //-------------------------
1974  // Kinematic/acceptance cuts
1975  //
1976  // Count the tracks in eta < 0.9 and different pT cuts
1977  Float_t ptTrack = fMomentum.Pt();
1978  if(TMath::Abs(track->Eta())< fTrackMultEtaCut)
1979  {
1980  for(Int_t iptCut = 0; iptCut < fTrackMultNPtCut; iptCut++ )
1981  {
1982  if(ptTrack > fTrackMultPtCut[iptCut])
1983  {
1984  fTrackMult [iptCut]++;
1985  fTrackSumPt[iptCut]+=ptTrack;
1986  }
1987  }
1988  }
1989 
1990  if(fCTSPtMin > ptTrack || fCTSPtMax < ptTrack) continue ;
1991 
1992  // Check effect of cuts on track BC
1993  if(fAccessTrackTOF && okTOF) SetTrackEventBCcut(trackBC+9);
1994 
1995  if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS)) continue;
1996 
1997  fhCTSTrackCutsPt[5]->Fill(track->Pt());
1998 
1999  // ------------------------------
2000  // Add selected tracks to array
2001  AliDebug(2,Form("Selected tracks pt %3.2f, phi %3.2f deg, eta %3.2f",
2002  fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2003 
2004  fCTSTracks->Add(track);
2005 
2006  // TODO, check if remove
2007  if (fMixedEvent) track->SetID(itrack);
2008 
2009  }// track loop
2010 
2011  if( fRecalculateVertexBC && (fVertexBC == 0 || fVertexBC == AliVTrack::kTOFBCNA))
2012  {
2013  if( bc0 ) fVertexBC = 0 ;
2014  else fVertexBC = AliVTrack::kTOFBCNA ;
2015  }
2016 
2017  AliDebug(1,Form("CTS entries %d, input tracks %d, pass status %d, multipliticy %d", fCTSTracks->GetEntriesFast(), nTracks, nstatus, fTrackMult[0]));//fCTSTracksNormalInputEntries);
2018 }
2019 
2020 //_______________________________________________________________________________
2038 //_______________________________________________________________________________
2040 {
2041  // Accept clusters with the proper label, only applicable for MC
2042  if ( clus->GetLabel() >= 0 ) // -1 corresponds to noisy MC
2043  {
2044  if ( !AcceptParticleMCLabel(clus->GetLabel()) ) return ;
2045  }
2046 
2047  // TODO, not sure if needed anymore
2048  Int_t vindex = 0 ;
2049  if (fMixedEvent)
2050  vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
2051 
2052  clus->GetMomentum(fMomentum, fVertex[vindex]);
2053 
2054  // No correction/cut applied yet
2055  fhEMCALClusterCutsE[0]->Fill(clus->E());
2056 
2057  // Get the maximum cell energy, its SM number and its col, row location, needed in
2058  // different places of this method, although not active by default, one can consider
2059  // deactivate this and only activate it when requiered.
2060  Int_t absIdMax= -1;
2061  Int_t iSupMod = -1;
2062  Int_t iphiMax = -1;
2063  Int_t ietaMax = -1;
2064  Bool_t shared = kFALSE;
2065  GetCaloUtils()->GetEMCALRecoUtils()->GetMaxEnergyCell(GetCaloUtils()->GetEMCALGeometry(),
2066  GetEMCALCells(),clus,absIdMax,iSupMod,
2067  ietaMax,iphiMax,shared);
2068 
2069  //if( (fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10 )
2070  AliDebug(2,Form("Input cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f, nCells %d, SM %d",
2071  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta(), clus->GetNCells(), iSupMod));
2072 
2073  //---------------------------
2074  // Embedding case
2075  // TO BE REVISED
2077  {
2078  if(clus->GetNLabels()==0 || clus->GetLabel() < 0) return;
2079  //else printf("Embedded cluster, %d, n label %d label %d \n",iclus,clus->GetNLabels(),clus->GetLabel());
2080  }
2081 
2082  //--------------------------------------
2083  // Apply some corrections in the cluster
2084  //
2086  {
2087  //Recalibrate the cluster energy
2089  {
2091 
2092  clus->SetE(energy);
2093  //printf("Recalibrated Energy %f\n",clus->E());
2094 
2097 
2098  } // recalculate E
2099 
2100  //Recalculate distance to bad channels, if new list of bad channels provided
2102 
2103  //Recalculate cluster position
2104  if(GetCaloUtils()->IsRecalculationOfClusterPositionOn())
2105  {
2107  //clus->GetPosition(pos);
2108  //printf("After Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
2109  }
2110 
2111  // Recalculate TOF
2112  if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn())
2113  {
2114  Double_t tof = clus->GetTOF();
2115  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2116 
2117  //additional L1 phase shift
2119  {
2120  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(iSupMod,fInputEvent->GetBunchCrossNumber(), tof, fCurrentParIndex);
2121  }
2122 
2123  clus->SetTOF(tof);
2124 
2125  }// Time recalibration
2126  }
2127 
2128  // Check effect of corrections
2129  fhEMCALClusterCutsE[1]->Fill(clus->E());
2130 
2131  //-----------------------------------------------------------------
2132  // Reject clusters with bad channels, close to borders and exotic
2133  //
2134  Bool_t goodCluster = GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,
2135  GetCaloUtils()->GetEMCALGeometry(),
2136  GetEMCALCells(),fInputEvent->GetBunchCrossNumber());
2137 
2138  if(!goodCluster)
2139  {
2140  //if( (fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10 )
2141  AliDebug(1,Form("Bad cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2142  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2143 
2144  return;
2145  }
2146 
2147  // Check effect of bad cluster removal
2148  fhEMCALClusterCutsE[2]->Fill(clus->E());
2149 
2150  //Float_t pos[3];
2151  //clus->GetPosition(pos);
2152  //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
2153 
2154  //-----------------------------------------------------
2155  // Correct non linearity, smear energy or scale per SM
2156  //
2158  {
2160 
2161  AliDebug(5,Form("Correct Non Lin: Old E %3.2f, New E %3.2f",
2162  fMomentum.E(),clus->E()));
2163  }
2164 
2165  // In case of MC analysis, to match resolution/calibration in real data
2166  // Not needed anymore, just leave for MC studies on systematics
2167  if( GetCaloUtils()->GetEMCALRecoUtils()->IsClusterEnergySmeared() )
2168  {
2169  Float_t rdmEnergy = GetCaloUtils()->GetEMCALRecoUtils()->SmearClusterEnergy(clus);
2170 
2171  AliDebug(5,Form("Smear energy: Old E %3.2f, New E %3.2f",clus->E(),rdmEnergy));
2172 
2173  clus->SetE(rdmEnergy);
2174  }
2175 
2176  // In case of uncalibrated data, or non final non linearity for MC
2177  // or calibration of MC for SMs behind TRD, apply a global scale factor per SM
2178  if ( fScaleEPerSM && iSupMod < 22 && iSupMod >=0)
2179  {
2180  Float_t scale = fScaleFactorPerSM[iSupMod];
2181 
2182  AliDebug(5,Form("Scale energy for SM %d: Old E %3.2f, scale factor %1.5f",iSupMod,clus->E(),scale));
2183 
2184  clus->SetE(clus->E()*scale);
2185  }
2186 
2187  clus->GetMomentum(fMomentum, fVertex[vindex]);
2188  fhEMCALClusterEtaPhi->Fill(fMomentum.Eta(),GetPhi(fMomentum.Phi()));
2189 
2190  // Check effect linearity correction, energy smearing
2191  fhEMCALClusterCutsE[3]->Fill(clus->E());
2192 
2193  // Check the event BC depending on EMCal clustr before final cuts
2194  Double_t tof = clus->GetTOF()*1e9;
2195 
2196  Int_t bc = TMath::Nint(tof/50) + 9;
2197  //printf("tof %2.2f, bc+5=%d\n",tof,bc);
2198 
2199  SetEMCalEventBC(bc);
2200 
2201  //--------------------------------------
2202  // Apply some kinematical/acceptance cuts
2203  //
2204  if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E())
2205  {
2206  AliDebug(2,Form("Cluster E out of range, %2.2f < %2.2f < %2.2f",fEMCALPtMin,clus->E(),fEMCALPtMax));
2207  return ;
2208  }
2209 
2210  // Select cluster fiducial region
2211  //
2212  Bool_t bEMCAL = kFALSE;
2213  Bool_t bDCAL = kFALSE;
2214  if(fCheckFidCut)
2215  {
2216  if(fFillEMCAL && fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) bEMCAL = kTRUE ;
2217  if(fFillDCAL && fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kDCAL )) bDCAL = kTRUE ;
2218  }
2219  else
2220  {
2221  bEMCAL = kTRUE;
2222  }
2223 
2224  //---------------------------------------------------------------------
2225  // Mask all cells in collumns facing ALICE thick material if requested
2226  //
2227  if ( GetCaloUtils()->GetNMaskCellColumns() )
2228  {
2229  AliDebug(2,Form("Masked collumn: cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2230  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2231 
2232  if ( GetCaloUtils()->MaskFrameCluster(iSupMod, ietaMax) )
2233  {
2234  AliDebug(2,"Mask cluster");
2235  return;
2236  }
2237  }
2238 
2239  // Check effect of energy and fiducial cuts
2240  if ( bEMCAL || bDCAL )
2241  {
2242  fhEMCALClusterCutsE[4]->Fill(clus->E());
2244  }
2245  else
2246  {
2247  AliDebug(2,"Cluster not on EMCal or DCal selected region");
2248  return ;
2249  }
2250 
2251  //----------------------------------------------------
2252  // Apply N cells cut
2253  //
2254  if(clus->GetNCells() <= fEMCALNCellsCut && fDataType != AliCaloTrackReader::kMC)
2255  {
2256  AliDebug(2,Form("Cluster with n cells %d < %d",clus->GetNCells(), fEMCALNCellsCut));
2257  return ;
2258  }
2259 
2260  // Check effect of n cells cut
2261  fhEMCALClusterCutsE[5]->Fill(clus->E());
2262 
2263  //----------------------------------------------------
2264  // Apply distance to bad channel cut
2265  //
2266  Double_t distBad = clus->GetDistanceToBadChannel() ; //Distance to bad channel
2267 
2268  if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2269 
2270  if(distBad < fEMCALBadChMinDist)
2271  {
2272  AliDebug(2, Form("Cluster close to bad, dist %2.2f < %2.2f",distBad,fEMCALBadChMinDist));
2273  return ;
2274  }
2275 
2276  // Check effect distance to bad channel cut
2277  fhEMCALClusterCutsE[6]->Fill(clus->E());
2278 
2279  //------------------------------------------
2280  // Apply time cut, count EMCal BC before cut
2281  //
2282  SetEMCalEventBCcut(bc);
2283 
2284  // Shift time in case of no calibration with rough factor
2285  Double_t tofShift = tof;
2286  //if(tof > 400) tofShift-=615;
2287  fhEMCALClusterTimeE->Fill(clus->E(),tofShift);
2288 
2289  if(!IsInTimeWindow(tof,clus->E()))
2290  {
2291  fNPileUpClusters++ ;
2292  if(fUseEMCALTimeCut)
2293  {
2294  AliDebug(2,Form("Out of time window E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f, time %e",
2295  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta(),tof));
2296 
2297  return ;
2298  }
2299  }
2300  else
2302 
2303  // Check effect of time cut
2304  fhEMCALClusterCutsE[7]->Fill(clus->E());
2305 
2306  //----------------------------------------
2307  // Apply cut on number of cells in different T-Card
2308  // than highest energy cell. At high E it should be more than 0
2309  // if not, sign of exotic cluster
2310 
2311  Int_t nDiff = 0, nSame = 0;
2312  Float_t eDiff = 0, eSame = 0;
2314  nDiff, nSame, eDiff, eSame,
2316  if ( nDiff == 0 && clus->E() > fEMCALHighEnergyNdiffCut )
2317  {
2318  printf("** Reader: Reject cluster with E = %2.1f (min %2.1f) and n cells in diff TCard = %d, for Ecell min = %1.2f; m02 %2.2f, ncells %d\n",
2319  clus->E(),fEMCALHighEnergyNdiffCut,nDiff,fEMCALMinCellEnNdiffCut,clus->GetM02(),clus->GetNCells());
2320  return;
2321  }
2322 
2323  fhEMCALClusterCutsE[8]->Fill(clus->E());
2324 
2325  //----------------------------------------------------
2326  // Smear the SS to try to match data and simulations,
2327  // do it only for simulations.
2328  //
2329  if ( fSmearShowerShape && clus->GetNCells() > 2 )
2330  {
2331  Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(clus, GetEMCALCells());
2332 // Int_t nMaxima = clus->GetNExMax(); // For Run2
2333 
2334  if ( nMaxima >= fSmearNLMMin && nMaxima <= fSmearNLMMax )
2335  {
2336  AliDebug(2,Form("Smear shower shape - Original: %2.4f", clus->GetM02()));
2338  {
2339  clus->SetM02( clus->GetM02() + fRandom.Landau(0, fSmearShowerShapeWidth) );
2340  }
2342  {
2343  if(iclus%3 == 0 && clus->GetM02() > 0.1) clus->SetM02( clus->GetM02() + fRandom.Landau(0.05, fSmearShowerShapeWidth) ); //fSmearShowerShapeWidth = 0.035
2344  }
2345  else if (fSmearingFunction == kNoSmearing)
2346  {
2347  clus->SetM02( clus->GetM02() );
2348  }
2349  //clus->SetM02( fRandom.Landau(clus->GetM02(), fSmearShowerShapeWidth) );
2350  AliDebug(2,Form("Width %2.4f Smeared : %2.4f", fSmearShowerShapeWidth,clus->GetM02()));
2351  }
2352  }
2353 
2354  //--------------------------------------------------------
2355  // Fill the corresponding array with the selected clusters
2356  // Usually just filling EMCal array with upper or lower clusters is enough,
2357  // but maybe we want to do EMCal-DCal correlations.
2358 
2359  //if((fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10)
2360  AliDebug(2,Form("Selected clusters (EMCAL%d, DCAL%d), E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2361  bEMCAL,bDCAL,fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2362 
2363 
2364  if (bEMCAL) fEMCALClusters->Add(clus);
2365  else if(bDCAL ) fDCALClusters ->Add(clus);
2366 
2367  // TODO, not sure if needed anymore
2368  if (fMixedEvent)
2369  clus->SetID(iclus) ;
2370 }
2371 
2372 //_______________________________________
2379 //_______________________________________
2381 {
2382  AliDebug(1,"Begin");
2383 
2384  // First recalibrate cells, time or energy
2385  // if(GetCaloUtils()->IsRecalibrationOn())
2386  // GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCells(GetCaloUtils()->GetEMCALGeometry(),
2387  // GetEMCALCells(),
2388  // fInputEvent->GetBunchCrossNumber());
2389 
2390  fNPileUpClusters = 0; // Init counter
2391  fNNonPileUpClusters = 0; // Init counter
2392  for(Int_t i = 0; i < 19; i++)
2393  {
2394  fEMCalBCEvent [i] = 0;
2395  fEMCalBCEventCut[i] = 0;
2396  }
2397 
2398  //Loop to select clusters in fiducial cut and fill container with aodClusters
2399  if(fEMCALClustersListName=="")
2400  {
2401  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2402  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2403  {
2404  AliVCluster * clus = 0;
2405  if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
2406  {
2407  if (clus->IsEMCAL())
2408  {
2409  FillInputEMCALAlgorithm(clus, iclus);
2410  }//EMCAL cluster
2411  }// cluster exists
2412  }// cluster loop
2413 
2414  //Recalculate track matching
2416 
2417  }//Get the clusters from the input event
2418  else
2419  {
2420  TClonesArray * clusterList = 0x0;
2421 
2422  if (fInputEvent->FindListObject(fEMCALClustersListName))
2423  {
2424  clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
2425  }
2426  else if(fOutputEvent)
2427  {
2428  clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
2429  }
2430 
2431  if(!clusterList)
2432  {
2433  AliWarning(Form("Wrong name of list with clusters? <%s>",fEMCALClustersListName.Data()));
2434  return;
2435  }
2436 
2437  Int_t nclusters = clusterList->GetEntriesFast();
2438  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2439  {
2440  AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
2441  //printf("E %f\n",clus->E());
2442  if (clus) FillInputEMCALAlgorithm(clus, iclus);
2443  else AliWarning("Null cluster in list!");
2444  }// cluster loop
2445 
2446  // Recalculate the pile-up time, in case long time clusters removed during clusterization
2447  //printf("Input event INIT : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
2448 
2449  fNPileUpClusters = 0; // Init counter
2450  fNNonPileUpClusters = 0; // Init counter
2451  for(Int_t i = 0; i < 19; i++)
2452  {
2453  fEMCalBCEvent [i] = 0;
2454  fEMCalBCEventCut[i] = 0;
2455  }
2456 
2457  for (Int_t iclus = 0; iclus < fInputEvent->GetNumberOfCaloClusters(); iclus++)
2458  {
2459  AliVCluster * clus = 0;
2460 
2461  if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
2462  {
2463  if (clus->IsEMCAL())
2464  {
2465 
2466  Float_t frac =-1;
2467  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
2468  Double_t tof = clus->GetTOF();
2469  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2470  //additional L1 phase shift
2472  {
2473  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absIdMax), fInputEvent->GetBunchCrossNumber(), tof, fCurrentParIndex);
2474  }
2475 
2476  tof*=1e9;
2477 
2478  //printf("Input event cluster : AbsIdMax %d, E %2.2f, time %2.2f \n", absIdMax,clus->E(),tof);
2479 
2480  //Reject clusters with bad channels, close to borders and exotic;
2481  if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber())) continue;
2482 
2483  Int_t bc = TMath::Nint(tof/50) + 9;
2484  SetEMCalEventBC(bc);
2485 
2486  if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E()) continue ;
2487 
2488  clus->GetMomentum(fMomentum, fVertex[0]);
2489 
2490  if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) return ;
2491 
2492  SetEMCalEventBCcut(bc);
2493 
2494  if(!IsInTimeWindow(tof,clus->E()))
2495  fNPileUpClusters++ ;
2496  else
2498 
2499  }
2500  }
2501  }
2502 
2503  //printf("Input event : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
2504 
2505  // Recalculate track matching, not necessary if already done in the reclusterization task.
2506  // in case it was not done ...
2508 
2509  }
2510 
2511  AliDebug(1,Form("EMCal selected clusters %d",
2512  fEMCALClusters->GetEntriesFast()));
2513  AliDebug(2,Form("\t n pile-up clusters %d, n non pile-up %d",
2515 }
2516 
2517 //_______________________________________
2519 //_______________________________________
2521 {
2522  AliDebug(1,"Begin");
2523 
2524  // Loop to select clusters in fiducial cut and fill container with aodClusters
2525  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2526  TString genName;
2527  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2528  {
2529  AliVCluster * clus = fInputEvent->GetCaloCluster(iclus) ;
2530  if ( !clus ) continue ;
2531 
2532  if ( !clus->IsPHOS() ) continue ;
2533 
2534  if(clus->GetLabel() >=0 ) // -1 corresponds to noisy MC
2535  {
2536  if ( !AcceptParticleMCLabel(clus->GetLabel()) ) continue ;
2537  }
2538 
2539  fhPHOSClusterCutsE[0]->Fill(clus->E());
2540 
2541  // Skip CPV input
2542  if( clus->GetType() == AliVCluster::kPHOSCharged ) continue ;
2543 
2544  fhPHOSClusterCutsE[1]->Fill(clus->E());
2545 
2546  //---------------------------------------------
2547  // Old feature, try to rely on PHOS tender
2548  //
2550  {
2551  // Recalibrate the cluster energy
2553  {
2554  Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
2555  clus->SetE(energy);
2556  }
2557  }
2558 
2559  //----------------------------------------------------------------------------------
2560  // Check if the cluster contains any bad channel and if close to calorimeter borders
2561  //
2562  // Old feature, try to rely on PHOS tender
2563  if( GetCaloUtils()->ClusterContainsBadChannel(kPHOS,clus->GetCellsAbsId(), clus->GetNCells()))
2564  continue;
2565 
2566  if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells()))
2567  continue;
2568 
2569  // TODO, add exotic cut???
2570 
2571  fhPHOSClusterCutsE[2]->Fill(clus->E());
2572 
2573  // TODO Dead code? remove?
2574  Int_t vindex = 0 ;
2575  if (fMixedEvent)
2576  vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
2577 
2578  clus->GetMomentum(fMomentum, fVertex[vindex]);
2579 
2580  //----------------------------------------------------------------------------------
2581  // Remove clusters close to borders
2582  //
2584  continue ;
2585 
2586  fhPHOSClusterCutsE[3]->Fill(clus->E());
2587 
2588  //----------------------------------------------------------------------------------
2589  // Remove clusters with too low energy
2590  //
2591  if (fPHOSPtMin > fMomentum.E() || fPHOSPtMax < fMomentum.E() )
2592  continue ;
2593 
2594  fhPHOSClusterCutsE[4]->Fill(clus->E());
2595 
2596  //----------------------------------------------------
2597  // Apply N cells cut
2598  //
2599  if(clus->GetNCells() <= fPHOSNCellsCut && fDataType != AliCaloTrackReader::kMC) return ;
2600 
2601  // Check effect of n cells cut
2602  fhPHOSClusterCutsE[5]->Fill(clus->E());
2603 
2604  //----------------------------------------------------
2605  // Apply distance to bad channel cut
2606  //
2607  Double_t distBad = clus->GetDistanceToBadChannel() ; //Distance to bad channel
2608 
2609  if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2610 
2611  if(distBad < fPHOSBadChMinDist) return ;
2612 
2613  // Check effect distance to bad channel cut
2614  fhPHOSClusterCutsE[6]->Fill(clus->E());
2615 
2616  // TODO, add time cut
2617 
2618  //----------------------------------------------------------------------------------
2619  // Add selected clusters to array
2620  //
2621  //if(fDebug > 2 && fMomentum.E() > 0.1)
2622  AliDebug(2,Form("Selected clusters E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2623  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2624 
2625  fPHOSClusters->Add(clus);
2626 
2627  // TODO Dead code? remove?
2628  if (fMixedEvent)
2629  clus->SetID(iclus) ;
2630 
2631  } // esd/aod cluster loop
2632 
2633  AliDebug(1,Form("PHOS selected clusters %d",fPHOSClusters->GetEntriesFast())) ;
2634 }
2635 
2636 //____________________________________________
2638 //____________________________________________
2640 {
2641  if(fEMCALCellsListName.Length() == 0)
2642  fEMCALCells = fInputEvent->GetEMCALCells();
2643  else
2644  fEMCALCells = (AliVCaloCells*) fInputEvent->FindListObject(fEMCALCellsListName);
2645 }
2646 
2647 //___________________________________________
2649 //___________________________________________
2651 {
2652  fPHOSCells = fInputEvent->GetPHOSCells();
2653 }
2654 
2655 //_______________________________________
2658 //_______________________________________
2660 {
2661  AliVVZERO* v0 = fInputEvent->GetVZEROData();
2662  //printf("Init V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
2663 
2664  if (v0)
2665  {
2666  AliESDVZERO* esdV0 = dynamic_cast<AliESDVZERO*> (v0);
2667  for (Int_t i = 0; i < 32; i++)
2668  {
2669  if(esdV0)
2670  {//Only available in ESDs
2671  fV0ADC[0] += (Int_t)esdV0->GetAdcV0C(i);
2672  fV0ADC[1] += (Int_t)esdV0->GetAdcV0A(i);
2673  }
2674 
2675  fV0Mul[0] += (Int_t)v0->GetMultiplicityV0C(i);
2676  fV0Mul[1] += (Int_t)v0->GetMultiplicityV0A(i);
2677  }
2678 
2679  AliDebug(1,Form("ADC (%d,%d), Multiplicity (%d,%d)",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]));
2680  }
2681  else
2682  {
2683  AliDebug(1,"Cannot retrieve V0 ESD! Run w/ null V0 charges");
2684  }
2685 }
2686 
2687 //_________________________________________________
2691 //_________________________________________________
2693 {
2694  AliDebug(2,"Begin");
2695 
2696  //
2697  //check if branch name is given
2698  if(!fInputNonStandardJetBranchName.Length())
2699  {
2700  fInputEvent->Print();
2701  AliFatal("No non-standard jet branch name specified. Specify among existing ones.");
2702  }
2703 
2704  fNonStandardJets = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fInputNonStandardJetBranchName.Data()));
2705 
2706  if(!fNonStandardJets)
2707  {
2708  //check if jet branch exist; exit if not
2709  fInputEvent->Print();
2710 
2711  AliFatal(Form("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fInputNonStandardJetBranchName.Data()));
2712  }
2713  else
2714  {
2715  AliDebug(1,Form("AOD input jets %d", fNonStandardJets->GetEntriesFast()));
2716  if(GetDebug()>3) fNonStandardJets->Print("");
2717  }
2718 }
2719 
2720 //_________________________________________________
2724 //_________________________________________________
2726 {
2727  AliDebug(1,"Begin");
2728  //
2729  //check if branch name is given
2730  if(!fInputBackgroundJetBranchName.Length())
2731  {
2732  fInputEvent->Print();
2733 
2734  AliFatal("No background jet branch name specified. Specify among existing ones.");
2735  }
2736 
2737  fBackgroundJets = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fInputBackgroundJetBranchName.Data()));
2738 
2739  if(!fBackgroundJets)
2740  {
2741  //check if jet branch exist; exit if not
2742  fInputEvent->Print();
2743 
2744  AliFatal(Form("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fInputBackgroundJetBranchName.Data()));
2745  }
2746  else
2747  {
2748  AliDebug(1,"FillInputBackgroundJets");
2749  if(GetDebug()>3) fBackgroundJets->Print("");
2750  }
2751 }
2752 
2753 //____________________________________________________________________
2759 //____________________________________________________________________
2761 {
2762  // init some variables
2763  Int_t trigtimes[30], globCol, globRow,ntimes, i;
2764  Int_t absId = -1; //[100];
2765  Int_t nPatch = 0;
2766 
2767  TArrayI patches(0);
2768 
2769  // get object pointer
2770  AliVCaloTrigger *caloTrigger = GetInputEvent()->GetCaloTrigger( "EMCAL" );
2771 
2772  if(!caloTrigger)
2773  {
2774  AliError("Trigger patches input (AliVCaloTrigger) not available in data!");
2775  return patches;
2776  }
2777 
2778  //printf("CaloTrigger Entries %d\n",caloTrigger->GetEntries() );
2779 
2780  // class is not empty
2781  if( caloTrigger->GetEntries() > 0 )
2782  {
2783  // must reset before usage, or the class will fail
2784  caloTrigger->Reset();
2785 
2786  // go throuth the trigger channels
2787  while( caloTrigger->Next() )
2788  {
2789  // get position in global 2x2 tower coordinates
2790  caloTrigger->GetPosition( globCol, globRow );
2791 
2792  //L0
2793  if(IsEventEMCALL0())
2794  {
2795  // get dimension of time arrays
2796  caloTrigger->GetNL0Times( ntimes );
2797 
2798  // no L0s in this channel
2799  // presence of the channel in the iterator still does not guarantee that L0 was produced!!
2800  if( ntimes < 1 )
2801  continue;
2802 
2803  // get timing array
2804  caloTrigger->GetL0Times( trigtimes );
2805  //printf("Get L0 patch : n times %d - trigger time window %d - %d\n",ntimes, tmin,tmax);
2806 
2807  // go through the array
2808  for( i = 0; i < ntimes; i++ )
2809  {
2810  // check if in cut - 8,9 shall be accepted in 2011
2811  if( trigtimes[i] >= tmin && trigtimes[i] <= tmax )
2812  {
2813  //printf("Accepted trigger time %d \n",trigtimes[i]);
2814  //if(nTrig > 99) continue;
2815  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(globCol,globRow, absId);
2816  //printf("pass the time cut globCol %d, globRow %d absId %d\n",globCol,globRow, absId);
2817  patches.Set(nPatch+1);
2818  patches.AddAt(absId,nPatch++);
2819  }
2820  } // trigger time array
2821  }//L0
2822  else if(IsEventEMCALL1()) // L1
2823  {
2824  Int_t bit = 0;
2825  caloTrigger->GetTriggerBits(bit);
2826 
2827  Int_t sum = 0;
2828  caloTrigger->GetL1TimeSum(sum);
2829  //fBitEGA-=2;
2830  Bool_t isEGA1 = ((bit >> fBitEGA ) & 0x1) && IsEventEMCALL1Gamma1() ;
2831  Bool_t isEGA2 = ((bit >> (fBitEGA+1)) & 0x1) && IsEventEMCALL1Gamma2() ;
2832  Bool_t isEJE1 = ((bit >> fBitEJE ) & 0x1) && IsEventEMCALL1Jet1 () ;
2833  Bool_t isEJE2 = ((bit >> (fBitEJE+1)) & 0x1) && IsEventEMCALL1Jet2 () ;
2834 
2835  //if((bit>> fBitEGA )&0x1) printf("Trig Bit %d - bit %d - EG1 %d - EG2 %d\n",fBitEGA ,bit,IsEventEMCALL1Gamma1(),IsEventEMCALL1Gamma2());
2836  //if((bit>>(fBitEGA+1))&0x1) printf("Trig Bit %d - bit %d - EG1 %d - EG2 %d\n",fBitEGA+1,bit,IsEventEMCALL1Gamma1(),IsEventEMCALL1Gamma2());
2837 
2838  if(!isEGA1 && !isEJE1 && !isEGA2 && !isEJE2) continue;
2839 
2840  Int_t patchsize = -1;
2841  if (isEGA1 || isEGA2) patchsize = 2;
2842  else if (isEJE1 || isEJE2) patchsize = 16;
2843 
2844  //printf("**** Get L1 Patch: Bit %x, sum %d, patchsize %d, EGA1 %d, EGA2 %d, EJE1 %d, EJE2 %d, EGA bit %d, EJE bit %d, Trigger Gamma %d, Trigger Jet %d\n",
2845  // bit,sum,patchsize,isEGA1,isEGA2,isEJE1,isEJE2,fBitEGA,fBitEJE,IsEventEMCALL1Gamma(),IsEventEMCALL1Jet());
2846 
2847 
2848  // add 2x2 (EGA) or 16x16 (EJE) patches
2849  for(Int_t irow=0; irow < patchsize; irow++)
2850  {
2851  for(Int_t icol=0; icol < patchsize; icol++)
2852  {
2853  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(globCol+icol,globRow+irow, absId);
2854  //printf("pass the time cut globCol %d, globRow %d absId %d\n",globCol,globRow, absId);
2855  patches.Set(nPatch+1);
2856  patches.AddAt(absId,nPatch++);
2857  }
2858  }
2859 
2860  } // L1
2861 
2862  } // trigger iterator
2863  } // go through triggers
2864 
2865  if(patches.GetSize()<=0) AliInfo(Form("No patch found! for triggers: %s and selected <%s>",
2867  //else printf(">>>>> N patches %d, test %d,first %d, last %d\n",patches.GetSize(), nPatch, patches.At(0), patches.At(patches.GetSize()-1));
2868 
2869  return patches;
2870 }
2871 
2872 //____________________________________________________________
2876 //____________________________________________________________
2878 {
2879  // Init info from previous event
2880  fTriggerClusterIndex = -1;
2881  fTriggerClusterId = -1;
2882  fTriggerClusterBC = -10000;
2883  fIsExoticEvent = kFALSE;
2884  fIsBadCellEvent = kFALSE;
2885  fIsBadMaxCellEvent = kFALSE;
2886 
2887  fIsTriggerMatch = kFALSE;
2888  fIsTriggerMatchOpenCut[0] = kFALSE;
2889  fIsTriggerMatchOpenCut[1] = kFALSE;
2890  fIsTriggerMatchOpenCut[2] = kFALSE;
2891 
2892  // Do only analysis for triggered events
2893  if(!IsEventEMCALL1() && !IsEventEMCALL0())
2894  {
2895  fTriggerClusterBC = 0;
2896  return;
2897  }
2898 
2899  //printf("***** Try to match trigger to cluster %d **** L0 %d, L1 %d\n",fTriggerPatchClusterMatch,IsEventEMCALL0(),IsEventEMCALL1());
2900 
2901  //Recover the list of clusters
2902  TClonesArray * clusterList = 0;
2903  if (fInputEvent->FindListObject(fEMCALClustersListName))
2904  {
2905  clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
2906  }
2907  else if(fOutputEvent)
2908  {
2909  clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
2910  }
2911 
2912  // Get number of clusters and of trigger patches
2913  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2914  if(clusterList)
2915  nclusters = clusterList->GetEntriesFast();
2916 
2917  Int_t nPatch = patches.GetSize();
2919 
2920  //Init some variables used in the cluster loop
2921  Float_t tofPatchMax = 100000;
2922  Float_t ePatchMax =-1;
2923 
2924  Float_t tofMax = 100000;
2925  Float_t eMax =-1;
2926 
2927  Int_t clusMax =-1;
2928  Int_t idclusMax =-1;
2929  Bool_t badClMax = kFALSE;
2930  Bool_t badCeMax = kFALSE;
2931  Bool_t exoMax = kFALSE;
2932 // Int_t absIdMaxTrig= -1;
2933  Int_t absIdMaxMax = -1;
2934 
2935  Int_t nOfHighECl = 0 ;
2936 
2937  //
2938  // Check what is the trigger threshold
2939  // set minimu energym of candidate for trigger cluster
2940  //
2942 
2943  Float_t triggerThreshold = fTriggerL1EventThreshold;
2944  if(IsEventEMCALL0()) triggerThreshold = fTriggerL0EventThreshold;
2945  Float_t minE = triggerThreshold / 2.;
2946 
2947  // This method is not really suitable for JET trigger
2948  // but in case, reduce the energy cut since we do not trigger on high energy particle
2949  if(IsEventEMCALL1Jet() || minE < 1) minE = 1;
2950 
2951  AliDebug(1,Form("IsL1Trigger %d, IsL1JetTrigger? %d, IsL0Trigger %d, L1 threshold %2.1f, L0 threshold %2.1f, Min cluster E %2.2f",IsEventEMCALL1Jet(), IsEventEMCALL1(), IsEventEMCALL0(), fTriggerL1EventThreshold,fTriggerL0EventThreshold,minE));
2952 
2953  //
2954  // Loop on the clusters, check if there is any that falls into one of the patches
2955  //
2956  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2957  {
2958  AliVCluster * clus = 0;
2959  if(clusterList) clus = (AliVCluster*) clusterList->At(iclus);
2960  else clus = fInputEvent->GetCaloCluster(iclus);
2961 
2962  if ( !clus ) continue ;
2963 
2964  if ( !clus->IsEMCAL() ) continue ;
2965 
2966  //Skip clusters with too low energy to be triggering
2967  if ( clus->E() < minE ) continue ;
2968 
2969  Float_t frac = -1;
2970  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fInputEvent->GetEMCALCells(), clus,frac);
2971 
2972  Bool_t badCluster = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),
2973  clus->GetCellsAbsId(),clus->GetNCells());
2974  UShort_t cellMax[] = {(UShort_t) absIdMax};
2975  Bool_t badCell = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),cellMax,1);
2976 
2977  // if cell is bad, it can happen that time calibration is not available,
2978  // when calculating if it is exotic, this can make it to be exotic by default
2979  // open it temporarily for this cluster
2980  if(badCell)
2982 
2983  Bool_t exotic = GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCluster(clus, fInputEvent->GetEMCALCells());
2984 
2985  // Set back the time cut on exotics
2986  if(badCell)
2988 
2989  // Energy threshold for exotic Ecross typically at 4 GeV,
2990  // for lower energy, check that there are more than 1 cell in the cluster
2991  if(!exotic && clus->GetNCells() < 2) exotic = kTRUE;
2992 
2993  Float_t energy = clus->E();
2994  Int_t idclus = clus->GetID();
2995 
2996  Double_t tof = clus->GetTOF();
2997  if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn() && fTriggerClusterTimeRecal){
2998  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2999  //additional L1 phase shift
3001  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absIdMax), fInputEvent->GetBunchCrossNumber(), tof, fCurrentParIndex);
3002  }
3003  }
3004  tof *=1.e9;
3005 
3006  //printf("cluster %d, ID %d, E %2.2f, tof %2.2f, AbsId max %d, exotic %d, bad Cluster %d, bad Cell %d\n",
3007  // iclus,idclus, energy,tof,absIdMax, exotic, badCluster,badCell);
3008 
3009  // Find the highest energy cluster, avobe trigger threshold
3010  // in the event in case no match to trigger is found
3011  if( energy > eMax )
3012  {
3013  tofMax = tof;
3014  eMax = energy;
3015  badClMax = badCluster;
3016  badCeMax = badCell;
3017  exoMax = exotic;
3018  clusMax = iclus;
3019  idclusMax = idclus;
3020  absIdMaxMax = absIdMax;
3021  }
3022 
3023  // count the good clusters in the event avobe the trigger threshold
3024  // to check the exotic events
3025  if(!badCluster && !exotic)
3026  nOfHighECl++;
3027 
3028  // Find match to trigger
3029  if(fTriggerPatchClusterMatch && nPatch>0)
3030  {
3031  for(Int_t iabsId =0; iabsId < nPatch; iabsId++)
3032  {
3033  Int_t absIDCell[4];
3034  GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patches.At(iabsId), absIDCell);
3035  //if(tof > 75 ) printf("E %2.2f TOF %2.2f Trigger patch %d, cells : %d, %d, %d, %d\n",
3036  // clus->E(),tof,patches.At(iabsId), absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
3037 
3038  for(Int_t ipatch = 0; ipatch < 4; ipatch++)
3039  {
3040  if(absIdMax == absIDCell[ipatch])
3041  {
3042  //printf("*** Patches : absId %d, E %2.1f, tof %f \n",absIdMax,clus->E(), tof);
3043  if(energy > ePatchMax)
3044  {
3045  tofPatchMax = tof;
3046  ePatchMax = energy;
3047  fIsBadCellEvent = badCluster;
3048  fIsBadMaxCellEvent = badCell;
3049  fIsExoticEvent = exotic;
3050  fTriggerClusterIndex = iclus;
3051  fTriggerClusterId = idclus;
3052  fIsTriggerMatch = kTRUE;
3053 // absIdMaxTrig = absIdMax;
3054  }
3055  }
3056  }// cell patch loop
3057  }// trigger patch loop
3058  } // Do trigger patch matching
3059 
3060  }// Cluster loop
3061 
3062  // If there was no match, assign as trigger
3063  // the highest energy cluster in the event
3064  if(!fIsTriggerMatch)
3065  {
3066  tofPatchMax = tofMax;
3067  ePatchMax = eMax;
3068  fIsBadCellEvent = badClMax;
3069  fIsBadMaxCellEvent = badCeMax;
3070  fIsExoticEvent = exoMax;
3071  fTriggerClusterIndex = clusMax;
3072  fTriggerClusterId = idclusMax;
3073  }
3074 
3075  Double_t tofPatchMaxUS = TMath::Abs(tofPatchMax);
3076 
3077  if (tofPatchMaxUS < 28 ) fTriggerClusterBC = 0 ;
3078  else if(tofPatchMaxUS < 75 ) fTriggerClusterBC = 1 ;
3079  else if(tofPatchMaxUS < 125) fTriggerClusterBC = 2 ;
3080  else if(tofPatchMaxUS < 175) fTriggerClusterBC = 3 ;
3081  else if(tofPatchMaxUS < 225) fTriggerClusterBC = 4 ;
3082  else if(tofPatchMaxUS < 275) fTriggerClusterBC = 5 ;
3083  else
3084  {
3085  //printf("AliCaloTrackReader::MatchTriggerCluster() - Large BC - tof %2.3f - Index %d\n",tofPatchMaxUS,fTriggerClusterIndex);
3087  else
3088  {
3089  fTriggerClusterIndex = -2;
3090  fTriggerClusterId = -2;
3091  }
3092  }
3093 
3094  if(tofPatchMax < 0) fTriggerClusterBC*=-1;
3095 
3096 
3097  // printf("AliCaloTrackReader::MatchTriggerCluster(TArrayI patches) - Trigger cluster: index %d, ID %d, E = %2.2f, tof = %2.2f (BC = %d), bad cluster? %d, bad cell? %d, exotic? %d, patch match? %d, n High E cluster %d, absId Max %d\n",
3098  // fTriggerClusterIndex, fTriggerClusterId,ePatchMax, tofPatchMax,
3099  // fTriggerClusterBC, fIsBadCellEvent,fIsBadMaxCellEvent,fIsExoticEvent, fIsTriggerMatch, nOfHighECl,absIdMaxMax);
3100  //
3101  // if(!fIsTriggerMatch) printf("\t highest energy cluster: index %d, ID %d, E = %2.2f, tof = %2.2f, bad cluster? %d, bad cell? %d, exotic? %d\n",
3102  // clusMax, idclusMax, eMax,tofMax, badClMax, badCeMax,exoMax);
3103 
3104  //Redo matching but open cuts
3105  if(!fIsTriggerMatch && fTriggerClusterId >= 0)
3106  {
3107  // Open time patch time
3108  TArrayI patchOpen = GetTriggerPatches(7,10);
3109 
3110  Int_t patchAbsIdOpenTime = -1;
3111  for(Int_t iabsId =0; iabsId < patchOpen.GetSize(); iabsId++)
3112  {
3113  Int_t absIDCell[4];
3114  patchAbsIdOpenTime = patchOpen.At(iabsId);
3115  GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patchAbsIdOpenTime, absIDCell);
3116  //if(tof > 75 ) printf("E %2.2f TOF %2.2f Trigger patch %d, cells : %d, %d, %d, %d\n",
3117  // clus->E(),tof,patches.At(iabsId), absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
3118 
3119  for(Int_t ipatch = 0; ipatch < 4; ipatch++)
3120  {
3121  if(absIdMaxMax == absIDCell[ipatch])
3122  {
3123  fIsTriggerMatchOpenCut[0] = kTRUE;
3124  break;
3125  }
3126  }// cell patch loop
3127  }// trigger patch loop
3128 
3129  // Check neighbour patches
3130  Int_t patchAbsId = -1;
3131  Int_t globalCol = -1;
3132  Int_t globalRow = -1;
3133  GetCaloUtils()->GetEMCALGeometry()->GetFastORIndexFromCellIndex(absIdMaxMax, patchAbsId);
3134  GetCaloUtils()->GetEMCALGeometry()->GetPositionInEMCALFromAbsFastORIndex(patchAbsId,globalCol,globalRow);
3135 
3136  // Check patches with strict time cut
3137  Int_t patchAbsIdNeigh = -1;
3138  for(Int_t icol = globalCol-1; icol <= globalCol+1; icol++)
3139  {
3140  if(icol < 0 || icol > 47) continue;
3141 
3142  for(Int_t irow = globalRow; irow <= globalRow+1; irow++)
3143  {
3144  if(irow < 0 || irow > 63) continue;
3145 
3146  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(icol, irow, patchAbsIdNeigh);
3147 
3148  if ( patchAbsIdNeigh < 0 ) continue;
3149 
3150  for(Int_t iabsId =0; iabsId < patches.GetSize(); iabsId++)
3151  {
3152  if(patchAbsIdNeigh == patches.At(iabsId))
3153  {
3154  fIsTriggerMatchOpenCut[1] = kTRUE;
3155  break;
3156  }
3157  }// trigger patch loop
3158 
3159  }// row
3160  }// col
3161 
3162  // Check patches with open time cut
3163  Int_t patchAbsIdNeighOpenTime = -1;
3164  for(Int_t icol = globalCol-1; icol <= globalCol+1; icol++)
3165  {
3166  if(icol < 0 || icol > 47) continue;
3167 
3168  for(Int_t irow = globalRow; irow <= globalRow+1; irow++)
3169  {
3170  if(irow < 0 || irow > 63) continue;
3171 
3172  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(icol, irow, patchAbsIdNeighOpenTime);
3173 
3174  if ( patchAbsIdNeighOpenTime < 0 ) continue;
3175 
3176  for(Int_t iabsId =0; iabsId < patchOpen.GetSize(); iabsId++)
3177  {
3178  if(patchAbsIdNeighOpenTime == patchOpen.At(iabsId))
3179  {
3180  fIsTriggerMatchOpenCut[2] = kTRUE;
3181  break;
3182  }
3183  }// trigger patch loop
3184 
3185  }// row
3186  }// col
3187 
3188  // printf("No match, new match: Open time %d-%d, open Neigh %d-%d, both open %d-%d\n",fIsTriggerMatchOpenCut[0],patchAbsIdOpenTime,
3189  // fIsTriggerMatchOpenCut[1],patchAbsIdNeigh,
3190  // fIsTriggerMatchOpenCut[2],patchAbsIdNeighOpenTime);
3191 
3192  patchOpen.Reset();
3193 
3194  }// No trigger match found
3195  //printf("Trigger BC %d, Id %d, Index %d\n",fTriggerClusterBC,fTriggerClusterId,fTriggerClusterIndex);
3196 }
3197 
3198 //_________________________________________________________
3202 //_________________________________________________________
3204 {
3206  {
3207  // get object pointer
3208  AliVCaloTrigger *caloTrigger = GetInputEvent()->GetCaloTrigger( "EMCAL" );
3209 
3210  if ( fBitEGA == 6 )
3211  {
3212  if (IsEventEMCALL1Gamma1()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(1);
3213  else if(IsEventEMCALL1Gamma2()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(3);
3214  else if(IsEventEMCALL1Jet1 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(0);
3215  else if(IsEventEMCALL1Jet2 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(2);
3216 
3217  // printf("L1 trigger Threshold Jet1 %f, Gamma1 %f, Jet2 %f, Gamma2 %f\n",
3218  // 0.07874*caloTrigger->GetL1Threshold(0),
3219  // 0.07874*caloTrigger->GetL1Threshold(1),
3220  // 0.07874*caloTrigger->GetL1Threshold(2),
3221  // 0.07874*caloTrigger->GetL1Threshold(3));
3222  }
3223  else
3224  {
3225  // Old AOD data format, in such case, order of thresholds different!!!
3226  if (IsEventEMCALL1Gamma1()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(0);
3227  else if(IsEventEMCALL1Gamma2()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(2);
3228  else if(IsEventEMCALL1Jet1 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(1);
3229  else if(IsEventEMCALL1Jet2 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(3);
3230 
3231  // printf("L1 trigger Threshold Jet1 %f, Gamma1 %f, Jet2 %f, Gamma2 %f\n",
3232  // 0.07874*caloTrigger->GetL1Threshold(1),
3233  // 0.07874*caloTrigger->GetL1Threshold(0),
3234  // 0.07874*caloTrigger->GetL1Threshold(3),
3235  // 0.07874*caloTrigger->GetL1Threshold(2));
3236  }
3237  }
3238 
3239  // Set L0 threshold, if not set by user
3241  {
3242  // Revise for periods > LHC11d
3243  Int_t runNumber = fInputEvent->GetRunNumber();
3244  if (runNumber < 146861) fTriggerL0EventThreshold = 3. ; // LHC11a
3245  else if(runNumber < 154000) fTriggerL0EventThreshold = 4. ; // LHC11b,c
3246  else if(runNumber < 165000) fTriggerL0EventThreshold = 5.5; // LHC11c,d,e
3247  else if(runNumber < 194000) fTriggerL0EventThreshold = 2 ; // LHC12
3248  else if(runNumber < 197400) fTriggerL0EventThreshold = 3 ; // LHC13def
3249  else if(runNumber < 197400) fTriggerL0EventThreshold = 2 ; // LHC13g
3250  else if(runNumber < 244300) fTriggerL0EventThreshold = 5 ; // LHC15 in, phys 1, 5 in phys2
3251  else if(runNumber < 266400) fTriggerL0EventThreshold = 2.5; // LHC16ir
3252  else fTriggerL0EventThreshold = 3.5; // LHC16s
3253  }
3254 }
3255 
3256 //________________________________________________________
3258 //________________________________________________________
3260 {
3261  if(! opt)
3262  return;
3263 
3264  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
3265  printf("Task name : %s\n", fTaskName.Data()) ;
3266  printf("Data type : %d\n", fDataType) ;
3267  printf("CTS Min pT : %2.1f GeV/c\n", fCTSPtMin) ;
3268  printf("EMCAL Min pT : %2.1f GeV/c\n", fEMCALPtMin) ;
3269  printf("PHOS Min pT : %2.1f GeV/c\n", fPHOSPtMin) ;
3270  printf("CTS Max pT : %2.1f GeV/c\n", fCTSPtMax) ;
3271  printf("EMCAL Max pT : %2.1f GeV/c\n", fEMCALPtMax) ;
3272  printf("PHOS Max pT : %2.1f GeV/c\n", fPHOSPtMax) ;
3273  printf("EMCAL Bad Dist > %2.1f \n" , fEMCALBadChMinDist) ;
3274  printf("PHOS Bad Dist > %2.1f \n" , fPHOSBadChMinDist) ;
3275  printf("EMCAL N cells > %d \n" , fEMCALNCellsCut) ;
3276  printf("PHOS N cells > %d \n" , fPHOSNCellsCut) ;
3277  printf("EMCAL Reject cluster N cells in diff = 0, for E>%2.2f and E cell > %1.2f \n", fEMCALHighEnergyNdiffCut,fEMCALMinCellEnNdiffCut) ;
3278  printf("EMCAL Time Cut: %3.1f < TOF < %3.1f\n", fEMCALTimeCutMin, fEMCALTimeCutMax);
3279  printf("Use CTS = %d\n", fFillCTS) ;
3280  printf("Use EMCAL = %d\n", fFillEMCAL) ;
3281  printf("Use DCAL = %d\n", fFillDCAL) ;
3282  printf("Use PHOS = %d\n", fFillPHOS) ;
3283  printf("Use EMCAL Cells = %d\n", fFillEMCALCells) ;
3284  printf("Use PHOS Cells = %d\n", fFillPHOSCells) ;
3285  printf("Track status = %d\n", (Int_t) fTrackStatus) ;
3286 
3287  printf("Track Mult Eta Cut = %2.2f\n", fTrackMultEtaCut) ;
3288 
3289  printf("Track Mult Pt Cuts:") ;
3290  for(Int_t icut = 0; icut < fTrackMultNPtCut; icut++) printf(" %2.2f GeV;",fTrackMultPtCut[icut]);
3291  printf(" \n") ;
3292 
3293  printf("Write delta AOD = %d\n", fWriteOutputDeltaAOD) ;
3294  printf("Recalculate Clusters = %d, E linearity = %d\n", fRecalculateClusters, fCorrectELinearity) ;
3295 
3296  printf("Use Triggers selected in SE base class %d; If not what Trigger Mask? %d; MB Trigger Mask for mixed %d \n",
3298 
3299  if ( fComparePtHardAndJetPt )
3300  printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
3301 
3303  printf("Compare cluster pt and pt hard to accept event, factor = %2.2f",fPtHardAndClusterPtFactor);
3304 
3305  if ( fRemoveLEDEvents > 0 )
3306  {
3307  printf("Remove LED events %d, %2.1f < Ecell < %1.2f:\n", fRemoveLEDEvents, fLEDMinCellEnergy, fLEDMaxCellEnergy );
3308  printf("\t SM - nCell >= %d - Sum E >= %2.0f; \n", fLEDHighNCellsCutSM, fLEDHighEnergyCutSM);
3309  printf("\t SM3: nCell <= %d - Sum E <= %2.0f \n", fLEDLowNCellsCutSM3, fLEDLowEnergyCutSM3);
3310 
3311  if ( fRemoveLEDStripEvents > 0 )
3312  {
3313  printf("Remove LED strip? %d, with n strip > %d\n: "
3314  "\t Full SM, nCell > %d, Sum E > %2.0f;\n "
3315  "\t 1/3 SM, nCell > %d, Sum E > %2.0f;\n "
3316  "\t SM3, nCell < %d, Sum E < %2.0f\n",
3321  }
3322  }
3323 
3324  printf("Delta AOD File Name = %s\n", fDeltaAODFileName.Data()) ;
3325  printf("Centrality: Class %s, Option %d, Bin [%d,%d] \n", fCentralityClass.Data(),fCentralityOpt,fCentralityBin[0], fCentralityBin[1]) ;
3326 
3327  printf(" \n") ;
3328 }
3329 
3330 //__________________________________________
3336 //__________________________________________
3338 {
3339  // For LHC11a
3340  // Count number of cells with energy larger than 0.1 in SM3, cut on this number
3341  Int_t ncellsSM[] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
3342  Float_t ecellsSM[] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
3343 
3344  // LHC11a case
3345  if ( fRemoveLEDEvents == 1 )
3346  {
3347  for(Int_t icell = 0; icell < fInputEvent->GetEMCALCells()->GetNumberOfCells(); icell++)
3348  {
3349  Int_t absID = fInputEvent->GetEMCALCells()->GetCellNumber(icell);
3350  Int_t sm = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
3351 
3352  if ( fInputEvent->GetEMCALCells()->GetAmplitude(icell) > 0.1 && sm == 3) ncellsSM[3]++;
3353  }
3354 
3355  Int_t ncellcut = 21;
3356  if ( GetFiredTriggerClasses().Contains("EMC") ) ncellcut = 35;
3357 
3358  if ( ncellsSM[3] >= ncellcut )
3359  {
3360  AliDebug(1,Form("Reject event with ncells in SM3 %d, cut %d, trig %s",
3361  ncellsSM[3],ncellcut,GetFiredTriggerClasses().Data()));
3362  return kTRUE;
3363  }
3364  }
3365  // Run2 and general case
3366  else if ( fRemoveLEDEvents > 1 )
3367  {
3368  for(Int_t icell = 0; icell < fInputEvent->GetEMCALCells()->GetNumberOfCells(); icell++)
3369  {
3370  Int_t absID = fInputEvent->GetEMCALCells()->GetCellNumber(icell);
3371  Int_t sm = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
3372  Float_t amp = fInputEvent->GetEMCALCells()->GetAmplitude(icell);
3373 
3374  if ( amp >= fLEDMinCellEnergy && amp <= fLEDMaxCellEnergy )
3375  {
3376  ncellsSM[sm]++;
3377  ecellsSM[sm]+=amp;
3378  }
3379  }
3380 
3381  for(Int_t ism = 0; ism < 20; ism++)
3382  fhEMCALNSumEnCellsPerSM->Fill(ncellsSM[ism],ecellsSM[ism]);
3383 
3384  // Run2
3385  if ( fRemoveLEDEvents == 2 )
3386  {
3387  // if there is some activity in SM3, accept the event
3388  if ( ncellsSM[3] <= fLEDLowNCellsCutSM3 || ecellsSM[3] <= fLEDLowEnergyCutSM3 )
3389  {
3390  for(Int_t ism = 0; ism < 20; ism++)
3391  {
3392  if ( ism == 3 ) continue;
3393 
3394  if ( ncellsSM[ism] >= fLEDHighNCellsCutSM )
3395  {
3396  printf("Reject event because of SM%d: ",ism);
3397  for(Int_t jsm = 0; jsm < 20; jsm++){
3398  if ( ncellsSM[jsm] > 0 ) printf("\t SM%d: ncells %d; sum E %3.1f \n",jsm,ncellsSM[jsm],ecellsSM[jsm]);}
3399  return kTRUE;
3400  }
3401 
3402  if ( ecellsSM[ism] >= fLEDHighEnergyCutSM )
3403  {
3404  printf("Reject event because of SM%d: ",ism);
3405  for(Int_t jsm = 0; jsm < 20; jsm++) {
3406  if ( ncellsSM[jsm] > 0 ) printf("\t SM%d: ncells %d; sum E %3.1f \n",jsm,ncellsSM[jsm],ecellsSM[jsm]);}
3407  return kTRUE;
3408  }
3409  }
3410  } // SM3 activity
3411 
3412  } // fRemoveLEDEvents == 2
3413 
3414  // General case without condition on SM3 low activity
3415  else //fRemoveLEDEvents > 2
3416  {
3417  for(Int_t ism = 0; ism < GetCaloUtils()->GetEMCALGeometry()->GetNumberOfSuperModules(); ism++)
3418  {
3419  if ( ncellsSM[ism] >= fLEDHighNCellsCutSM )
3420  {
3421  printf("Reject event because of SM%d: ",ism);
3422  for(Int_t jsm = 0; jsm < 20; jsm++){
3423  if ( ncellsSM[jsm] > 0 ) printf("\t SM%d: ncells %d; sum E %3.1f \n",jsm,ncellsSM[jsm],ecellsSM[jsm]);}
3424  return kTRUE;
3425  }
3426 
3427  if ( ecellsSM[ism] >= fLEDHighEnergyCutSM )
3428  {
3429  printf("Reject event because of SM%d: ",ism);
3430  for(Int_t jsm = 0; jsm < 20; jsm++) {
3431  if ( ncellsSM[jsm] > 0 ) printf("\t SM%d: ncells %d; sum E %3.1f \n",jsm,ncellsSM[jsm],ecellsSM[jsm]);}
3432  return kTRUE;
3433  }
3434  } // SM loop
3435  } //fRemoveLEDEvents > 2
3436 
3437  for(Int_t ism = 0; ism < 20; ism++)
3438  fhEMCALNSumEnCellsPerSMAfter->Fill(ncellsSM[ism],ecellsSM[ism]);
3439 
3440  } // fRemoveLEDEvents > 1
3441 
3442  // Check activity inside all the strips
3443  // (24 strips, 2x24 cells in full SM, 2x8 cellsin 1/3 SM),
3444  // n cells (max 48) and sum of cells energy
3445  if ( fRemoveLEDStripEvents )
3446  {
3447  Float_t amp1 = 0., amp2 = 0. ;
3448  Int_t absId1 = -1, absId2 = -1 ;
3449  Int_t eventNStripActiveSM[20];
3450  Float_t enCellsStrip[20][24];
3451  Int_t nCellsStrip[20][24];
3452 
3453  for (Int_t ism = 0; ism < 20; ism++)
3454  {
3455  eventNStripActiveSM[ism] = 0;
3456 
3457  for (Int_t ieta = 0; ieta < 48; ieta=ieta+2)
3458  {
3459  enCellsStrip[ism][ieta/2] = 0.;
3460  nCellsStrip [ism][ieta/2] = 0 ;
3461 
3462  for (Int_t iphi = 0; iphi < 24; iphi++)
3463  {
3464  absId1 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(ism, iphi, ieta);
3465  if ( absId1 < 0 || absId1 > 17664 ) continue;
3466 
3467  absId2 = GetCaloUtils()->GetEMCALGeometry()->GetAbsCellIdFromCellIndexes(ism, iphi, ieta+1);
3468  if ( absId2 < 0 || absId2 > 17664 ) continue;
3469 
3470  amp1 = fInputEvent->GetEMCALCells()->GetCellAmplitude(absId1);
3471  amp2 = fInputEvent->GetEMCALCells()->GetCellAmplitude(absId2);
3472 
3473  if ( amp1 > fLEDMinCellEnergy && amp1 < fLEDMaxCellEnergy )
3474  {
3475  nCellsStrip[ism][ieta/2]++;
3476  enCellsStrip[ism][ieta/2]+=amp1;
3477 // printf("Reader: cell1 %d amp %f; SM %d, strip %d, n cells %d, sum E %f \n",
3478 // absId1, amp1,
3479 // ism,ieta,
3480 // nCellsStrip[ism][ieta/2],
3481 // enCellsStrip[ism][ieta/2]
3482 // );
3483  }
3484 
3485  if ( amp2 > fLEDMinCellEnergy && amp2 < fLEDMaxCellEnergy )
3486  {
3487  nCellsStrip[ism][ieta/2]++;
3488  enCellsStrip[ism][ieta/2]+=amp2;
3489 // printf("Reader: cell2 %d amp %f; SM %d, strip %d, n cells %d, sum E %f \n",
3490 // absId2, amp2,
3491 // ism,ieta,
3492 // nCellsStrip[ism][ieta/2],
3493 // enCellsStrip[ism][ieta/2]
3494 // );
3495  }
3496  }// iphi
3497 
3498  fhEMCALNSumEnCellsPerStrip->Fill(nCellsStrip[ism][ieta/2],enCellsStrip[ism][ieta/2]);
3499 
3500  } // ieta
3501  } // ism
3502 
3503  // Count per event over event cut
3504  // Low activity on SM3 for emin = 0.5
3505  Bool_t bSM3StripsLowActivity = kTRUE;
3506  for (Int_t ieta = 0; ieta < 24; ieta++)
3507  {
3508  if ( enCellsStrip[3][ieta] > fLEDLowEnergyCutSM3Strip ||
3509  nCellsStrip[3][ieta] > fLEDLowNCellsCutSM3Strip )
3510  bSM3StripsLowActivity = kFALSE;
3511  }
3512 
3513  // Count number or active strips, depending on cuts
3514  //
3515  Int_t nStrips = 0;
3516 
3517  if ( bSM3StripsLowActivity )
3518  {
3519  Int_t maxNCells = 24;
3520  Float_t maxECells = 80;
3521  for (Int_t ism = 0; ism < 20; ism++)
3522  {
3523  if ( ism == 3 ) continue ;
3524 
3525  maxNCells = fLEDHighNCellsCutStrip[0];
3526  maxECells = fLEDHighEnergyCutStrip[0];
3527  if ( ism == 10 || ism == 11 ||
3528  ism == 18 || ism == 19 )
3529  {
3530  maxNCells = fLEDHighNCellsCutStrip[1];
3531  maxECells = fLEDHighEnergyCutStrip[1];
3532  }
3533 
3534  for (Int_t ieta = 0; ieta < 24; ieta++)
3535  {
3536 // if ( nCellsStrip[ism][ieta] > 0 || enCellsStrip[ism][ieta] > 0 )
3537 // printf("Reader: SM %d, strip %d, n cells %d < %d , sum E %f < %f \n",
3538 // ism,ieta,
3539 // nCellsStrip[ism][ieta],maxNCells,
3540 // enCellsStrip[ism][ieta],maxECells
3541 // );
3542 
3543  if( enCellsStrip[ism][ieta] >= maxECells ||
3544  nCellsStrip[ism][ieta] >= maxNCells )
3545  {
3546  nStrips++;
3547  }
3548  } // ieta
3549  } // ism
3550  } // bSM03
3551 
3552  //printf("Reader: Number of strips %d\n",nStrips);
3553 
3554  if ( nStrips > fLEDEventMaxNumberOfStrips ) return kTRUE;
3555 
3556  for (Int_t ism = 0; ism < 20; ism++)
3557  {
3558  if ( fRemoveLEDEvents > 1 )
3559  fhEMCALNSumEnCellsPerSMAfterStripCut->Fill(ncellsSM[ism],ecellsSM[ism]);
3560 
3561  for (Int_t ieta = 0; ieta < 24; ieta++)
3562  {
3563  fhEMCALNSumEnCellsPerStripAfter->Fill(nCellsStrip[ism][ieta],enCellsStrip[ism][ieta]);
3564  }
3565  }
3566 
3567  } // remove strip LED events
3568 
3569  return kFALSE;
3570 }
3571 
3572 //_________________________________________________________
3576 //_________________________________________________________
3578 {
3579  if(label < 0) return ;
3580 
3581  AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fInputEvent) ;
3582  if(!evt) return ;
3583 
3584  TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
3585  if(!arr) return ;
3586 
3587  if(label < arr->GetEntriesFast())
3588  {
3589  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
3590  if(!particle) return ;
3591 
3592  if(label == particle->Label()) return ; // label already OK
3593  //else printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
3594  }
3595  //else printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
3596 
3597  // loop on the particles list and check if there is one with the same label
3598  for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
3599  {
3600  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
3601  if(!particle) continue ;
3602 
3603  if(label == particle->Label())
3604  {
3605  label = ind;
3606  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
3607  return;
3608  }
3609  }
3610 
3611  label = -1;
3612 
3613  //printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label not found set to -1 \n");
3614 }
3615 
3616 //___________________________________
3618 //___________________________________
3620 {
3621  if(fCTSTracks) fCTSTracks -> Clear();
3622  if(fEMCALClusters) fEMCALClusters -> Clear("C");
3623  if(fPHOSClusters) fPHOSClusters -> Clear("C");
3624 
3625  fV0ADC[0] = 0; fV0ADC[1] = 0;
3626  fV0Mul[0] = 0; fV0Mul[1] = 0;
3627 
3628  if(fNonStandardJets) fNonStandardJets -> Clear("C");
3629  if(fBackgroundJets) fBackgroundJets -> Clear("C");
3630  //fBackgroundJets->Reset();
3631 }
3632 
3633 //___________________________________________
3637 //___________________________________________
3639 {
3640  fEventTrigMinBias = kFALSE;
3641  fEventTrigCentral = kFALSE;
3642  fEventTrigSemiCentral = kFALSE;
3643  fEventTrigEMCALL0 = kFALSE;
3644  fEventTrigEMCALL1Gamma1 = kFALSE;
3645  fEventTrigEMCALL1Gamma2 = kFALSE;
3646  fEventTrigEMCALL1Jet1 = kFALSE;
3647  fEventTrigEMCALL1Jet2 = kFALSE;
3648 
3649  AliDebug(1,Form("Select trigger mask bit %d - Trigger Event %s - Select <%s>",
3651 
3652  if(fEventTriggerMask <=0 )// in case no mask set
3653  {
3654  // EMC triggered event? Which type?
3655  if( GetFiredTriggerClasses().Contains("-B-") || GetFiredTriggerClasses().Contains("-S-") || GetFiredTriggerClasses().Contains("-I-") )
3656  {
3657  if ( GetFiredTriggerClasses().Contains("EGA" ) ||
3658  GetFiredTriggerClasses().Contains("EG1" ) )
3659  {
3660  fEventTrigEMCALL1Gamma1 = kTRUE;
3661  if( GetFiredTriggerClasses().Contains("EG1" ) && !fFiredTriggerClassName.Contains("EG1") ) fEventTrigEMCALL1Gamma1 = kFALSE;
3662  }
3663  else if( GetFiredTriggerClasses().Contains("EG2" ) )
3664  {
3665  fEventTrigEMCALL1Gamma2 = kTRUE;
3666  if( !fFiredTriggerClassName.Contains("EG2") ) fEventTrigEMCALL1Gamma2 = kFALSE;
3667  }
3668  else if( GetFiredTriggerClasses().Contains("EJE" ) ||
3669  GetFiredTriggerClasses().Contains("EJ1" ) )
3670  {
3671  fEventTrigEMCALL1Jet1 = kTRUE;
3672  if( GetFiredTriggerClasses().Contains("EJ1" ) && !fFiredTriggerClassName.Contains("EJ1") )
3673  fEventTrigEMCALL1Jet1 = kFALSE;
3674  }
3675  else if( GetFiredTriggerClasses().Contains("EJ2" ) )
3676  {
3677  fEventTrigEMCALL1Jet2 = kTRUE;
3678  if( !fFiredTriggerClassName.Contains("EJ2") ) fEventTrigEMCALL1Jet2 = kFALSE;
3679  }
3680  else if( GetFiredTriggerClasses().Contains("CEMC") &&
3681  !GetFiredTriggerClasses().Contains("EGA" ) &&
3682  !GetFiredTriggerClasses().Contains("EJE" ) &&
3683  !GetFiredTriggerClasses().Contains("EG1" ) &&
3684  !GetFiredTriggerClasses().Contains("EJ1" ) &&
3685  !GetFiredTriggerClasses().Contains("EG2" ) &&
3686  !GetFiredTriggerClasses().Contains("EJ2" ) )
3687  {
3688  fEventTrigEMCALL0 = kTRUE;
3689  }
3690 
3691  //Min bias event trigger?
3692  if (GetFiredTriggerClasses().Contains("CCENT_R2-B-NOPF-ALLNOTRD"))
3693  {
3694  fEventTrigCentral = kTRUE;
3695  }
3696  else if(GetFiredTriggerClasses().Contains("CSEMI_R1-B-NOPF-ALLNOTRD"))
3697  {
3698  fEventTrigSemiCentral = kTRUE;
3699  }
3700  else if((GetFiredTriggerClasses().Contains("CINT") || GetFiredTriggerClasses().Contains("CPBI2_B1") ) &&
3701  GetFiredTriggerClasses().Contains("-NOPF-ALLNOTRD") )
3702  {
3703  fEventTrigMinBias = kTRUE;
3704  }
3705  }
3706  }
3707  else
3708  {
3709  // EMC L1 Gamma
3710  if ( fEventTriggerMask & AliVEvent::kEMCEGA )
3711  {
3712  //printf("EGA trigger bit\n");
3713  if (GetFiredTriggerClasses().Contains("EG"))
3714  {
3715  if (GetFiredTriggerClasses().Contains("EGA")) fEventTrigEMCALL1Gamma1 = kTRUE;
3716  else
3717  {
3718  if(GetFiredTriggerClasses().Contains("EG1")) fEventTrigEMCALL1Gamma1 = kTRUE;
3719  if(GetFiredTriggerClasses().Contains("EG2")) fEventTrigEMCALL1Gamma2 = kTRUE;
3720  }
3721  }
3722  }
3723  // EMC L1 Jet
3724  else if( fEventTriggerMask & AliVEvent::kEMCEJE )
3725  {
3726  //printf("EGA trigger bit\n");
3727  if (GetFiredTriggerClasses().Contains("EJ"))
3728  {
3729  if (GetFiredTriggerClasses().Contains("EJE")) fEventTrigEMCALL1Jet1 = kTRUE;
3730  else
3731  {
3732  if(GetFiredTriggerClasses().Contains("EJ1")) fEventTrigEMCALL1Jet1 = kTRUE;
3733  if(GetFiredTriggerClasses().Contains("EJ2")) fEventTrigEMCALL1Jet2 = kTRUE;
3734  }
3735  }
3736  }
3737  // EMC L0
3738  else if((fEventTriggerMask & AliVEvent::kEMC7) ||
3739  (fEventTriggerMask & AliVEvent::kEMC1) )
3740  {
3741  //printf("L0 trigger bit\n");
3742  fEventTrigEMCALL0 = kTRUE;
3743  }
3744  // Min Bias Pb-Pb
3745  else if( fEventTriggerMask & AliVEvent::kCentral )
3746  {
3747  //printf("MB semi central trigger bit\n");
3748  fEventTrigSemiCentral = kTRUE;
3749  }
3750  // Min Bias Pb-Pb
3751  else if( fEventTriggerMask & AliVEvent::kSemiCentral )
3752  {
3753  //printf("MB central trigger bit\n");
3754  fEventTrigCentral = kTRUE;
3755  }
3756  // Min Bias pp, PbPb, pPb
3757  else if((fEventTriggerMask & AliVEvent::kMB ) ||
3758  (fEventTriggerMask & AliVEvent::kINT7) ||
3759  (fEventTriggerMask & AliVEvent::kINT8) ||
3760  (fEventTriggerMask & AliVEvent::kAnyINT) )
3761  {
3762  //printf("MB trigger bit\n");
3763  fEventTrigMinBias = kTRUE;
3764  }
3765  }
3766 
3767  AliDebug(1,Form("Event bits: \n \t MB %d, Cen %d, Sem %d, L0 %d, L1G1 %d, L1G2 %d, L1J1 %d, L1J2 %d",
3771 
3772  // L1 trigger bit
3773  if( fBitEGA == 0 && fBitEJE == 0 )
3774  {
3775  // Init the trigger bit once, correct depending on AliESD(AOD)CaloTrigger header version
3776 
3777  // Simpler way to do it ...
3778 // if( fInputEvent->GetRunNumber() < 172000 )
3779 // reader->SetEventTriggerL1Bit(4,5); // current LHC11 data
3780 // else
3781 // reader->SetEventTriggerL1Bit(6,8); // LHC12-13 data
3782 
3783  // Old values
3784  fBitEGA = 4;
3785  fBitEJE = 5;
3786 
3787  TFile* file = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
3788 
3789  const TList *clist = file->GetStreamerInfoCache();
3790 
3791  if(clist)
3792  {
3793  TStreamerInfo *cinfo = (TStreamerInfo*)clist->FindObject("AliESDCaloTrigger");
3794  Int_t verid = 5; // newer ESD header version
3795  if(!cinfo)
3796  {
3797  cinfo = (TStreamerInfo*)clist->FindObject("AliAODCaloTrigger");
3798  verid = 3; // newer AOD header version
3799  }
3800 
3801  if(cinfo)
3802  {
3803  Int_t classversionid = cinfo->GetClassVersion();
3804  //printf("********* Header class version %d *********** \n",classversionid);
3805 
3806  if (classversionid >= verid)
3807  {
3808  fBitEGA = 6;
3809  fBitEJE = 8;
3810  }
3811  } else AliInfo("AliCaloTrackReader()::SetEventTriggerBit() - Streamer info for trigger class not available, bit not changed");
3812  } else AliInfo("AliCaloTrackReader::SetEventTriggerBit() - Streamer list not available!, bit not changed");
3813 
3814  } // set once the EJE, EGA trigger bit
3815 }
3816 
3817 //____________________________________________________________
3820 //____________________________________________________________
3821 void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
3822 {
3823  fInputEvent = input;
3824  fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ;
3825  if (fMixedEvent)
3826  fNMixedEvent = fMixedEvent->GetNumberOfEvents() ;
3827 
3828  //Delete previous vertex
3829  if(fVertex)
3830  {
3831  for (Int_t i = 0; i < fNMixedEvent; i++)
3832  {
3833  delete [] fVertex[i] ;
3834  }
3835  delete [] fVertex ;
3836  }
3837 
3838  fVertex = new Double_t*[fNMixedEvent] ;
3839  for (Int_t i = 0; i < fNMixedEvent; i++)
3840  {
3841  fVertex[i] = new Double_t[3] ;
3842  fVertex[i][0] = 0.0 ;
3843  fVertex[i][1] = 0.0 ;
3844  fVertex[i][2] = 0.0 ;
3845  }
3846 }
3847 
3848 
Bool_t IsPileUpFromSPD() const
TH2F * fhEMCALNSumEnCellsPerStrip
! Control histogram of LED events on strips rejection, after LED SM rejection
virtual Int_t GetDebug() const
Float_t fLEDHighEnergyCutStrip[2]
SM strip is too active if energy above this value, likely LED event. [0] Full SM, [1] 1/3 SM...
Double_t fEventWeight
Weight assigned to the event when filling histograms.
Bool_t fUseTrackDCACut
Do DCA selection.
TArrayI GetTriggerPatches(Int_t tmin, Int_t tmax)
virtual void FillInputVZERO()
virtual AliMCEvent * GetMC() const
Int_t fV0ADC[2]
Integrated V0 signal.
Bool_t IsL1PhaseInTimeRecalibrationOn() const
Bool_t fComparePtHardAndClusterPt
In MonteCarlo, jet events, reject events with too large cluster energy.
AliAnaWeights * fWeightUtils
Pointer to AliAnaWeights.
Double_t fTimeStampEventCTPBCCorrMax
Maximum value of time stamp corrected by CTP in run.
AliCalorimeterUtils * GetCaloUtils() const
TH2F * fhEMCALClusterEtaPhiFidCut
! Control histogram on EMCAL clusters acceptance, after fiducial cuts
Float_t fTimeStampEventFracMin
Minimum value of time stamp fraction event.
double Double_t
Definition: External.C:58
virtual void FillInputNonStandardJets()
Double_t fEMCALTimeCutMax
Remove clusters/cells with time larger than this value, in ns.
Int_t fBitEJE
Trigger bit on VCaloTrigger for EJE.
TLorentzVector fMomentum
! Temporal TLorentzVector container, avoid declaration of TLorentzVectors per event.
Double_t fTimeStampRunMin
Minimum value of time stamp in run.
Bool_t fDoPileUpEventRejection
Select pile-up events by SPD.
void SetCentrality(Float_t cen)
Definition: AliAnaWeights.h:61
TObjArray * fPHOSClusters
Temporal array with PHOS CaloClusters.
TString GetGeneratorNameAndIndex(Int_t index, Int_t &genIndex) const
Double_t fTrackDCACut[3]
Remove tracks with DCA larger than cut, parameters of function stored here.
Definition: External.C:236
Bool_t fUseEventsWithPrimaryVertex
Select events with primary vertex.
virtual AliVCaloCells * GetPHOSCells() const
Bool_t fIsBadCellEvent
Bad cell triggered event flag, any cell in cluster is bad.
Bool_t IsCorrectionOfClusterEnergyOn() const
Bool_t fEventTrigEMCALL1Gamma1
Event is L1-Gamma, threshold 1 on its name, it should correspond kEMCEGA.
virtual AliHeader * GetHeader() const
Bool_t IsEventEMCALL1Jet1() const
Bool_t IsEventEMCALL1Gamma1() const
Int_t fRemoveLEDEvents
Remove events where LED was wrongly firing - only EMCAL LHC11a for this equal to 1, generalized to any SM for larger.
Int_t fLEDEventMaxNumberOfStrips
Cut on events with a number of too active strips.
AliEMCALRecoUtils * GetEMCALRecoUtils() const
void GetEnergyAndNumberOfCellsInTCard(AliVCluster *clus, Int_t absIdMax, AliVCaloCells *cells, Int_t &nDiff, Int_t &nSame, Float_t &eDiff, Float_t &eSame, Float_t emin=0.)
Int_t fSmearNLMMin
Do smearing for clusters with at least this value.
TString fEventPlaneMethod
Name of event plane method, by default "Q".
Int_t fTrackBCEventCut[19]
Fill one entry per event if there is a track in a given BC, depend on track pT, acceptance cut...
Bool_t fDoVertexBCEventSelection
Select events with vertex on BC=0 or -100.
TH2F * fhEMCALClusterTimeE
! Control histogram on EMCAL timing
AliMixedEvent * fMixedEvent
! Mixed event object. This class is not the owner.
virtual AliVEvent * GetInputEvent() const
Float_t fEMCALHighEnergyNdiffCut
Minimum energy for which the cut on n diff T-Card = 0 is applied.
Bool_t IsPileUpFromSPDAndNotEMCal() const
Check if event is from pile-up determined by SPD and not by EMCal.
Bool_t IsEventEMCALL1() const
Float_t fPtHardAndClusterPtFactor
Factor between ptHard and cluster pT to reject/accept event.
Bool_t fAcceptFastCluster
Accept events from fast cluster, exclude these events for LHC11a.
energy
Definition: HFPtSpectrum.C:45
Double_t bz
static AliGenPythiaEventHeader * GetPythiaEventHeader(AliMCEvent *mcevent, TString selecHeaderName, TString &genName, TString &processName, Int_t &process, Int_t &firstParticle, Int_t &pythiaVersion)
Bool_t fSmearShowerShape
Smear shower shape (use in MC).
AliVEvent * fInputEvent
! pointer to esd or aod input.
Calculate the weight to the event to be applied when filling histograms.
Definition: AliAnaWeights.h:34
virtual void SetInputEvent(AliVEvent *input)
Bool_t IsExoticCluster(const AliVCluster *cluster, AliVCaloCells *cells, Int_t bc=0)
Double_t fTrackTimeCutMax
Remove tracks with time larger than this value, in ns.
Int_t fTriggerClusterBC
Event triggered by a cluster in BC -5 0 to 5.
Float_t fLEDMinCellEnergy
Count or sum cells energy above this value to determine if event had LEDs.
virtual AliMultSelection * GetMultSelCen() const
TString fEMCALClustersListName
Alternative list of clusters produced elsewhere and not from InputEvent.
Float_t fPHOSBadChMinDist
Minimal distance to bad channel to accept cluster in PHOS, cm.
ULong64_t GetGlobalIDPar(Short_t par)
Bool_t fWriteOutputDeltaAOD
Write the created delta AOD objects into file.
Int_t fV0Mul[2]
Integrated V0 Multiplicity.
void RecalculateClusterDistanceToBadChannel(AliVCaloCells *cells, AliVCluster *clu)
Float_t fEMCALPtMin
pT Threshold on emcal clusters.
Bool_t fTriggerClusterTimeRecal
In case cluster already calibrated, do not try to recalibrate even if recalib on in AliEMCALRecoUtils...
Double_t ** fVertex
! Vertex array 3 dim for each mixed event buffer.
TH1F * fhCTSTrackCutsPt[6]
! Control histogram on the different CTS tracks selection cuts, pT
Int_t GetDebug() const
Definition: AliAnaWeights.h:48
UInt_t fEventTriggerMask
Select this triggerered event.
virtual Int_t GetEventCentrality() const
virtual AliCentrality * GetCentrality() const
AliVCaloCells * fPHOSCells
! Temporal array with PHOS AliVCaloCells.
TList * fOutputContainer
! Output container with cut control histograms.
Bool_t IsEventEMCALL1Gamma2() const
virtual AliMixedEvent * GetMixedEvent() const
Int_t fEnergyHistogramNbins
Binning of the control histograms, min and max window.
Bool_t fFillInputBackgroundJetBranch
Flag to use data from background jets.
void SetTrackEventBC(Int_t bc)
TList * fAODBranchList
List with AOD branches created and needed in analysis.
void RemapMCLabelForAODs(Int_t &label)
Bool_t fEventTrigSemiCentral
Event is AliVEvent::kSemiCentral on its name, it should correspond to PbPb.
TRandom3 fRandom
! Random generator.
Bool_t AcceptDCA(Float_t pt, Float_t dca)
Bool_t IsPileUpFromNotSPDAndNotEMCal() const
Check if event not from pile-up determined neither by SPD nor by EMCal.
Bool_t fFillCTS
Use data from CTS.
virtual AliGenEventHeader * GetGenEventHeader() const
virtual void InitParameters()
Initialize the parameters with default.
virtual void GetVertex(Double_t v[3]) const
Float_t fEMCALMinCellEnNdiffCut
Minimum energy of cells used counting n diff in T-Card.
Bool_t IsRecalibrationOn() const
Int_t fTrackMult[10]
Track multiplicity, count for different pT cuts.
virtual Bool_t FillInputEvent(Int_t iEntry, const char *currentFileName)
Float_t fScaleFactorPerSM[22]
Scale factor depending on SM number to be applied to cluster energy.
Float_t fPHOSPtMin
pT Threshold on phos clusters.
Bool_t fSelectEmbeddedClusters
Use only simulated clusters that come from embedding.
Bool_t fEventTrigMinBias
Event is min bias on its name, it should correspond to AliVEvent::kMB, AliVEvent::kAnyInt.
Float_t fEMCALParamTimeCutMin[4]
Remove clusters/cells with time smaller than parametrized value, in ns.
Int_t fEventNumber
Event number.
void RecalculateClusterPID(AliVCluster *clu)
TString fTaskName
Name of task that executes the analysis.
Bool_t fRemoveBadTriggerEvents
Remove triggered events because trigger was exotic, bad, or out of BC.
Float_t fTrackMultPtCut[10]
Track multiplicity and sum pt cuts list.
Bool_t IsPileUpFromSPDOrEMCal() const
Check if event is from pile-up determined by SPD or EMCal.
Bool_t IsPileUpFromEMCalAndNotSPD() const
Check if event is from pile-up determined by EMCal, not by SPD.
Bool_t fEventTrigEMCALL1Gamma2
Event is L1-Gamma, threshold 2 on its name, it should correspond kEMCEGA.
Float_t fLEDLowEnergyCutSM3Strip
SM3 strip low activity if energy below this value, check activity on other SM for LED event (Run2) ...
Bool_t fEventTrigEMCALL1Jet1
Event is L1-Gamma, threshold 1 on its name, it should correspond kEMCEGA.
Double_t fTimeStampEventCTPBCCorrMin
Minimum value of time stamp corrected by CTP in run.
Bool_t fFillEMCALCells
Use data from EMCAL.
AliGenPythiaEventHeader * fGenPythiaEventHeader
! Event header casted to pythia
Float_t fEnergyHistogramLimit[2]
Binning of the control histograms, number of bins.
Int_t fTriggerClusterId
Id of trigger cluster (cluster->GetID()).
Int_t fPHOSNCellsCut
Accept for the analysis PHOS clusters with more than fNCellsCut cells.
Float_t fCTSPtMax
pT Threshold on charged particles.
Int_t GetCocktailGeneratorAndIndex(Int_t index, TString &nameGen) const
Double_t fEMCALParamTimeCutMax[4]
Remove clusters/cells with time larger than parametrized value, in ns.
virtual void FillInputEMCALCells()
Connects the array with EMCAL cells and the pointer.
Bool_t fFillPHOSCells
Use data from PHOS.
virtual AliVCaloCells * GetEMCALCells() const
TH2F * fhEMCALClusterEtaPhi
! Control histogram on EMCAL clusters acceptance, before fiducial cuts
AliEMCALGeometry * GetEMCALGeometry() const
Bool_t fIsBadMaxCellEvent
Bad cell triggered event flag, only max energy cell is bad.
int Int_t
Definition: External.C:63
Int_t fLEDHighNCellsCutSM
SM is too active if n cells above this value, likely LED event.
void RecalculateClusterShowerShapeParameters(AliVCaloCells *cells, AliVCluster *clu)
Bool_t fSelectSPDHitTracks
Ensure that track hits SPD layers.
TString fFiredTriggerClassName
Name of trigger event type used to do the analysis.
virtual TClonesArray * GetAODMCParticles() const
Definition: External.C:204
TString GetFiredTriggerClasses() const
TH2F * fhEMCALNSumEnCellsPerStripAfter
! Control histogram of LED events on strips rejection, after strip LED and SM rejection ...
unsigned int UInt_t
Definition: External.C:33
TObjArray * fEMCALClusters
Temporal array with EMCAL CaloClusters.
Bool_t IsEventEMCALL1Jet() const
AliGenEventHeader * fGenEventHeader
! Event header
TClonesArray * fBackgroundJets
! Background jets.
Float_t fPHOSPtMax
pT Threshold on phos clusters.
Bool_t fAccessTrackTOF
Access the track TOF, in case of problems when accessing GetTOFBunchCrossing.
float Float_t
Definition: External.C:68
virtual Bool_t CheckForPrimaryVertex() const
TH2F * fhEMCALNSumEnCellsPerSMAfter
! Control histogram of LED events rejection, after cut
Int_t fEMCalBCEvent[19]
Fill one entry per event if there is a cluster in a given BC.
TString fInputBackgroundJetBranchName
Name of background jet branch.
virtual AliEventplane * GetEventPlane() const
Bool_t fEventTriggerAtSE
Select triggered event at SE base task or here.
Int_t fNMCGenerToAccept
Number of MC generators that should not be included in analysis.
Float_t fEMCALPtMax
pT Threshold on emcal clusters.
Int_t fEventType
Set the event species: 7 physics, 0 MC, 8 LED (not useful now)
Float_t GetExoticCellDiffTimeCut() const
void RecalculateClusterTrackMatching(AliVEvent *event, TObjArray *clusterArray=0x0, AliMCEvent *mc=0x0)
Bool_t IsInFiducialCut(Float_t eta, Float_t phi, Int_t det) const
void SetExoticCellDiffTimeCut(Float_t dt)
virtual Bool_t ComparePtHardAndClusterPt(Int_t process, TString processName)
Bool_t fEventTrigCentral
Event is AliVEvent::kCentral on its name, it should correspond to PbPb.
Float_t fLEDLowEnergyCutSM3
SM3 low activity if energy below this value, check activity on other SM for LED event (Run2) ...
void RecalibrateCellTimeL1Phase(Int_t iSM, Int_t bc, Double_t &time, Short_t par=0) const
virtual TObjString * GetListOfParameters()
Save parameters used for analysis in a string.
Float_t fTrackMultEtaCut
Track multiplicity eta cut.
TString fMCGenerToAccept[5]
List with name of generators that should not be included.
Int_t fCentralityBin[2]
Minimum and maximum value of the centrality for the analysis.
Short_t fCurrentParIndex
! temporal PAR number based on event global to get L1 phase correction in PAR runs ...
virtual void FillInputPHOSCells()
Connects the array with PHOS cells and the pointer.
Base class for event, clusters and tracks filtering and preparation for the analysis.
TObjArray * fDCALClusters
Temporal array with DCAL CaloClusters, not needed in the normal case, use just EMCal array with DCal ...
Bool_t fUseTrackTimeCut
Do time cut selection.
Int_t GetDebug() const
virtual TString GetEventPlaneMethod() const
virtual Int_t GetTrackID(AliVTrack *track)
Bool_t fScaleEPerSM
Scale cluster energy by a constant factor, depending on SM.
TArrayI fAcceptEventsWithBit
Accept events if trigger bit is on.
unsigned long ULong_t
Definition: External.C:38
virtual Bool_t ComparePtHardAndJetPt(Int_t process, TString processName)
virtual void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
void MatchTriggerCluster(TArrayI patches)
Bool_t fTimeStampEventSelect
Select events within a fraction of data taking time.
TH2F * fhEMCALNSumEnCellsPerSM
! Control histogram of LED events rejection
Double_t fTimeStampRunMax
Maximum value of time stamp in run.
Double_t fEMCALTimeCutMin
Remove clusters/cells with time smaller than this value, in ns.
Bool_t fEventTrigEMCALL1Jet2
Event is L1-Gamma, threshold 2 on its name, it should correspond kEMCEGA.
void SetEMCalEventBCcut(Int_t bc)
Int_t fNPileUpClustersCut
Cut to select event as pile-up.
Int_t fDebug
Debugging level.
Bool_t fCheckFidCut
Do analysis for clusters in defined region.
virtual ~AliCaloTrackReader()
Destructor.
Bool_t fEventTrigEMCALL0
Event is EMCal L0 on its name, it should correspond to AliVEvent::kEMC7, AliVEvent::kEMC1.
Int_t GetNumberOfLocalMaxima(AliVCluster *cluster, AliVCaloCells *cells)
Find the number of local maxima in cluster.
Bool_t IsWeightSettingOn() const
Definition: AliAnaWeights.h:78
Int_t fTrackBCEvent[19]
Fill one entry per event if there is a track in a given BC.
Int_t fLEDHighNCellsCutStrip[2]
SM strip is too active if n cells above this value, likely LED event. [0] Full SM, [1] 1/3 SM.
virtual void FillInputEMCALAlgorithm(AliVCluster *clus, Int_t iclus)
Bool_t fComparePtHardAndJetPt
In MonteCarlo, jet events, reject fake events with wrong jet energy.
Int_t fNNonPileUpClusters
Number of clusters with time below 20 ns.
Int_t fNMixedEvent
Number of events in mixed event 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)
Int_t fCentralityOpt
Option for the returned value of the centrality, possible options 5, 10, 100.
Float_t fZvtxCut
Cut on vertex position.
short Short_t
Definition: External.C:23
TH2F * fhEMCALNSumEnCellsPerSMAfterStripCut
! Control histogram of LED events rejection, after LED strip rejection
void RecalculateClusterPosition(AliVCaloCells *cells, AliVCluster *clu)
AliCaloTrackReader()
Constructor. Initialize parameters.
Bool_t fTriggerPatchClusterMatch
Search for the trigger patch and check if associated cluster was the trigger.
Bool_t fUseParamTimeCut
Use simple or parametrized time cut.
Float_t fLEDMaxCellEnergy
Count or sum cells energy below this value to determine if event had LEDs.
Int_t GetNMaskCellColumns() const
AliFiducialCut * fFiducialCut
Acceptance cuts.
Bool_t fRejectEMCalTriggerEventsWith2Tresholds
Reject events EG2 also triggered by EG1 or EJ2 also triggered by EJ1.
Bool_t IsGoodCluster(AliVCluster *cluster, const AliEMCALGeometry *geom, AliVCaloCells *cells, Int_t bc=-1)
TH1F * fhPHOSClusterCutsE[7]
! Control histogram on the different PHOS cluster selection cuts, E
void GetMaxEnergyCell(const AliEMCALGeometry *geom, AliVCaloCells *cells, const AliVCluster *clu, Int_t &absId, Int_t &iSupMod, Int_t &ieta, Int_t &iphi, Bool_t &shared)
virtual AliMCAnalysisUtils * GetMCAnalysisUtils()
virtual void FillInputEMCAL()
Bool_t fFillEMCAL
Use data from EMCAL.
Bool_t fFillPHOS
Use data from PHOS.
AliVCaloCells * fEMCALCells
! Temporal array with EMCAL AliVCaloCells.
virtual AliAODMCHeader * GetAODMCHeader() const
Int_t fEMCALNCellsCut
Accept for the analysis EMCAL clusters with more than fNCellsCut cells.
void CorrectClusterEnergy(AliVCluster *cl)
Correct cluster energy non linearity.
Float_t fCTSPtMin
pT Threshold on charged particles.
Bool_t fMultWithEventSel
Embedded event selection in multiplicity task activated.
Int_t fLEDLowNCellsCutSM3Strip
SM3 strip low activity if n cells below this value, check activity on other SM LED event (Run2) ...
Int_t fSmearNLMMax
Do smearing for clusters with at maximum this value.
Int_t fRemoveLEDStripEvents
Remove events where an LED strip or more was wrongly firing - only EMCAL.
Float_t fTriggerL0EventThreshold
L0 Threshold to look for triggered events, set outside.
Float_t fTimeStampEventFracMax
Maximum value of time stamp fraction event.
virtual void ResetLists()
Reset lists, called in AliAnaCaloTrackCorrMaker.
Int_t fNPileUpClusters
Number of clusters with time avobe 20 ns.
Float_t fLEDHighEnergyCutSM
SM is too active if energy above this value, likely LED event.
AliAODEvent * fOutputEvent
! pointer to aod output.
Bool_t fRecalculateClusters
Correct clusters, recalculate them if recalibration parameters is given.
Bool_t fDoRejectNoTrackEvents
Reject events with no selected tracks in event.
Bool_t IsEventEMCALL1Jet2() const
Float_t fEMCALBadChMinDist
Minimal distance to bad channel to accept cluster in EMCal, cell units.
Bool_t fTimeStampEventCTPBCCorrExclude
Activate event selection within a range of data taking time CTP corrected. ESD only.
Bool_t fIsTriggerMatchOpenCut[3]
Could not match the event to a trigger patch?, retry opening cuts.
TString fEMCALCellsListName
Alternative list of cells produced elsewhere and not from InputEvent.
UInt_t fMixEventTriggerMask
Select this triggerered event for mixing, tipically kMB or kAnyINT.
Float_t RadToDeg(Float_t rad) const
TFile * file
TList with histograms for a given trigger.
virtual void FillVertexArray()
Bool_t IsInTimeWindow(Double_t tof, Float_t energy) const
Bool_t fRemoveUnMatchedTriggers
Analyze events where trigger patch and cluster where found or not.
void SetCurrentParNumber(Short_t par)
virtual Double_t GetEventPlaneAngle() const
Int_t fLEDLowNCellsCutSM3
SM3 low activity if n cells below this value, check activity on other SM LED event (Run2) ...
void SetPythiaEventHeader(AliGenPythiaEventHeader *py)
Definition: AliAnaWeights.h:83
TObjArray * fCTSTracks
Temporal array with tracks.
Int_t fTriggerPatchTimeWindow[2]
Trigger patch selection window.
unsigned short UShort_t
Definition: External.C:28
Bool_t AcceptParticleMCLabel(Int_t mcLabel) const
AliMCEvent * fMC
! Monte Carlo Event Handler.
TH1I * fhNEventsAfterCut
! Each bin represents number of events resulting after a given selection cut: vertex, trigger, ...
Float_t fPtHardAndJetPtFactor
Factor between ptHard and jet pT to reject/accept event.
TString fMCGenerEventHeaderToAccept
Accept events that contain at least this event header name.
const char Option_t
Definition: External.C:48
TH1F * fhEMCALClusterCutsE[9]
! Control histogram on the different EMCal cluster selection cuts, E
TString fDeltaAODFileName
Delta AOD file name.
Int_t GetVertexBC() const
Float_t fTrackSumPt[10]
Track sum pT, count for different pT cuts.
Bool_t IsPileUpFromEMCal() const
Check if event is from pile-up determined by EMCal.
virtual Double_t GetWeight()
Bool_t fFillInputNonStandardJetBranch
Flag to use data from non standard jets.
Bool_t IsParRun() const
TClonesArray * fNonStandardJets
! Temporal array with jets.
Int_t fEMCalBCEventCut[19]
Fill one entry per event if there is a cluster in a given BC, depend on cluster E, acceptance cut.
Int_t fTriggerClusterIndex
Index in clusters array of trigger cluster.
Bool_t fTriggerL1EventThresholdFix
L1 Threshold is fix and set outside.
bool Bool_t
Definition: External.C:53
Float_t fSmearShowerShapeWidth
Smear shower shape landau function "width" (use in MC).
Int_t fBitEGA
Trigger bit on VCaloTrigger for EGA.
Int_t fTrackMultNPtCut
Track multiplicty, number of pt cuts.
AliMCAnalysisUtils * fMCUtils
MonteCarlo Analysis utils. Initialized in SetMC()
virtual void FillInputBackgroundJets()
Bool_t fDoV0ANDEventSelection
Select events depending on V0AND.
Bool_t ClusterContainsBadChannel(const AliEMCALGeometry *geom, const UShort_t *cellList, Int_t nCells)
Class with analysis utils for simulations.
Bool_t fCorrectELinearity
Correct cluster linearity, always on.
virtual TList * GetCreateControlHistograms()
Int_t fVertexBC
Vertex BC.
Bool_t fIsExoticEvent
Exotic trigger event flag.
TString fInputNonStandardJetBranchName
Name of non standard jet branch.
TArrayI fRejectEventsWithBit
Reject events if trigger bit is on.
Bool_t fUseEMCALTimeCut
Do time cut selection.
Float_t GetPhi(Float_t phi) const
Shift phi angle in case of negative value 360 degrees. Example TLorenzVector::Phi defined in -pi to p...
Int_t fDataType
Select MC: Kinematics, Data: ESD/AOD, MCData: Both.
virtual Bool_t SelectTrack(AliVTrack *, Double_t *)
Int_t GetMaxEnergyCell(AliVCaloCells *cells, AliVCluster *clu, Float_t &fraction) const
For a given CaloCluster, it gets the absId of the cell with maximum energy deposit.
Bool_t fFillDCAL
Use data from DCAL, not needed in the normal case, use just EMCal array with DCal limits...
void SetTrackEventBCcut(Int_t bc)
Bool_t fIsTriggerMatch
Could match the event to a trigger patch?
Float_t fTriggerL1EventThreshold
L1 Threshold to look for triggered events, set in data.
TString fCentralityClass
Name of selected centrality class.
Int_t fMCGenerIndexToAccept[5]
List with index of generators that should not be included.
void RecalibrateCellTime(Int_t absId, Int_t bc, Double_t &time, Bool_t isLGon=kFALSE) const
ULong_t fTrackStatus
Track selection bit, select tracks refitted in TPC, ITS ...
Bool_t IsPileUpFromSPDAndEMCal() const
Check if event is from pile-up determined by SPD and EMCal.
virtual void FillInputPHOS()
Fill the array with PHOS filtered clusters.
virtual TObjArray * GetCTSTracks() const
Bool_t fUseAliCentrality
Select as centrality estimator AliCentrality (Run1) or AliMultSelection (Run1 and Run2) ...
Int_t fSmearingFunction
Choice of smearing function. 0 no smearing. 1 smearing from Gustavo (Landau center at 0)...
Bool_t fRecalculateVertexBC
Recalculate vertex BC from tracks pointing to vertex.
Float_t RecalibrateClusterEnergy(AliVCluster *cluster, AliVCaloCells *cells)
Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that co...
void SetEMCalEventBC(Int_t bc)
Bool_t reject
Float_t SmearClusterEnergy(const AliVCluster *clu)
Bool_t CheckCellFiducialRegion(AliVCluster *cluster, AliVCaloCells *cells) const
Bool_t IsEventEMCALL0() const