AliPhysics  b752f14 (b752f14)
 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 
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(kFALSE),
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",
1862  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
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()) return ;
1989 
1990  // Select cluster fiducial region
1991  //
1992  Bool_t bEMCAL = kFALSE;
1993  Bool_t bDCAL = kFALSE;
1994  if(fCheckFidCut)
1995  {
1996  if(fFillEMCAL && fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) bEMCAL = kTRUE ;
1997  if(fFillDCAL && fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kDCAL )) bDCAL = kTRUE ;
1998  }
1999  else
2000  {
2001  bEMCAL = kTRUE;
2002  }
2003 
2004  //---------------------------------------------------------------------
2005  // Mask all cells in collumns facing ALICE thick material if requested
2006  //
2008  {
2009  Int_t absId = -1;
2010  Int_t iSupMod = -1;
2011  Int_t iphi = -1;
2012  Int_t ieta = -1;
2013  Bool_t shared = kFALSE;
2014  GetCaloUtils()->GetEMCALRecoUtils()->GetMaxEnergyCell(GetCaloUtils()->GetEMCALGeometry(), GetEMCALCells(),clus,absId,iSupMod,ieta,iphi,shared);
2015 
2016  AliDebug(2,Form("Masked collumn: cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2017  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2018 
2019 
2020  if(GetCaloUtils()->MaskFrameCluster(iSupMod, ieta)) return;
2021  }
2022 
2023  // Check effect of energy and fiducial cuts
2024  fhEMCALClusterCutsE[4]->Fill(clus->E());
2025 
2026  //----------------------------------------------------
2027  // Apply N cells cut
2028  //
2029  if(clus->GetNCells() <= fEMCALNCellsCut && fDataType != AliCaloTrackReader::kMC) return ;
2030 
2031  // Check effect of n cells cut
2032  fhEMCALClusterCutsE[5]->Fill(clus->E());
2033 
2034  //----------------------------------------------------
2035  // Apply distance to bad channel cut
2036  //
2037  Double_t distBad = clus->GetDistanceToBadChannel() ; //Distance to bad channel
2038 
2039  if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2040 
2041  if(distBad < fEMCALBadChMinDist) return ;
2042 
2043  // Check effect distance to bad channel cut
2044  fhEMCALClusterCutsE[6]->Fill(clus->E());
2045 
2046  //------------------------------------------
2047  // Apply time cut, count EMCal BC before cut
2048  //
2049  SetEMCalEventBCcut(bc);
2050 
2051  // Shift time in case of no calibration with rough factor
2052  Double_t tofShift = tof;
2053  if(tof > 400) tofShift-=615;
2054  fhEMCALClusterTimeE->Fill(clus->E(),tofShift);
2055 
2056  if(!IsInTimeWindow(tof,clus->E()))
2057  {
2058  fNPileUpClusters++ ;
2059  if(fUseEMCALTimeCut)
2060  {
2061  AliDebug(2,Form("Out of time window E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f, time %e",
2062  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta(),tof));
2063 
2064  return ;
2065  }
2066  }
2067  else
2069 
2070  // Check effect of time cut
2071  fhEMCALClusterCutsE[7]->Fill(clus->E());
2072 
2073 
2074  //----------------------------------------------------
2075  // Smear the SS to try to match data and simulations,
2076  // do it only for simulations.
2077  //
2078  Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(clus, GetEMCALCells());
2079  // Int_t nMaxima = clus->GetNExMax(); // For Run2
2080  if( fSmearShowerShape && clus->GetNCells() > 2 &&
2081  nMaxima >= fSmearNLMMin && nMaxima <= fSmearNLMMax )
2082  {
2083  AliDebug(2,Form("Smear shower shape - Original: %2.4f\n", clus->GetM02()));
2085  {
2086  clus->SetM02( clus->GetM02() + fRandom.Landau(0, fSmearShowerShapeWidth) );
2087  }
2089  {
2090  if(iclus%3 == 0 && clus->GetM02() > 0.1) clus->SetM02( clus->GetM02() + fRandom.Landau(0.05, fSmearShowerShapeWidth) ); //fSmearShowerShapeWidth = 0.035
2091  }
2092  else if (fSmearingFunction == kNoSmearing)
2093  {
2094  clus->SetM02( clus->GetM02() );
2095  }
2096  //clus->SetM02( fRandom.Landau(clus->GetM02(), fSmearShowerShapeWidth) );
2097  AliDebug(2,Form("Width %2.4f Smeared : %2.4f\n", fSmearShowerShapeWidth,clus->GetM02()));
2098  }
2099 
2100  //--------------------------------------------------------
2101  // Fill the corresponding array with the selected clusters
2102  // Usually just filling EMCal array with upper or lower clusters is enough,
2103  // but maybe we want to do EMCal-DCal correlations.
2104 
2105  //if((fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10)
2106  AliDebug(2,Form("Selected clusters (EMCAL%d, DCAL%d), E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2107  bEMCAL,bDCAL,fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2108 
2109 
2110  if (bEMCAL) fEMCALClusters->Add(clus);
2111  else if(bDCAL ) fDCALClusters ->Add(clus);
2112 
2113  // TODO, not sure if needed anymore
2114  if (fMixedEvent)
2115  clus->SetID(iclus) ;
2116 }
2117 
2118 //_______________________________________
2125 //_______________________________________
2127 {
2128  AliDebug(1,"Begin");
2129 
2130  // First recalibrate cells, time or energy
2131  // if(GetCaloUtils()->IsRecalibrationOn())
2132  // GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCells(GetCaloUtils()->GetEMCALGeometry(),
2133  // GetEMCALCells(),
2134  // fInputEvent->GetBunchCrossNumber());
2135 
2136  fNPileUpClusters = 0; // Init counter
2137  fNNonPileUpClusters = 0; // Init counter
2138  for(Int_t i = 0; i < 19; i++)
2139  {
2140  fEMCalBCEvent [i] = 0;
2141  fEMCalBCEventCut[i] = 0;
2142  }
2143 
2144  //Loop to select clusters in fiducial cut and fill container with aodClusters
2145  if(fEMCALClustersListName=="")
2146  {
2147  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2148  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2149  {
2150  AliVCluster * clus = 0;
2151  if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
2152  {
2153  if (clus->IsEMCAL())
2154  {
2155  FillInputEMCALAlgorithm(clus, iclus);
2156  }//EMCAL cluster
2157  }// cluster exists
2158  }// cluster loop
2159 
2160  //Recalculate track matching
2162 
2163  }//Get the clusters from the input event
2164  else
2165  {
2166  TClonesArray * clusterList = 0x0;
2167 
2168  if (fInputEvent->FindListObject(fEMCALClustersListName))
2169  {
2170  clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
2171  }
2172  else if(fOutputEvent)
2173  {
2174  clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
2175  }
2176 
2177  if(!clusterList)
2178  {
2179  AliWarning(Form("Wrong name of list with clusters? <%s>",fEMCALClustersListName.Data()));
2180  return;
2181  }
2182 
2183  Int_t nclusters = clusterList->GetEntriesFast();
2184  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2185  {
2186  AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
2187  //printf("E %f\n",clus->E());
2188  if (clus) FillInputEMCALAlgorithm(clus, iclus);
2189  else AliWarning("Null cluster in list!");
2190  }// cluster loop
2191 
2192  // Recalculate the pile-up time, in case long time clusters removed during clusterization
2193  //printf("Input event INIT : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
2194 
2195  fNPileUpClusters = 0; // Init counter
2196  fNNonPileUpClusters = 0; // Init counter
2197  for(Int_t i = 0; i < 19; i++)
2198  {
2199  fEMCalBCEvent [i] = 0;
2200  fEMCalBCEventCut[i] = 0;
2201  }
2202 
2203  for (Int_t iclus = 0; iclus < fInputEvent->GetNumberOfCaloClusters(); iclus++)
2204  {
2205  AliVCluster * clus = 0;
2206 
2207  if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
2208  {
2209  if (clus->IsEMCAL())
2210  {
2211 
2212  Float_t frac =-1;
2213  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
2214  Double_t tof = clus->GetTOF();
2215  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2216  //additional L1 phase shift
2217  if(GetCaloUtils()->GetEMCALRecoUtils()->IsL1PhaseInTimeRecalibrationOn())
2218  {
2219  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absIdMax), fInputEvent->GetBunchCrossNumber(), tof);
2220  }
2221 
2222  tof*=1e9;
2223 
2224  //printf("Input event cluster : AbsIdMax %d, E %2.2f, time %2.2f \n", absIdMax,clus->E(),tof);
2225 
2226  //Reject clusters with bad channels, close to borders and exotic;
2227  if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber())) continue;
2228 
2229  Int_t bc = TMath::Nint(tof/50) + 9;
2230  SetEMCalEventBC(bc);
2231 
2232  if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E()) continue ;
2233 
2234  clus->GetMomentum(fMomentum, fVertex[0]);
2235 
2236  if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) return ;
2237 
2238  SetEMCalEventBCcut(bc);
2239 
2240  if(!IsInTimeWindow(tof,clus->E()))
2241  fNPileUpClusters++ ;
2242  else
2244 
2245  }
2246  }
2247  }
2248 
2249  //printf("Input event : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
2250 
2251  // Recalculate track matching, not necessary if already done in the reclusterization task.
2252  // in case it was not done ...
2254 
2255  }
2256 
2257  AliDebug(1,Form("EMCal selected clusters %d",
2258  fEMCALClusters->GetEntriesFast()));
2259  AliDebug(2,Form("\t n pile-up clusters %d, n non pile-up %d",
2261 }
2262 
2263 //_______________________________________
2265 //_______________________________________
2267 {
2268  AliDebug(1,"Begin");
2269 
2270  // Loop to select clusters in fiducial cut and fill container with aodClusters
2271  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2272  TString genName;
2273  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2274  {
2275  AliVCluster * clus = fInputEvent->GetCaloCluster(iclus) ;
2276  if ( !clus ) continue ;
2277 
2278  if ( !clus->IsPHOS() ) continue ;
2279 
2280  if(clus->GetLabel() >=0 ) // -1 corresponds to noisy MC
2281  {
2282  if ( !AcceptParticleMCLabel(clus->GetLabel()) ) continue ;
2283  }
2284 
2285  fhPHOSClusterCutsE[0]->Fill(clus->E());
2286 
2287  // Skip CPV input
2288  if( clus->GetType() == AliVCluster::kPHOSCharged ) continue ;
2289 
2290  fhPHOSClusterCutsE[1]->Fill(clus->E());
2291 
2292  //---------------------------------------------
2293  // Old feature, try to rely on PHOS tender
2294  //
2296  {
2297  // Recalibrate the cluster energy
2299  {
2300  Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
2301  clus->SetE(energy);
2302  }
2303  }
2304 
2305  //----------------------------------------------------------------------------------
2306  // Check if the cluster contains any bad channel and if close to calorimeter borders
2307  //
2308  // Old feature, try to rely on PHOS tender
2309  if( GetCaloUtils()->ClusterContainsBadChannel(kPHOS,clus->GetCellsAbsId(), clus->GetNCells()))
2310  continue;
2311 
2312  if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells()))
2313  continue;
2314 
2315  // TODO, add exotic cut???
2316 
2317  fhPHOSClusterCutsE[2]->Fill(clus->E());
2318 
2319  // TODO Dead code? remove?
2320  Int_t vindex = 0 ;
2321  if (fMixedEvent)
2322  vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
2323 
2324  clus->GetMomentum(fMomentum, fVertex[vindex]);
2325 
2326  //----------------------------------------------------------------------------------
2327  // Remove clusters close to borders
2328  //
2330  continue ;
2331 
2332  fhPHOSClusterCutsE[3]->Fill(clus->E());
2333 
2334  //----------------------------------------------------------------------------------
2335  // Remove clusters with too low energy
2336  //
2337  if (fPHOSPtMin > fMomentum.E() || fPHOSPtMax < fMomentum.E() )
2338  continue ;
2339 
2340  fhPHOSClusterCutsE[4]->Fill(clus->E());
2341 
2342  //----------------------------------------------------
2343  // Apply N cells cut
2344  //
2345  if(clus->GetNCells() <= fPHOSNCellsCut && fDataType != AliCaloTrackReader::kMC) return ;
2346 
2347  // Check effect of n cells cut
2348  fhPHOSClusterCutsE[5]->Fill(clus->E());
2349 
2350  //----------------------------------------------------
2351  // Apply distance to bad channel cut
2352  //
2353  Double_t distBad = clus->GetDistanceToBadChannel() ; //Distance to bad channel
2354 
2355  if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2356 
2357  if(distBad < fPHOSBadChMinDist) return ;
2358 
2359  // Check effect distance to bad channel cut
2360  fhPHOSClusterCutsE[6]->Fill(clus->E());
2361 
2362  // TODO, add time cut
2363 
2364  //----------------------------------------------------------------------------------
2365  // Add selected clusters to array
2366  //
2367  //if(fDebug > 2 && fMomentum.E() > 0.1)
2368  AliDebug(2,Form("Selected clusters E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2369  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2370 
2371  fPHOSClusters->Add(clus);
2372 
2373  // TODO Dead code? remove?
2374  if (fMixedEvent)
2375  clus->SetID(iclus) ;
2376 
2377  } // esd/aod cluster loop
2378 
2379  AliDebug(1,Form("PHOS selected clusters %d",fPHOSClusters->GetEntriesFast())) ;
2380 }
2381 
2382 //____________________________________________
2384 //____________________________________________
2386 {
2387  fEMCALCells = fInputEvent->GetEMCALCells();
2388 }
2389 
2390 //___________________________________________
2392 //___________________________________________
2394 {
2395  fPHOSCells = fInputEvent->GetPHOSCells();
2396 }
2397 
2398 //_______________________________________
2401 //_______________________________________
2403 {
2404  AliVVZERO* v0 = fInputEvent->GetVZEROData();
2405  //printf("Init V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
2406 
2407  if (v0)
2408  {
2409  AliESDVZERO* esdV0 = dynamic_cast<AliESDVZERO*> (v0);
2410  for (Int_t i = 0; i < 32; i++)
2411  {
2412  if(esdV0)
2413  {//Only available in ESDs
2414  fV0ADC[0] += (Int_t)esdV0->GetAdcV0C(i);
2415  fV0ADC[1] += (Int_t)esdV0->GetAdcV0A(i);
2416  }
2417 
2418  fV0Mul[0] += (Int_t)v0->GetMultiplicityV0C(i);
2419  fV0Mul[1] += (Int_t)v0->GetMultiplicityV0A(i);
2420  }
2421 
2422  AliDebug(1,Form("ADC (%d,%d), Multiplicity (%d,%d)",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]));
2423  }
2424  else
2425  {
2426  AliDebug(1,"Cannot retrieve V0 ESD! Run w/ null V0 charges");
2427  }
2428 }
2429 
2430 //_________________________________________________
2434 //_________________________________________________
2436 {
2437  AliDebug(2,"Begin");
2438 
2439  //
2440  //check if branch name is given
2441  if(!fInputNonStandardJetBranchName.Length())
2442  {
2443  fInputEvent->Print();
2444  AliFatal("No non-standard jet branch name specified. Specify among existing ones.");
2445  }
2446 
2447  fNonStandardJets = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fInputNonStandardJetBranchName.Data()));
2448 
2449  if(!fNonStandardJets)
2450  {
2451  //check if jet branch exist; exit if not
2452  fInputEvent->Print();
2453 
2454  AliFatal(Form("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fInputNonStandardJetBranchName.Data()));
2455  }
2456  else
2457  {
2458  AliDebug(1,Form("AOD input jets %d", fNonStandardJets->GetEntriesFast()));
2459  }
2460 }
2461 
2462 //_________________________________________________
2466 //_________________________________________________
2468 {
2469  AliDebug(1,"Begin");
2470  //
2471  //check if branch name is given
2472  if(!fInputBackgroundJetBranchName.Length())
2473  {
2474  fInputEvent->Print();
2475 
2476  AliFatal("No background jet branch name specified. Specify among existing ones.");
2477  }
2478 
2479  fBackgroundJets = (AliAODJetEventBackground*)(fInputEvent->FindListObject(fInputBackgroundJetBranchName.Data()));
2480 
2481  if(!fBackgroundJets)
2482  {
2483  //check if jet branch exist; exit if not
2484  fInputEvent->Print();
2485 
2486  AliFatal(Form("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fInputBackgroundJetBranchName.Data()));
2487  }
2488  else
2489  {
2490  AliDebug(1,"FillInputBackgroundJets");
2491  fBackgroundJets->Print("");
2492  }
2493 }
2494 
2495 //____________________________________________________________________
2501 //____________________________________________________________________
2503 {
2504  // init some variables
2505  Int_t trigtimes[30], globCol, globRow,ntimes, i;
2506  Int_t absId = -1; //[100];
2507  Int_t nPatch = 0;
2508 
2509  TArrayI patches(0);
2510 
2511  // get object pointer
2512  AliVCaloTrigger *caloTrigger = GetInputEvent()->GetCaloTrigger( "EMCAL" );
2513 
2514  if(!caloTrigger)
2515  {
2516  AliError("Trigger patches input (AliVCaloTrigger) not available in data!");
2517  return patches;
2518  }
2519 
2520  //printf("CaloTrigger Entries %d\n",caloTrigger->GetEntries() );
2521 
2522  // class is not empty
2523  if( caloTrigger->GetEntries() > 0 )
2524  {
2525  // must reset before usage, or the class will fail
2526  caloTrigger->Reset();
2527 
2528  // go throuth the trigger channels
2529  while( caloTrigger->Next() )
2530  {
2531  // get position in global 2x2 tower coordinates
2532  caloTrigger->GetPosition( globCol, globRow );
2533 
2534  //L0
2535  if(IsEventEMCALL0())
2536  {
2537  // get dimension of time arrays
2538  caloTrigger->GetNL0Times( ntimes );
2539 
2540  // no L0s in this channel
2541  // presence of the channel in the iterator still does not guarantee that L0 was produced!!
2542  if( ntimes < 1 )
2543  continue;
2544 
2545  // get timing array
2546  caloTrigger->GetL0Times( trigtimes );
2547  //printf("Get L0 patch : n times %d - trigger time window %d - %d\n",ntimes, tmin,tmax);
2548 
2549  // go through the array
2550  for( i = 0; i < ntimes; i++ )
2551  {
2552  // check if in cut - 8,9 shall be accepted in 2011
2553  if( trigtimes[i] >= tmin && trigtimes[i] <= tmax )
2554  {
2555  //printf("Accepted trigger time %d \n",trigtimes[i]);
2556  //if(nTrig > 99) continue;
2557  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(globCol,globRow, absId);
2558  //printf("pass the time cut globCol %d, globRow %d absId %d\n",globCol,globRow, absId);
2559  patches.Set(nPatch+1);
2560  patches.AddAt(absId,nPatch++);
2561  }
2562  } // trigger time array
2563  }//L0
2564  else if(IsEventEMCALL1()) // L1
2565  {
2566  Int_t bit = 0;
2567  caloTrigger->GetTriggerBits(bit);
2568 
2569  Int_t sum = 0;
2570  caloTrigger->GetL1TimeSum(sum);
2571  //fBitEGA-=2;
2572  Bool_t isEGA1 = ((bit >> fBitEGA ) & 0x1) && IsEventEMCALL1Gamma1() ;
2573  Bool_t isEGA2 = ((bit >> (fBitEGA+1)) & 0x1) && IsEventEMCALL1Gamma2() ;
2574  Bool_t isEJE1 = ((bit >> fBitEJE ) & 0x1) && IsEventEMCALL1Jet1 () ;
2575  Bool_t isEJE2 = ((bit >> (fBitEJE+1)) & 0x1) && IsEventEMCALL1Jet2 () ;
2576 
2577  //if((bit>> fBitEGA )&0x1) printf("Trig Bit %d - bit %d - EG1 %d - EG2 %d\n",fBitEGA ,bit,IsEventEMCALL1Gamma1(),IsEventEMCALL1Gamma2());
2578  //if((bit>>(fBitEGA+1))&0x1) printf("Trig Bit %d - bit %d - EG1 %d - EG2 %d\n",fBitEGA+1,bit,IsEventEMCALL1Gamma1(),IsEventEMCALL1Gamma2());
2579 
2580  if(!isEGA1 && !isEJE1 && !isEGA2 && !isEJE2) continue;
2581 
2582  Int_t patchsize = -1;
2583  if (isEGA1 || isEGA2) patchsize = 2;
2584  else if (isEJE1 || isEJE2) patchsize = 16;
2585 
2586  //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",
2587  // bit,sum,patchsize,isEGA1,isEGA2,isEJE1,isEJE2,fBitEGA,fBitEJE,IsEventEMCALL1Gamma(),IsEventEMCALL1Jet());
2588 
2589 
2590  // add 2x2 (EGA) or 16x16 (EJE) patches
2591  for(Int_t irow=0; irow < patchsize; irow++)
2592  {
2593  for(Int_t icol=0; icol < patchsize; icol++)
2594  {
2595  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(globCol+icol,globRow+irow, absId);
2596  //printf("pass the time cut globCol %d, globRow %d absId %d\n",globCol,globRow, absId);
2597  patches.Set(nPatch+1);
2598  patches.AddAt(absId,nPatch++);
2599  }
2600  }
2601 
2602  } // L1
2603 
2604  } // trigger iterator
2605  } // go through triggers
2606 
2607  if(patches.GetSize()<=0) AliInfo(Form("No patch found! for triggers: %s and selected <%s>",
2609  //else printf(">>>>> N patches %d, test %d,first %d, last %d\n",patches.GetSize(), nPatch, patches.At(0), patches.At(patches.GetSize()-1));
2610 
2611  return patches;
2612 }
2613 
2614 //____________________________________________________________
2618 //____________________________________________________________
2620 {
2621  // Init info from previous event
2622  fTriggerClusterIndex = -1;
2623  fTriggerClusterId = -1;
2624  fTriggerClusterBC = -10000;
2625  fIsExoticEvent = kFALSE;
2626  fIsBadCellEvent = kFALSE;
2627  fIsBadMaxCellEvent = kFALSE;
2628 
2629  fIsTriggerMatch = kFALSE;
2630  fIsTriggerMatchOpenCut[0] = kFALSE;
2631  fIsTriggerMatchOpenCut[1] = kFALSE;
2632  fIsTriggerMatchOpenCut[2] = kFALSE;
2633 
2634  // Do only analysis for triggered events
2635  if(!IsEventEMCALL1() && !IsEventEMCALL0())
2636  {
2637  fTriggerClusterBC = 0;
2638  return;
2639  }
2640 
2641  //printf("***** Try to match trigger to cluster %d **** L0 %d, L1 %d\n",fTriggerPatchClusterMatch,IsEventEMCALL0(),IsEventEMCALL1());
2642 
2643  //Recover the list of clusters
2644  TClonesArray * clusterList = 0;
2645  if (fInputEvent->FindListObject(fEMCALClustersListName))
2646  {
2647  clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
2648  }
2649  else if(fOutputEvent)
2650  {
2651  clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
2652  }
2653 
2654  // Get number of clusters and of trigger patches
2655  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2656  if(clusterList)
2657  nclusters = clusterList->GetEntriesFast();
2658 
2659  Int_t nPatch = patches.GetSize();
2660  Float_t exoDiffTime = GetCaloUtils()->GetEMCALRecoUtils()->GetExoticCellDiffTimeCut();
2661 
2662  //Init some variables used in the cluster loop
2663  Float_t tofPatchMax = 100000;
2664  Float_t ePatchMax =-1;
2665 
2666  Float_t tofMax = 100000;
2667  Float_t eMax =-1;
2668 
2669  Int_t clusMax =-1;
2670  Int_t idclusMax =-1;
2671  Bool_t badClMax = kFALSE;
2672  Bool_t badCeMax = kFALSE;
2673  Bool_t exoMax = kFALSE;
2674 // Int_t absIdMaxTrig= -1;
2675  Int_t absIdMaxMax = -1;
2676 
2677  Int_t nOfHighECl = 0 ;
2678 
2679  //
2680  // Check what is the trigger threshold
2681  // set minimu energym of candidate for trigger cluster
2682  //
2684 
2685  Float_t triggerThreshold = fTriggerL1EventThreshold;
2686  if(IsEventEMCALL0()) triggerThreshold = fTriggerL0EventThreshold;
2687  Float_t minE = triggerThreshold / 2.;
2688 
2689  // This method is not really suitable for JET trigger
2690  // but in case, reduce the energy cut since we do not trigger on high energy particle
2691  if(IsEventEMCALL1Jet() || minE < 1) minE = 1;
2692 
2693  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));
2694 
2695  //
2696  // Loop on the clusters, check if there is any that falls into one of the patches
2697  //
2698  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2699  {
2700  AliVCluster * clus = 0;
2701  if(clusterList) clus = (AliVCluster*) clusterList->At(iclus);
2702  else clus = fInputEvent->GetCaloCluster(iclus);
2703 
2704  if ( !clus ) continue ;
2705 
2706  if ( !clus->IsEMCAL() ) continue ;
2707 
2708  //Skip clusters with too low energy to be triggering
2709  if ( clus->E() < minE ) continue ;
2710 
2711  Float_t frac = -1;
2712  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fInputEvent->GetEMCALCells(), clus,frac);
2713 
2714  Bool_t badCluster = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),
2715  clus->GetCellsAbsId(),clus->GetNCells());
2716  UShort_t cellMax[] = {(UShort_t) absIdMax};
2717  Bool_t badCell = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),cellMax,1);
2718 
2719  // if cell is bad, it can happen that time calibration is not available,
2720  // when calculating if it is exotic, this can make it to be exotic by default
2721  // open it temporarily for this cluster
2722  if(badCell)
2723  GetCaloUtils()->GetEMCALRecoUtils()->SetExoticCellDiffTimeCut(10000000);
2724 
2725  Bool_t exotic = GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCluster(clus, fInputEvent->GetEMCALCells());
2726 
2727  // Set back the time cut on exotics
2728  if(badCell)
2729  GetCaloUtils()->GetEMCALRecoUtils()->SetExoticCellDiffTimeCut(exoDiffTime);
2730 
2731  // Energy threshold for exotic Ecross typically at 4 GeV,
2732  // for lower energy, check that there are more than 1 cell in the cluster
2733  if(!exotic && clus->GetNCells() < 2) exotic = kTRUE;
2734 
2735  Float_t energy = clus->E();
2736  Int_t idclus = clus->GetID();
2737 
2738  Double_t tof = clus->GetTOF();
2739  if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn() && fTriggerClusterTimeRecal){
2740  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2741  //additional L1 phase shift
2742  if(GetCaloUtils()->GetEMCALRecoUtils()->IsL1PhaseInTimeRecalibrationOn()){
2743  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absIdMax), fInputEvent->GetBunchCrossNumber(), tof);
2744  }
2745  }
2746  tof *=1.e9;
2747 
2748  //printf("cluster %d, ID %d, E %2.2f, tof %2.2f, AbsId max %d, exotic %d, bad Cluster %d, bad Cell %d\n",
2749  // iclus,idclus, energy,tof,absIdMax, exotic, badCluster,badCell);
2750 
2751  // Find the highest energy cluster, avobe trigger threshold
2752  // in the event in case no match to trigger is found
2753  if( energy > eMax )
2754  {
2755  tofMax = tof;
2756  eMax = energy;
2757  badClMax = badCluster;
2758  badCeMax = badCell;
2759  exoMax = exotic;
2760  clusMax = iclus;
2761  idclusMax = idclus;
2762  absIdMaxMax = absIdMax;
2763  }
2764 
2765  // count the good clusters in the event avobe the trigger threshold
2766  // to check the exotic events
2767  if(!badCluster && !exotic)
2768  nOfHighECl++;
2769 
2770  // Find match to trigger
2771  if(fTriggerPatchClusterMatch && nPatch>0)
2772  {
2773  for(Int_t iabsId =0; iabsId < nPatch; iabsId++)
2774  {
2775  Int_t absIDCell[4];
2776  GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patches.At(iabsId), absIDCell);
2777  //if(tof > 75 ) printf("E %2.2f TOF %2.2f Trigger patch %d, cells : %d, %d, %d, %d\n",
2778  // clus->E(),tof,patches.At(iabsId), absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
2779 
2780  for(Int_t ipatch = 0; ipatch < 4; ipatch++)
2781  {
2782  if(absIdMax == absIDCell[ipatch])
2783  {
2784  //printf("*** Patches : absId %d, E %2.1f, tof %f \n",absIdMax,clus->E(), tof);
2785  if(energy > ePatchMax)
2786  {
2787  tofPatchMax = tof;
2788  ePatchMax = energy;
2789  fIsBadCellEvent = badCluster;
2790  fIsBadMaxCellEvent = badCell;
2791  fIsExoticEvent = exotic;
2792  fTriggerClusterIndex = iclus;
2793  fTriggerClusterId = idclus;
2794  fIsTriggerMatch = kTRUE;
2795 // absIdMaxTrig = absIdMax;
2796  }
2797  }
2798  }// cell patch loop
2799  }// trigger patch loop
2800  } // Do trigger patch matching
2801 
2802  }// Cluster loop
2803 
2804  // If there was no match, assign as trigger
2805  // the highest energy cluster in the event
2806  if(!fIsTriggerMatch)
2807  {
2808  tofPatchMax = tofMax;
2809  ePatchMax = eMax;
2810  fIsBadCellEvent = badClMax;
2811  fIsBadMaxCellEvent = badCeMax;
2812  fIsExoticEvent = exoMax;
2813  fTriggerClusterIndex = clusMax;
2814  fTriggerClusterId = idclusMax;
2815  }
2816 
2817  Double_t tofPatchMaxUS = TMath::Abs(tofPatchMax);
2818 
2819  if (tofPatchMaxUS < 28 ) fTriggerClusterBC = 0 ;
2820  else if(tofPatchMaxUS < 75 ) fTriggerClusterBC = 1 ;
2821  else if(tofPatchMaxUS < 125) fTriggerClusterBC = 2 ;
2822  else if(tofPatchMaxUS < 175) fTriggerClusterBC = 3 ;
2823  else if(tofPatchMaxUS < 225) fTriggerClusterBC = 4 ;
2824  else if(tofPatchMaxUS < 275) fTriggerClusterBC = 5 ;
2825  else
2826  {
2827  //printf("AliCaloTrackReader::MatchTriggerCluster() - Large BC - tof %2.3f - Index %d\n",tofPatchMaxUS,fTriggerClusterIndex);
2829  else
2830  {
2831  fTriggerClusterIndex = -2;
2832  fTriggerClusterId = -2;
2833  }
2834  }
2835 
2836  if(tofPatchMax < 0) fTriggerClusterBC*=-1;
2837 
2838 
2839  // 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",
2840  // fTriggerClusterIndex, fTriggerClusterId,ePatchMax, tofPatchMax,
2841  // fTriggerClusterBC, fIsBadCellEvent,fIsBadMaxCellEvent,fIsExoticEvent, fIsTriggerMatch, nOfHighECl,absIdMaxMax);
2842  //
2843  // 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",
2844  // clusMax, idclusMax, eMax,tofMax, badClMax, badCeMax,exoMax);
2845 
2846  //Redo matching but open cuts
2847  if(!fIsTriggerMatch && fTriggerClusterId >= 0)
2848  {
2849  // Open time patch time
2850  TArrayI patchOpen = GetTriggerPatches(7,10);
2851 
2852  Int_t patchAbsIdOpenTime = -1;
2853  for(Int_t iabsId =0; iabsId < patchOpen.GetSize(); iabsId++)
2854  {
2855  Int_t absIDCell[4];
2856  patchAbsIdOpenTime = patchOpen.At(iabsId);
2857  GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patchAbsIdOpenTime, absIDCell);
2858  //if(tof > 75 ) printf("E %2.2f TOF %2.2f Trigger patch %d, cells : %d, %d, %d, %d\n",
2859  // clus->E(),tof,patches.At(iabsId), absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
2860 
2861  for(Int_t ipatch = 0; ipatch < 4; ipatch++)
2862  {
2863  if(absIdMaxMax == absIDCell[ipatch])
2864  {
2865  fIsTriggerMatchOpenCut[0] = kTRUE;
2866  break;
2867  }
2868  }// cell patch loop
2869  }// trigger patch loop
2870 
2871  // Check neighbour patches
2872  Int_t patchAbsId = -1;
2873  Int_t globalCol = -1;
2874  Int_t globalRow = -1;
2875  GetCaloUtils()->GetEMCALGeometry()->GetFastORIndexFromCellIndex(absIdMaxMax, patchAbsId);
2876  GetCaloUtils()->GetEMCALGeometry()->GetPositionInEMCALFromAbsFastORIndex(patchAbsId,globalCol,globalRow);
2877 
2878  // Check patches with strict time cut
2879  Int_t patchAbsIdNeigh = -1;
2880  for(Int_t icol = globalCol-1; icol <= globalCol+1; icol++)
2881  {
2882  if(icol < 0 || icol > 47) continue;
2883 
2884  for(Int_t irow = globalRow; irow <= globalRow+1; irow++)
2885  {
2886  if(irow < 0 || irow > 63) continue;
2887 
2888  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(icol, irow, patchAbsIdNeigh);
2889 
2890  if ( patchAbsIdNeigh < 0 ) continue;
2891 
2892  for(Int_t iabsId =0; iabsId < patches.GetSize(); iabsId++)
2893  {
2894  if(patchAbsIdNeigh == patches.At(iabsId))
2895  {
2896  fIsTriggerMatchOpenCut[1] = kTRUE;
2897  break;
2898  }
2899  }// trigger patch loop
2900 
2901  }// row
2902  }// col
2903 
2904  // Check patches with open time cut
2905  Int_t patchAbsIdNeighOpenTime = -1;
2906  for(Int_t icol = globalCol-1; icol <= globalCol+1; icol++)
2907  {
2908  if(icol < 0 || icol > 47) continue;
2909 
2910  for(Int_t irow = globalRow; irow <= globalRow+1; irow++)
2911  {
2912  if(irow < 0 || irow > 63) continue;
2913 
2914  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(icol, irow, patchAbsIdNeighOpenTime);
2915 
2916  if ( patchAbsIdNeighOpenTime < 0 ) continue;
2917 
2918  for(Int_t iabsId =0; iabsId < patchOpen.GetSize(); iabsId++)
2919  {
2920  if(patchAbsIdNeighOpenTime == patchOpen.At(iabsId))
2921  {
2922  fIsTriggerMatchOpenCut[2] = kTRUE;
2923  break;
2924  }
2925  }// trigger patch loop
2926 
2927  }// row
2928  }// col
2929 
2930  // printf("No match, new match: Open time %d-%d, open Neigh %d-%d, both open %d-%d\n",fIsTriggerMatchOpenCut[0],patchAbsIdOpenTime,
2931  // fIsTriggerMatchOpenCut[1],patchAbsIdNeigh,
2932  // fIsTriggerMatchOpenCut[2],patchAbsIdNeighOpenTime);
2933 
2934  patchOpen.Reset();
2935 
2936  }// No trigger match found
2937  //printf("Trigger BC %d, Id %d, Index %d\n",fTriggerClusterBC,fTriggerClusterId,fTriggerClusterIndex);
2938 }
2939 
2940 //_________________________________________________________
2944 //_________________________________________________________
2946 {
2948  {
2949  // get object pointer
2950  AliVCaloTrigger *caloTrigger = GetInputEvent()->GetCaloTrigger( "EMCAL" );
2951 
2952  if ( fBitEGA == 6 )
2953  {
2954  if (IsEventEMCALL1Gamma1()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(1);
2955  else if(IsEventEMCALL1Gamma2()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(3);
2956  else if(IsEventEMCALL1Jet1 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(0);
2957  else if(IsEventEMCALL1Jet2 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(2);
2958 
2959  // printf("L1 trigger Threshold Jet1 %f, Gamma1 %f, Jet2 %f, Gamma2 %f\n",
2960  // 0.07874*caloTrigger->GetL1Threshold(0),
2961  // 0.07874*caloTrigger->GetL1Threshold(1),
2962  // 0.07874*caloTrigger->GetL1Threshold(2),
2963  // 0.07874*caloTrigger->GetL1Threshold(3));
2964  }
2965  else
2966  {
2967  // Old AOD data format, in such case, order of thresholds different!!!
2968  if (IsEventEMCALL1Gamma1()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(0);
2969  else if(IsEventEMCALL1Gamma2()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(2);
2970  else if(IsEventEMCALL1Jet1 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(1);
2971  else if(IsEventEMCALL1Jet2 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(3);
2972 
2973  // printf("L1 trigger Threshold Jet1 %f, Gamma1 %f, Jet2 %f, Gamma2 %f\n",
2974  // 0.07874*caloTrigger->GetL1Threshold(1),
2975  // 0.07874*caloTrigger->GetL1Threshold(0),
2976  // 0.07874*caloTrigger->GetL1Threshold(3),
2977  // 0.07874*caloTrigger->GetL1Threshold(2));
2978  }
2979  }
2980 
2981  // Set L0 threshold, if not set by user
2983  {
2984  // Revise for periods > LHC11d
2985  Int_t runNumber = fInputEvent->GetRunNumber();
2986  if (runNumber < 146861) fTriggerL0EventThreshold = 3. ; // LHC11a
2987  else if(runNumber < 154000) fTriggerL0EventThreshold = 4. ; // LHC11b,c
2988  else if(runNumber < 165000) fTriggerL0EventThreshold = 5.5; // LHC11c,d,e
2989  else if(runNumber < 194000) fTriggerL0EventThreshold = 2 ; // LHC12
2990  else if(runNumber < 197400) fTriggerL0EventThreshold = 3 ; // LHC13def
2991  else if(runNumber < 197400) fTriggerL0EventThreshold = 2 ; // LHC13g
2992  else if(runNumber < 244300) fTriggerL0EventThreshold = 5 ; // LHC15 in, phys 1, 5 in phys2
2993  else if(runNumber < 266400) fTriggerL0EventThreshold = 2.5; // LHC16ir
2994  else fTriggerL0EventThreshold = 3.5; // LHC16s
2995  }
2996 }
2997 
2998 //________________________________________________________
3000 //________________________________________________________
3001 void AliCaloTrackReader::Print(const Option_t * opt) const
3002 {
3003  if(! opt)
3004  return;
3005 
3006  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
3007  printf("Task name : %s\n", fTaskName.Data()) ;
3008  printf("Data type : %d\n", fDataType) ;
3009  printf("CTS Min pT : %2.1f GeV/c\n", fCTSPtMin) ;
3010  printf("EMCAL Min pT : %2.1f GeV/c\n", fEMCALPtMin) ;
3011  printf("PHOS Min pT : %2.1f GeV/c\n", fPHOSPtMin) ;
3012  printf("CTS Max pT : %2.1f GeV/c\n", fCTSPtMax) ;
3013  printf("EMCAL Max pT : %2.1f GeV/c\n", fEMCALPtMax) ;
3014  printf("PHOS Max pT : %2.1f GeV/c\n", fPHOSPtMax) ;
3015  printf("EMCAL Bad Dist > %2.1f \n" , fEMCALBadChMinDist) ;
3016  printf("PHOS Bad Dist > %2.1f \n" , fPHOSBadChMinDist) ;
3017  printf("EMCAL N cells > %d \n" , fEMCALNCellsCut) ;
3018  printf("PHOS N cells > %d \n" , fPHOSNCellsCut) ;
3019  printf("EMCAL Time Cut: %3.1f < TOF < %3.1f\n", fEMCALTimeCutMin, fEMCALTimeCutMax);
3020  printf("Use CTS = %d\n", fFillCTS) ;
3021  printf("Use EMCAL = %d\n", fFillEMCAL) ;
3022  printf("Use DCAL = %d\n", fFillDCAL) ;
3023  printf("Use PHOS = %d\n", fFillPHOS) ;
3024  printf("Use EMCAL Cells = %d\n", fFillEMCALCells) ;
3025  printf("Use PHOS Cells = %d\n", fFillPHOSCells) ;
3026  printf("Track status = %d\n", (Int_t) fTrackStatus) ;
3027 
3028  printf("Track Mult Eta Cut = %2.2f\n", fTrackMultEtaCut) ;
3029 
3030  printf("Track Mult Pt Cuts:") ;
3031  for(Int_t icut = 0; icut < fTrackMultNPtCut; icut++) printf(" %2.2f GeV;",fTrackMultPtCut[icut]);
3032  printf(" \n") ;
3033 
3034  printf("Write delta AOD = %d\n", fWriteOutputDeltaAOD) ;
3035  printf("Recalculate Clusters = %d, E linearity = %d\n", fRecalculateClusters, fCorrectELinearity) ;
3036 
3037  printf("Use Triggers selected in SE base class %d; If not what Trigger Mask? %d; MB Trigger Mask for mixed %d \n",
3039 
3041  printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
3042 
3044  printf("Compare cluster pt and pt hard to accept event, factor = %2.2f",fPtHardAndClusterPtFactor);
3045 
3046  printf("Delta AOD File Name = %s\n", fDeltaAODFileName.Data()) ;
3047  printf("Centrality: Class %s, Option %d, Bin [%d,%d] \n", fCentralityClass.Data(),fCentralityOpt,fCentralityBin[0], fCentralityBin[1]) ;
3048 
3049  printf(" \n") ;
3050 }
3051 
3052 //__________________________________________
3056 //__________________________________________
3058 {
3059  // Count number of cells with energy larger than 0.1 in SM3, cut on this number
3060  Int_t ncellsSM3 = 0;
3061  for(Int_t icell = 0; icell < fInputEvent->GetEMCALCells()->GetNumberOfCells(); icell++)
3062  {
3063  Int_t absID = fInputEvent->GetEMCALCells()->GetCellNumber(icell);
3064  Int_t sm = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
3065  if(fInputEvent->GetEMCALCells()->GetAmplitude(icell) > 0.1 && sm==3) ncellsSM3++;
3066  }
3067 
3068  Int_t ncellcut = 21;
3069  if(GetFiredTriggerClasses().Contains("EMC")) ncellcut = 35;
3070 
3071  if(ncellsSM3 >= ncellcut)
3072  {
3073  AliDebug(1,Form("Reject event with ncells in SM3 %d, cut %d, trig %s",
3074  ncellsSM3,ncellcut,GetFiredTriggerClasses().Data()));
3075  return kTRUE;
3076  }
3077 
3078  return kFALSE;
3079 }
3080 
3081 //_________________________________________________________
3085 //_________________________________________________________
3087 {
3088  if(label < 0) return ;
3089 
3090  AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fInputEvent) ;
3091  if(!evt) return ;
3092 
3093  TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
3094  if(!arr) return ;
3095 
3096  if(label < arr->GetEntriesFast())
3097  {
3098  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
3099  if(!particle) return ;
3100 
3101  if(label == particle->Label()) return ; // label already OK
3102  //else printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
3103  }
3104  //else printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
3105 
3106  // loop on the particles list and check if there is one with the same label
3107  for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
3108  {
3109  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
3110  if(!particle) continue ;
3111 
3112  if(label == particle->Label())
3113  {
3114  label = ind;
3115  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
3116  return;
3117  }
3118  }
3119 
3120  label = -1;
3121 
3122  //printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label not found set to -1 \n");
3123 }
3124 
3125 //___________________________________
3127 //___________________________________
3129 {
3130  if(fCTSTracks) fCTSTracks -> Clear();
3131  if(fEMCALClusters) fEMCALClusters -> Clear("C");
3132  if(fPHOSClusters) fPHOSClusters -> Clear("C");
3133 
3134  fV0ADC[0] = 0; fV0ADC[1] = 0;
3135  fV0Mul[0] = 0; fV0Mul[1] = 0;
3136 
3137  if(fNonStandardJets) fNonStandardJets -> Clear("C");
3138  fBackgroundJets->Reset();
3139 }
3140 
3141 //___________________________________________
3145 //___________________________________________
3147 {
3148  fEventTrigMinBias = kFALSE;
3149  fEventTrigCentral = kFALSE;
3150  fEventTrigSemiCentral = kFALSE;
3151  fEventTrigEMCALL0 = kFALSE;
3152  fEventTrigEMCALL1Gamma1 = kFALSE;
3153  fEventTrigEMCALL1Gamma2 = kFALSE;
3154  fEventTrigEMCALL1Jet1 = kFALSE;
3155  fEventTrigEMCALL1Jet2 = kFALSE;
3156 
3157  AliDebug(1,Form("Select trigger mask bit %d - Trigger Event %s - Select <%s>",
3159 
3160  if(fEventTriggerMask <=0 )// in case no mask set
3161  {
3162  // EMC triggered event? Which type?
3163  if( GetFiredTriggerClasses().Contains("-B-") || GetFiredTriggerClasses().Contains("-S-") || GetFiredTriggerClasses().Contains("-I-") )
3164  {
3165  if ( GetFiredTriggerClasses().Contains("EGA" ) ||
3166  GetFiredTriggerClasses().Contains("EG1" ) )
3167  {
3168  fEventTrigEMCALL1Gamma1 = kTRUE;
3169  if( GetFiredTriggerClasses().Contains("EG1" ) && !fFiredTriggerClassName.Contains("EG1") ) fEventTrigEMCALL1Gamma1 = kFALSE;
3170  }
3171  else if( GetFiredTriggerClasses().Contains("EG2" ) )
3172  {
3173  fEventTrigEMCALL1Gamma2 = kTRUE;
3174  if( !fFiredTriggerClassName.Contains("EG2") ) fEventTrigEMCALL1Gamma2 = kFALSE;
3175  }
3176  else if( GetFiredTriggerClasses().Contains("EJE" ) ||
3177  GetFiredTriggerClasses().Contains("EJ1" ) )
3178  {
3179  fEventTrigEMCALL1Jet1 = kTRUE;
3180  if( GetFiredTriggerClasses().Contains("EJ1" ) && !fFiredTriggerClassName.Contains("EJ1") )
3181  fEventTrigEMCALL1Jet1 = kFALSE;
3182  }
3183  else if( GetFiredTriggerClasses().Contains("EJ2" ) )
3184  {
3185  fEventTrigEMCALL1Jet2 = kTRUE;
3186  if( !fFiredTriggerClassName.Contains("EJ2") ) fEventTrigEMCALL1Jet2 = kFALSE;
3187  }
3188  else if( GetFiredTriggerClasses().Contains("CEMC") &&
3189  !GetFiredTriggerClasses().Contains("EGA" ) &&
3190  !GetFiredTriggerClasses().Contains("EJE" ) &&
3191  !GetFiredTriggerClasses().Contains("EG1" ) &&
3192  !GetFiredTriggerClasses().Contains("EJ1" ) &&
3193  !GetFiredTriggerClasses().Contains("EG2" ) &&
3194  !GetFiredTriggerClasses().Contains("EJ2" ) )
3195  {
3196  fEventTrigEMCALL0 = kTRUE;
3197  }
3198 
3199  //Min bias event trigger?
3200  if (GetFiredTriggerClasses().Contains("CCENT_R2-B-NOPF-ALLNOTRD"))
3201  {
3202  fEventTrigCentral = kTRUE;
3203  }
3204  else if(GetFiredTriggerClasses().Contains("CSEMI_R1-B-NOPF-ALLNOTRD"))
3205  {
3206  fEventTrigSemiCentral = kTRUE;
3207  }
3208  else if((GetFiredTriggerClasses().Contains("CINT") || GetFiredTriggerClasses().Contains("CPBI2_B1") ) &&
3209  GetFiredTriggerClasses().Contains("-NOPF-ALLNOTRD") )
3210  {
3211  fEventTrigMinBias = kTRUE;
3212  }
3213  }
3214  }
3215  else
3216  {
3217  // EMC L1 Gamma
3218  if ( fEventTriggerMask & AliVEvent::kEMCEGA )
3219  {
3220  //printf("EGA trigger bit\n");
3221  if (GetFiredTriggerClasses().Contains("EG"))
3222  {
3223  if (GetFiredTriggerClasses().Contains("EGA")) fEventTrigEMCALL1Gamma1 = kTRUE;
3224  else
3225  {
3226  if(GetFiredTriggerClasses().Contains("EG1")) fEventTrigEMCALL1Gamma1 = kTRUE;
3227  if(GetFiredTriggerClasses().Contains("EG2")) fEventTrigEMCALL1Gamma2 = kTRUE;
3228  }
3229  }
3230  }
3231  // EMC L1 Jet
3232  else if( fEventTriggerMask & AliVEvent::kEMCEJE )
3233  {
3234  //printf("EGA trigger bit\n");
3235  if (GetFiredTriggerClasses().Contains("EJ"))
3236  {
3237  if (GetFiredTriggerClasses().Contains("EJE")) fEventTrigEMCALL1Jet1 = kTRUE;
3238  else
3239  {
3240  if(GetFiredTriggerClasses().Contains("EJ1")) fEventTrigEMCALL1Jet1 = kTRUE;
3241  if(GetFiredTriggerClasses().Contains("EJ2")) fEventTrigEMCALL1Jet2 = kTRUE;
3242  }
3243  }
3244  }
3245  // EMC L0
3246  else if((fEventTriggerMask & AliVEvent::kEMC7) ||
3247  (fEventTriggerMask & AliVEvent::kEMC1) )
3248  {
3249  //printf("L0 trigger bit\n");
3250  fEventTrigEMCALL0 = kTRUE;
3251  }
3252  // Min Bias Pb-Pb
3253  else if( fEventTriggerMask & AliVEvent::kCentral )
3254  {
3255  //printf("MB semi central trigger bit\n");
3256  fEventTrigSemiCentral = kTRUE;
3257  }
3258  // Min Bias Pb-Pb
3259  else if( fEventTriggerMask & AliVEvent::kSemiCentral )
3260  {
3261  //printf("MB central trigger bit\n");
3262  fEventTrigCentral = kTRUE;
3263  }
3264  // Min Bias pp, PbPb, pPb
3265  else if((fEventTriggerMask & AliVEvent::kMB ) ||
3266  (fEventTriggerMask & AliVEvent::kINT7) ||
3267  (fEventTriggerMask & AliVEvent::kINT8) ||
3268  (fEventTriggerMask & AliVEvent::kAnyINT) )
3269  {
3270  //printf("MB trigger bit\n");
3271  fEventTrigMinBias = kTRUE;
3272  }
3273  }
3274 
3275  AliDebug(1,Form("Event bits: \n \t MB %d, Cen %d, Sem %d, L0 %d, L1G1 %d, L1G2 %d, L1J1 %d, L1J2 %d",
3279 
3280  // L1 trigger bit
3281  if( fBitEGA == 0 && fBitEJE == 0 )
3282  {
3283  // Init the trigger bit once, correct depending on AliESD(AOD)CaloTrigger header version
3284 
3285  // Simpler way to do it ...
3286 // if( fInputEvent->GetRunNumber() < 172000 )
3287 // reader->SetEventTriggerL1Bit(4,5); // current LHC11 data
3288 // else
3289 // reader->SetEventTriggerL1Bit(6,8); // LHC12-13 data
3290 
3291  // Old values
3292  fBitEGA = 4;
3293  fBitEJE = 5;
3294 
3295  TFile* file = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
3296 
3297  const TList *clist = file->GetStreamerInfoCache();
3298 
3299  if(clist)
3300  {
3301  TStreamerInfo *cinfo = (TStreamerInfo*)clist->FindObject("AliESDCaloTrigger");
3302  Int_t verid = 5; // newer ESD header version
3303  if(!cinfo)
3304  {
3305  cinfo = (TStreamerInfo*)clist->FindObject("AliAODCaloTrigger");
3306  verid = 3; // newer AOD header version
3307  }
3308 
3309  if(cinfo)
3310  {
3311  Int_t classversionid = cinfo->GetClassVersion();
3312  //printf("********* Header class version %d *********** \n",classversionid);
3313 
3314  if (classversionid >= verid)
3315  {
3316  fBitEGA = 6;
3317  fBitEJE = 8;
3318  }
3319  } else AliInfo("AliCaloTrackReader()::SetEventTriggerBit() - Streamer info for trigger class not available, bit not changed");
3320  } else AliInfo("AliCaloTrackReader::SetEventTriggerBit() - Streamer list not available!, bit not changed");
3321 
3322  } // set once the EJE, EGA trigger bit
3323 }
3324 
3325 //____________________________________________________________
3328 //____________________________________________________________
3329 void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
3330 {
3331  fInputEvent = input;
3332  fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ;
3333  if (fMixedEvent)
3334  fNMixedEvent = fMixedEvent->GetNumberOfEvents() ;
3335 
3336  //Delete previous vertex
3337  if(fVertex)
3338  {
3339  for (Int_t i = 0; i < fNMixedEvent; i++)
3340  {
3341  delete [] fVertex[i] ;
3342  }
3343  delete [] fVertex ;
3344  }
3345 
3346  fVertex = new Double_t*[fNMixedEvent] ;
3347  for (Int_t i = 0; i < fNMixedEvent; i++)
3348  {
3349  fVertex[i] = new Double_t[3] ;
3350  fVertex[i][0] = 0.0 ;
3351  fVertex[i][1] = 0.0 ;
3352  fVertex[i][2] = 0.0 ;
3353  }
3354 }
3355 
3356 
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
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
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
Bool_t fRemoveLEDEvents
Remove events where LED was wrongly firing - EMCAL LHC11a.
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.
ClassImp(AliAnalysisTaskDeltaPt) AliAnalysisTaskDeltaPt
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