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