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