AliPhysics  d497547 (d497547)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
48 // ---- Jets ----
49 #include "AliAODJet.h"
50 #include "AliAODJetEventBackground.h"
51 
53 ClassImp(AliCaloTrackReader) ;
55 
56 //________________________________________
58 //________________________________________
60 TObject(), fEventNumber(-1), //fCurrentFileName(""),
61 fDataType(0), fDebug(0),
62 fFiducialCut(0x0), fCheckFidCut(kFALSE),
63 fComparePtHardAndJetPt(0), fPtHardAndJetPtFactor(0),
64 fComparePtHardAndClusterPt(0),fPtHardAndClusterPtFactor(0),
65 fCTSPtMin(0), fEMCALPtMin(0), fPHOSPtMin(0),
66 fCTSPtMax(0), fEMCALPtMax(0), fPHOSPtMax(0),
67 fEMCALBadChMinDist(0), fPHOSBadChMinDist (0),
68 fEMCALNCellsCut(0), fPHOSNCellsCut(0),
69 fUseEMCALTimeCut(1), fUseParamTimeCut(0),
70 fUseTrackTimeCut(0), fAccessTrackTOF(0),
71 fEMCALTimeCutMin(-10000), fEMCALTimeCutMax(10000),
72 fEMCALParamTimeCutMin(), fEMCALParamTimeCutMax(),
73 fTrackTimeCutMin(-10000), fTrackTimeCutMax(10000),
74 fUseTrackDCACut(0),
75 fAODBranchList(0x0),
76 fCTSTracks(0x0), fEMCALClusters(0x0),
77 fDCALClusters(0x0), fPHOSClusters(0x0),
78 fEMCALCells(0x0), fPHOSCells(0x0),
79 fInputEvent(0x0), fOutputEvent(0x0),fMC(0x0),
80 fFillCTS(0), fFillEMCAL(0),
81 fFillDCAL(0), fFillPHOS(0),
82 fFillEMCALCells(0), fFillPHOSCells(0),
83 fRecalculateClusters(kFALSE),fCorrectELinearity(kTRUE),
84 fSelectEmbeddedClusters(kFALSE),
85 fSmearShowerShape(0), fSmearShowerShapeWidth(0), fRandom(),
86 fSmearingFunction(0), fSmearNLMMin(0), fSmearNLMMax(0),
87 fTrackStatus(0), fSelectSPDHitTracks(0),
88 fTrackMultNPtCut(0), fTrackMultEtaCut(0.9),
89 fDeltaAODFileName(""), fFiredTriggerClassName(""),
90 
91 fEventTriggerMask(0), fMixEventTriggerMask(0), fEventTriggerAtSE(0),
92 fEventTrigMinBias(0), fEventTrigCentral(0),
93 fEventTrigSemiCentral(0), fEventTrigEMCALL0(0),
94 fEventTrigEMCALL1Gamma1(0), fEventTrigEMCALL1Gamma2(0),
95 fEventTrigEMCALL1Jet1(0), fEventTrigEMCALL1Jet2(0),
96 fBitEGA(0), fBitEJE(0),
97 
98 fEventType(-1),
99 fTaskName(""), fCaloUtils(0x0),
100 fWeightUtils(0x0), fEventWeight(1),
101 fMixedEvent(NULL), fNMixedEvent(0), fVertex(NULL),
102 fListMixedTracksEvents(), fListMixedCaloEvents(),
103 fLastMixedTracksEvent(-1), fLastMixedCaloEvent(-1),
104 fWriteOutputDeltaAOD(kFALSE),
105 fEMCALClustersListName(""), fZvtxCut(0.),
106 fAcceptFastCluster(kFALSE), fRemoveLEDEvents(0),
107 //Trigger rejection
108 fRemoveBadTriggerEvents(0), fTriggerPatchClusterMatch(0),
109 fTriggerPatchTimeWindow(), fTriggerL0EventThreshold(0),
110 fTriggerL1EventThreshold(0), fTriggerL1EventThresholdFix(0),
111 fTriggerClusterBC(0), fTriggerClusterIndex(0), fTriggerClusterId(0),
112 fIsExoticEvent(0), fIsBadCellEvent(0), fIsBadMaxCellEvent(0),
113 fIsTriggerMatch(0), fIsTriggerMatchOpenCut(),
114 fTriggerClusterTimeRecal(kTRUE), fRemoveUnMatchedTriggers(kTRUE),
115 fDoPileUpEventRejection(kFALSE), fDoV0ANDEventSelection(kFALSE),
116 fDoVertexBCEventSelection(kFALSE),
117 fDoRejectNoTrackEvents(kFALSE),
118 fUseEventsWithPrimaryVertex(kFALSE),
119 //fTriggerAnalysis (0x0),
120 fTimeStampEventSelect(0),
121 fTimeStampEventFracMin(0), fTimeStampEventFracMax(0),
122 fTimeStampRunMin(0), fTimeStampRunMax(0),
123 fNPileUpClusters(-1), fNNonPileUpClusters(-1), fNPileUpClustersCut(3),
124 fVertexBC(-200), fRecalculateVertexBC(0),
125 fUseAliCentrality(0), fCentralityClass(""), fCentralityOpt(0),
126 fEventPlaneMethod(""),
127 fFillInputNonStandardJetBranch(kFALSE),
128 fNonStandardJets(new TClonesArray("AliAODJet",100)), fInputNonStandardJetBranchName("jets"),
129 fFillInputBackgroundJetBranch(kFALSE),
130 fBackgroundJets(0x0),fInputBackgroundJetBranchName("jets"),
131 fAcceptEventsWithBit(0), fRejectEventsWithBit(0), fRejectEMCalTriggerEventsWith2Tresholds(0),
132 fMomentum(), fOutputContainer(0x0),
133 fhEMCALClusterEtaPhi(0), fhEMCALClusterTimeE(0),
134 fEnergyHistogramNbins(0),
135 fhNEventsAfterCut(0), fNMCGenerToAccept(0), fMCGenerEventHeaderToAccept("")
136 {
137  for(Int_t i = 0; i < 8; i++) fhEMCALClusterCutsE [i]= 0x0 ;
138  for(Int_t i = 0; i < 7; i++) fhPHOSClusterCutsE [i]= 0x0 ;
139  for(Int_t i = 0; i < 6; i++) fhCTSTrackCutsPt [i]= 0x0 ;
140  for(Int_t j = 0; j < 5; j++) { fMCGenerToAccept [j] = ""; fMCGenerIndexToAccept[j] = -1; }
141 
142  InitParameters();
143 }
144 
145 //_______________________________________
147 //_______________________________________
149 {
150  DeletePointers();
151 }
152 
153 //_______________________________________
156 //_______________________________________
158 {
159  delete fFiducialCut ;
160 
161  if(fAODBranchList)
162  {
163  fAODBranchList->Delete();
164  delete fAODBranchList ;
165  }
166 
167  if(fCTSTracks)
168  {
169  if(fDataType!=kMC)fCTSTracks->Clear() ;
170  else fCTSTracks->Delete() ;
171  delete fCTSTracks ;
172  }
173 
174  if(fEMCALClusters)
175  {
176  if(fDataType!=kMC)fEMCALClusters->Clear("C") ;
177  else fEMCALClusters->Delete() ;
178  delete fEMCALClusters ;
179  }
180 
181  if(fDCALClusters)
182  {
183  if(fDataType!=kMC)fDCALClusters->Clear("C") ;
184  else fDCALClusters->Delete() ;
185  delete fDCALClusters ;
186  }
187 
188  if(fPHOSClusters)
189  {
190  if(fDataType!=kMC)fPHOSClusters->Clear("C") ;
191  else fPHOSClusters->Delete() ;
192  delete fPHOSClusters ;
193  }
194 
195  if(fVertex)
196  {
197  for (Int_t i = 0; i < fNMixedEvent; i++)
198  {
199  delete [] fVertex[i] ;
200 
201  }
202  delete [] fVertex ;
203  }
204 
205  //delete fTriggerAnalysis;
206 
207  if(fNonStandardJets)
208  {
209  if(fDataType!=kMC) fNonStandardJets->Clear("C") ;
210  else fNonStandardJets->Delete() ;
211  delete fNonStandardJets ;
212  }
213  delete fBackgroundJets ;
214 
215  fRejectEventsWithBit.Reset();
216  fAcceptEventsWithBit.Reset();
217 
218  if ( fWeightUtils ) delete fWeightUtils ;
219 
220  // Pointers not owned, done by the analysis frame
221  // if(fInputEvent) delete fInputEvent ;
222  // if(fOutputEvent) delete fOutputEvent ;
223  // if(fMC) delete fMC ;
224  // Pointer not owned, deleted by maker
225  // if (fCaloUtils) delete fCaloUtils ;
226 }
227 
228 //____________________________________________________________
232 //____________________________________________________________
234 {
235  Float_t cut = fTrackDCACut[0]+fTrackDCACut[1]/TMath::Power(pt,fTrackDCACut[2]);
236 
237  if(TMath::Abs(dca) < cut)
238  return kTRUE;
239  else
240  return kFALSE;
241 }
242 
243 //_____________________________________________________
246 //_____________________________________________________
248 {
249  Int_t nAccept = fAcceptEventsWithBit.GetSize();
250 
251  //printf("N accept %d\n", nAccept);
252 
253  if( nAccept <= 0 )
254  return kTRUE ; // accept the event
255 
256  UInt_t trigFired = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
257 
258  for(Int_t ibit = 0; ibit < nAccept; ibit++)
259  {
260  Bool_t accept = (trigFired & fAcceptEventsWithBit.At(ibit));
261 
262  //printf("accept %d, ibit %d, bit %d \n",accept, ibit,fAcceptEventsWithBit.At(ibit));
263  if(accept) return kTRUE ; // accept the event
264  }
265 
266  return kFALSE ; // reject the event
267 }
268 
269 //_____________________________________________________
274 //_____________________________________________________
276 {
277  if( !fMC || fNMCGenerToAccept <= 0 ) return kTRUE;
278 
279  TString genName;
280  Int_t genIndex;
281  genIndex = GetCocktailGeneratorAndIndex(mcLabel, genName);
282  //fMC->GetCocktailGenerator(mcLabel,genName);
283 
284  Bool_t generOK = kFALSE;
285  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
286  {
287  if ( fMCGenerToAccept[ig].Contains(genName) ) generOK = kTRUE;
288 
289  if ( generOK && fMCGenerIndexToAccept[ig] >= 0 && fMCGenerToAccept[ig] != genIndex) generOK = kFALSE;
290  }
291 
292  if ( !generOK ) AliDebug(1, Form("skip label %d, gen %s",mcLabel,genName.Data()) );
293 
294  return generOK;
295 }
296 
297 //_____________________________________________________________________
305 //_____________________________________________________________________
307 {
308  //method that gives the generator for a given particle with label index (or that of the corresponding primary)
309  AliVParticle* mcpart0 = (AliVParticle*) GetMC()->GetTrack(index);
310  Int_t genIndex = -1;
311 
312  if(!mcpart0)
313  {
314  printf("AliMCEvent-BREAK: No valid AliMCParticle at label %i\n",index);
315  return -1;
316  }
317 
318  nameGen = GetGeneratorNameAndIndex(index,genIndex);
319 
320  if(nameGen.Contains("nococktailheader") ) return -1;
321 
322  Int_t lab=index;
323 
324  while(nameGen.IsWhitespace())
325  {
326  AliVParticle* mcpart = (AliVParticle*) GetMC()->GetTrack(lab);
327 
328  if(!mcpart)
329  {
330  printf("AliMCEvent-BREAK: No valid AliMCParticle at label %i\n",lab);
331  break;
332  }
333 
334  Int_t mother=0;
335  mother = mcpart->GetMother();
336 
337  if(mother<0)
338  {
339  printf("AliMCEvent - BREAK: Reached primary particle without valid mother\n");
340  break;
341  }
342 
343  AliVParticle* mcmom = (AliVParticle*) GetMC()->GetTrack(mother);
344  if(!mcmom)
345  {
346  printf("AliMCEvent-BREAK: No valid AliMCParticle mother at label %i\n",mother);
347  break;
348  }
349 
350  lab=mother;
351 
352  nameGen = GetGeneratorNameAndIndex(mother,genIndex);
353  }
354 
355  return genIndex;
356 }
357 //_____________________________________________________________________
365 //_____________________________________________________________________
367 {
368  Int_t nsumpart = GetMC()->GetNumberOfPrimaries();
369 
370  genIndex = -1;
371 
372  TList* lh = GetMC()->GetCocktailList();
373  if(!lh)
374  {
375  TString noheader="nococktailheader";
376  return noheader;
377  }
378 
379  Int_t nh = lh->GetEntries();
380 
381  for (Int_t i = nh-1; i >= 0; i--)
382  {
383  AliGenEventHeader* gh = (AliGenEventHeader*)lh->At(i);
384 
385  TString genname = gh->GetName();
386 
387  Int_t npart=gh->NProduced();
388 
389  if (i == 0) npart = nsumpart;
390 
391  if(index < nsumpart && index >= (nsumpart-npart))
392  {
393  genIndex = i ;
394  return genname;
395  }
396 
397  nsumpart-=npart;
398  }
399 
400  TString empty="";
401  return empty;
402 }
403 
404 
405 //_____________________________________________________
408 //_____________________________________________________
410 {
411  Int_t nReject = fRejectEventsWithBit.GetSize();
412 
413  //printf("N reject %d\n", nReject);
414 
415  if( nReject <= 0 )
416  return kTRUE ; // accept the event
417 
418  UInt_t trigFired = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
419 
420  for(Int_t ibit = 0; ibit < nReject; ibit++)
421  {
422  Bool_t reject = (trigFired & fRejectEventsWithBit.At(ibit));
423 
424  //printf("reject %d, ibit %d, bit %d \n",reject, ibit,fRejectEventsWithBit.At(ibit));
425  if(reject) return kFALSE ; // reject the event
426  }
427 
428  return kTRUE ; // accept the event
429 }
430 
431 //_____________________________________________
435 //_____________________________________________
437 {
438  //-----------------------------------------------------------
439  // Reject events depending on the event species type
440  //-----------------------------------------------------------
441 
442  // Event types:
443  // kStartOfRun = 1, // START_OF_RUN
444  // kEndOfRun = 2, // END_OF_RUN
445  // kStartOfRunFiles = 3, // START_OF_RUN_FILES
446  // kEndOfRunFiles = 4, // END_OF_RUN_FILES
447  // kStartOfBurst = 5, // START_OF_BURST
448  // kEndOfBurst = 6, // END_OF_BURST
449  // kPhysicsEvent = 7, // PHYSICS_EVENT
450  // kCalibrationEvent = 8, // CALIBRATION_EVENT
451  // kFormatError = 9, // EVENT_FORMAT_ERROR
452  // kStartOfData = 10, // START_OF_DATA
453  // kEndOfData = 11, // END_OF_DATA
454  // kSystemSoftwareTriggerEvent = 12, // SYSTEM_SOFTWARE_TRIGGER_EVENT
455  // kDetectorSoftwareTriggerEvent = 13 // DETECTOR_SOFTWARE_TRIGGER_EVENT
456 
457  Int_t eventType = 0;
458  if(fInputEvent->GetHeader())
459  eventType = ((AliVHeader*)fInputEvent->GetHeader())->GetEventType();
460 
461  AliDebug(3,Form("Event type %d",eventType));
462 
463  // Select only Physics events in data, eventType = 7,
464  // usually done by physics selection
465  // MC not set, eventType =0, I think, so do not apply a value to fEventType by default
466  // LED events have eventType = 8, implemented a selection cut in the past
467  // but not really useful for EMCal data since LED are not reconstructed (unless wrongly assumed as physics)
468 
469  if ( fEventType >= 0 && eventType != fEventType ) return kFALSE ;
470 
471  AliDebug(1,"Pass event species selection");
472 
473  fhNEventsAfterCut->Fill(1.5);
474 
475  //-----------------------------------------------------------------
476  // In case of mixing analysis, select here the trigger of the event
477  //-----------------------------------------------------------------
478 
479  UInt_t isTrigger = kFALSE;
480  UInt_t isMB = kFALSE;
481 
482  if(!fEventTriggerAtSE)
483  {
484  // In case of mixing analysis, accept MB events, not only Trigger
485  // Track and cluster arrays filled for MB in order to create the pool in the corresponding analysis
486  // via de method in the base class FillMixedEventPool()
487 
488  AliAnalysisManager *manager = AliAnalysisManager::GetAnalysisManager();
489  AliInputEventHandler *inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
490 
491  if(!inputHandler) return kFALSE ; // to content coverity
492 
493  isTrigger = inputHandler->IsEventSelected() & fEventTriggerMask;
494  isMB = inputHandler->IsEventSelected() & fMixEventTriggerMask;
495 
496  if(!isTrigger && !isMB) return kFALSE;
497 
498  //printf("Selected triggered event : %s\n",GetFiredTriggerClasses().Data());
499  AliDebug(1,"Pass uninteresting triggered events rejection in case of mixing analysis");
500 
501  fhNEventsAfterCut->Fill(2.5);
502  }
503 
504  //-----------------------------------------------------------
505  // Reject events depending on the trigger name
506  // Careful!, if a special MB event string is selected but the option
507  // to select events via the mask in the reader is done, it will not
508  // be taken into account.
509  //-----------------------------------------------------------
510 
511  AliDebug(1,Form("FiredTriggerClass <%s>, selected class <%s>, compare name %d",
514 
515  if ( fFiredTriggerClassName != "" && !isMB )
516  {
518  return kFALSE;
519 
520  AliDebug(1,"Accepted triggered event");
521 
522  fhNEventsAfterCut->Fill(3.5);
523  }
524 
525  //-------------------------------------------------------------------------------------
526  // Reject or accept events depending on the trigger bit
527  //-------------------------------------------------------------------------------------
528 
531 
532  //printf("AliCaloTrackReader::FillInputEvent() - Accept event? %d, Reject event %d? \n",okA,okR);
533 
534  if(!okA || !okR) return kFALSE;
535 
536  AliDebug(1,"Pass event bit rejection");
537 
538  fhNEventsAfterCut->Fill(4.5);
539 
540  //----------------------------------------------------------------------
541  // Do not count events that were likely triggered by an exotic cluster
542  // or out BC cluster
543  //----------------------------------------------------------------------
544 
545  // Set a bit with the event kind, MB, L0, L1 ...
547 
548  // In case of Mixing, avoid checking the triggers in the min bias events
549  if(!fEventTriggerAtSE && (isMB && !isTrigger)) return kTRUE;
550 
552  {
554  {
555  // Reject triggered events when there is coincidence on both EMCal trigger thresholds,
556  // but the requested trigger is the low trigger threshold
557  if(IsEventEMCALL1Jet1 () && IsEventEMCALL1Jet2 () && fFiredTriggerClassName.Contains("EJ2")) return kFALSE;
558  if(IsEventEMCALL1Gamma1() && IsEventEMCALL1Gamma2() && fFiredTriggerClassName.Contains("EG2")) return kFALSE;
559  }
560 
561  //Get Patches that triggered
563 
564  MatchTriggerCluster(patches);
565 
566  patches.Reset();
567 
568  // If requested, remove badly triggeed events, but only when the EMCal trigger bit is set
570  {
571  AliDebug(1,Form("ACCEPT triggered event? \n exotic? %d - bad cell %d - bad Max cell %d - BC %d - Matched %d\n",
573 
574  if (fIsExoticEvent) return kFALSE;
575  else if(fIsBadCellEvent) return kFALSE;
576  else if(fRemoveUnMatchedTriggers && !fIsTriggerMatch) return kFALSE ;
577  else if(fTriggerClusterBC != 0) return kFALSE;
578 
579  AliDebug(1,Form("\t *** YES for %s",GetFiredTriggerClasses().Data()));
580  }
581 
582  AliDebug(1,"Pass EMCal triggered event rejection \n");
583 
584  fhNEventsAfterCut->Fill(5.5);
585  }
586 
587  //-------------------------------------------------------------------------------------
588  // Select events only fired by a certain trigger configuration if it is provided
589  //-------------------------------------------------------------------------------------
590 
591  if (GetFiredTriggerClasses().Contains("FAST") && !GetFiredTriggerClasses().Contains("ALL") && !fAcceptFastCluster)
592  {
593  AliDebug(1,Form("Do not count events from fast cluster, trigger name %s\n",fFiredTriggerClassName.Data()));
594  return kFALSE;
595 
596  fhNEventsAfterCut->Fill(6.5);
597  }
598 
599  //-------------------------------------------------------------------------------------
600  // Reject event if large clusters with large energy
601  // Use only for LHC11a data for the moment, and if input is clusterizer V1 or V1+unfolding
602  // If clusterzer NxN or V2 it does not help
603  //-------------------------------------------------------------------------------------
604 
605  //Int_t run = fInputEvent->GetRunNumber();
606 
607  if ( fRemoveLEDEvents )
608  {
610 
611  if(reject) return kFALSE;
612 
613  AliDebug(1,"Pass LED event rejection");
614 
615  fhNEventsAfterCut->Fill(7.5);
616  } // Remove LED events
617 
618  // All selection criteria passed, accept the event
619  return kTRUE ;
620 }
621 
622 //________________________________________________
627 //________________________________________________
629 {
630  //printf("AliCaloTrackReader::ComparePtHardAndJetPt() - GenHeaderName : %s\n",GetGenEventHeader()->ClassName());
631 
632  if ( !GetGenEventHeader() )
633  {
634  AliError("Skip event, event header is not available!");
635  return kFALSE;
636  }
637 
638  if ( !strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader") )
639  {
640  AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
641 
642  // Do this check only for gamma-jet or jet-jet productions
643  Int_t process = pygeh->ProcessType();
644  //printf("Process %d \n", process);
645  if ( process != 14 && process != 18 && process != 29 && process != 114 && process != 115 && // prompt gamma
646  process != 11 && process != 12 && process != 13 && process != 28 && process != 53 && // jet-jet
647  process != 68 && process != 96 // jet-jet
648  ) return kTRUE ;
649  //printf("\t yes \n");
650 
651  Int_t nTriggerJets = pygeh->NTriggerJets();
652  Float_t ptHard = pygeh->GetPtHard();
653  TParticle * jet = 0;
654 
655  AliDebug(1,Form("Njets: %d, pT Hard %f",nTriggerJets, ptHard));
656 
657  Float_t tmpjet[]={0,0,0,0};
658  for(Int_t ijet = 0; ijet< nTriggerJets; ijet++)
659  {
660  pygeh->TriggerJet(ijet, tmpjet);
661  jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
662 
663  AliDebug(1,Form("jet %d; pycell jet pT %f",ijet, jet->Pt()));
664 
665  //Compare jet pT and pt Hard
666  if(jet->Pt() > fPtHardAndJetPtFactor * ptHard)
667  {
668  AliInfo(Form("Reject jet event with : process %d, pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
669  process, ptHard, jet->Pt(), fPtHardAndJetPtFactor));
670  return kFALSE;
671  }
672  }
673 
674  if(jet) delete jet;
675  }
676 
677  return kTRUE ;
678 }
679 
680 //____________________________________________________
685 //____________________________________________________
687 {
688  if ( !GetGenEventHeader() )
689  {
690  AliError("Skip event, event header is not available!");
691  return kFALSE;
692  }
693 
694  if ( !strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader") )
695  {
696  AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
697 
698  // Do this check only for gamma-jet productions
699  Int_t process = pygeh->ProcessType();
700  //printf("Process %d \n", process);
701  if ( process != 14 && process != 18 && process != 29 && process != 114 && process != 115 // && // prompt gamma
702  //process != 11 && process != 12 && process != 13 && process != 28 && process != 53 && // jet-jet
703  //process != 68 && process != 96 // jet-jet
704  ) return kTRUE ;
705  //printf("\t yes \n");
706 
707  Float_t ptHard = pygeh->GetPtHard();
708 
709  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
710  for (Int_t iclus = 0; iclus < nclusters; iclus++)
711  {
712  AliVCluster * clus = fInputEvent->GetCaloCluster(iclus) ;
713  Float_t ecluster = clus->E();
714 
715  if(ecluster > fPtHardAndClusterPtFactor*ptHard)
716  {
717  AliInfo(Form("Reject : process %d, ecluster %2.2f, calo %d, factor %2.2f, ptHard %f",
718  process, ecluster,clus->GetType(),fPtHardAndClusterPtFactor,ptHard));
719 
720  return kFALSE;
721  }
722  }
723 
724  }
725 
726  return kTRUE ;
727 }
728 
729 //___________________________________________________
732 //___________________________________________________
734 {
735  fhNEventsAfterCut = new TH1I("hNEventsAfterCut", "Number of analyzed events", 20, 0, 20) ;
736  //fhNEventsAfterCut->SetXTitle("Selection");
737  fhNEventsAfterCut->SetYTitle("# events");
738  fhNEventsAfterCut->GetXaxis()->SetBinLabel(1 ,"1=Input");
739  fhNEventsAfterCut->GetXaxis()->SetBinLabel(2 ,"2=Event Type");
740  fhNEventsAfterCut->GetXaxis()->SetBinLabel(3 ,"3=Mixing Event");
741  fhNEventsAfterCut->GetXaxis()->SetBinLabel(4 ,"4=Trigger string");
742  fhNEventsAfterCut->GetXaxis()->SetBinLabel(5 ,"5=Trigger Bit");
743  fhNEventsAfterCut->GetXaxis()->SetBinLabel(6 ,"6=Good EMC Trigger");
744  fhNEventsAfterCut->GetXaxis()->SetBinLabel(7 ,"7=!Fast Cluster");
745  fhNEventsAfterCut->GetXaxis()->SetBinLabel(8 ,"8=!LED");
746  fhNEventsAfterCut->GetXaxis()->SetBinLabel(9 ,"9=Time stamp");
747  fhNEventsAfterCut->GetXaxis()->SetBinLabel(10,"10=Primary vertex");
748  fhNEventsAfterCut->GetXaxis()->SetBinLabel(11,"11=Null 3 vertex");
749  fhNEventsAfterCut->GetXaxis()->SetBinLabel(12,"12=Z vertex window");
750  fhNEventsAfterCut->GetXaxis()->SetBinLabel(13,"13=Pile-up");
751  fhNEventsAfterCut->GetXaxis()->SetBinLabel(14,"14=V0AND");
752  fhNEventsAfterCut->GetXaxis()->SetBinLabel(15,"15=Centrality");
753  fhNEventsAfterCut->GetXaxis()->SetBinLabel(16,"16=GenHeader");
754  fhNEventsAfterCut->GetXaxis()->SetBinLabel(17,"17=PtHard-Jet");
755  fhNEventsAfterCut->GetXaxis()->SetBinLabel(18,"18=PtHard-Cluster");
756  fhNEventsAfterCut->GetXaxis()->SetBinLabel(19,"19=N Track>0");
757  fhNEventsAfterCut->GetXaxis()->SetBinLabel(20,"20=TOF BC");
759 
760  if(fFillEMCAL)
761  {
762  for(Int_t i = 0; i < 8; i++)
763  {
764  TString names[] =
765  { "NoCut", "Corrected", "GoodCluster", "NonLinearity",
766  "EnergyAndFidutial", "NCells", "BadDist", "Time" } ;
767 
768  fhEMCALClusterCutsE[i] = new TH1F(Form("hEMCALReaderClusterCuts_%d_%s",i,names[i].Data()),
769  Form("EMCal %d, %s",i,names[i].Data()),
771  fhEMCALClusterCutsE[i]->SetYTitle("# clusters");
772  fhEMCALClusterCutsE[i]->SetXTitle("#it{E} (GeV)");
774  }
775 
777  ("hEMCALReaderTimeE","time vs #it{E} after cuts (if no calib, shifted -615 ns)", 100,0,100,400,-400,400);
778  fhEMCALClusterTimeE->SetXTitle("#it{E} (GeV)");
779  fhEMCALClusterTimeE->SetYTitle("#it{time} (ns)");
781 
783  ("hEMCALReaderEtaPhi","#eta vs #varphi",40,-2, 2,50, 0,10);
784  // Very open limits to check problems
785  fhEMCALClusterEtaPhi->SetXTitle("#eta");
786  fhEMCALClusterEtaPhi->SetYTitle("#varphi (rad)");
788  }
789 
790  if(fFillPHOS)
791  {
792  for(Int_t i = 0; i < 7; i++)
793  {
794  TString names[] = {"NoCut", "ExcludeCPV", "BorderCut", "FiducialCut", "EnergyCut", "NCells", "BadDist"};
795 
796  fhPHOSClusterCutsE[i] = new TH1F(Form("hPHOSReaderClusterCuts_%d_%s",i,names[i].Data()),
797  Form("PHOS Cut %d, %s",i,names[i].Data()),
799  fhPHOSClusterCutsE[i]->SetYTitle("# clusters");
800  fhPHOSClusterCutsE[i]->SetXTitle("#it{E} (GeV)");
802  }
803  }
804 
805  if(fFillCTS)
806  {
807  for(Int_t i = 0; i < 6; i++)
808  {
809  TString names[] = {"NoCut", "Status", "ESD_AOD", "TOF", "DCA","PtAcceptanceMult"};
810 
811  fhCTSTrackCutsPt[i] = new TH1F(Form("hCTSReaderClusterCuts_%d_%s",i,names[i].Data()),
812  Form("CTS Cut %d, %s",i,names[i].Data()),
814  fhCTSTrackCutsPt[i]->SetYTitle("# tracks");
815  fhCTSTrackCutsPt[i]->SetXTitle("#it{p}_{T} (GeV)");
817  }
818  }
819 
820  return fOutputContainer ;
821 }
822 
823 //_____________________________________________________
825 //_____________________________________________________
827 {
828  TString parList ; //this will be list of parameters used for this analysis.
829 
830  const Int_t buffersize = 255;
831  char onePar[buffersize] ;
832 
833  snprintf(onePar,buffersize,"--- Reader ---:") ;
834  parList+=onePar ;
835  snprintf(onePar,buffersize,"Data type: %d; zvertex cut %2.2f; EMC cluster name: <%s> ",fDataType, fZvtxCut, fEMCALClustersListName.Data()) ;
836  parList+=onePar ;
837  snprintf(onePar,buffersize,"Use detector: EMC %d, DCA %d, PHOS %d, CTS %d, EMCcells %d, PHOScells %d ; ",
839  snprintf(onePar,buffersize,"E-pT window: EMC (%2.1f,%2.1f), PHOS (%2.1f,%2.1f), CTS (%2.1f,%2.1f); ",
841  parList+=onePar ;
842  snprintf(onePar,buffersize,"Dist to bad channel: EMC > %2.1f, PHOS > %2.1f; ",fEMCALBadChMinDist,fPHOSBadChMinDist) ;
843  parList+=onePar ;
844  snprintf(onePar,buffersize,"N cells: EMC > %d, PHOS > %d; ",fEMCALNCellsCut,fPHOSNCellsCut) ;
845  parList+=onePar ;
846  snprintf(onePar,buffersize,"EMC time cut single window (%2.2f,%2.2f); ",fEMCALTimeCutMin,fEMCALTimeCutMax) ;
847  parList+=onePar ;
848  snprintf(onePar,buffersize,"Check: calo fid cut %d; ",fCheckFidCut) ;
849  parList+=onePar ;
850  snprintf(onePar,buffersize,"Track: status %d, SPD hit %d; ",(Int_t) fTrackStatus, fSelectSPDHitTracks) ;
851  parList+=onePar ;
852  snprintf(onePar,buffersize,"multip. eta cut %1.1f; npt cuts %d;",fTrackMultEtaCut, fTrackMultNPtCut) ;
853  parList+=onePar ;
854 
855  if(fUseTrackDCACut)
856  {
857  snprintf(onePar,buffersize,"DCA cut ON, param (%2.4f,%2.4f,%2.4f); ",fTrackDCACut[0],fTrackDCACut[1],fTrackDCACut[2]) ;
858  parList+=onePar ;
859  }
860 
861  snprintf(onePar,buffersize,"Recalculate Clusters = %d, E linearity = %d; ",fRecalculateClusters, fCorrectELinearity) ;
862  parList+=onePar ;
863 
864  snprintf(onePar,buffersize,"SE trigger sel. %d, not? trigger Mask? %d, MB Trigger Mask for mixed %d; ",
866  parList+=onePar ;
867 
868  snprintf(onePar,buffersize,"Select fired trigger %s; Remove Bad trigger event %d, unmatched %d; Accept fastcluster %d, Reject LED %d ",
870  parList+=onePar ;
871 
873  {
874  snprintf(onePar,buffersize,"Accept only labels from: ");
875  parList+=onePar ;
876  for(Int_t i = 0; i< fNMCGenerToAccept; i++)
877  parList+=(fMCGenerToAccept[i]+" ") ;
878 
879  snprintf(onePar,buffersize,"; ");
880  parList+=onePar ;
881  }
882 
884  {
885  snprintf(onePar,buffersize,"EMC M02 smear ON, function %d, param %2.4f, %d<=NLM<=%d; ",
887  parList+=onePar ;
888  }
889 
891  {
892  snprintf(onePar,buffersize,"jet pt / pt hard < %2.1f; ",fPtHardAndJetPtFactor);
893  parList+=onePar ;
894  }
895 
897  {
898  snprintf(onePar,buffersize,"cluster pt / pt hard < %2.2f",fPtHardAndClusterPtFactor);
899  parList+=onePar ;
900  }
901 
902  snprintf(onePar,buffersize,"Centrality: Class %s, Option %d, Bin [%d,%d]; New centrality %d; Event plane method %s; ",
904  parList+=onePar ;
905 
906  return new TObjString(parList) ;
907 }
908 
909 
910 //______________________________________________
912 //______________________________________________
914 {
915  if(fMC)
916  {
917  return fMC->Header();
918  }
919  else
920  {
921  AliInfo("Header is not available");
922  return 0x0 ;
923  }
924 }
925 
926 //_________________________________________________________
929 //_________________________________________________________
931 {
932  AliInfo("Input are not AODs");
933 
934  return NULL ;
935 }
936 
937 //________________________________________________________
940 //________________________________________________________
941 AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader() const
942 {
943  AliInfo("Input are not AODs");
944 
945  return NULL ;
946 }
947 
948 //___________________________________________________________
956 //___________________________________________________________
957 Int_t AliCaloTrackReader::GetVertexBC(const AliVVertex * vtx)
958 {
959  Int_t vertexBC=vtx->GetBC();
960 
961  if(!fRecalculateVertexBC) return vertexBC;
962 
963  // Value not available, recalculate it.
964 
965  Double_t bz = fInputEvent->GetMagneticField();
966  Bool_t bc0 = kFALSE;
967  Int_t ntr = GetCTSTracks()->GetEntriesFast();
968  //printf("N Tracks %d\n",ntr);
969 
970  for(Int_t i = 0 ; i < ntr ; i++)
971  {
972  AliVTrack * track = (AliVTrack*) (GetCTSTracks()->At(i));
973 
974  //Check if has TOF info, if not skip
975  ULong_t status = track->GetStatus();
976  Bool_t okTOF = (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ;
977  vertexBC = track->GetTOFBunchCrossing(bz);
978  Float_t pt = track->Pt();
979 
980  if(!okTOF) continue;
981 
982  // Get DCA x, y
983  Double_t dca[2] = {1e6,1e6};
984  Double_t covar[3] = {1e6,1e6,1e6};
985  track->PropagateToDCA(vtx,bz,100.,dca,covar);
986 
987  if(AcceptDCA(pt,dca[0]))
988  {
989  if (vertexBC !=0 && fVertexBC != AliVTrack::kTOFBCNA) return vertexBC;
990  else if(vertexBC == 0) bc0 = kTRUE;
991  }
992  }
993 
994  if( bc0 ) vertexBC = 0 ;
995  else vertexBC = AliVTrack::kTOFBCNA ;
996 
997  return vertexBC;
998 }
999 
1000 //_____________________________
1007 //_____________________________
1009 {
1010  return track->GetID();
1011 }
1012 
1013 //_____________________________
1016 //_____________________________
1018 {
1019  // Activate debug level in AliAnaWeights
1020  if( fWeightUtils->GetDebug() >= 0 )
1021  (AliAnalysisManager::GetAnalysisManager())->AddClassDebug(fWeightUtils->ClassName(), fWeightUtils->GetDebug());
1022 }
1023 
1024 //_______________________________________
1026 //_______________________________________
1028 {
1029  fDataType = kESD ;
1030  fCTSPtMin = 0.1 ;
1031  fEMCALPtMin = 0.1 ;
1032  fPHOSPtMin = 0.1 ;
1033  fCTSPtMax = 1000. ;
1034  fEMCALPtMax = 1000. ;
1035  fPHOSPtMax = 1000. ;
1036 
1037  fEMCALBadChMinDist = 0; // open, 2; // standard
1038  fPHOSBadChMinDist = 0; // open, 2; // standard
1039 
1040  fEMCALNCellsCut = 0; // open, 1; // standard
1041  fPHOSNCellsCut = 0; // open, 2; // standard
1042 
1043  //Track DCA cuts
1044  // dca_xy cut = 0.0105+0.0350/TMath::Power(pt,1.1);
1045  fTrackDCACut[0] = 0.0105;
1046  fTrackDCACut[1] = 0.0350;
1047  fTrackDCACut[2] = 1.1;
1048 
1049  //Do not filter the detectors input by default.
1050  fFillEMCAL = kFALSE;
1051  fFillDCAL = kFALSE;
1052  fFillPHOS = kFALSE;
1053  fFillCTS = kFALSE;
1054  fFillEMCALCells = kFALSE;
1055  fFillPHOSCells = kFALSE;
1056 
1057  fDeltaAODFileName = "deltaAODPartCorr.root";
1059  fEventTriggerMask = AliVEvent::kAny;
1060  fMixEventTriggerMask = AliVEvent::kAnyINT;
1061  fEventTriggerAtSE = kTRUE; // Use only events that pass event selection at SE base class
1062 
1063  fAcceptFastCluster = kTRUE;
1064  fEventType = -1;
1065 
1066  //We want tracks fitted in the detectors:
1067  //fTrackStatus=AliESDtrack::kTPCrefit;
1068  //fTrackStatus|=AliESDtrack::kITSrefit;
1069  fTrackStatus = 0;
1070 
1071  fV0ADC[0] = 0; fV0ADC[1] = 0;
1072  fV0Mul[0] = 0; fV0Mul[1] = 0;
1073 
1074  fZvtxCut = 10.;
1075 
1076  fNMixedEvent = 1;
1077 
1078  fPtHardAndJetPtFactor = 4. ;
1080 
1081  //Centrality
1082  fUseAliCentrality = kFALSE;
1083  fCentralityClass = "V0M";
1084  fCentralityOpt = 100;
1085  fCentralityBin[0] = fCentralityBin[1]=-1;
1086 
1087  fEventPlaneMethod = "V0";
1088 
1089  // Allocate memory (not sure this is the right place)
1090  fCTSTracks = new TObjArray();
1091  fEMCALClusters = new TObjArray();
1092  fDCALClusters = new TObjArray();
1093  fPHOSClusters = new TObjArray();
1094  //fTriggerAnalysis = new AliTriggerAnalysis;
1095  fAODBranchList = new TList ;
1096  fOutputContainer = new TList ;
1097 
1098  fEnergyHistogramNbins = 200;
1099  fEnergyHistogramLimit[0] = 0 ;
1100  fEnergyHistogramLimit[1] = 100;
1101 
1102  fPileUpParamSPD[0] = 3 ; fPileUpParamSPD[1] = 0.8 ;
1103  fPileUpParamSPD[2] = 3.0 ; fPileUpParamSPD[3] = 2.0 ; fPileUpParamSPD[4] = 5.0;
1104 
1105  // Parametrized time cut (LHC11d)
1108 
1109  // Parametrized time cut (LHC11c)
1110  //fEMCALParamTimeCutMin[0] =-5; fEMCALParamTimeCutMin[1] =-1 ; fEMCALParamTimeCutMin[2] = 1.87; fEMCALParamTimeCutMin[3] = 0.4;
1111  //fEMCALParamTimeCutMax[0] = 3.5; fEMCALParamTimeCutMax[1] = 50; fEMCALParamTimeCutMax[2] = 0.15; fEMCALParamTimeCutMax[3] = 1.6;
1112 
1113  fTimeStampRunMin = -1;
1114  fTimeStampRunMax = 1e12;
1117 
1118  for(Int_t i = 0; i < 19; i++)
1119  {
1120  fEMCalBCEvent [i] = 0;
1121  fEMCalBCEventCut[i] = 0;
1122  fTrackBCEvent [i] = 0;
1123  fTrackBCEventCut[i] = 0;
1124  }
1125 
1126  // Trigger match-rejection
1127  fTriggerPatchTimeWindow[0] = 8;
1128  fTriggerPatchTimeWindow[1] = 9;
1129 
1130  fTriggerClusterBC = -10000 ;
1133  fTriggerClusterIndex = -1;
1134  fTriggerClusterId = -1;
1135 
1136  //Jets
1139  if(!fNonStandardJets) fNonStandardJets = new TClonesArray("AliAODJet",100);
1142  if(!fBackgroundJets) fBackgroundJets = new AliAODJetEventBackground();
1143 
1144  fSmearShowerShapeWidth = 0.005;
1145  fSmearNLMMin = 1;
1146  fSmearNLMMax = 1;
1147 
1148  fWeightUtils = new AliAnaWeights() ;
1149  fEventWeight = 1 ;
1150 
1151  fTrackMultNPtCut = 8;
1152  fTrackMultPtCut[0] = 0.15; fTrackMultPtCut[1] = 0.5; fTrackMultPtCut[2] = 1.0;
1153  fTrackMultPtCut[3] = 2.0 ; fTrackMultPtCut[4] = 4.0; fTrackMultPtCut[5] = 6.0;
1154  fTrackMultPtCut[6] = 8.0 ; fTrackMultPtCut[7] = 10.;
1155  fTrackMultPtCut[8] = 15.0; fTrackMultPtCut[9] = 20.;
1156 }
1157 
1158 //__________________________________________________________________________
1161 //__________________________________________________________________________
1163 {
1164  // Parametrized cut depending on E
1165  if(fUseParamTimeCut)
1166  {
1169  //printf("tof %f, minCut %f, maxCut %f\n",tof,minCut,maxCut);
1170  if( tof < minCut || tof > maxCut ) return kFALSE ;
1171  }
1172 
1173  //In any case, the time should to be larger than the fixed window ...
1174  if( tof < fEMCALTimeCutMin || tof > fEMCALTimeCutMax ) return kFALSE ;
1175 
1176  return kTRUE ;
1177 }
1178 
1179 //________________________________________________
1182 //________________________________________________
1184 {
1185  return fInputEvent->IsPileupFromSPD((Int_t) fPileUpParamSPD[0] , fPileUpParamSPD[1] ,
1186  fPileUpParamSPD[2] , fPileUpParamSPD[3] , fPileUpParamSPD[4] );
1187  //printf("Param : %d, %2.2f, %2.2f, %2.2f, %2.2f\n",(Int_t) fPileUpParamSPD[0], fPileUpParamSPD[1], fPileUpParamSPD[2], fPileUpParamSPD[3], fPileUpParamSPD[4]);
1188 }
1189 
1190 //__________________________________________________
1192 //__________________________________________________
1194 {
1195  if(fNPileUpClusters > fNPileUpClustersCut) return kTRUE ;
1196  else return kFALSE;
1197 }
1198 
1199 //________________________________________________________
1201 //________________________________________________________
1203 {
1204  if( IsPileUpFromSPD() && IsPileUpFromEMCal()) return kTRUE ;
1205  else return kFALSE;
1206 }
1207 
1208 //_______________________________________________________
1210 //_______________________________________________________
1212 {
1213  if( IsPileUpFromSPD() || IsPileUpFromEMCal()) return kTRUE ;
1214  else return kFALSE;
1215 }
1216 
1217 //___________________________________________________________
1219 //___________________________________________________________
1221 {
1222  if( IsPileUpFromSPD() && !IsPileUpFromEMCal()) return kTRUE ;
1223  else return kFALSE;
1224 }
1225 
1226 //___________________________________________________________
1228 //___________________________________________________________
1230 {
1231  if( !IsPileUpFromSPD() && IsPileUpFromEMCal()) return kTRUE ;
1232  else return kFALSE;
1233 }
1234 
1235 //______________________________________________________________
1237 //______________________________________________________________
1239 {
1240  if( !IsPileUpFromSPD() && !IsPileUpFromEMCal()) return kTRUE ;
1241  else return kFALSE;
1242 }
1243 
1244 //___________________________________________________________________________________
1252 //___________________________________________________________________________________
1253 Bool_t AliCaloTrackReader::FillInputEvent(Int_t iEntry, const char * /*curFileName*/)
1254 {
1255  fEventNumber = iEntry;
1256  fTriggerClusterIndex = -1;
1257  fTriggerClusterId = -1;
1258  fIsTriggerMatch = kFALSE;
1259  fTriggerClusterBC = -10000;
1260  fIsExoticEvent = kFALSE;
1261  fIsBadCellEvent = kFALSE;
1262  fIsBadMaxCellEvent = kFALSE;
1263 
1264  fIsTriggerMatchOpenCut[0] = kFALSE ;
1265  fIsTriggerMatchOpenCut[1] = kFALSE ;
1266  fIsTriggerMatchOpenCut[2] = kFALSE ;
1267 
1268  //fCurrentFileName = TString(currentFileName);
1269  if(!fInputEvent)
1270  {
1271  AliInfo("Input event not available, skip event analysis");
1272  return kFALSE;
1273  }
1274 
1275  fhNEventsAfterCut->Fill(0.5);
1276 
1277  //-----------------------------------------------
1278  // Select the event depending on the trigger type
1279  // and other event characteristics
1280  // like the goodness of the EMCal trigger
1281  //-----------------------------------------------
1282 
1283  Bool_t accept = CheckEventTriggers();
1284  if(!accept) return kFALSE;
1285 
1286  AliDebug(1,"Pass Event trigger selection");
1287 
1288 
1289  //------------------------------------------------------
1290  // Event rejection depending on time stamp
1291  //------------------------------------------------------
1292 
1294  {
1295  AliESDEvent* esd = dynamic_cast<AliESDEvent*> (fInputEvent);
1296  if(esd)
1297  {
1298  Int_t timeStamp = esd->GetTimeStamp();
1299  Float_t timeStampFrac = 1.*(timeStamp-fTimeStampRunMin) / (fTimeStampRunMax-fTimeStampRunMin);
1300 
1301  //printf("stamp0 %d, max0 %d, frac %f\n", timeStamp-fTimeStampRunMin,fTimeStampRunMax-fTimeStampRunMin, timeStampFrac);
1302 
1303  if(timeStampFrac < fTimeStampEventFracMin || timeStampFrac > fTimeStampEventFracMax) return kFALSE;
1304  }
1305 
1306  AliDebug(1,"Pass Time Stamp rejection");
1307 
1308  fhNEventsAfterCut->Fill(8.5);
1309  }
1310 
1311  //------------------------------------------------------
1312  // Event rejection depending on vertex, pileup, v0and
1313  //------------------------------------------------------
1314 
1315  FillVertexArray();
1316 
1318  {
1319  if( !CheckForPrimaryVertex() ) return kFALSE; // algorithm in ESD/AOD Readers
1320 
1321  fhNEventsAfterCut->Fill(9.5);
1322 
1323  if( TMath::Abs(fVertex[0][0] ) < 1.e-6 &&
1324  TMath::Abs(fVertex[0][1] ) < 1.e-6 &&
1325  TMath::Abs(fVertex[0][2] ) < 1.e-6 ) return kFALSE;
1326 
1327  AliDebug(1,"Pass primary vertex/null rejection");
1328 
1329  fhNEventsAfterCut->Fill(10.5);
1330  }
1331 
1332  //Reject events with Z vertex too large, only for SE analysis, if not, cut on the analysis code
1333  if(!GetMixedEvent() && TMath::Abs(fVertex[0][2]) > fZvtxCut) return kFALSE;
1334 
1335  fhNEventsAfterCut->Fill(11.5);
1336 
1337  AliDebug(1,"Pass z vertex rejection");
1338 
1339 
1340  //printf("Reader : IsPileUp %d, Multi %d\n",IsPileUpFromSPD(),fInputEvent->IsPileupFromSPDInMultBins());
1341 
1343  {
1344  // Do not analyze events with pileup
1345  Bool_t bPileup = IsPileUpFromSPD();
1346  //IsPileupFromSPDInMultBins() // method to try
1347  //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]);
1348  if(bPileup) return kFALSE;
1349 
1350  AliDebug(1,"Pass Pile-Up event rejection");
1351 
1352  fhNEventsAfterCut->Fill(12.5);
1353  }
1354 
1356  {
1357  AliVVZERO* v0 = fInputEvent->GetVZEROData();
1358 
1359  Bool_t bV0AND = ((v0->GetV0ADecision()==1) && (v0->GetV0CDecision()==1));
1360  //bV0AND = fTriggerAnalysis->IsOfflineTriggerFired((AliESDEvent*)fInputEvent, AliTriggerAnalysis::kV0AND);
1361  //printf("V0AND event? %d\n",bV0AND);
1362 
1363  if(!bV0AND)
1364  {
1365  AliDebug(1,"Reject event by V0AND");
1366  return kFALSE;
1367  }
1368 
1369  AliDebug(1,"Pass V0AND event rejection");
1370 
1371  fhNEventsAfterCut->Fill(13.5);
1372  }
1373 
1374  //------------------------------------------------------
1375  // Check if there is a centrality value, PbPb analysis,
1376  // and if a centrality bin selection is requested
1377  //------------------------------------------------------
1378 
1379  //If we need a centrality bin, we select only those events in the corresponding bin.
1380  Int_t cen = -1;
1381  if ( fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
1382  {
1383  cen = GetEventCentrality();
1384 
1385  if(cen > fCentralityBin[1] || cen <= fCentralityBin[0]) return kFALSE; //reject events out of bin.
1386 
1387  AliDebug(1,"Pass centrality rejection");
1388 
1389  fhNEventsAfterCut->Fill(14.5);
1390  }
1391 
1392  //----------------------------------------------------------------
1393  // Reject the event if the event header name is not
1394  // the one requested among the possible generators.
1395  // Needed in case of cocktail MC generation with multiple options.
1396  //----------------------------------------------------------------
1398  {
1399  if(!GetGenEventHeader()) return kFALSE;
1400 
1401  AliDebug(1,"Pass Event header selection");
1402 
1403  fhNEventsAfterCut->Fill(15.5);
1404  }
1405 
1406  //---------------------------------------------------------------------------
1407  // In case of analysis of events with jets, skip those with jet pt > 5 pt hard
1408  // To be used on for MC data in pT hard bins
1409  //---------------------------------------------------------------------------
1410 
1412  {
1413  if(!ComparePtHardAndJetPt()) return kFALSE ;
1414 
1415  AliDebug(1,"Pass Pt Hard - Jet rejection");
1416 
1417  fhNEventsAfterCut->Fill(16.5);
1418  }
1419 
1421  {
1422  if(!ComparePtHardAndClusterPt()) return kFALSE ;
1423 
1424  AliDebug(1,"Pass Pt Hard - Cluster rejection");
1425 
1426  fhNEventsAfterCut->Fill(17.5);
1427  }
1428 
1429  //------------------------------------------------------------------
1430  // Recover the weight assigned to the event, if provided
1431  // right now only for pT-hard bins and centrality depedent weights
1432  //------------------------------------------------------------------
1434  {
1436 
1437  fWeightUtils->SetPythiaEventHeader(((AliGenPythiaEventHeader*)GetGenEventHeader()));
1438 
1440  }
1441 
1442  //-------------------------------------------------------
1443  // Get the main vertex BC, in case not available
1444  // it is calculated in FillCTS checking the BC of tracks
1445  //------------------------------------------------------
1446  fVertexBC = fInputEvent->GetPrimaryVertex()->GetBC();
1447 
1448  //-----------------------------------------------
1449  // Fill the arrays with cluster/tracks/cells data
1450  //-----------------------------------------------
1451 
1452  if(fFillCTS)
1453  {
1454  FillInputCTS();
1455  //Accept events with at least one track
1456  if(fTrackMult[0] == 0 && fDoRejectNoTrackEvents) return kFALSE ;
1457 
1458  AliDebug(1,"Pass rejection of null track events");
1459 
1460  fhNEventsAfterCut->Fill(18.5);
1461  }
1462 
1464  {
1465  if(fVertexBC != 0 && fVertexBC != AliVTrack::kTOFBCNA) return kFALSE ;
1466 
1467  AliDebug(1,"Pass rejection of events with vertex at BC!=0");
1468 
1469  fhNEventsAfterCut->Fill(19.5);
1470  }
1471 
1472  if(fFillEMCALCells)
1474 
1475  if(fFillPHOSCells)
1477 
1478  if(fFillEMCAL || fFillDCAL)
1479  FillInputEMCAL();
1480 
1481  if(fFillPHOS)
1482  FillInputPHOS();
1483 
1484  FillInputVZERO();
1485 
1486  //one specified jet branch
1491 
1492  AliDebug(1,"Event accepted for analysis");
1493 
1494  return kTRUE ;
1495 }
1496 
1497 //__________________________________________________
1500 //__________________________________________________
1502 {
1503  if(fUseAliCentrality)
1504  {
1505  if ( !GetCentrality() ) return -1;
1506 
1507  AliDebug(1,Form("Cent. Percentile: V0M %2.2f, CL0 %2.2f, CL1 %2.2f; selected class %s",
1508  GetCentrality()->GetCentralityPercentile("V0M"),
1509  GetCentrality()->GetCentralityPercentile("CL0"),
1510  GetCentrality()->GetCentralityPercentile("CL1"),
1511  fCentralityClass.Data()));
1512 
1513  if (fCentralityOpt == 100) return (Int_t) GetCentrality()->GetCentralityPercentile(fCentralityClass); // 100 bins max
1514  else if(fCentralityOpt == 10) return GetCentrality()->GetCentralityClass10(fCentralityClass);// 10 bins max
1515  else if(fCentralityOpt == 20) return GetCentrality()->GetCentralityClass5(fCentralityClass); // 20 bins max
1516  else
1517  {
1518  AliInfo(Form("Unknown centrality option %d, use 10, 20 or 100\n",fCentralityOpt));
1519  return -1;
1520  }
1521  }
1522  else
1523  {
1524  if ( !GetMultSelCen() ) return -1;
1525 
1526  AliDebug(1,Form("Mult. Percentile: V0M %2.2f, CL0 %2.2f, CL1 %2.2f; selected class %s",
1527  GetMultSelCen()->GetMultiplicityPercentile("V0M",1),
1528  GetMultSelCen()->GetMultiplicityPercentile("CL0",1),
1529  GetMultSelCen()->GetMultiplicityPercentile("CL1",1),
1530  fCentralityClass.Data()));
1531 
1532  return GetMultSelCen()->GetMultiplicityPercentile(fCentralityClass, kTRUE); // returns centrality only for events used in calibration
1533 
1534  // equivalent to
1535  //GetMultSelCen()->GetMultiplicityPercentile("V0M", kFALSE); // returns centrality for any event
1536  //Int_t qual = GetMultSelCen()->GetEvSelCode(); if (qual ! = 0) cent = qual;
1537  }
1538 }
1539 
1540 //_____________________________________________________
1543 //_____________________________________________________
1545 {
1546  if( !GetEventPlane() ) return -1000;
1547 
1548  Float_t ep = GetEventPlane()->GetEventplane(GetEventPlaneMethod(), GetInputEvent());
1549 
1550  if(GetEventPlaneMethod()=="Q" && (ep < 0 || ep > TMath::Pi()))
1551  {
1552  AliDebug(1,Form("Bad EP for <Q> method : %f\n",ep));
1553  return -1000;
1554  }
1555  else if(GetEventPlaneMethod().Contains("V0") )
1556  {
1557  if((ep > TMath::Pi()/2 || ep < -TMath::Pi()/2))
1558  {
1559  AliDebug(1,Form("Bad EP for <%s> method : %f\n",GetEventPlaneMethod().Data(), ep));
1560  return -1000;
1561  }
1562 
1563  ep+=TMath::Pi()/2; // put same range as for <Q> method
1564  }
1565 
1566  AliDebug(3,Form("Event plane angle %f",ep));
1567 
1568 // if(fDebug > 0 )
1569 // {
1570 // if (ep > TMath::Pi()) printf("AliCaloTrackReader::GetEventPlaneAngle() - Too large angle = %f\n",ep);
1571 // else if(ep < 0 ) printf("AliCaloTrackReader::GetEventPlaneAngle() - Negative angle = %f\n" ,ep);
1572 // }
1573 
1574  return ep;
1575 }
1576 
1577 //__________________________________________________________
1579 //__________________________________________________________
1581 {
1582  vertex[0] = fVertex[0][0];
1583  vertex[1] = fVertex[0][1];
1584  vertex[2] = fVertex[0][2];
1585 }
1586 
1587 //__________________________________________________________________________
1589 //__________________________________________________________________________
1590 void AliCaloTrackReader::GetVertex(Double_t vertex[3], Int_t evtIndex) const
1591 {
1592  vertex[0] = fVertex[evtIndex][0];
1593  vertex[1] = fVertex[evtIndex][1];
1594  vertex[2] = fVertex[evtIndex][2];
1595 }
1596 
1597 //________________________________________
1600 //________________________________________
1602 {
1603  // Delete previous vertex
1604  if(fVertex)
1605  {
1606  for (Int_t i = 0; i < fNMixedEvent; i++)
1607  {
1608  delete [] fVertex[i] ;
1609  }
1610  delete [] fVertex ;
1611  }
1612 
1613  fVertex = new Double_t*[fNMixedEvent] ;
1614  for (Int_t i = 0; i < fNMixedEvent; i++)
1615  {
1616  fVertex[i] = new Double_t[3] ;
1617  fVertex[i][0] = 0.0 ;
1618  fVertex[i][1] = 0.0 ;
1619  fVertex[i][2] = 0.0 ;
1620  }
1621 
1622  if (!fMixedEvent)
1623  { // Single event analysis
1624  if(fDataType!=kMC)
1625  {
1626 
1627  if(fInputEvent->GetPrimaryVertex())
1628  {
1629  fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]);
1630  }
1631  else
1632  {
1633  AliWarning("NULL primary vertex");
1634  fVertex[0][0]=0.; fVertex[0][1]=0.; fVertex[0][2]=0.;
1635  }//Primary vertex pointer do not exist
1636 
1637  } else
1638  {// MC read event
1639  fVertex[0][0]=0.; fVertex[0][1]=0.; fVertex[0][2]=0.;
1640  }
1641 
1642  AliDebug(1,Form("Single Event Vertex : %f,%f,%f\n",fVertex[0][0],fVertex[0][1],fVertex[0][2]));
1643 
1644  } else
1645  { // MultiEvent analysis
1646  for (Int_t iev = 0; iev < fNMixedEvent; iev++)
1647  {
1648  if (fMixedEvent->GetVertexOfEvent(iev))
1649  fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
1650  else
1651  AliWarning("No vertex found");
1652 
1653  AliDebug(1,Form("Multi Event %d Vertex : %f,%f,%f",iev,fVertex[iev][0],fVertex[iev][1],fVertex[iev][2]));
1654  }
1655  }
1656 }
1657 
1658 //_____________________________________
1664 //_____________________________________
1666 {
1667  AliDebug(1,"Begin");
1668 
1669  Double_t pTrack[3] = {0,0,0};
1670 
1671  Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
1672  Int_t nstatus = 0;
1673  Double_t bz = GetInputEvent()->GetMagneticField();
1674 
1675  for(Int_t i = 0; i < 19; i++)
1676  {
1677  fTrackBCEvent [i] = 0;
1678  fTrackBCEventCut[i] = 0;
1679  }
1680 
1681  for(Int_t iptCut = 0; iptCut < fTrackMultNPtCut; iptCut++ )
1682  {
1683  fTrackMult [iptCut] = 0;
1684  fTrackSumPt[iptCut] = 0;
1685  }
1686 
1687  Bool_t bc0 = kFALSE;
1688  if(fRecalculateVertexBC) fVertexBC = AliVTrack::kTOFBCNA;
1689 
1690  for (Int_t itrack = 0; itrack < nTracks; itrack++)
1691  {
1692  AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
1693 
1694  if ( !AcceptParticleMCLabel( TMath::Abs(track->GetLabel()) ) ) continue ;
1695 
1696  fhCTSTrackCutsPt[0]->Fill(track->Pt());
1697 
1698  //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
1699  ULong_t status = track->GetStatus();
1700 
1701  if (fTrackStatus && !((status & fTrackStatus) == fTrackStatus))
1702  continue ;
1703 
1704  fhCTSTrackCutsPt[1]->Fill(track->Pt());
1705 
1706  nstatus++;
1707 
1708  //-------------------------
1709  // Select the tracks depending on cuts of AOD or ESD
1710  if(!SelectTrack(track, pTrack)) continue ;
1711 
1712  fhCTSTrackCutsPt[2]->Fill(track->Pt());
1713 
1714  //-------------------------
1715  // TOF cuts
1716  Bool_t okTOF = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
1717  Double_t tof = -1000;
1718  Int_t trackBC = -1000 ;
1719 
1720  if(fAccessTrackTOF)
1721  {
1722  if(okTOF)
1723  {
1724  trackBC = track->GetTOFBunchCrossing(bz);
1725  SetTrackEventBC(trackBC+9);
1726 
1727  tof = track->GetTOFsignal()*1e-3;
1728 
1729  //After selecting tracks with small DCA, pointing to vertex, set vertex BC depeding on tracks BC
1731  {
1732  if (trackBC != 0 && trackBC != AliVTrack::kTOFBCNA) fVertexBC = trackBC;
1733  else if(trackBC == 0) bc0 = kTRUE;
1734  }
1735 
1736  //In any case, the time should to be larger than the fixed window ...
1737  if( fUseTrackTimeCut && (trackBC !=0 || tof < fTrackTimeCutMin || tof > fTrackTimeCutMax) )
1738  {
1739  //printf("Remove track time %f and bc = %d\n",tof,trackBC);
1740  continue ;
1741  }
1742  //else printf("Accept track time %f and bc = %d\n",tof,trackBC);
1743  }
1744  }
1745 
1746  fhCTSTrackCutsPt[3]->Fill(track->Pt());
1747 
1748  //---------------------
1749  // DCA cuts
1750  //
1751  fMomentum.SetPxPyPzE(pTrack[0],pTrack[1],pTrack[2],0);
1752 
1753  if(fUseTrackDCACut)
1754  {
1755  Float_t dcaTPC =-999;
1756  //In case of AODs, TPC tracks cannot be propagated back to primary vertex,
1757  if( fDataType == kAOD ) dcaTPC = ((AliAODTrack*) track)->DCA();
1758 
1759  //normal way to get the dca, cut on dca_xy
1760  if(dcaTPC==-999)
1761  {
1762  Double_t dca[2] = {1e6,1e6};
1763  Double_t covar[3] = {1e6,1e6,1e6};
1764  Bool_t okDCA = track->PropagateToDCA(fInputEvent->GetPrimaryVertex(),bz,100.,dca,covar);
1765  if( okDCA) okDCA = AcceptDCA(fMomentum.Pt(),dca[0]);
1766  if(!okDCA)
1767  {
1768  //printf("AliCaloTrackReader::FillInputCTS() - Reject track pt %2.2f, dca_xy %2.4f\n",fMomentum.Pt(),dca[0]);
1769  continue ;
1770  }
1771  }
1772  }// DCA cuts
1773 
1774  fhCTSTrackCutsPt[4]->Fill(track->Pt());
1775 
1776  //-------------------------
1777  // Kinematic/acceptance cuts
1778  //
1779  // Count the tracks in eta < 0.9 and different pT cuts
1780  Float_t ptTrack = fMomentum.Pt();
1781  if(TMath::Abs(track->Eta())< fTrackMultEtaCut)
1782  {
1783  for(Int_t iptCut = 0; iptCut < fTrackMultNPtCut; iptCut++ )
1784  {
1785  if(ptTrack > fTrackMultPtCut[iptCut])
1786  {
1787  fTrackMult [iptCut]++;
1788  fTrackSumPt[iptCut]+=ptTrack;
1789  }
1790  }
1791  }
1792 
1793  if(fCTSPtMin > ptTrack || fCTSPtMax < ptTrack) continue ;
1794 
1795  // Check effect of cuts on track BC
1796  if(fAccessTrackTOF && okTOF) SetTrackEventBCcut(trackBC+9);
1797 
1798  if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS)) continue;
1799 
1800  fhCTSTrackCutsPt[5]->Fill(track->Pt());
1801 
1802  // ------------------------------
1803  // Add selected tracks to array
1804  AliDebug(2,Form("Selected tracks pt %3.2f, phi %3.2f deg, eta %3.2f",
1805  fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
1806 
1807  fCTSTracks->Add(track);
1808 
1809  // TODO, check if remove
1810  if (fMixedEvent) track->SetID(itrack);
1811 
1812  }// track loop
1813 
1814  if( fRecalculateVertexBC && (fVertexBC == 0 || fVertexBC == AliVTrack::kTOFBCNA))
1815  {
1816  if( bc0 ) fVertexBC = 0 ;
1817  else fVertexBC = AliVTrack::kTOFBCNA ;
1818  }
1819 
1820  AliDebug(1,Form("CTS entries %d, input tracks %d, pass status %d, multipliticy %d", fCTSTracks->GetEntriesFast(), nTracks, nstatus, fTrackMult[0]));//fCTSTracksNormalInputEntries);
1821 }
1822 
1823 //_______________________________________________________________________________
1841 //_______________________________________________________________________________
1843 {
1844  // Accept clusters with the proper label, only applicable for MC
1845  if ( clus->GetLabel() >= 0 ) // -1 corresponds to noisy MC
1846  {
1847  if ( !AcceptParticleMCLabel(clus->GetLabel()) ) return ;
1848  }
1849 
1850  // TODO, not sure if needed anymore
1851  Int_t vindex = 0 ;
1852  if (fMixedEvent)
1853  vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
1854 
1855  clus->GetMomentum(fMomentum, fVertex[vindex]);
1856 
1857  // No correction/cut applied yet
1858  fhEMCALClusterCutsE[0]->Fill(clus->E());
1859 
1860  //if( (fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10 )
1861  AliDebug(2,Form("Input cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f, nCells %d",
1862  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta(), clus->GetNCells()));
1863 
1864  //---------------------------
1865  // Embedding case
1866  // TO BE REVISED
1868  {
1869  if(clus->GetNLabels()==0 || clus->GetLabel() < 0) return;
1870  //else printf("Embedded cluster, %d, n label %d label %d \n",iclus,clus->GetNLabels(),clus->GetLabel());
1871  }
1872 
1873  //--------------------------------------
1874  // Apply some corrections in the cluster
1875  //
1877  {
1878  //Recalibrate the cluster energy
1880  {
1882 
1883  clus->SetE(energy);
1884  //printf("Recalibrated Energy %f\n",clus->E());
1885 
1888 
1889  } // recalculate E
1890 
1891  //Recalculate distance to bad channels, if new list of bad channels provided
1893 
1894  //Recalculate cluster position
1895  if(GetCaloUtils()->IsRecalculationOfClusterPositionOn())
1896  {
1898  //clus->GetPosition(pos);
1899  //printf("After Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1900  }
1901 
1902  // Recalculate TOF
1903  if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1904  {
1905  Double_t tof = clus->GetTOF();
1906  Float_t frac =-1;
1907  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
1908 
1909  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1910 
1911  //additional L1 phase shift
1912  if(GetCaloUtils()->GetEMCALRecoUtils()->IsL1PhaseInTimeRecalibrationOn())
1913  {
1914  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absIdMax),
1915  fInputEvent->GetBunchCrossNumber(), tof);
1916  }
1917 
1918  clus->SetTOF(tof);
1919 
1920  }// Time recalibration
1921  }
1922 
1923  // Check effect of corrections
1924  fhEMCALClusterCutsE[1]->Fill(clus->E());
1925 
1926  //-----------------------------------------------------------------
1927  // Reject clusters with bad channels, close to borders and exotic
1928  //
1929  Bool_t goodCluster = GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,
1930  GetCaloUtils()->GetEMCALGeometry(),
1931  GetEMCALCells(),fInputEvent->GetBunchCrossNumber());
1932 
1933  if(!goodCluster)
1934  {
1935  //if( (fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10 )
1936  AliDebug(1,Form("Bad cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
1937  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
1938 
1939  return;
1940  }
1941 
1942  // Check effect of bad cluster removal
1943  fhEMCALClusterCutsE[2]->Fill(clus->E());
1944 
1945  //Float_t pos[3];
1946  //clus->GetPosition(pos);
1947  //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1948 
1949  //--------------------------------------
1950  // Correct non linearity or smear energy
1951  //
1953  {
1955 
1956  AliDebug(5,Form("Correct Non Lin: Old E %3.2f, New E %3.2f",
1957  fMomentum.E(),clus->E()));
1958 
1959  // In case of MC analysis, to match resolution/calibration in real data
1960  // Not needed anymore, just leave for MC studies on systematics
1961  if( GetCaloUtils()->GetEMCALRecoUtils()->IsClusterEnergySmeared() )
1962  {
1963  Float_t rdmEnergy = GetCaloUtils()->GetEMCALRecoUtils()->SmearClusterEnergy(clus);
1964 
1965  AliDebug(5,Form("Smear energy: Old E %3.2f, New E %3.2f",clus->E(),rdmEnergy));
1966 
1967  clus->SetE(rdmEnergy);
1968  }
1969  }
1970 
1971  clus->GetMomentum(fMomentum, fVertex[vindex]);
1972  fhEMCALClusterEtaPhi->Fill(fMomentum.Eta(),GetPhi(fMomentum.Phi()));
1973 
1974  // Check effect linearity correction, energy smearing
1975  fhEMCALClusterCutsE[3]->Fill(clus->E());
1976 
1977  // Check the event BC depending on EMCal clustr before final cuts
1978  Double_t tof = clus->GetTOF()*1e9;
1979 
1980  Int_t bc = TMath::Nint(tof/50) + 9;
1981  //printf("tof %2.2f, bc+5=%d\n",tof,bc);
1982 
1983  SetEMCalEventBC(bc);
1984 
1985  //--------------------------------------
1986  // Apply some kinematical/acceptance cuts
1987  //
1988  if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E())
1989  {
1990  AliDebug(2,Form("Cluster E out of range, %2.2f < %2.2f < %2.2f",fEMCALPtMin,clus->E(),fEMCALPtMax));
1991  return ;
1992  }
1993 
1994  // Select cluster fiducial region
1995  //
1996  Bool_t bEMCAL = kFALSE;
1997  Bool_t bDCAL = kFALSE;
1998  if(fCheckFidCut)
1999  {
2000  if(fFillEMCAL && fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) bEMCAL = kTRUE ;
2001  if(fFillDCAL && fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kDCAL )) bDCAL = kTRUE ;
2002  }
2003  else
2004  {
2005  bEMCAL = kTRUE;
2006  }
2007 
2008  //---------------------------------------------------------------------
2009  // Mask all cells in collumns facing ALICE thick material if requested
2010  //
2012  {
2013  Int_t absId = -1;
2014  Int_t iSupMod = -1;
2015  Int_t iphi = -1;
2016  Int_t ieta = -1;
2017  Bool_t shared = kFALSE;
2018  GetCaloUtils()->GetEMCALRecoUtils()->GetMaxEnergyCell(GetCaloUtils()->GetEMCALGeometry(), GetEMCALCells(),clus,absId,iSupMod,ieta,iphi,shared);
2019 
2020  AliDebug(2,Form("Masked collumn: cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2021  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2022 
2023 
2024  if(GetCaloUtils()->MaskFrameCluster(iSupMod, ieta))
2025  {
2026  AliDebug(2,"Mask cluster");
2027  return;
2028  }
2029  }
2030 
2031  // Check effect of energy and fiducial cuts
2032  fhEMCALClusterCutsE[4]->Fill(clus->E());
2033 
2034  //----------------------------------------------------
2035  // Apply N cells cut
2036  //
2037  if(clus->GetNCells() <= fEMCALNCellsCut && fDataType != AliCaloTrackReader::kMC)
2038  {
2039  AliDebug(2,Form("Cluster with n cells %d < %d",clus->GetNCells(), fEMCALNCellsCut));
2040  return ;
2041  }
2042 
2043  // Check effect of n cells cut
2044  fhEMCALClusterCutsE[5]->Fill(clus->E());
2045 
2046  //----------------------------------------------------
2047  // Apply distance to bad channel cut
2048  //
2049  Double_t distBad = clus->GetDistanceToBadChannel() ; //Distance to bad channel
2050 
2051  if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2052 
2053  if(distBad < fEMCALBadChMinDist)
2054  {
2055  AliDebug(2, Form("Cluster close to bad, dist %2.2f < %2.2f",distBad,fEMCALBadChMinDist));
2056  return ;
2057  }
2058 
2059  // Check effect distance to bad channel cut
2060  fhEMCALClusterCutsE[6]->Fill(clus->E());
2061 
2062  //------------------------------------------
2063  // Apply time cut, count EMCal BC before cut
2064  //
2065  SetEMCalEventBCcut(bc);
2066 
2067  // Shift time in case of no calibration with rough factor
2068  Double_t tofShift = tof;
2069  if(tof > 400) tofShift-=615;
2070  fhEMCALClusterTimeE->Fill(clus->E(),tofShift);
2071 
2072  if(!IsInTimeWindow(tof,clus->E()))
2073  {
2074  fNPileUpClusters++ ;
2075  if(fUseEMCALTimeCut)
2076  {
2077  AliDebug(2,Form("Out of time window E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f, time %e",
2078  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta(),tof));
2079 
2080  return ;
2081  }
2082  }
2083  else
2085 
2086  // Check effect of time cut
2087  fhEMCALClusterCutsE[7]->Fill(clus->E());
2088 
2089 
2090  //----------------------------------------------------
2091  // Smear the SS to try to match data and simulations,
2092  // do it only for simulations.
2093  //
2094  if ( fSmearShowerShape && clus->GetNCells() > 2 )
2095  {
2096  Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(clus, GetEMCALCells());
2097 // Int_t nMaxima = clus->GetNExMax(); // For Run2
2098 
2099  if ( nMaxima >= fSmearNLMMin && nMaxima <= fSmearNLMMax )
2100  {
2101  AliDebug(2,Form("Smear shower shape - Original: %2.4f\n", clus->GetM02()));
2103  {
2104  clus->SetM02( clus->GetM02() + fRandom.Landau(0, fSmearShowerShapeWidth) );
2105  }
2107  {
2108  if(iclus%3 == 0 && clus->GetM02() > 0.1) clus->SetM02( clus->GetM02() + fRandom.Landau(0.05, fSmearShowerShapeWidth) ); //fSmearShowerShapeWidth = 0.035
2109  }
2110  else if (fSmearingFunction == kNoSmearing)
2111  {
2112  clus->SetM02( clus->GetM02() );
2113  }
2114  //clus->SetM02( fRandom.Landau(clus->GetM02(), fSmearShowerShapeWidth) );
2115  AliDebug(2,Form("Width %2.4f Smeared : %2.4f\n", fSmearShowerShapeWidth,clus->GetM02()));
2116  }
2117  }
2118 
2119  //--------------------------------------------------------
2120  // Fill the corresponding array with the selected clusters
2121  // Usually just filling EMCal array with upper or lower clusters is enough,
2122  // but maybe we want to do EMCal-DCal correlations.
2123 
2124  //if((fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10)
2125  AliDebug(2,Form("Selected clusters (EMCAL%d, DCAL%d), E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2126  bEMCAL,bDCAL,fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2127 
2128 
2129  if (bEMCAL) fEMCALClusters->Add(clus);
2130  else if(bDCAL ) fDCALClusters ->Add(clus);
2131 
2132  // TODO, not sure if needed anymore
2133  if (fMixedEvent)
2134  clus->SetID(iclus) ;
2135 }
2136 
2137 //_______________________________________
2144 //_______________________________________
2146 {
2147  AliDebug(1,"Begin");
2148 
2149  // First recalibrate cells, time or energy
2150  // if(GetCaloUtils()->IsRecalibrationOn())
2151  // GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCells(GetCaloUtils()->GetEMCALGeometry(),
2152  // GetEMCALCells(),
2153  // fInputEvent->GetBunchCrossNumber());
2154 
2155  fNPileUpClusters = 0; // Init counter
2156  fNNonPileUpClusters = 0; // Init counter
2157  for(Int_t i = 0; i < 19; i++)
2158  {
2159  fEMCalBCEvent [i] = 0;
2160  fEMCalBCEventCut[i] = 0;
2161  }
2162 
2163  //Loop to select clusters in fiducial cut and fill container with aodClusters
2164  if(fEMCALClustersListName=="")
2165  {
2166  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2167  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2168  {
2169  AliVCluster * clus = 0;
2170  if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
2171  {
2172  if (clus->IsEMCAL())
2173  {
2174  FillInputEMCALAlgorithm(clus, iclus);
2175  }//EMCAL cluster
2176  }// cluster exists
2177  }// cluster loop
2178 
2179  //Recalculate track matching
2181 
2182  }//Get the clusters from the input event
2183  else
2184  {
2185  TClonesArray * clusterList = 0x0;
2186 
2187  if (fInputEvent->FindListObject(fEMCALClustersListName))
2188  {
2189  clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
2190  }
2191  else if(fOutputEvent)
2192  {
2193  clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
2194  }
2195 
2196  if(!clusterList)
2197  {
2198  AliWarning(Form("Wrong name of list with clusters? <%s>",fEMCALClustersListName.Data()));
2199  return;
2200  }
2201 
2202  Int_t nclusters = clusterList->GetEntriesFast();
2203  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2204  {
2205  AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
2206  //printf("E %f\n",clus->E());
2207  if (clus) FillInputEMCALAlgorithm(clus, iclus);
2208  else AliWarning("Null cluster in list!");
2209  }// cluster loop
2210 
2211  // Recalculate the pile-up time, in case long time clusters removed during clusterization
2212  //printf("Input event INIT : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
2213 
2214  fNPileUpClusters = 0; // Init counter
2215  fNNonPileUpClusters = 0; // Init counter
2216  for(Int_t i = 0; i < 19; i++)
2217  {
2218  fEMCalBCEvent [i] = 0;
2219  fEMCalBCEventCut[i] = 0;
2220  }
2221 
2222  for (Int_t iclus = 0; iclus < fInputEvent->GetNumberOfCaloClusters(); iclus++)
2223  {
2224  AliVCluster * clus = 0;
2225 
2226  if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
2227  {
2228  if (clus->IsEMCAL())
2229  {
2230 
2231  Float_t frac =-1;
2232  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
2233  Double_t tof = clus->GetTOF();
2234  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2235  //additional L1 phase shift
2236  if(GetCaloUtils()->GetEMCALRecoUtils()->IsL1PhaseInTimeRecalibrationOn())
2237  {
2238  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absIdMax), fInputEvent->GetBunchCrossNumber(), tof);
2239  }
2240 
2241  tof*=1e9;
2242 
2243  //printf("Input event cluster : AbsIdMax %d, E %2.2f, time %2.2f \n", absIdMax,clus->E(),tof);
2244 
2245  //Reject clusters with bad channels, close to borders and exotic;
2246  if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber())) continue;
2247 
2248  Int_t bc = TMath::Nint(tof/50) + 9;
2249  SetEMCalEventBC(bc);
2250 
2251  if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E()) continue ;
2252 
2253  clus->GetMomentum(fMomentum, fVertex[0]);
2254 
2255  if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) return ;
2256 
2257  SetEMCalEventBCcut(bc);
2258 
2259  if(!IsInTimeWindow(tof,clus->E()))
2260  fNPileUpClusters++ ;
2261  else
2263 
2264  }
2265  }
2266  }
2267 
2268  //printf("Input event : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
2269 
2270  // Recalculate track matching, not necessary if already done in the reclusterization task.
2271  // in case it was not done ...
2273 
2274  }
2275 
2276  AliDebug(1,Form("EMCal selected clusters %d",
2277  fEMCALClusters->GetEntriesFast()));
2278  AliDebug(2,Form("\t n pile-up clusters %d, n non pile-up %d",
2280 }
2281 
2282 //_______________________________________
2284 //_______________________________________
2286 {
2287  AliDebug(1,"Begin");
2288 
2289  // Loop to select clusters in fiducial cut and fill container with aodClusters
2290  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2291  TString genName;
2292  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2293  {
2294  AliVCluster * clus = fInputEvent->GetCaloCluster(iclus) ;
2295  if ( !clus ) continue ;
2296 
2297  if ( !clus->IsPHOS() ) continue ;
2298 
2299  if(clus->GetLabel() >=0 ) // -1 corresponds to noisy MC
2300  {
2301  if ( !AcceptParticleMCLabel(clus->GetLabel()) ) continue ;
2302  }
2303 
2304  fhPHOSClusterCutsE[0]->Fill(clus->E());
2305 
2306  // Skip CPV input
2307  if( clus->GetType() == AliVCluster::kPHOSCharged ) continue ;
2308 
2309  fhPHOSClusterCutsE[1]->Fill(clus->E());
2310 
2311  //---------------------------------------------
2312  // Old feature, try to rely on PHOS tender
2313  //
2315  {
2316  // Recalibrate the cluster energy
2318  {
2319  Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
2320  clus->SetE(energy);
2321  }
2322  }
2323 
2324  //----------------------------------------------------------------------------------
2325  // Check if the cluster contains any bad channel and if close to calorimeter borders
2326  //
2327  // Old feature, try to rely on PHOS tender
2328  if( GetCaloUtils()->ClusterContainsBadChannel(kPHOS,clus->GetCellsAbsId(), clus->GetNCells()))
2329  continue;
2330 
2331  if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells()))
2332  continue;
2333 
2334  // TODO, add exotic cut???
2335 
2336  fhPHOSClusterCutsE[2]->Fill(clus->E());
2337 
2338  // TODO Dead code? remove?
2339  Int_t vindex = 0 ;
2340  if (fMixedEvent)
2341  vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
2342 
2343  clus->GetMomentum(fMomentum, fVertex[vindex]);
2344 
2345  //----------------------------------------------------------------------------------
2346  // Remove clusters close to borders
2347  //
2349  continue ;
2350 
2351  fhPHOSClusterCutsE[3]->Fill(clus->E());
2352 
2353  //----------------------------------------------------------------------------------
2354  // Remove clusters with too low energy
2355  //
2356  if (fPHOSPtMin > fMomentum.E() || fPHOSPtMax < fMomentum.E() )
2357  continue ;
2358 
2359  fhPHOSClusterCutsE[4]->Fill(clus->E());
2360 
2361  //----------------------------------------------------
2362  // Apply N cells cut
2363  //
2364  if(clus->GetNCells() <= fPHOSNCellsCut && fDataType != AliCaloTrackReader::kMC) return ;
2365 
2366  // Check effect of n cells cut
2367  fhPHOSClusterCutsE[5]->Fill(clus->E());
2368 
2369  //----------------------------------------------------
2370  // Apply distance to bad channel cut
2371  //
2372  Double_t distBad = clus->GetDistanceToBadChannel() ; //Distance to bad channel
2373 
2374  if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2375 
2376  if(distBad < fPHOSBadChMinDist) return ;
2377 
2378  // Check effect distance to bad channel cut
2379  fhPHOSClusterCutsE[6]->Fill(clus->E());
2380 
2381  // TODO, add time cut
2382 
2383  //----------------------------------------------------------------------------------
2384  // Add selected clusters to array
2385  //
2386  //if(fDebug > 2 && fMomentum.E() > 0.1)
2387  AliDebug(2,Form("Selected clusters E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2388  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2389 
2390  fPHOSClusters->Add(clus);
2391 
2392  // TODO Dead code? remove?
2393  if (fMixedEvent)
2394  clus->SetID(iclus) ;
2395 
2396  } // esd/aod cluster loop
2397 
2398  AliDebug(1,Form("PHOS selected clusters %d",fPHOSClusters->GetEntriesFast())) ;
2399 }
2400 
2401 //____________________________________________
2403 //____________________________________________
2405 {
2406  fEMCALCells = fInputEvent->GetEMCALCells();
2407 }
2408 
2409 //___________________________________________
2411 //___________________________________________
2413 {
2414  fPHOSCells = fInputEvent->GetPHOSCells();
2415 }
2416 
2417 //_______________________________________
2420 //_______________________________________
2422 {
2423  AliVVZERO* v0 = fInputEvent->GetVZEROData();
2424  //printf("Init V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
2425 
2426  if (v0)
2427  {
2428  AliESDVZERO* esdV0 = dynamic_cast<AliESDVZERO*> (v0);
2429  for (Int_t i = 0; i < 32; i++)
2430  {
2431  if(esdV0)
2432  {//Only available in ESDs
2433  fV0ADC[0] += (Int_t)esdV0->GetAdcV0C(i);
2434  fV0ADC[1] += (Int_t)esdV0->GetAdcV0A(i);
2435  }
2436 
2437  fV0Mul[0] += (Int_t)v0->GetMultiplicityV0C(i);
2438  fV0Mul[1] += (Int_t)v0->GetMultiplicityV0A(i);
2439  }
2440 
2441  AliDebug(1,Form("ADC (%d,%d), Multiplicity (%d,%d)",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]));
2442  }
2443  else
2444  {
2445  AliDebug(1,"Cannot retrieve V0 ESD! Run w/ null V0 charges");
2446  }
2447 }
2448 
2449 //_________________________________________________
2453 //_________________________________________________
2455 {
2456  AliDebug(2,"Begin");
2457 
2458  //
2459  //check if branch name is given
2460  if(!fInputNonStandardJetBranchName.Length())
2461  {
2462  fInputEvent->Print();
2463  AliFatal("No non-standard jet branch name specified. Specify among existing ones.");
2464  }
2465 
2466  fNonStandardJets = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fInputNonStandardJetBranchName.Data()));
2467 
2468  if(!fNonStandardJets)
2469  {
2470  //check if jet branch exist; exit if not
2471  fInputEvent->Print();
2472 
2473  AliFatal(Form("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fInputNonStandardJetBranchName.Data()));
2474  }
2475  else
2476  {
2477  AliDebug(1,Form("AOD input jets %d", fNonStandardJets->GetEntriesFast()));
2478  }
2479 }
2480 
2481 //_________________________________________________
2485 //_________________________________________________
2487 {
2488  AliDebug(1,"Begin");
2489  //
2490  //check if branch name is given
2491  if(!fInputBackgroundJetBranchName.Length())
2492  {
2493  fInputEvent->Print();
2494 
2495  AliFatal("No background jet branch name specified. Specify among existing ones.");
2496  }
2497 
2498  fBackgroundJets = (AliAODJetEventBackground*)(fInputEvent->FindListObject(fInputBackgroundJetBranchName.Data()));
2499 
2500  if(!fBackgroundJets)
2501  {
2502  //check if jet branch exist; exit if not
2503  fInputEvent->Print();
2504 
2505  AliFatal(Form("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fInputBackgroundJetBranchName.Data()));
2506  }
2507  else
2508  {
2509  AliDebug(1,"FillInputBackgroundJets");
2510  fBackgroundJets->Print("");
2511  }
2512 }
2513 
2514 //____________________________________________________________________
2520 //____________________________________________________________________
2522 {
2523  // init some variables
2524  Int_t trigtimes[30], globCol, globRow,ntimes, i;
2525  Int_t absId = -1; //[100];
2526  Int_t nPatch = 0;
2527 
2528  TArrayI patches(0);
2529 
2530  // get object pointer
2531  AliVCaloTrigger *caloTrigger = GetInputEvent()->GetCaloTrigger( "EMCAL" );
2532 
2533  if(!caloTrigger)
2534  {
2535  AliError("Trigger patches input (AliVCaloTrigger) not available in data!");
2536  return patches;
2537  }
2538 
2539  //printf("CaloTrigger Entries %d\n",caloTrigger->GetEntries() );
2540 
2541  // class is not empty
2542  if( caloTrigger->GetEntries() > 0 )
2543  {
2544  // must reset before usage, or the class will fail
2545  caloTrigger->Reset();
2546 
2547  // go throuth the trigger channels
2548  while( caloTrigger->Next() )
2549  {
2550  // get position in global 2x2 tower coordinates
2551  caloTrigger->GetPosition( globCol, globRow );
2552 
2553  //L0
2554  if(IsEventEMCALL0())
2555  {
2556  // get dimension of time arrays
2557  caloTrigger->GetNL0Times( ntimes );
2558 
2559  // no L0s in this channel
2560  // presence of the channel in the iterator still does not guarantee that L0 was produced!!
2561  if( ntimes < 1 )
2562  continue;
2563 
2564  // get timing array
2565  caloTrigger->GetL0Times( trigtimes );
2566  //printf("Get L0 patch : n times %d - trigger time window %d - %d\n",ntimes, tmin,tmax);
2567 
2568  // go through the array
2569  for( i = 0; i < ntimes; i++ )
2570  {
2571  // check if in cut - 8,9 shall be accepted in 2011
2572  if( trigtimes[i] >= tmin && trigtimes[i] <= tmax )
2573  {
2574  //printf("Accepted trigger time %d \n",trigtimes[i]);
2575  //if(nTrig > 99) continue;
2576  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(globCol,globRow, absId);
2577  //printf("pass the time cut globCol %d, globRow %d absId %d\n",globCol,globRow, absId);
2578  patches.Set(nPatch+1);
2579  patches.AddAt(absId,nPatch++);
2580  }
2581  } // trigger time array
2582  }//L0
2583  else if(IsEventEMCALL1()) // L1
2584  {
2585  Int_t bit = 0;
2586  caloTrigger->GetTriggerBits(bit);
2587 
2588  Int_t sum = 0;
2589  caloTrigger->GetL1TimeSum(sum);
2590  //fBitEGA-=2;
2591  Bool_t isEGA1 = ((bit >> fBitEGA ) & 0x1) && IsEventEMCALL1Gamma1() ;
2592  Bool_t isEGA2 = ((bit >> (fBitEGA+1)) & 0x1) && IsEventEMCALL1Gamma2() ;
2593  Bool_t isEJE1 = ((bit >> fBitEJE ) & 0x1) && IsEventEMCALL1Jet1 () ;
2594  Bool_t isEJE2 = ((bit >> (fBitEJE+1)) & 0x1) && IsEventEMCALL1Jet2 () ;
2595 
2596  //if((bit>> fBitEGA )&0x1) printf("Trig Bit %d - bit %d - EG1 %d - EG2 %d\n",fBitEGA ,bit,IsEventEMCALL1Gamma1(),IsEventEMCALL1Gamma2());
2597  //if((bit>>(fBitEGA+1))&0x1) printf("Trig Bit %d - bit %d - EG1 %d - EG2 %d\n",fBitEGA+1,bit,IsEventEMCALL1Gamma1(),IsEventEMCALL1Gamma2());
2598 
2599  if(!isEGA1 && !isEJE1 && !isEGA2 && !isEJE2) continue;
2600 
2601  Int_t patchsize = -1;
2602  if (isEGA1 || isEGA2) patchsize = 2;
2603  else if (isEJE1 || isEJE2) patchsize = 16;
2604 
2605  //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",
2606  // bit,sum,patchsize,isEGA1,isEGA2,isEJE1,isEJE2,fBitEGA,fBitEJE,IsEventEMCALL1Gamma(),IsEventEMCALL1Jet());
2607 
2608 
2609  // add 2x2 (EGA) or 16x16 (EJE) patches
2610  for(Int_t irow=0; irow < patchsize; irow++)
2611  {
2612  for(Int_t icol=0; icol < patchsize; icol++)
2613  {
2614  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(globCol+icol,globRow+irow, absId);
2615  //printf("pass the time cut globCol %d, globRow %d absId %d\n",globCol,globRow, absId);
2616  patches.Set(nPatch+1);
2617  patches.AddAt(absId,nPatch++);
2618  }
2619  }
2620 
2621  } // L1
2622 
2623  } // trigger iterator
2624  } // go through triggers
2625 
2626  if(patches.GetSize()<=0) AliInfo(Form("No patch found! for triggers: %s and selected <%s>",
2628  //else printf(">>>>> N patches %d, test %d,first %d, last %d\n",patches.GetSize(), nPatch, patches.At(0), patches.At(patches.GetSize()-1));
2629 
2630  return patches;
2631 }
2632 
2633 //____________________________________________________________
2637 //____________________________________________________________
2639 {
2640  // Init info from previous event
2641  fTriggerClusterIndex = -1;
2642  fTriggerClusterId = -1;
2643  fTriggerClusterBC = -10000;
2644  fIsExoticEvent = kFALSE;
2645  fIsBadCellEvent = kFALSE;
2646  fIsBadMaxCellEvent = kFALSE;
2647 
2648  fIsTriggerMatch = kFALSE;
2649  fIsTriggerMatchOpenCut[0] = kFALSE;
2650  fIsTriggerMatchOpenCut[1] = kFALSE;
2651  fIsTriggerMatchOpenCut[2] = kFALSE;
2652 
2653  // Do only analysis for triggered events
2654  if(!IsEventEMCALL1() && !IsEventEMCALL0())
2655  {
2656  fTriggerClusterBC = 0;
2657  return;
2658  }
2659 
2660  //printf("***** Try to match trigger to cluster %d **** L0 %d, L1 %d\n",fTriggerPatchClusterMatch,IsEventEMCALL0(),IsEventEMCALL1());
2661 
2662  //Recover the list of clusters
2663  TClonesArray * clusterList = 0;
2664  if (fInputEvent->FindListObject(fEMCALClustersListName))
2665  {
2666  clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
2667  }
2668  else if(fOutputEvent)
2669  {
2670  clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
2671  }
2672 
2673  // Get number of clusters and of trigger patches
2674  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2675  if(clusterList)
2676  nclusters = clusterList->GetEntriesFast();
2677 
2678  Int_t nPatch = patches.GetSize();
2679  Float_t exoDiffTime = GetCaloUtils()->GetEMCALRecoUtils()->GetExoticCellDiffTimeCut();
2680 
2681  //Init some variables used in the cluster loop
2682  Float_t tofPatchMax = 100000;
2683  Float_t ePatchMax =-1;
2684 
2685  Float_t tofMax = 100000;
2686  Float_t eMax =-1;
2687 
2688  Int_t clusMax =-1;
2689  Int_t idclusMax =-1;
2690  Bool_t badClMax = kFALSE;
2691  Bool_t badCeMax = kFALSE;
2692  Bool_t exoMax = kFALSE;
2693 // Int_t absIdMaxTrig= -1;
2694  Int_t absIdMaxMax = -1;
2695 
2696  Int_t nOfHighECl = 0 ;
2697 
2698  //
2699  // Check what is the trigger threshold
2700  // set minimu energym of candidate for trigger cluster
2701  //
2703 
2704  Float_t triggerThreshold = fTriggerL1EventThreshold;
2705  if(IsEventEMCALL0()) triggerThreshold = fTriggerL0EventThreshold;
2706  Float_t minE = triggerThreshold / 2.;
2707 
2708  // This method is not really suitable for JET trigger
2709  // but in case, reduce the energy cut since we do not trigger on high energy particle
2710  if(IsEventEMCALL1Jet() || minE < 1) minE = 1;
2711 
2712  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));
2713 
2714  //
2715  // Loop on the clusters, check if there is any that falls into one of the patches
2716  //
2717  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2718  {
2719  AliVCluster * clus = 0;
2720  if(clusterList) clus = (AliVCluster*) clusterList->At(iclus);
2721  else clus = fInputEvent->GetCaloCluster(iclus);
2722 
2723  if ( !clus ) continue ;
2724 
2725  if ( !clus->IsEMCAL() ) continue ;
2726 
2727  //Skip clusters with too low energy to be triggering
2728  if ( clus->E() < minE ) continue ;
2729 
2730  Float_t frac = -1;
2731  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fInputEvent->GetEMCALCells(), clus,frac);
2732 
2733  Bool_t badCluster = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),
2734  clus->GetCellsAbsId(),clus->GetNCells());
2735  UShort_t cellMax[] = {(UShort_t) absIdMax};
2736  Bool_t badCell = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),cellMax,1);
2737 
2738  // if cell is bad, it can happen that time calibration is not available,
2739  // when calculating if it is exotic, this can make it to be exotic by default
2740  // open it temporarily for this cluster
2741  if(badCell)
2742  GetCaloUtils()->GetEMCALRecoUtils()->SetExoticCellDiffTimeCut(10000000);
2743 
2744  Bool_t exotic = GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCluster(clus, fInputEvent->GetEMCALCells());
2745 
2746  // Set back the time cut on exotics
2747  if(badCell)
2748  GetCaloUtils()->GetEMCALRecoUtils()->SetExoticCellDiffTimeCut(exoDiffTime);
2749 
2750  // Energy threshold for exotic Ecross typically at 4 GeV,
2751  // for lower energy, check that there are more than 1 cell in the cluster
2752  if(!exotic && clus->GetNCells() < 2) exotic = kTRUE;
2753 
2754  Float_t energy = clus->E();
2755  Int_t idclus = clus->GetID();
2756 
2757  Double_t tof = clus->GetTOF();
2758  if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn() && fTriggerClusterTimeRecal){
2759  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2760  //additional L1 phase shift
2761  if(GetCaloUtils()->GetEMCALRecoUtils()->IsL1PhaseInTimeRecalibrationOn()){
2762  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absIdMax), fInputEvent->GetBunchCrossNumber(), tof);
2763  }
2764  }
2765  tof *=1.e9;
2766 
2767  //printf("cluster %d, ID %d, E %2.2f, tof %2.2f, AbsId max %d, exotic %d, bad Cluster %d, bad Cell %d\n",
2768  // iclus,idclus, energy,tof,absIdMax, exotic, badCluster,badCell);
2769 
2770  // Find the highest energy cluster, avobe trigger threshold
2771  // in the event in case no match to trigger is found
2772  if( energy > eMax )
2773  {
2774  tofMax = tof;
2775  eMax = energy;
2776  badClMax = badCluster;
2777  badCeMax = badCell;
2778  exoMax = exotic;
2779  clusMax = iclus;
2780  idclusMax = idclus;
2781  absIdMaxMax = absIdMax;
2782  }
2783 
2784  // count the good clusters in the event avobe the trigger threshold
2785  // to check the exotic events
2786  if(!badCluster && !exotic)
2787  nOfHighECl++;
2788 
2789  // Find match to trigger
2790  if(fTriggerPatchClusterMatch && nPatch>0)
2791  {
2792  for(Int_t iabsId =0; iabsId < nPatch; iabsId++)
2793  {
2794  Int_t absIDCell[4];
2795  GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patches.At(iabsId), absIDCell);
2796  //if(tof > 75 ) printf("E %2.2f TOF %2.2f Trigger patch %d, cells : %d, %d, %d, %d\n",
2797  // clus->E(),tof,patches.At(iabsId), absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
2798 
2799  for(Int_t ipatch = 0; ipatch < 4; ipatch++)
2800  {
2801  if(absIdMax == absIDCell[ipatch])
2802  {
2803  //printf("*** Patches : absId %d, E %2.1f, tof %f \n",absIdMax,clus->E(), tof);
2804  if(energy > ePatchMax)
2805  {
2806  tofPatchMax = tof;
2807  ePatchMax = energy;
2808  fIsBadCellEvent = badCluster;
2809  fIsBadMaxCellEvent = badCell;
2810  fIsExoticEvent = exotic;
2811  fTriggerClusterIndex = iclus;
2812  fTriggerClusterId = idclus;
2813  fIsTriggerMatch = kTRUE;
2814 // absIdMaxTrig = absIdMax;
2815  }
2816  }
2817  }// cell patch loop
2818  }// trigger patch loop
2819  } // Do trigger patch matching
2820 
2821  }// Cluster loop
2822 
2823  // If there was no match, assign as trigger
2824  // the highest energy cluster in the event
2825  if(!fIsTriggerMatch)
2826  {
2827  tofPatchMax = tofMax;
2828  ePatchMax = eMax;
2829  fIsBadCellEvent = badClMax;
2830  fIsBadMaxCellEvent = badCeMax;
2831  fIsExoticEvent = exoMax;
2832  fTriggerClusterIndex = clusMax;
2833  fTriggerClusterId = idclusMax;
2834  }
2835 
2836  Double_t tofPatchMaxUS = TMath::Abs(tofPatchMax);
2837 
2838  if (tofPatchMaxUS < 28 ) fTriggerClusterBC = 0 ;
2839  else if(tofPatchMaxUS < 75 ) fTriggerClusterBC = 1 ;
2840  else if(tofPatchMaxUS < 125) fTriggerClusterBC = 2 ;
2841  else if(tofPatchMaxUS < 175) fTriggerClusterBC = 3 ;
2842  else if(tofPatchMaxUS < 225) fTriggerClusterBC = 4 ;
2843  else if(tofPatchMaxUS < 275) fTriggerClusterBC = 5 ;
2844  else
2845  {
2846  //printf("AliCaloTrackReader::MatchTriggerCluster() - Large BC - tof %2.3f - Index %d\n",tofPatchMaxUS,fTriggerClusterIndex);
2848  else
2849  {
2850  fTriggerClusterIndex = -2;
2851  fTriggerClusterId = -2;
2852  }
2853  }
2854 
2855  if(tofPatchMax < 0) fTriggerClusterBC*=-1;
2856 
2857 
2858  // 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",
2859  // fTriggerClusterIndex, fTriggerClusterId,ePatchMax, tofPatchMax,
2860  // fTriggerClusterBC, fIsBadCellEvent,fIsBadMaxCellEvent,fIsExoticEvent, fIsTriggerMatch, nOfHighECl,absIdMaxMax);
2861  //
2862  // 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",
2863  // clusMax, idclusMax, eMax,tofMax, badClMax, badCeMax,exoMax);
2864 
2865  //Redo matching but open cuts
2866  if(!fIsTriggerMatch && fTriggerClusterId >= 0)
2867  {
2868  // Open time patch time
2869  TArrayI patchOpen = GetTriggerPatches(7,10);
2870 
2871  Int_t patchAbsIdOpenTime = -1;
2872  for(Int_t iabsId =0; iabsId < patchOpen.GetSize(); iabsId++)
2873  {
2874  Int_t absIDCell[4];
2875  patchAbsIdOpenTime = patchOpen.At(iabsId);
2876  GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patchAbsIdOpenTime, absIDCell);
2877  //if(tof > 75 ) printf("E %2.2f TOF %2.2f Trigger patch %d, cells : %d, %d, %d, %d\n",
2878  // clus->E(),tof,patches.At(iabsId), absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
2879 
2880  for(Int_t ipatch = 0; ipatch < 4; ipatch++)
2881  {
2882  if(absIdMaxMax == absIDCell[ipatch])
2883  {
2884  fIsTriggerMatchOpenCut[0] = kTRUE;
2885  break;
2886  }
2887  }// cell patch loop
2888  }// trigger patch loop
2889 
2890  // Check neighbour patches
2891  Int_t patchAbsId = -1;
2892  Int_t globalCol = -1;
2893  Int_t globalRow = -1;
2894  GetCaloUtils()->GetEMCALGeometry()->GetFastORIndexFromCellIndex(absIdMaxMax, patchAbsId);
2895  GetCaloUtils()->GetEMCALGeometry()->GetPositionInEMCALFromAbsFastORIndex(patchAbsId,globalCol,globalRow);
2896 
2897  // Check patches with strict time cut
2898  Int_t patchAbsIdNeigh = -1;
2899  for(Int_t icol = globalCol-1; icol <= globalCol+1; icol++)
2900  {
2901  if(icol < 0 || icol > 47) continue;
2902 
2903  for(Int_t irow = globalRow; irow <= globalRow+1; irow++)
2904  {
2905  if(irow < 0 || irow > 63) continue;
2906 
2907  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(icol, irow, patchAbsIdNeigh);
2908 
2909  if ( patchAbsIdNeigh < 0 ) continue;
2910 
2911  for(Int_t iabsId =0; iabsId < patches.GetSize(); iabsId++)
2912  {
2913  if(patchAbsIdNeigh == patches.At(iabsId))
2914  {
2915  fIsTriggerMatchOpenCut[1] = kTRUE;
2916  break;
2917  }
2918  }// trigger patch loop
2919 
2920  }// row
2921  }// col
2922 
2923  // Check patches with open time cut
2924  Int_t patchAbsIdNeighOpenTime = -1;
2925  for(Int_t icol = globalCol-1; icol <= globalCol+1; icol++)
2926  {
2927  if(icol < 0 || icol > 47) continue;
2928 
2929  for(Int_t irow = globalRow; irow <= globalRow+1; irow++)
2930  {
2931  if(irow < 0 || irow > 63) continue;
2932 
2933  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(icol, irow, patchAbsIdNeighOpenTime);
2934 
2935  if ( patchAbsIdNeighOpenTime < 0 ) continue;
2936 
2937  for(Int_t iabsId =0; iabsId < patchOpen.GetSize(); iabsId++)
2938  {
2939  if(patchAbsIdNeighOpenTime == patchOpen.At(iabsId))
2940  {
2941  fIsTriggerMatchOpenCut[2] = kTRUE;
2942  break;
2943  }
2944  }// trigger patch loop
2945 
2946  }// row
2947  }// col
2948 
2949  // printf("No match, new match: Open time %d-%d, open Neigh %d-%d, both open %d-%d\n",fIsTriggerMatchOpenCut[0],patchAbsIdOpenTime,
2950  // fIsTriggerMatchOpenCut[1],patchAbsIdNeigh,
2951  // fIsTriggerMatchOpenCut[2],patchAbsIdNeighOpenTime);
2952 
2953  patchOpen.Reset();
2954 
2955  }// No trigger match found
2956  //printf("Trigger BC %d, Id %d, Index %d\n",fTriggerClusterBC,fTriggerClusterId,fTriggerClusterIndex);
2957 }
2958 
2959 //_________________________________________________________
2963 //_________________________________________________________
2965 {
2967  {
2968  // get object pointer
2969  AliVCaloTrigger *caloTrigger = GetInputEvent()->GetCaloTrigger( "EMCAL" );
2970 
2971  if ( fBitEGA == 6 )
2972  {
2973  if (IsEventEMCALL1Gamma1()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(1);
2974  else if(IsEventEMCALL1Gamma2()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(3);
2975  else if(IsEventEMCALL1Jet1 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(0);
2976  else if(IsEventEMCALL1Jet2 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(2);
2977 
2978  // printf("L1 trigger Threshold Jet1 %f, Gamma1 %f, Jet2 %f, Gamma2 %f\n",
2979  // 0.07874*caloTrigger->GetL1Threshold(0),
2980  // 0.07874*caloTrigger->GetL1Threshold(1),
2981  // 0.07874*caloTrigger->GetL1Threshold(2),
2982  // 0.07874*caloTrigger->GetL1Threshold(3));
2983  }
2984  else
2985  {
2986  // Old AOD data format, in such case, order of thresholds different!!!
2987  if (IsEventEMCALL1Gamma1()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(0);
2988  else if(IsEventEMCALL1Gamma2()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(2);
2989  else if(IsEventEMCALL1Jet1 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(1);
2990  else if(IsEventEMCALL1Jet2 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(3);
2991 
2992  // printf("L1 trigger Threshold Jet1 %f, Gamma1 %f, Jet2 %f, Gamma2 %f\n",
2993  // 0.07874*caloTrigger->GetL1Threshold(1),
2994  // 0.07874*caloTrigger->GetL1Threshold(0),
2995  // 0.07874*caloTrigger->GetL1Threshold(3),
2996  // 0.07874*caloTrigger->GetL1Threshold(2));
2997  }
2998  }
2999 
3000  // Set L0 threshold, if not set by user
3002  {
3003  // Revise for periods > LHC11d
3004  Int_t runNumber = fInputEvent->GetRunNumber();
3005  if (runNumber < 146861) fTriggerL0EventThreshold = 3. ; // LHC11a
3006  else if(runNumber < 154000) fTriggerL0EventThreshold = 4. ; // LHC11b,c
3007  else if(runNumber < 165000) fTriggerL0EventThreshold = 5.5; // LHC11c,d,e
3008  else if(runNumber < 194000) fTriggerL0EventThreshold = 2 ; // LHC12
3009  else if(runNumber < 197400) fTriggerL0EventThreshold = 3 ; // LHC13def
3010  else if(runNumber < 197400) fTriggerL0EventThreshold = 2 ; // LHC13g
3011  else if(runNumber < 244300) fTriggerL0EventThreshold = 5 ; // LHC15 in, phys 1, 5 in phys2
3012  else if(runNumber < 266400) fTriggerL0EventThreshold = 2.5; // LHC16ir
3013  else fTriggerL0EventThreshold = 3.5; // LHC16s
3014  }
3015 }
3016 
3017 //________________________________________________________
3019 //________________________________________________________
3020 void AliCaloTrackReader::Print(const Option_t * opt) const
3021 {
3022  if(! opt)
3023  return;
3024 
3025  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
3026  printf("Task name : %s\n", fTaskName.Data()) ;
3027  printf("Data type : %d\n", fDataType) ;
3028  printf("CTS Min pT : %2.1f GeV/c\n", fCTSPtMin) ;
3029  printf("EMCAL Min pT : %2.1f GeV/c\n", fEMCALPtMin) ;
3030  printf("PHOS Min pT : %2.1f GeV/c\n", fPHOSPtMin) ;
3031  printf("CTS Max pT : %2.1f GeV/c\n", fCTSPtMax) ;
3032  printf("EMCAL Max pT : %2.1f GeV/c\n", fEMCALPtMax) ;
3033  printf("PHOS Max pT : %2.1f GeV/c\n", fPHOSPtMax) ;
3034  printf("EMCAL Bad Dist > %2.1f \n" , fEMCALBadChMinDist) ;
3035  printf("PHOS Bad Dist > %2.1f \n" , fPHOSBadChMinDist) ;
3036  printf("EMCAL N cells > %d \n" , fEMCALNCellsCut) ;
3037  printf("PHOS N cells > %d \n" , fPHOSNCellsCut) ;
3038  printf("EMCAL Time Cut: %3.1f < TOF < %3.1f\n", fEMCALTimeCutMin, fEMCALTimeCutMax);
3039  printf("Use CTS = %d\n", fFillCTS) ;
3040  printf("Use EMCAL = %d\n", fFillEMCAL) ;
3041  printf("Use DCAL = %d\n", fFillDCAL) ;
3042  printf("Use PHOS = %d\n", fFillPHOS) ;
3043  printf("Use EMCAL Cells = %d\n", fFillEMCALCells) ;
3044  printf("Use PHOS Cells = %d\n", fFillPHOSCells) ;
3045  printf("Track status = %d\n", (Int_t) fTrackStatus) ;
3046 
3047  printf("Track Mult Eta Cut = %2.2f\n", fTrackMultEtaCut) ;
3048 
3049  printf("Track Mult Pt Cuts:") ;
3050  for(Int_t icut = 0; icut < fTrackMultNPtCut; icut++) printf(" %2.2f GeV;",fTrackMultPtCut[icut]);
3051  printf(" \n") ;
3052 
3053  printf("Write delta AOD = %d\n", fWriteOutputDeltaAOD) ;
3054  printf("Recalculate Clusters = %d, E linearity = %d\n", fRecalculateClusters, fCorrectELinearity) ;
3055 
3056  printf("Use Triggers selected in SE base class %d; If not what Trigger Mask? %d; MB Trigger Mask for mixed %d \n",
3058 
3060  printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
3061 
3063  printf("Compare cluster pt and pt hard to accept event, factor = %2.2f",fPtHardAndClusterPtFactor);
3064 
3065  printf("Delta AOD File Name = %s\n", fDeltaAODFileName.Data()) ;
3066  printf("Centrality: Class %s, Option %d, Bin [%d,%d] \n", fCentralityClass.Data(),fCentralityOpt,fCentralityBin[0], fCentralityBin[1]) ;
3067 
3068  printf(" \n") ;
3069 }
3070 
3071 //__________________________________________
3077 //__________________________________________
3079 {
3080  // For LHC11a
3081  // Count number of cells with energy larger than 0.1 in SM3, cut on this number
3082  if ( fRemoveLEDEvents == 1 )
3083  {
3084  Int_t ncellsSM3 = 0;
3085  for(Int_t icell = 0; icell < fInputEvent->GetEMCALCells()->GetNumberOfCells(); icell++)
3086  {
3087  Int_t absID = fInputEvent->GetEMCALCells()->GetCellNumber(icell);
3088  Int_t sm = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
3089  if(fInputEvent->GetEMCALCells()->GetAmplitude(icell) > 0.1 && sm==3) ncellsSM3++;
3090  }
3091 
3092  Int_t ncellcut = 21;
3093  if(GetFiredTriggerClasses().Contains("EMC")) ncellcut = 35;
3094 
3095  if(ncellsSM3 >= ncellcut)
3096  {
3097  AliDebug(1,Form("Reject event with ncells in SM3 %d, cut %d, trig %s",
3098  ncellsSM3,ncellcut,GetFiredTriggerClasses().Data()));
3099  return kTRUE;
3100  }
3101  }
3102  // For testing
3103  // Count number of cells with energy larger than 0.1 any SM, cut on this number
3104  else if ( fRemoveLEDEvents > 1 )
3105  {
3106  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};
3107 
3108  for(Int_t icell = 0; icell < fInputEvent->GetEMCALCells()->GetNumberOfCells(); icell++)
3109  {
3110  Int_t absID = fInputEvent->GetEMCALCells()->GetCellNumber(icell);
3111  Int_t sm = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
3112  if(fInputEvent->GetEMCALCells()->GetAmplitude(icell) > 0.1) ncellsSM[sm]++;
3113  }
3114 
3115  Int_t ncellcut = 21;
3116  if(GetFiredTriggerClasses().Contains("EMC")) ncellcut = 35;
3117 
3118  for(Int_t ism = 0; ism < GetCaloUtils()->GetEMCALGeometry()->GetNumberOfSuperModules(); ism++)
3119  {
3120  if(ncellsSM[ism] >= ncellcut)
3121  {
3122  AliDebug(1,Form("Reject event with ncells in SM%d %d, cut %d, trig %s",
3123  ism,ncellsSM[ism],ncellcut,GetFiredTriggerClasses().Data()));
3124 
3125  printf("Reject event with ncells in SM%d %d, cut %d, trig %s\n",
3126  ism,ncellsSM[ism],ncellcut,GetFiredTriggerClasses().Data());
3127  return kTRUE;
3128  }
3129  }
3130  }
3131 
3132  return kFALSE;
3133 }
3134 
3135 //_________________________________________________________
3139 //_________________________________________________________
3141 {
3142  if(label < 0) return ;
3143 
3144  AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fInputEvent) ;
3145  if(!evt) return ;
3146 
3147  TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
3148  if(!arr) return ;
3149 
3150  if(label < arr->GetEntriesFast())
3151  {
3152  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
3153  if(!particle) return ;
3154 
3155  if(label == particle->Label()) return ; // label already OK
3156  //else printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
3157  }
3158  //else printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
3159 
3160  // loop on the particles list and check if there is one with the same label
3161  for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
3162  {
3163  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
3164  if(!particle) continue ;
3165 
3166  if(label == particle->Label())
3167  {
3168  label = ind;
3169  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
3170  return;
3171  }
3172  }
3173 
3174  label = -1;
3175 
3176  //printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label not found set to -1 \n");
3177 }
3178 
3179 //___________________________________
3181 //___________________________________
3183 {
3184  if(fCTSTracks) fCTSTracks -> Clear();
3185  if(fEMCALClusters) fEMCALClusters -> Clear("C");
3186  if(fPHOSClusters) fPHOSClusters -> Clear("C");
3187 
3188  fV0ADC[0] = 0; fV0ADC[1] = 0;
3189  fV0Mul[0] = 0; fV0Mul[1] = 0;
3190 
3191  if(fNonStandardJets) fNonStandardJets -> Clear("C");
3192  fBackgroundJets->Reset();
3193 }
3194 
3195 //___________________________________________
3199 //___________________________________________
3201 {
3202  fEventTrigMinBias = kFALSE;
3203  fEventTrigCentral = kFALSE;
3204  fEventTrigSemiCentral = kFALSE;
3205  fEventTrigEMCALL0 = kFALSE;
3206  fEventTrigEMCALL1Gamma1 = kFALSE;
3207  fEventTrigEMCALL1Gamma2 = kFALSE;
3208  fEventTrigEMCALL1Jet1 = kFALSE;
3209  fEventTrigEMCALL1Jet2 = kFALSE;
3210 
3211  AliDebug(1,Form("Select trigger mask bit %d - Trigger Event %s - Select <%s>",
3213 
3214  if(fEventTriggerMask <=0 )// in case no mask set
3215  {
3216  // EMC triggered event? Which type?
3217  if( GetFiredTriggerClasses().Contains("-B-") || GetFiredTriggerClasses().Contains("-S-") || GetFiredTriggerClasses().Contains("-I-") )
3218  {
3219  if ( GetFiredTriggerClasses().Contains("EGA" ) ||
3220  GetFiredTriggerClasses().Contains("EG1" ) )
3221  {
3222  fEventTrigEMCALL1Gamma1 = kTRUE;
3223  if( GetFiredTriggerClasses().Contains("EG1" ) && !fFiredTriggerClassName.Contains("EG1") ) fEventTrigEMCALL1Gamma1 = kFALSE;
3224  }
3225  else if( GetFiredTriggerClasses().Contains("EG2" ) )
3226  {
3227  fEventTrigEMCALL1Gamma2 = kTRUE;
3228  if( !fFiredTriggerClassName.Contains("EG2") ) fEventTrigEMCALL1Gamma2 = kFALSE;
3229  }
3230  else if( GetFiredTriggerClasses().Contains("EJE" ) ||
3231  GetFiredTriggerClasses().Contains("EJ1" ) )
3232  {
3233  fEventTrigEMCALL1Jet1 = kTRUE;
3234  if( GetFiredTriggerClasses().Contains("EJ1" ) && !fFiredTriggerClassName.Contains("EJ1") )
3235  fEventTrigEMCALL1Jet1 = kFALSE;
3236  }
3237  else if( GetFiredTriggerClasses().Contains("EJ2" ) )
3238  {
3239  fEventTrigEMCALL1Jet2 = kTRUE;
3240  if( !fFiredTriggerClassName.Contains("EJ2") ) fEventTrigEMCALL1Jet2 = kFALSE;
3241  }
3242  else if( GetFiredTriggerClasses().Contains("CEMC") &&
3243  !GetFiredTriggerClasses().Contains("EGA" ) &&
3244  !GetFiredTriggerClasses().Contains("EJE" ) &&
3245  !GetFiredTriggerClasses().Contains("EG1" ) &&
3246  !GetFiredTriggerClasses().Contains("EJ1" ) &&
3247  !GetFiredTriggerClasses().Contains("EG2" ) &&
3248  !GetFiredTriggerClasses().Contains("EJ2" ) )
3249  {
3250  fEventTrigEMCALL0 = kTRUE;
3251  }
3252 
3253  //Min bias event trigger?
3254  if (GetFiredTriggerClasses().Contains("CCENT_R2-B-NOPF-ALLNOTRD"))
3255  {
3256  fEventTrigCentral = kTRUE;
3257  }
3258  else if(GetFiredTriggerClasses().Contains("CSEMI_R1-B-NOPF-ALLNOTRD"))
3259  {
3260  fEventTrigSemiCentral = kTRUE;
3261  }
3262  else if((GetFiredTriggerClasses().Contains("CINT") || GetFiredTriggerClasses().Contains("CPBI2_B1") ) &&
3263  GetFiredTriggerClasses().Contains("-NOPF-ALLNOTRD") )
3264  {
3265  fEventTrigMinBias = kTRUE;
3266  }
3267  }
3268  }
3269  else
3270  {
3271  // EMC L1 Gamma
3272  if ( fEventTriggerMask & AliVEvent::kEMCEGA )
3273  {
3274  //printf("EGA trigger bit\n");
3275  if (GetFiredTriggerClasses().Contains("EG"))
3276  {
3277  if (GetFiredTriggerClasses().Contains("EGA")) fEventTrigEMCALL1Gamma1 = kTRUE;
3278  else
3279  {
3280  if(GetFiredTriggerClasses().Contains("EG1")) fEventTrigEMCALL1Gamma1 = kTRUE;
3281  if(GetFiredTriggerClasses().Contains("EG2")) fEventTrigEMCALL1Gamma2 = kTRUE;
3282  }
3283  }
3284  }
3285  // EMC L1 Jet
3286  else if( fEventTriggerMask & AliVEvent::kEMCEJE )
3287  {
3288  //printf("EGA trigger bit\n");
3289  if (GetFiredTriggerClasses().Contains("EJ"))
3290  {
3291  if (GetFiredTriggerClasses().Contains("EJE")) fEventTrigEMCALL1Jet1 = kTRUE;
3292  else
3293  {
3294  if(GetFiredTriggerClasses().Contains("EJ1")) fEventTrigEMCALL1Jet1 = kTRUE;
3295  if(GetFiredTriggerClasses().Contains("EJ2")) fEventTrigEMCALL1Jet2 = kTRUE;
3296  }
3297  }
3298  }
3299  // EMC L0
3300  else if((fEventTriggerMask & AliVEvent::kEMC7) ||
3301  (fEventTriggerMask & AliVEvent::kEMC1) )
3302  {
3303  //printf("L0 trigger bit\n");
3304  fEventTrigEMCALL0 = kTRUE;
3305  }
3306  // Min Bias Pb-Pb
3307  else if( fEventTriggerMask & AliVEvent::kCentral )
3308  {
3309  //printf("MB semi central trigger bit\n");
3310  fEventTrigSemiCentral = kTRUE;
3311  }
3312  // Min Bias Pb-Pb
3313  else if( fEventTriggerMask & AliVEvent::kSemiCentral )
3314  {
3315  //printf("MB central trigger bit\n");
3316  fEventTrigCentral = kTRUE;
3317  }
3318  // Min Bias pp, PbPb, pPb
3319  else if((fEventTriggerMask & AliVEvent::kMB ) ||
3320  (fEventTriggerMask & AliVEvent::kINT7) ||
3321  (fEventTriggerMask & AliVEvent::kINT8) ||
3322  (fEventTriggerMask & AliVEvent::kAnyINT) )
3323  {
3324  //printf("MB trigger bit\n");
3325  fEventTrigMinBias = kTRUE;
3326  }
3327  }
3328 
3329  AliDebug(1,Form("Event bits: \n \t MB %d, Cen %d, Sem %d, L0 %d, L1G1 %d, L1G2 %d, L1J1 %d, L1J2 %d",
3333 
3334  // L1 trigger bit
3335  if( fBitEGA == 0 && fBitEJE == 0 )
3336  {
3337  // Init the trigger bit once, correct depending on AliESD(AOD)CaloTrigger header version
3338 
3339  // Simpler way to do it ...
3340 // if( fInputEvent->GetRunNumber() < 172000 )
3341 // reader->SetEventTriggerL1Bit(4,5); // current LHC11 data
3342 // else
3343 // reader->SetEventTriggerL1Bit(6,8); // LHC12-13 data
3344 
3345  // Old values
3346  fBitEGA = 4;
3347  fBitEJE = 5;
3348 
3349  TFile* file = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
3350 
3351  const TList *clist = file->GetStreamerInfoCache();
3352 
3353  if(clist)
3354  {
3355  TStreamerInfo *cinfo = (TStreamerInfo*)clist->FindObject("AliESDCaloTrigger");
3356  Int_t verid = 5; // newer ESD header version
3357  if(!cinfo)
3358  {
3359  cinfo = (TStreamerInfo*)clist->FindObject("AliAODCaloTrigger");
3360  verid = 3; // newer AOD header version
3361  }
3362 
3363  if(cinfo)
3364  {
3365  Int_t classversionid = cinfo->GetClassVersion();
3366  //printf("********* Header class version %d *********** \n",classversionid);
3367 
3368  if (classversionid >= verid)
3369  {
3370  fBitEGA = 6;
3371  fBitEJE = 8;
3372  }
3373  } else AliInfo("AliCaloTrackReader()::SetEventTriggerBit() - Streamer info for trigger class not available, bit not changed");
3374  } else AliInfo("AliCaloTrackReader::SetEventTriggerBit() - Streamer list not available!, bit not changed");
3375 
3376  } // set once the EJE, EGA trigger bit
3377 }
3378 
3379 //____________________________________________________________
3382 //____________________________________________________________
3383 void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
3384 {
3385  fInputEvent = input;
3386  fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ;
3387  if (fMixedEvent)
3388  fNMixedEvent = fMixedEvent->GetNumberOfEvents() ;
3389 
3390  //Delete previous vertex
3391  if(fVertex)
3392  {
3393  for (Int_t i = 0; i < fNMixedEvent; i++)
3394  {
3395  delete [] fVertex[i] ;
3396  }
3397  delete [] fVertex ;
3398  }
3399 
3400  fVertex = new Double_t*[fNMixedEvent] ;
3401  for (Int_t i = 0; i < fNMixedEvent; i++)
3402  {
3403  fVertex[i] = new Double_t[3] ;
3404  fVertex[i][0] = 0.0 ;
3405  fVertex[i][1] = 0.0 ;
3406  fVertex[i][2] = 0.0 ;
3407  }
3408 }
3409 
3410 
Bool_t IsPileUpFromSPD() const
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 fComparePtHardAndClusterPt
In MonteCarlo, jet events, reject events with too large cluster energy.
AliAnaWeights * fWeightUtils
Pointer to AliAnaWeights.
AliCalorimeterUtils * GetCaloUtils() const
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.
AliEMCALRecoUtils * GetEMCALRecoUtils() const
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
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:44
Double_t bz
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)
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.
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.
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)
virtual Bool_t ComparePtHardAndClusterPt()
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
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 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.
Bool_t fEventTrigEMCALL1Jet1
Event is L1-Gamma, threshold 1 on its name, it should correspond kEMCEGA.
Bool_t fFillEMCALCells
Use data from EMCAL.
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
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
unsigned int UInt_t
Definition: External.C:33
TObjArray * fEMCALClusters
Temporal array with EMCAL CaloClusters.
Bool_t IsEventEMCALL1Jet() const
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
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)
virtual Bool_t ComparePtHardAndJetPt()
void RecalculateClusterTrackMatching(AliVEvent *event, TObjArray *clusterArray=0x0, AliMCEvent *mc=0x0)
Bool_t IsInFiducialCut(Float_t eta, Float_t phi, Int_t det) const
Bool_t fEventTrigCentral
Event is AliVEvent::kCentral on its name, it should correspond to PbPb.
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.
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.
virtual TString GetEventPlaneMethod() const
virtual Int_t GetTrackID(AliVTrack *track)
TArrayI fAcceptEventsWithBit
Accept events if trigger bit is on.
unsigned long ULong_t
Definition: External.C:38
TH1F * fhEMCALClusterCutsE[8]
! Control histogram on the different EMCal cluster selection cuts, E
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.
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.
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.
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.
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.
Int_t GetNMaskCellColumns() const
AliFiducialCut * fFiducialCut
Acceptance cuts.
Bool_t fRejectEMCalTriggerEventsWith2Tresholds
Reject events EG2 also triggered by EG1 or EJ2 also triggered by EJ1.
TH1F * fhPHOSClusterCutsE[7]
! Control histogram on the different PHOS cluster selection cuts, E
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.
Int_t fSmearNLMMax
Do smearing for clusters with at maximum this value.
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.
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 fIsTriggerMatchOpenCut[3]
Could not match the event to a trigger patch?, retry opening cuts.
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.
virtual Double_t GetEventPlaneAngle() const
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
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.
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.
virtual void FillInputBackgroundJets()
Bool_t fDoV0ANDEventSelection
Select events depending on V0AND.
AliAODJetEventBackground * fBackgroundJets
! Background jets.
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.
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
Bool_t CheckCellFiducialRegion(AliVCluster *cluster, AliVCaloCells *cells) const
Bool_t IsEventEMCALL0() const