AliPhysics  63e47e1 (63e47e1)
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",
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  AliDebug(1,Form("Centrality %d in [%d,%d]?", cen, fCentralityBin[0], fCentralityBin[1]));
1421 
1422  if ( cen >= fCentralityBin[1] || cen < fCentralityBin[0] ) return kFALSE; //reject events out of bin.
1423 
1424  AliDebug(1,"Pass centrality rejection");
1425 
1426  fhNEventsAfterCut->Fill(14.5);
1427  }
1428 
1429  //----------------------------------------------------------------
1430  // MC events selections
1431  //----------------------------------------------------------------
1432  if ( GetMC() )
1433  {
1434  //----------------------------------------------------------------
1435  // Get the event headers
1436  //----------------------------------------------------------------
1437 
1438  // Main header
1439  // Init it first to 0 to tell the method to recover it.
1440  fGenEventHeader = 0;
1442 
1443  if ( fGenEventHeader )
1444  {
1445  AliDebug(1,Form("Selected event header class <%s>, name <%s>; cocktail %p",
1446  fGenEventHeader->ClassName(),
1447  fGenEventHeader->GetName(),
1448  GetMC()->GetCocktailList()));
1449  }
1450 
1451  //----------------------------------------------------------------
1452  // Reject the event if the event header name is not
1453  // the one requested among the possible generators.
1454  // Needed in case of cocktail MC generation with multiple options.
1455  //----------------------------------------------------------------
1457  {
1458  if ( !fGenEventHeader ) return kFALSE;
1459 
1460  AliDebug(1,"Pass Event header selection");
1461 
1462  fhNEventsAfterCut->Fill(15.5);
1463  }
1464 
1465  // Pythia header
1466  TString pyGenName = "";
1467  TString pyProcessName = "";
1468  Int_t pyProcess = 0;
1469  Int_t pyFirstGenPart = 0;
1470  Int_t pythiaVersion = 0;
1471 
1472  // Init it first to 0 to tell the method to recover it.
1475  pyGenName,pyProcessName,pyProcess,pyFirstGenPart,pythiaVersion);
1476 
1478  {
1479  AliDebug(2,Form("Pythia v%d name <%s>, process %d <%s>, first generated particle %d",
1480  pythiaVersion, pyGenName.Data(), pyProcess, pyProcessName.Data(), pyFirstGenPart));
1481 
1482  //---------------------------------------------------------------------------
1483  // In case of analysis of events with jets, skip those with jet pt > 5 pt hard
1484  // To be used on for MC data in pT hard bins
1485  //---------------------------------------------------------------------------
1486 
1488  {
1489  if ( !ComparePtHardAndJetPt(pyProcess, pyProcessName) ) return kFALSE ;
1490 
1491  AliDebug(1,"Pass Pt Hard - Jet rejection");
1492 
1493  fhNEventsAfterCut->Fill(16.5);
1494  }
1495 
1497  {
1498  if ( !ComparePtHardAndClusterPt(pyProcess, pyProcessName) ) return kFALSE ;
1499 
1500  AliDebug(1,"Pass Pt Hard - Cluster rejection");
1501 
1502  fhNEventsAfterCut->Fill(17.5);
1503  }
1504  } // pythia header
1505  } // MC
1506 
1507  //------------------------------------------------------------------
1508  // Recover the weight assigned to the event, if provided
1509  // right now only for pT-hard bins and centrality dependent weights
1510  //------------------------------------------------------------------
1512  {
1514 
1516 
1518  }
1519 
1520  //-------------------------------------------------------
1521  // Get the main vertex BC, in case not available
1522  // it is calculated in FillCTS checking the BC of tracks
1523  //------------------------------------------------------
1524  fVertexBC = fInputEvent->GetPrimaryVertex()->GetBC();
1525 
1526  //-----------------------------------------------
1527  // Fill the arrays with cluster/tracks/cells data
1528  //-----------------------------------------------
1529 
1530  if(fFillCTS)
1531  {
1532  FillInputCTS();
1533  //Accept events with at least one track
1534  if(fTrackMult[0] == 0 && fDoRejectNoTrackEvents) return kFALSE ;
1535 
1536  AliDebug(1,"Pass rejection of null track events");
1537 
1538  fhNEventsAfterCut->Fill(18.5);
1539  }
1540 
1542  {
1543  if(fVertexBC != 0 && fVertexBC != AliVTrack::kTOFBCNA) return kFALSE ;
1544 
1545  AliDebug(1,"Pass rejection of events with vertex at BC!=0");
1546 
1547  fhNEventsAfterCut->Fill(19.5);
1548  }
1549 
1550  if(fFillEMCALCells)
1552 
1553  if(fFillPHOSCells)
1555 
1556  if(fFillEMCAL || fFillDCAL)
1557  FillInputEMCAL();
1558 
1559  if(fFillPHOS)
1560  FillInputPHOS();
1561 
1562  FillInputVZERO();
1563 
1564  //one specified jet branch
1569 
1570  AliDebug(1,"Event accepted for analysis");
1571 
1572  return kTRUE ;
1573 }
1574 
1575 //__________________________________________________
1578 //__________________________________________________
1580 {
1581  if(fUseAliCentrality)
1582  {
1583  if ( !GetCentrality() ) return -1;
1584 
1585  AliDebug(1,Form("Cent. Percentile: V0M %2.2f, CL0 %2.2f, CL1 %2.2f; selected class %s",
1586  GetCentrality()->GetCentralityPercentile("V0M"),
1587  GetCentrality()->GetCentralityPercentile("CL0"),
1588  GetCentrality()->GetCentralityPercentile("CL1"),
1589  fCentralityClass.Data()));
1590 
1591  if (fCentralityOpt == 100) return (Int_t) GetCentrality()->GetCentralityPercentile(fCentralityClass); // 100 bins max
1592  else if(fCentralityOpt == 10) return GetCentrality()->GetCentralityClass10(fCentralityClass);// 10 bins max
1593  else if(fCentralityOpt == 20) return GetCentrality()->GetCentralityClass5(fCentralityClass); // 20 bins max
1594  else
1595  {
1596  AliInfo(Form("Unknown centrality option %d, use 10, 20 or 100",fCentralityOpt));
1597  return -1;
1598  }
1599  }
1600  else
1601  {
1602  if ( !GetMultSelCen() ) return -1;
1603 
1604  AliDebug(1,Form("Mult. Percentile: V0M %2.2f, CL0 %2.2f, CL1 %2.2f; selected class %s",
1605  GetMultSelCen()->GetMultiplicityPercentile("V0M",1),
1606  GetMultSelCen()->GetMultiplicityPercentile("CL0",1),
1607  GetMultSelCen()->GetMultiplicityPercentile("CL1",1),
1608  fCentralityClass.Data()));
1609 
1610  return GetMultSelCen()->GetMultiplicityPercentile(fCentralityClass, kTRUE); // returns centrality only for events used in calibration
1611 
1612  // equivalent to
1613  //GetMultSelCen()->GetMultiplicityPercentile("V0M", kFALSE); // returns centrality for any event
1614  //Int_t qual = GetMultSelCen()->GetEvSelCode(); if (qual ! = 0) cent = qual;
1615  }
1616 }
1617 
1618 //_____________________________________________________
1621 //_____________________________________________________
1623 {
1624  if( !GetEventPlane() ) return -1000;
1625 
1626  Float_t ep = GetEventPlane()->GetEventplane(GetEventPlaneMethod(), GetInputEvent());
1627 
1628  if(GetEventPlaneMethod()=="Q" && (ep < 0 || ep > TMath::Pi()))
1629  {
1630  AliDebug(1,Form("Bad EP for <Q> method : %f",ep));
1631  return -1000;
1632  }
1633  else if(GetEventPlaneMethod().Contains("V0") )
1634  {
1635  if((ep > TMath::Pi()/2 || ep < -TMath::Pi()/2))
1636  {
1637  AliDebug(1,Form("Bad EP for <%s> method : %f",GetEventPlaneMethod().Data(), ep));
1638  return -1000;
1639  }
1640 
1641  ep+=TMath::Pi()/2; // put same range as for <Q> method
1642  }
1643 
1644  AliDebug(3,Form("Event plane angle %f",ep));
1645 
1646 // if(fDebug > 0 )
1647 // {
1648 // if (ep > TMath::Pi()) printf("AliCaloTrackReader::GetEventPlaneAngle() - Too large angle = %f\n",ep);
1649 // else if(ep < 0 ) printf("AliCaloTrackReader::GetEventPlaneAngle() - Negative angle = %f\n" ,ep);
1650 // }
1651 
1652  return ep;
1653 }
1654 
1655 //__________________________________________________________
1657 //__________________________________________________________
1659 {
1660  vertex[0] = fVertex[0][0];
1661  vertex[1] = fVertex[0][1];
1662  vertex[2] = fVertex[0][2];
1663 }
1664 
1665 //__________________________________________________________________________
1667 //__________________________________________________________________________
1668 void AliCaloTrackReader::GetVertex(Double_t vertex[3], Int_t evtIndex) const
1669 {
1670  vertex[0] = fVertex[evtIndex][0];
1671  vertex[1] = fVertex[evtIndex][1];
1672  vertex[2] = fVertex[evtIndex][2];
1673 }
1674 
1675 //________________________________________
1678 //________________________________________
1680 {
1681  // Delete previous vertex
1682  if(fVertex)
1683  {
1684  for (Int_t i = 0; i < fNMixedEvent; i++)
1685  {
1686  delete [] fVertex[i] ;
1687  }
1688  delete [] fVertex ;
1689  }
1690 
1691  fVertex = new Double_t*[fNMixedEvent] ;
1692  for (Int_t i = 0; i < fNMixedEvent; i++)
1693  {
1694  fVertex[i] = new Double_t[3] ;
1695  fVertex[i][0] = 0.0 ;
1696  fVertex[i][1] = 0.0 ;
1697  fVertex[i][2] = 0.0 ;
1698  }
1699 
1700  if (!fMixedEvent)
1701  { // Single event analysis
1702  if(fDataType!=kMC)
1703  {
1704 
1705  if(fInputEvent->GetPrimaryVertex())
1706  {
1707  fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]);
1708  }
1709  else
1710  {
1711  AliWarning("NULL primary vertex");
1712  fVertex[0][0]=0.; fVertex[0][1]=0.; fVertex[0][2]=0.;
1713  }//Primary vertex pointer do not exist
1714 
1715  } else
1716  {// MC read event
1717  fVertex[0][0]=0.; fVertex[0][1]=0.; fVertex[0][2]=0.;
1718  }
1719 
1720  AliDebug(1,Form("Single Event Vertex : %f,%f,%f",fVertex[0][0],fVertex[0][1],fVertex[0][2]));
1721 
1722  } else
1723  { // MultiEvent analysis
1724  for (Int_t iev = 0; iev < fNMixedEvent; iev++)
1725  {
1726  if (fMixedEvent->GetVertexOfEvent(iev))
1727  fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
1728  else
1729  AliWarning("No vertex found");
1730 
1731  AliDebug(1,Form("Multi Event %d Vertex : %f,%f,%f",iev,fVertex[iev][0],fVertex[iev][1],fVertex[iev][2]));
1732  }
1733  }
1734 }
1735 
1736 //_____________________________________
1742 //_____________________________________
1744 {
1745  AliDebug(1,"Begin");
1746 
1747  Double_t pTrack[3] = {0,0,0};
1748 
1749  Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
1750  Int_t nstatus = 0;
1751  Double_t bz = GetInputEvent()->GetMagneticField();
1752 
1753  for(Int_t i = 0; i < 19; i++)
1754  {
1755  fTrackBCEvent [i] = 0;
1756  fTrackBCEventCut[i] = 0;
1757  }
1758 
1759  for(Int_t iptCut = 0; iptCut < fTrackMultNPtCut; iptCut++ )
1760  {
1761  fTrackMult [iptCut] = 0;
1762  fTrackSumPt[iptCut] = 0;
1763  }
1764 
1765  Bool_t bc0 = kFALSE;
1766  if(fRecalculateVertexBC) fVertexBC = AliVTrack::kTOFBCNA;
1767 
1768  for (Int_t itrack = 0; itrack < nTracks; itrack++)
1769  {
1770  AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
1771 
1772  if ( !AcceptParticleMCLabel( TMath::Abs(track->GetLabel()) ) ) continue ;
1773 
1774  fhCTSTrackCutsPt[0]->Fill(track->Pt());
1775 
1776  //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
1777  ULong_t status = track->GetStatus();
1778 
1779  if (fTrackStatus && !((status & fTrackStatus) == fTrackStatus))
1780  continue ;
1781 
1782  fhCTSTrackCutsPt[1]->Fill(track->Pt());
1783 
1784  nstatus++;
1785 
1786  //-------------------------
1787  // Select the tracks depending on cuts of AOD or ESD
1788  if(!SelectTrack(track, pTrack)) continue ;
1789 
1790  fhCTSTrackCutsPt[2]->Fill(track->Pt());
1791 
1792  //-------------------------
1793  // TOF cuts
1794  Bool_t okTOF = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
1795  Double_t tof = -1000;
1796  Int_t trackBC = -1000 ;
1797 
1798  if(fAccessTrackTOF)
1799  {
1800  if(okTOF)
1801  {
1802  trackBC = track->GetTOFBunchCrossing(bz);
1803  SetTrackEventBC(trackBC+9);
1804 
1805  tof = track->GetTOFsignal()*1e-3;
1806 
1807  //After selecting tracks with small DCA, pointing to vertex, set vertex BC depeding on tracks BC
1809  {
1810  if (trackBC != 0 && trackBC != AliVTrack::kTOFBCNA) fVertexBC = trackBC;
1811  else if(trackBC == 0) bc0 = kTRUE;
1812  }
1813 
1814  //In any case, the time should to be larger than the fixed window ...
1815  if( fUseTrackTimeCut && (trackBC !=0 || tof < fTrackTimeCutMin || tof > fTrackTimeCutMax) )
1816  {
1817  //printf("Remove track time %f and bc = %d\n",tof,trackBC);
1818  continue ;
1819  }
1820  //else printf("Accept track time %f and bc = %d\n",tof,trackBC);
1821  }
1822  }
1823 
1824  fhCTSTrackCutsPt[3]->Fill(track->Pt());
1825 
1826  //---------------------
1827  // DCA cuts
1828  //
1829  fMomentum.SetPxPyPzE(pTrack[0],pTrack[1],pTrack[2],0);
1830 
1831  if(fUseTrackDCACut)
1832  {
1833  Float_t dcaTPC =-999;
1834  //In case of AODs, TPC tracks cannot be propagated back to primary vertex,
1835  if( fDataType == kAOD ) dcaTPC = ((AliAODTrack*) track)->DCA();
1836 
1837  //normal way to get the dca, cut on dca_xy
1838  if(dcaTPC==-999)
1839  {
1840  Double_t dca[2] = {1e6,1e6};
1841  Double_t covar[3] = {1e6,1e6,1e6};
1842  Bool_t okDCA = track->PropagateToDCA(fInputEvent->GetPrimaryVertex(),bz,100.,dca,covar);
1843  if( okDCA) okDCA = AcceptDCA(fMomentum.Pt(),dca[0]);
1844  if(!okDCA)
1845  {
1846  //printf("AliCaloTrackReader::FillInputCTS() - Reject track pt %2.2f, dca_xy %2.4f\n",fMomentum.Pt(),dca[0]);
1847  continue ;
1848  }
1849  }
1850  }// DCA cuts
1851 
1852  fhCTSTrackCutsPt[4]->Fill(track->Pt());
1853 
1854  //-------------------------
1855  // Kinematic/acceptance cuts
1856  //
1857  // Count the tracks in eta < 0.9 and different pT cuts
1858  Float_t ptTrack = fMomentum.Pt();
1859  if(TMath::Abs(track->Eta())< fTrackMultEtaCut)
1860  {
1861  for(Int_t iptCut = 0; iptCut < fTrackMultNPtCut; iptCut++ )
1862  {
1863  if(ptTrack > fTrackMultPtCut[iptCut])
1864  {
1865  fTrackMult [iptCut]++;
1866  fTrackSumPt[iptCut]+=ptTrack;
1867  }
1868  }
1869  }
1870 
1871  if(fCTSPtMin > ptTrack || fCTSPtMax < ptTrack) continue ;
1872 
1873  // Check effect of cuts on track BC
1874  if(fAccessTrackTOF && okTOF) SetTrackEventBCcut(trackBC+9);
1875 
1876  if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS)) continue;
1877 
1878  fhCTSTrackCutsPt[5]->Fill(track->Pt());
1879 
1880  // ------------------------------
1881  // Add selected tracks to array
1882  AliDebug(2,Form("Selected tracks pt %3.2f, phi %3.2f deg, eta %3.2f",
1883  fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
1884 
1885  fCTSTracks->Add(track);
1886 
1887  // TODO, check if remove
1888  if (fMixedEvent) track->SetID(itrack);
1889 
1890  }// track loop
1891 
1892  if( fRecalculateVertexBC && (fVertexBC == 0 || fVertexBC == AliVTrack::kTOFBCNA))
1893  {
1894  if( bc0 ) fVertexBC = 0 ;
1895  else fVertexBC = AliVTrack::kTOFBCNA ;
1896  }
1897 
1898  AliDebug(1,Form("CTS entries %d, input tracks %d, pass status %d, multipliticy %d", fCTSTracks->GetEntriesFast(), nTracks, nstatus, fTrackMult[0]));//fCTSTracksNormalInputEntries);
1899 }
1900 
1901 //_______________________________________________________________________________
1919 //_______________________________________________________________________________
1921 {
1922  // Accept clusters with the proper label, only applicable for MC
1923  if ( clus->GetLabel() >= 0 ) // -1 corresponds to noisy MC
1924  {
1925  if ( !AcceptParticleMCLabel(clus->GetLabel()) ) return ;
1926  }
1927 
1928  // TODO, not sure if needed anymore
1929  Int_t vindex = 0 ;
1930  if (fMixedEvent)
1931  vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
1932 
1933  clus->GetMomentum(fMomentum, fVertex[vindex]);
1934 
1935  // No correction/cut applied yet
1936  fhEMCALClusterCutsE[0]->Fill(clus->E());
1937 
1938  // Get the maximum cell energy, its SM number and its col, row location, needed in
1939  // different places of this method, although not active by default, one can consider
1940  // deactivate this and only activate it when requiered.
1941  Int_t absIdMax= -1;
1942  Int_t iSupMod = -1;
1943  Int_t iphiMax = -1;
1944  Int_t ietaMax = -1;
1945  Bool_t shared = kFALSE;
1946  GetCaloUtils()->GetEMCALRecoUtils()->GetMaxEnergyCell(GetCaloUtils()->GetEMCALGeometry(),
1947  GetEMCALCells(),clus,absIdMax,iSupMod,
1948  ietaMax,iphiMax,shared);
1949 
1950  //if( (fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10 )
1951  AliDebug(2,Form("Input cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f, nCells %d, SM %d",
1952  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta(), clus->GetNCells(), iSupMod));
1953 
1954  //---------------------------
1955  // Embedding case
1956  // TO BE REVISED
1958  {
1959  if(clus->GetNLabels()==0 || clus->GetLabel() < 0) return;
1960  //else printf("Embedded cluster, %d, n label %d label %d \n",iclus,clus->GetNLabels(),clus->GetLabel());
1961  }
1962 
1963  //--------------------------------------
1964  // Apply some corrections in the cluster
1965  //
1967  {
1968  //Recalibrate the cluster energy
1970  {
1972 
1973  clus->SetE(energy);
1974  //printf("Recalibrated Energy %f\n",clus->E());
1975 
1978 
1979  } // recalculate E
1980 
1981  //Recalculate distance to bad channels, if new list of bad channels provided
1983 
1984  //Recalculate cluster position
1985  if(GetCaloUtils()->IsRecalculationOfClusterPositionOn())
1986  {
1988  //clus->GetPosition(pos);
1989  //printf("After Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1990  }
1991 
1992  // Recalculate TOF
1993  if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1994  {
1995  Double_t tof = clus->GetTOF();
1996  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1997 
1998  //additional L1 phase shift
2000  {
2001  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(iSupMod,fInputEvent->GetBunchCrossNumber(), tof);
2002  }
2003 
2004  clus->SetTOF(tof);
2005 
2006  }// Time recalibration
2007  }
2008 
2009  // Check effect of corrections
2010  fhEMCALClusterCutsE[1]->Fill(clus->E());
2011 
2012  //-----------------------------------------------------------------
2013  // Reject clusters with bad channels, close to borders and exotic
2014  //
2015  Bool_t goodCluster = GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,
2016  GetCaloUtils()->GetEMCALGeometry(),
2017  GetEMCALCells(),fInputEvent->GetBunchCrossNumber());
2018 
2019  if(!goodCluster)
2020  {
2021  //if( (fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10 )
2022  AliDebug(1,Form("Bad cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2023  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2024 
2025  return;
2026  }
2027 
2028  // Check effect of bad cluster removal
2029  fhEMCALClusterCutsE[2]->Fill(clus->E());
2030 
2031  //Float_t pos[3];
2032  //clus->GetPosition(pos);
2033  //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
2034 
2035  //-----------------------------------------------------
2036  // Correct non linearity, smear energy or scale per SM
2037  //
2039  {
2041 
2042  AliDebug(5,Form("Correct Non Lin: Old E %3.2f, New E %3.2f",
2043  fMomentum.E(),clus->E()));
2044  }
2045 
2046  // In case of MC analysis, to match resolution/calibration in real data
2047  // Not needed anymore, just leave for MC studies on systematics
2048  if( GetCaloUtils()->GetEMCALRecoUtils()->IsClusterEnergySmeared() )
2049  {
2050  Float_t rdmEnergy = GetCaloUtils()->GetEMCALRecoUtils()->SmearClusterEnergy(clus);
2051 
2052  AliDebug(5,Form("Smear energy: Old E %3.2f, New E %3.2f",clus->E(),rdmEnergy));
2053 
2054  clus->SetE(rdmEnergy);
2055  }
2056 
2057  // In case of uncalibrated data, or non final non linearity for MC
2058  // or calibration of MC for SMs behind TRD, apply a global scale factor per SM
2059  if ( fScaleEPerSM && iSupMod < 22 && iSupMod >=0)
2060  {
2061  Float_t scale = fScaleFactorPerSM[iSupMod];
2062 
2063  AliDebug(5,Form("Scale energy for SM %d: Old E %3.2f, scale factor %1.5f",iSupMod,clus->E(),scale));
2064 
2065  clus->SetE(clus->E()*scale);
2066  }
2067 
2068  clus->GetMomentum(fMomentum, fVertex[vindex]);
2069  fhEMCALClusterEtaPhi->Fill(fMomentum.Eta(),GetPhi(fMomentum.Phi()));
2070 
2071  // Check effect linearity correction, energy smearing
2072  fhEMCALClusterCutsE[3]->Fill(clus->E());
2073 
2074  // Check the event BC depending on EMCal clustr before final cuts
2075  Double_t tof = clus->GetTOF()*1e9;
2076 
2077  Int_t bc = TMath::Nint(tof/50) + 9;
2078  //printf("tof %2.2f, bc+5=%d\n",tof,bc);
2079 
2080  SetEMCalEventBC(bc);
2081 
2082  //--------------------------------------
2083  // Apply some kinematical/acceptance cuts
2084  //
2085  if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E())
2086  {
2087  AliDebug(2,Form("Cluster E out of range, %2.2f < %2.2f < %2.2f",fEMCALPtMin,clus->E(),fEMCALPtMax));
2088  return ;
2089  }
2090 
2091  // Select cluster fiducial region
2092  //
2093  Bool_t bEMCAL = kFALSE;
2094  Bool_t bDCAL = kFALSE;
2095  if(fCheckFidCut)
2096  {
2097  if(fFillEMCAL && fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) bEMCAL = kTRUE ;
2098  if(fFillDCAL && fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kDCAL )) bDCAL = kTRUE ;
2099  }
2100  else
2101  {
2102  bEMCAL = kTRUE;
2103  }
2104 
2105  //---------------------------------------------------------------------
2106  // Mask all cells in collumns facing ALICE thick material if requested
2107  //
2108  if ( GetCaloUtils()->GetNMaskCellColumns() )
2109  {
2110  AliDebug(2,Form("Masked collumn: cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2111  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2112 
2113  if ( GetCaloUtils()->MaskFrameCluster(iSupMod, ietaMax) )
2114  {
2115  AliDebug(2,"Mask cluster");
2116  return;
2117  }
2118  }
2119 
2120  // Check effect of energy and fiducial cuts
2121  fhEMCALClusterCutsE[4]->Fill(clus->E());
2122 
2123  if ( bEMCAL || bDCAL ) fhEMCALClusterEtaPhiFidCut->Fill(fMomentum.Eta(),GetPhi(fMomentum.Phi()));
2124 
2125  //----------------------------------------------------
2126  // Apply N cells cut
2127  //
2128  if(clus->GetNCells() <= fEMCALNCellsCut && fDataType != AliCaloTrackReader::kMC)
2129  {
2130  AliDebug(2,Form("Cluster with n cells %d < %d",clus->GetNCells(), fEMCALNCellsCut));
2131  return ;
2132  }
2133 
2134  // Check effect of n cells cut
2135  fhEMCALClusterCutsE[5]->Fill(clus->E());
2136 
2137  //----------------------------------------------------
2138  // Apply distance to bad channel cut
2139  //
2140  Double_t distBad = clus->GetDistanceToBadChannel() ; //Distance to bad channel
2141 
2142  if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2143 
2144  if(distBad < fEMCALBadChMinDist)
2145  {
2146  AliDebug(2, Form("Cluster close to bad, dist %2.2f < %2.2f",distBad,fEMCALBadChMinDist));
2147  return ;
2148  }
2149 
2150  // Check effect distance to bad channel cut
2151  fhEMCALClusterCutsE[6]->Fill(clus->E());
2152 
2153  //------------------------------------------
2154  // Apply time cut, count EMCal BC before cut
2155  //
2156  SetEMCalEventBCcut(bc);
2157 
2158  // Shift time in case of no calibration with rough factor
2159  Double_t tofShift = tof;
2160  if(tof > 400) tofShift-=615;
2161  fhEMCALClusterTimeE->Fill(clus->E(),tofShift);
2162 
2163  if(!IsInTimeWindow(tof,clus->E()))
2164  {
2165  fNPileUpClusters++ ;
2166  if(fUseEMCALTimeCut)
2167  {
2168  AliDebug(2,Form("Out of time window E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f, time %e",
2169  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta(),tof));
2170 
2171  return ;
2172  }
2173  }
2174  else
2176 
2177  // Check effect of time cut
2178  fhEMCALClusterCutsE[7]->Fill(clus->E());
2179 
2180  //----------------------------------------------------
2181  // Smear the SS to try to match data and simulations,
2182  // do it only for simulations.
2183  //
2184  if ( fSmearShowerShape && clus->GetNCells() > 2 )
2185  {
2186  Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(clus, GetEMCALCells());
2187 // Int_t nMaxima = clus->GetNExMax(); // For Run2
2188 
2189  if ( nMaxima >= fSmearNLMMin && nMaxima <= fSmearNLMMax )
2190  {
2191  AliDebug(2,Form("Smear shower shape - Original: %2.4f", clus->GetM02()));
2193  {
2194  clus->SetM02( clus->GetM02() + fRandom.Landau(0, fSmearShowerShapeWidth) );
2195  }
2197  {
2198  if(iclus%3 == 0 && clus->GetM02() > 0.1) clus->SetM02( clus->GetM02() + fRandom.Landau(0.05, fSmearShowerShapeWidth) ); //fSmearShowerShapeWidth = 0.035
2199  }
2200  else if (fSmearingFunction == kNoSmearing)
2201  {
2202  clus->SetM02( clus->GetM02() );
2203  }
2204  //clus->SetM02( fRandom.Landau(clus->GetM02(), fSmearShowerShapeWidth) );
2205  AliDebug(2,Form("Width %2.4f Smeared : %2.4f", fSmearShowerShapeWidth,clus->GetM02()));
2206  }
2207  }
2208 
2209  //--------------------------------------------------------
2210  // Fill the corresponding array with the selected clusters
2211  // Usually just filling EMCal array with upper or lower clusters is enough,
2212  // but maybe we want to do EMCal-DCal correlations.
2213 
2214  //if((fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10)
2215  AliDebug(2,Form("Selected clusters (EMCAL%d, DCAL%d), E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2216  bEMCAL,bDCAL,fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2217 
2218 
2219  if (bEMCAL) fEMCALClusters->Add(clus);
2220  else if(bDCAL ) fDCALClusters ->Add(clus);
2221 
2222  // TODO, not sure if needed anymore
2223  if (fMixedEvent)
2224  clus->SetID(iclus) ;
2225 }
2226 
2227 //_______________________________________
2234 //_______________________________________
2236 {
2237  AliDebug(1,"Begin");
2238 
2239  // First recalibrate cells, time or energy
2240  // if(GetCaloUtils()->IsRecalibrationOn())
2241  // GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCells(GetCaloUtils()->GetEMCALGeometry(),
2242  // GetEMCALCells(),
2243  // fInputEvent->GetBunchCrossNumber());
2244 
2245  fNPileUpClusters = 0; // Init counter
2246  fNNonPileUpClusters = 0; // Init counter
2247  for(Int_t i = 0; i < 19; i++)
2248  {
2249  fEMCalBCEvent [i] = 0;
2250  fEMCalBCEventCut[i] = 0;
2251  }
2252 
2253  //Loop to select clusters in fiducial cut and fill container with aodClusters
2254  if(fEMCALClustersListName=="")
2255  {
2256  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2257  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2258  {
2259  AliVCluster * clus = 0;
2260  if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
2261  {
2262  if (clus->IsEMCAL())
2263  {
2264  FillInputEMCALAlgorithm(clus, iclus);
2265  }//EMCAL cluster
2266  }// cluster exists
2267  }// cluster loop
2268 
2269  //Recalculate track matching
2271 
2272  }//Get the clusters from the input event
2273  else
2274  {
2275  TClonesArray * clusterList = 0x0;
2276 
2277  if (fInputEvent->FindListObject(fEMCALClustersListName))
2278  {
2279  clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
2280  }
2281  else if(fOutputEvent)
2282  {
2283  clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
2284  }
2285 
2286  if(!clusterList)
2287  {
2288  AliWarning(Form("Wrong name of list with clusters? <%s>",fEMCALClustersListName.Data()));
2289  return;
2290  }
2291 
2292  Int_t nclusters = clusterList->GetEntriesFast();
2293  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2294  {
2295  AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
2296  //printf("E %f\n",clus->E());
2297  if (clus) FillInputEMCALAlgorithm(clus, iclus);
2298  else AliWarning("Null cluster in list!");
2299  }// cluster loop
2300 
2301  // Recalculate the pile-up time, in case long time clusters removed during clusterization
2302  //printf("Input event INIT : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
2303 
2304  fNPileUpClusters = 0; // Init counter
2305  fNNonPileUpClusters = 0; // Init counter
2306  for(Int_t i = 0; i < 19; i++)
2307  {
2308  fEMCalBCEvent [i] = 0;
2309  fEMCalBCEventCut[i] = 0;
2310  }
2311 
2312  for (Int_t iclus = 0; iclus < fInputEvent->GetNumberOfCaloClusters(); iclus++)
2313  {
2314  AliVCluster * clus = 0;
2315 
2316  if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
2317  {
2318  if (clus->IsEMCAL())
2319  {
2320 
2321  Float_t frac =-1;
2322  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
2323  Double_t tof = clus->GetTOF();
2324  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2325  //additional L1 phase shift
2327  {
2328  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absIdMax), fInputEvent->GetBunchCrossNumber(), tof);
2329  }
2330 
2331  tof*=1e9;
2332 
2333  //printf("Input event cluster : AbsIdMax %d, E %2.2f, time %2.2f \n", absIdMax,clus->E(),tof);
2334 
2335  //Reject clusters with bad channels, close to borders and exotic;
2336  if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber())) continue;
2337 
2338  Int_t bc = TMath::Nint(tof/50) + 9;
2339  SetEMCalEventBC(bc);
2340 
2341  if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E()) continue ;
2342 
2343  clus->GetMomentum(fMomentum, fVertex[0]);
2344 
2345  if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) return ;
2346 
2347  SetEMCalEventBCcut(bc);
2348 
2349  if(!IsInTimeWindow(tof,clus->E()))
2350  fNPileUpClusters++ ;
2351  else
2353 
2354  }
2355  }
2356  }
2357 
2358  //printf("Input event : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
2359 
2360  // Recalculate track matching, not necessary if already done in the reclusterization task.
2361  // in case it was not done ...
2363 
2364  }
2365 
2366  AliDebug(1,Form("EMCal selected clusters %d",
2367  fEMCALClusters->GetEntriesFast()));
2368  AliDebug(2,Form("\t n pile-up clusters %d, n non pile-up %d",
2370 }
2371 
2372 //_______________________________________
2374 //_______________________________________
2376 {
2377  AliDebug(1,"Begin");
2378 
2379  // Loop to select clusters in fiducial cut and fill container with aodClusters
2380  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2381  TString genName;
2382  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2383  {
2384  AliVCluster * clus = fInputEvent->GetCaloCluster(iclus) ;
2385  if ( !clus ) continue ;
2386 
2387  if ( !clus->IsPHOS() ) continue ;
2388 
2389  if(clus->GetLabel() >=0 ) // -1 corresponds to noisy MC
2390  {
2391  if ( !AcceptParticleMCLabel(clus->GetLabel()) ) continue ;
2392  }
2393 
2394  fhPHOSClusterCutsE[0]->Fill(clus->E());
2395 
2396  // Skip CPV input
2397  if( clus->GetType() == AliVCluster::kPHOSCharged ) continue ;
2398 
2399  fhPHOSClusterCutsE[1]->Fill(clus->E());
2400 
2401  //---------------------------------------------
2402  // Old feature, try to rely on PHOS tender
2403  //
2405  {
2406  // Recalibrate the cluster energy
2408  {
2409  Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
2410  clus->SetE(energy);
2411  }
2412  }
2413 
2414  //----------------------------------------------------------------------------------
2415  // Check if the cluster contains any bad channel and if close to calorimeter borders
2416  //
2417  // Old feature, try to rely on PHOS tender
2418  if( GetCaloUtils()->ClusterContainsBadChannel(kPHOS,clus->GetCellsAbsId(), clus->GetNCells()))
2419  continue;
2420 
2421  if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells()))
2422  continue;
2423 
2424  // TODO, add exotic cut???
2425 
2426  fhPHOSClusterCutsE[2]->Fill(clus->E());
2427 
2428  // TODO Dead code? remove?
2429  Int_t vindex = 0 ;
2430  if (fMixedEvent)
2431  vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
2432 
2433  clus->GetMomentum(fMomentum, fVertex[vindex]);
2434 
2435  //----------------------------------------------------------------------------------
2436  // Remove clusters close to borders
2437  //
2439  continue ;
2440 
2441  fhPHOSClusterCutsE[3]->Fill(clus->E());
2442 
2443  //----------------------------------------------------------------------------------
2444  // Remove clusters with too low energy
2445  //
2446  if (fPHOSPtMin > fMomentum.E() || fPHOSPtMax < fMomentum.E() )
2447  continue ;
2448 
2449  fhPHOSClusterCutsE[4]->Fill(clus->E());
2450 
2451  //----------------------------------------------------
2452  // Apply N cells cut
2453  //
2454  if(clus->GetNCells() <= fPHOSNCellsCut && fDataType != AliCaloTrackReader::kMC) return ;
2455 
2456  // Check effect of n cells cut
2457  fhPHOSClusterCutsE[5]->Fill(clus->E());
2458 
2459  //----------------------------------------------------
2460  // Apply distance to bad channel cut
2461  //
2462  Double_t distBad = clus->GetDistanceToBadChannel() ; //Distance to bad channel
2463 
2464  if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2465 
2466  if(distBad < fPHOSBadChMinDist) return ;
2467 
2468  // Check effect distance to bad channel cut
2469  fhPHOSClusterCutsE[6]->Fill(clus->E());
2470 
2471  // TODO, add time cut
2472 
2473  //----------------------------------------------------------------------------------
2474  // Add selected clusters to array
2475  //
2476  //if(fDebug > 2 && fMomentum.E() > 0.1)
2477  AliDebug(2,Form("Selected clusters E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2478  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2479 
2480  fPHOSClusters->Add(clus);
2481 
2482  // TODO Dead code? remove?
2483  if (fMixedEvent)
2484  clus->SetID(iclus) ;
2485 
2486  } // esd/aod cluster loop
2487 
2488  AliDebug(1,Form("PHOS selected clusters %d",fPHOSClusters->GetEntriesFast())) ;
2489 }
2490 
2491 //____________________________________________
2493 //____________________________________________
2495 {
2496  if(fEMCALCellsListName.Length() == 0)
2497  fEMCALCells = fInputEvent->GetEMCALCells();
2498  else
2499  fEMCALCells = (AliVCaloCells*) fInputEvent->FindListObject(fEMCALCellsListName);
2500 }
2501 
2502 //___________________________________________
2504 //___________________________________________
2506 {
2507  fPHOSCells = fInputEvent->GetPHOSCells();
2508 }
2509 
2510 //_______________________________________
2513 //_______________________________________
2515 {
2516  AliVVZERO* v0 = fInputEvent->GetVZEROData();
2517  //printf("Init V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
2518 
2519  if (v0)
2520  {
2521  AliESDVZERO* esdV0 = dynamic_cast<AliESDVZERO*> (v0);
2522  for (Int_t i = 0; i < 32; i++)
2523  {
2524  if(esdV0)
2525  {//Only available in ESDs
2526  fV0ADC[0] += (Int_t)esdV0->GetAdcV0C(i);
2527  fV0ADC[1] += (Int_t)esdV0->GetAdcV0A(i);
2528  }
2529 
2530  fV0Mul[0] += (Int_t)v0->GetMultiplicityV0C(i);
2531  fV0Mul[1] += (Int_t)v0->GetMultiplicityV0A(i);
2532  }
2533 
2534  AliDebug(1,Form("ADC (%d,%d), Multiplicity (%d,%d)",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]));
2535  }
2536  else
2537  {
2538  AliDebug(1,"Cannot retrieve V0 ESD! Run w/ null V0 charges");
2539  }
2540 }
2541 
2542 //_________________________________________________
2546 //_________________________________________________
2548 {
2549  AliDebug(2,"Begin");
2550 
2551  //
2552  //check if branch name is given
2553  if(!fInputNonStandardJetBranchName.Length())
2554  {
2555  fInputEvent->Print();
2556  AliFatal("No non-standard jet branch name specified. Specify among existing ones.");
2557  }
2558 
2559  fNonStandardJets = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fInputNonStandardJetBranchName.Data()));
2560 
2561  if(!fNonStandardJets)
2562  {
2563  //check if jet branch exist; exit if not
2564  fInputEvent->Print();
2565 
2566  AliFatal(Form("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fInputNonStandardJetBranchName.Data()));
2567  }
2568  else
2569  {
2570  AliDebug(1,Form("AOD input jets %d", fNonStandardJets->GetEntriesFast()));
2571  }
2572 }
2573 
2574 //_________________________________________________
2578 //_________________________________________________
2580 {
2581  AliDebug(1,"Begin");
2582  //
2583  //check if branch name is given
2584  if(!fInputBackgroundJetBranchName.Length())
2585  {
2586  fInputEvent->Print();
2587 
2588  AliFatal("No background jet branch name specified. Specify among existing ones.");
2589  }
2590 
2591  fBackgroundJets = (AliAODJetEventBackground*)(fInputEvent->FindListObject(fInputBackgroundJetBranchName.Data()));
2592 
2593  if(!fBackgroundJets)
2594  {
2595  //check if jet branch exist; exit if not
2596  fInputEvent->Print();
2597 
2598  AliFatal(Form("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fInputBackgroundJetBranchName.Data()));
2599  }
2600  else
2601  {
2602  AliDebug(1,"FillInputBackgroundJets");
2603  fBackgroundJets->Print("");
2604  }
2605 }
2606 
2607 //____________________________________________________________________
2613 //____________________________________________________________________
2615 {
2616  // init some variables
2617  Int_t trigtimes[30], globCol, globRow,ntimes, i;
2618  Int_t absId = -1; //[100];
2619  Int_t nPatch = 0;
2620 
2621  TArrayI patches(0);
2622 
2623  // get object pointer
2624  AliVCaloTrigger *caloTrigger = GetInputEvent()->GetCaloTrigger( "EMCAL" );
2625 
2626  if(!caloTrigger)
2627  {
2628  AliError("Trigger patches input (AliVCaloTrigger) not available in data!");
2629  return patches;
2630  }
2631 
2632  //printf("CaloTrigger Entries %d\n",caloTrigger->GetEntries() );
2633 
2634  // class is not empty
2635  if( caloTrigger->GetEntries() > 0 )
2636  {
2637  // must reset before usage, or the class will fail
2638  caloTrigger->Reset();
2639 
2640  // go throuth the trigger channels
2641  while( caloTrigger->Next() )
2642  {
2643  // get position in global 2x2 tower coordinates
2644  caloTrigger->GetPosition( globCol, globRow );
2645 
2646  //L0
2647  if(IsEventEMCALL0())
2648  {
2649  // get dimension of time arrays
2650  caloTrigger->GetNL0Times( ntimes );
2651 
2652  // no L0s in this channel
2653  // presence of the channel in the iterator still does not guarantee that L0 was produced!!
2654  if( ntimes < 1 )
2655  continue;
2656 
2657  // get timing array
2658  caloTrigger->GetL0Times( trigtimes );
2659  //printf("Get L0 patch : n times %d - trigger time window %d - %d\n",ntimes, tmin,tmax);
2660 
2661  // go through the array
2662  for( i = 0; i < ntimes; i++ )
2663  {
2664  // check if in cut - 8,9 shall be accepted in 2011
2665  if( trigtimes[i] >= tmin && trigtimes[i] <= tmax )
2666  {
2667  //printf("Accepted trigger time %d \n",trigtimes[i]);
2668  //if(nTrig > 99) continue;
2669  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(globCol,globRow, absId);
2670  //printf("pass the time cut globCol %d, globRow %d absId %d\n",globCol,globRow, absId);
2671  patches.Set(nPatch+1);
2672  patches.AddAt(absId,nPatch++);
2673  }
2674  } // trigger time array
2675  }//L0
2676  else if(IsEventEMCALL1()) // L1
2677  {
2678  Int_t bit = 0;
2679  caloTrigger->GetTriggerBits(bit);
2680 
2681  Int_t sum = 0;
2682  caloTrigger->GetL1TimeSum(sum);
2683  //fBitEGA-=2;
2684  Bool_t isEGA1 = ((bit >> fBitEGA ) & 0x1) && IsEventEMCALL1Gamma1() ;
2685  Bool_t isEGA2 = ((bit >> (fBitEGA+1)) & 0x1) && IsEventEMCALL1Gamma2() ;
2686  Bool_t isEJE1 = ((bit >> fBitEJE ) & 0x1) && IsEventEMCALL1Jet1 () ;
2687  Bool_t isEJE2 = ((bit >> (fBitEJE+1)) & 0x1) && IsEventEMCALL1Jet2 () ;
2688 
2689  //if((bit>> fBitEGA )&0x1) printf("Trig Bit %d - bit %d - EG1 %d - EG2 %d\n",fBitEGA ,bit,IsEventEMCALL1Gamma1(),IsEventEMCALL1Gamma2());
2690  //if((bit>>(fBitEGA+1))&0x1) printf("Trig Bit %d - bit %d - EG1 %d - EG2 %d\n",fBitEGA+1,bit,IsEventEMCALL1Gamma1(),IsEventEMCALL1Gamma2());
2691 
2692  if(!isEGA1 && !isEJE1 && !isEGA2 && !isEJE2) continue;
2693 
2694  Int_t patchsize = -1;
2695  if (isEGA1 || isEGA2) patchsize = 2;
2696  else if (isEJE1 || isEJE2) patchsize = 16;
2697 
2698  //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",
2699  // bit,sum,patchsize,isEGA1,isEGA2,isEJE1,isEJE2,fBitEGA,fBitEJE,IsEventEMCALL1Gamma(),IsEventEMCALL1Jet());
2700 
2701 
2702  // add 2x2 (EGA) or 16x16 (EJE) patches
2703  for(Int_t irow=0; irow < patchsize; irow++)
2704  {
2705  for(Int_t icol=0; icol < patchsize; icol++)
2706  {
2707  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(globCol+icol,globRow+irow, absId);
2708  //printf("pass the time cut globCol %d, globRow %d absId %d\n",globCol,globRow, absId);
2709  patches.Set(nPatch+1);
2710  patches.AddAt(absId,nPatch++);
2711  }
2712  }
2713 
2714  } // L1
2715 
2716  } // trigger iterator
2717  } // go through triggers
2718 
2719  if(patches.GetSize()<=0) AliInfo(Form("No patch found! for triggers: %s and selected <%s>",
2721  //else printf(">>>>> N patches %d, test %d,first %d, last %d\n",patches.GetSize(), nPatch, patches.At(0), patches.At(patches.GetSize()-1));
2722 
2723  return patches;
2724 }
2725 
2726 //____________________________________________________________
2730 //____________________________________________________________
2732 {
2733  // Init info from previous event
2734  fTriggerClusterIndex = -1;
2735  fTriggerClusterId = -1;
2736  fTriggerClusterBC = -10000;
2737  fIsExoticEvent = kFALSE;
2738  fIsBadCellEvent = kFALSE;
2739  fIsBadMaxCellEvent = kFALSE;
2740 
2741  fIsTriggerMatch = kFALSE;
2742  fIsTriggerMatchOpenCut[0] = kFALSE;
2743  fIsTriggerMatchOpenCut[1] = kFALSE;
2744  fIsTriggerMatchOpenCut[2] = kFALSE;
2745 
2746  // Do only analysis for triggered events
2747  if(!IsEventEMCALL1() && !IsEventEMCALL0())
2748  {
2749  fTriggerClusterBC = 0;
2750  return;
2751  }
2752 
2753  //printf("***** Try to match trigger to cluster %d **** L0 %d, L1 %d\n",fTriggerPatchClusterMatch,IsEventEMCALL0(),IsEventEMCALL1());
2754 
2755  //Recover the list of clusters
2756  TClonesArray * clusterList = 0;
2757  if (fInputEvent->FindListObject(fEMCALClustersListName))
2758  {
2759  clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
2760  }
2761  else if(fOutputEvent)
2762  {
2763  clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
2764  }
2765 
2766  // Get number of clusters and of trigger patches
2767  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2768  if(clusterList)
2769  nclusters = clusterList->GetEntriesFast();
2770 
2771  Int_t nPatch = patches.GetSize();
2773 
2774  //Init some variables used in the cluster loop
2775  Float_t tofPatchMax = 100000;
2776  Float_t ePatchMax =-1;
2777 
2778  Float_t tofMax = 100000;
2779  Float_t eMax =-1;
2780 
2781  Int_t clusMax =-1;
2782  Int_t idclusMax =-1;
2783  Bool_t badClMax = kFALSE;
2784  Bool_t badCeMax = kFALSE;
2785  Bool_t exoMax = kFALSE;
2786 // Int_t absIdMaxTrig= -1;
2787  Int_t absIdMaxMax = -1;
2788 
2789  Int_t nOfHighECl = 0 ;
2790 
2791  //
2792  // Check what is the trigger threshold
2793  // set minimu energym of candidate for trigger cluster
2794  //
2796 
2797  Float_t triggerThreshold = fTriggerL1EventThreshold;
2798  if(IsEventEMCALL0()) triggerThreshold = fTriggerL0EventThreshold;
2799  Float_t minE = triggerThreshold / 2.;
2800 
2801  // This method is not really suitable for JET trigger
2802  // but in case, reduce the energy cut since we do not trigger on high energy particle
2803  if(IsEventEMCALL1Jet() || minE < 1) minE = 1;
2804 
2805  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));
2806 
2807  //
2808  // Loop on the clusters, check if there is any that falls into one of the patches
2809  //
2810  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2811  {
2812  AliVCluster * clus = 0;
2813  if(clusterList) clus = (AliVCluster*) clusterList->At(iclus);
2814  else clus = fInputEvent->GetCaloCluster(iclus);
2815 
2816  if ( !clus ) continue ;
2817 
2818  if ( !clus->IsEMCAL() ) continue ;
2819 
2820  //Skip clusters with too low energy to be triggering
2821  if ( clus->E() < minE ) continue ;
2822 
2823  Float_t frac = -1;
2824  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fInputEvent->GetEMCALCells(), clus,frac);
2825 
2826  Bool_t badCluster = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),
2827  clus->GetCellsAbsId(),clus->GetNCells());
2828  UShort_t cellMax[] = {(UShort_t) absIdMax};
2829  Bool_t badCell = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),cellMax,1);
2830 
2831  // if cell is bad, it can happen that time calibration is not available,
2832  // when calculating if it is exotic, this can make it to be exotic by default
2833  // open it temporarily for this cluster
2834  if(badCell)
2836 
2837  Bool_t exotic = GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCluster(clus, fInputEvent->GetEMCALCells());
2838 
2839  // Set back the time cut on exotics
2840  if(badCell)
2842 
2843  // Energy threshold for exotic Ecross typically at 4 GeV,
2844  // for lower energy, check that there are more than 1 cell in the cluster
2845  if(!exotic && clus->GetNCells() < 2) exotic = kTRUE;
2846 
2847  Float_t energy = clus->E();
2848  Int_t idclus = clus->GetID();
2849 
2850  Double_t tof = clus->GetTOF();
2851  if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn() && fTriggerClusterTimeRecal){
2852  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2853  //additional L1 phase shift
2855  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absIdMax), fInputEvent->GetBunchCrossNumber(), tof);
2856  }
2857  }
2858  tof *=1.e9;
2859 
2860  //printf("cluster %d, ID %d, E %2.2f, tof %2.2f, AbsId max %d, exotic %d, bad Cluster %d, bad Cell %d\n",
2861  // iclus,idclus, energy,tof,absIdMax, exotic, badCluster,badCell);
2862 
2863  // Find the highest energy cluster, avobe trigger threshold
2864  // in the event in case no match to trigger is found
2865  if( energy > eMax )
2866  {
2867  tofMax = tof;
2868  eMax = energy;
2869  badClMax = badCluster;
2870  badCeMax = badCell;
2871  exoMax = exotic;
2872  clusMax = iclus;
2873  idclusMax = idclus;
2874  absIdMaxMax = absIdMax;
2875  }
2876 
2877  // count the good clusters in the event avobe the trigger threshold
2878  // to check the exotic events
2879  if(!badCluster && !exotic)
2880  nOfHighECl++;
2881 
2882  // Find match to trigger
2883  if(fTriggerPatchClusterMatch && nPatch>0)
2884  {
2885  for(Int_t iabsId =0; iabsId < nPatch; iabsId++)
2886  {
2887  Int_t absIDCell[4];
2888  GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patches.At(iabsId), absIDCell);
2889  //if(tof > 75 ) printf("E %2.2f TOF %2.2f Trigger patch %d, cells : %d, %d, %d, %d\n",
2890  // clus->E(),tof,patches.At(iabsId), absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
2891 
2892  for(Int_t ipatch = 0; ipatch < 4; ipatch++)
2893  {
2894  if(absIdMax == absIDCell[ipatch])
2895  {
2896  //printf("*** Patches : absId %d, E %2.1f, tof %f \n",absIdMax,clus->E(), tof);
2897  if(energy > ePatchMax)
2898  {
2899  tofPatchMax = tof;
2900  ePatchMax = energy;
2901  fIsBadCellEvent = badCluster;
2902  fIsBadMaxCellEvent = badCell;
2903  fIsExoticEvent = exotic;
2904  fTriggerClusterIndex = iclus;
2905  fTriggerClusterId = idclus;
2906  fIsTriggerMatch = kTRUE;
2907 // absIdMaxTrig = absIdMax;
2908  }
2909  }
2910  }// cell patch loop
2911  }// trigger patch loop
2912  } // Do trigger patch matching
2913 
2914  }// Cluster loop
2915 
2916  // If there was no match, assign as trigger
2917  // the highest energy cluster in the event
2918  if(!fIsTriggerMatch)
2919  {
2920  tofPatchMax = tofMax;
2921  ePatchMax = eMax;
2922  fIsBadCellEvent = badClMax;
2923  fIsBadMaxCellEvent = badCeMax;
2924  fIsExoticEvent = exoMax;
2925  fTriggerClusterIndex = clusMax;
2926  fTriggerClusterId = idclusMax;
2927  }
2928 
2929  Double_t tofPatchMaxUS = TMath::Abs(tofPatchMax);
2930 
2931  if (tofPatchMaxUS < 28 ) fTriggerClusterBC = 0 ;
2932  else if(tofPatchMaxUS < 75 ) fTriggerClusterBC = 1 ;
2933  else if(tofPatchMaxUS < 125) fTriggerClusterBC = 2 ;
2934  else if(tofPatchMaxUS < 175) fTriggerClusterBC = 3 ;
2935  else if(tofPatchMaxUS < 225) fTriggerClusterBC = 4 ;
2936  else if(tofPatchMaxUS < 275) fTriggerClusterBC = 5 ;
2937  else
2938  {
2939  //printf("AliCaloTrackReader::MatchTriggerCluster() - Large BC - tof %2.3f - Index %d\n",tofPatchMaxUS,fTriggerClusterIndex);
2941  else
2942  {
2943  fTriggerClusterIndex = -2;
2944  fTriggerClusterId = -2;
2945  }
2946  }
2947 
2948  if(tofPatchMax < 0) fTriggerClusterBC*=-1;
2949 
2950 
2951  // 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",
2952  // fTriggerClusterIndex, fTriggerClusterId,ePatchMax, tofPatchMax,
2953  // fTriggerClusterBC, fIsBadCellEvent,fIsBadMaxCellEvent,fIsExoticEvent, fIsTriggerMatch, nOfHighECl,absIdMaxMax);
2954  //
2955  // 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",
2956  // clusMax, idclusMax, eMax,tofMax, badClMax, badCeMax,exoMax);
2957 
2958  //Redo matching but open cuts
2959  if(!fIsTriggerMatch && fTriggerClusterId >= 0)
2960  {
2961  // Open time patch time
2962  TArrayI patchOpen = GetTriggerPatches(7,10);
2963 
2964  Int_t patchAbsIdOpenTime = -1;
2965  for(Int_t iabsId =0; iabsId < patchOpen.GetSize(); iabsId++)
2966  {
2967  Int_t absIDCell[4];
2968  patchAbsIdOpenTime = patchOpen.At(iabsId);
2969  GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patchAbsIdOpenTime, absIDCell);
2970  //if(tof > 75 ) printf("E %2.2f TOF %2.2f Trigger patch %d, cells : %d, %d, %d, %d\n",
2971  // clus->E(),tof,patches.At(iabsId), absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
2972 
2973  for(Int_t ipatch = 0; ipatch < 4; ipatch++)
2974  {
2975  if(absIdMaxMax == absIDCell[ipatch])
2976  {
2977  fIsTriggerMatchOpenCut[0] = kTRUE;
2978  break;
2979  }
2980  }// cell patch loop
2981  }// trigger patch loop
2982 
2983  // Check neighbour patches
2984  Int_t patchAbsId = -1;
2985  Int_t globalCol = -1;
2986  Int_t globalRow = -1;
2987  GetCaloUtils()->GetEMCALGeometry()->GetFastORIndexFromCellIndex(absIdMaxMax, patchAbsId);
2988  GetCaloUtils()->GetEMCALGeometry()->GetPositionInEMCALFromAbsFastORIndex(patchAbsId,globalCol,globalRow);
2989 
2990  // Check patches with strict time cut
2991  Int_t patchAbsIdNeigh = -1;
2992  for(Int_t icol = globalCol-1; icol <= globalCol+1; icol++)
2993  {
2994  if(icol < 0 || icol > 47) continue;
2995 
2996  for(Int_t irow = globalRow; irow <= globalRow+1; irow++)
2997  {
2998  if(irow < 0 || irow > 63) continue;
2999 
3000  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(icol, irow, patchAbsIdNeigh);
3001 
3002  if ( patchAbsIdNeigh < 0 ) continue;
3003 
3004  for(Int_t iabsId =0; iabsId < patches.GetSize(); iabsId++)
3005  {
3006  if(patchAbsIdNeigh == patches.At(iabsId))
3007  {
3008  fIsTriggerMatchOpenCut[1] = kTRUE;
3009  break;
3010  }
3011  }// trigger patch loop
3012 
3013  }// row
3014  }// col
3015 
3016  // Check patches with open time cut
3017  Int_t patchAbsIdNeighOpenTime = -1;
3018  for(Int_t icol = globalCol-1; icol <= globalCol+1; icol++)
3019  {
3020  if(icol < 0 || icol > 47) continue;
3021 
3022  for(Int_t irow = globalRow; irow <= globalRow+1; irow++)
3023  {
3024  if(irow < 0 || irow > 63) continue;
3025 
3026  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(icol, irow, patchAbsIdNeighOpenTime);
3027 
3028  if ( patchAbsIdNeighOpenTime < 0 ) continue;
3029 
3030  for(Int_t iabsId =0; iabsId < patchOpen.GetSize(); iabsId++)
3031  {
3032  if(patchAbsIdNeighOpenTime == patchOpen.At(iabsId))
3033  {
3034  fIsTriggerMatchOpenCut[2] = kTRUE;
3035  break;
3036  }
3037  }// trigger patch loop
3038 
3039  }// row
3040  }// col
3041 
3042  // printf("No match, new match: Open time %d-%d, open Neigh %d-%d, both open %d-%d\n",fIsTriggerMatchOpenCut[0],patchAbsIdOpenTime,
3043  // fIsTriggerMatchOpenCut[1],patchAbsIdNeigh,
3044  // fIsTriggerMatchOpenCut[2],patchAbsIdNeighOpenTime);
3045 
3046  patchOpen.Reset();
3047 
3048  }// No trigger match found
3049  //printf("Trigger BC %d, Id %d, Index %d\n",fTriggerClusterBC,fTriggerClusterId,fTriggerClusterIndex);
3050 }
3051 
3052 //_________________________________________________________
3056 //_________________________________________________________
3058 {
3060  {
3061  // get object pointer
3062  AliVCaloTrigger *caloTrigger = GetInputEvent()->GetCaloTrigger( "EMCAL" );
3063 
3064  if ( fBitEGA == 6 )
3065  {
3066  if (IsEventEMCALL1Gamma1()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(1);
3067  else if(IsEventEMCALL1Gamma2()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(3);
3068  else if(IsEventEMCALL1Jet1 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(0);
3069  else if(IsEventEMCALL1Jet2 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(2);
3070 
3071  // printf("L1 trigger Threshold Jet1 %f, Gamma1 %f, Jet2 %f, Gamma2 %f\n",
3072  // 0.07874*caloTrigger->GetL1Threshold(0),
3073  // 0.07874*caloTrigger->GetL1Threshold(1),
3074  // 0.07874*caloTrigger->GetL1Threshold(2),
3075  // 0.07874*caloTrigger->GetL1Threshold(3));
3076  }
3077  else
3078  {
3079  // Old AOD data format, in such case, order of thresholds different!!!
3080  if (IsEventEMCALL1Gamma1()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(0);
3081  else if(IsEventEMCALL1Gamma2()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(2);
3082  else if(IsEventEMCALL1Jet1 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(1);
3083  else if(IsEventEMCALL1Jet2 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(3);
3084 
3085  // printf("L1 trigger Threshold Jet1 %f, Gamma1 %f, Jet2 %f, Gamma2 %f\n",
3086  // 0.07874*caloTrigger->GetL1Threshold(1),
3087  // 0.07874*caloTrigger->GetL1Threshold(0),
3088  // 0.07874*caloTrigger->GetL1Threshold(3),
3089  // 0.07874*caloTrigger->GetL1Threshold(2));
3090  }
3091  }
3092 
3093  // Set L0 threshold, if not set by user
3095  {
3096  // Revise for periods > LHC11d
3097  Int_t runNumber = fInputEvent->GetRunNumber();
3098  if (runNumber < 146861) fTriggerL0EventThreshold = 3. ; // LHC11a
3099  else if(runNumber < 154000) fTriggerL0EventThreshold = 4. ; // LHC11b,c
3100  else if(runNumber < 165000) fTriggerL0EventThreshold = 5.5; // LHC11c,d,e
3101  else if(runNumber < 194000) fTriggerL0EventThreshold = 2 ; // LHC12
3102  else if(runNumber < 197400) fTriggerL0EventThreshold = 3 ; // LHC13def
3103  else if(runNumber < 197400) fTriggerL0EventThreshold = 2 ; // LHC13g
3104  else if(runNumber < 244300) fTriggerL0EventThreshold = 5 ; // LHC15 in, phys 1, 5 in phys2
3105  else if(runNumber < 266400) fTriggerL0EventThreshold = 2.5; // LHC16ir
3106  else fTriggerL0EventThreshold = 3.5; // LHC16s
3107  }
3108 }
3109 
3110 //________________________________________________________
3112 //________________________________________________________
3114 {
3115  if(! opt)
3116  return;
3117 
3118  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
3119  printf("Task name : %s\n", fTaskName.Data()) ;
3120  printf("Data type : %d\n", fDataType) ;
3121  printf("CTS Min pT : %2.1f GeV/c\n", fCTSPtMin) ;
3122  printf("EMCAL Min pT : %2.1f GeV/c\n", fEMCALPtMin) ;
3123  printf("PHOS Min pT : %2.1f GeV/c\n", fPHOSPtMin) ;
3124  printf("CTS Max pT : %2.1f GeV/c\n", fCTSPtMax) ;
3125  printf("EMCAL Max pT : %2.1f GeV/c\n", fEMCALPtMax) ;
3126  printf("PHOS Max pT : %2.1f GeV/c\n", fPHOSPtMax) ;
3127  printf("EMCAL Bad Dist > %2.1f \n" , fEMCALBadChMinDist) ;
3128  printf("PHOS Bad Dist > %2.1f \n" , fPHOSBadChMinDist) ;
3129  printf("EMCAL N cells > %d \n" , fEMCALNCellsCut) ;
3130  printf("PHOS N cells > %d \n" , fPHOSNCellsCut) ;
3131  printf("EMCAL Time Cut: %3.1f < TOF < %3.1f\n", fEMCALTimeCutMin, fEMCALTimeCutMax);
3132  printf("Use CTS = %d\n", fFillCTS) ;
3133  printf("Use EMCAL = %d\n", fFillEMCAL) ;
3134  printf("Use DCAL = %d\n", fFillDCAL) ;
3135  printf("Use PHOS = %d\n", fFillPHOS) ;
3136  printf("Use EMCAL Cells = %d\n", fFillEMCALCells) ;
3137  printf("Use PHOS Cells = %d\n", fFillPHOSCells) ;
3138  printf("Track status = %d\n", (Int_t) fTrackStatus) ;
3139 
3140  printf("Track Mult Eta Cut = %2.2f\n", fTrackMultEtaCut) ;
3141 
3142  printf("Track Mult Pt Cuts:") ;
3143  for(Int_t icut = 0; icut < fTrackMultNPtCut; icut++) printf(" %2.2f GeV;",fTrackMultPtCut[icut]);
3144  printf(" \n") ;
3145 
3146  printf("Write delta AOD = %d\n", fWriteOutputDeltaAOD) ;
3147  printf("Recalculate Clusters = %d, E linearity = %d\n", fRecalculateClusters, fCorrectELinearity) ;
3148 
3149  printf("Use Triggers selected in SE base class %d; If not what Trigger Mask? %d; MB Trigger Mask for mixed %d \n",
3151 
3153  printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
3154 
3156  printf("Compare cluster pt and pt hard to accept event, factor = %2.2f",fPtHardAndClusterPtFactor);
3157 
3158  printf("Delta AOD File Name = %s\n", fDeltaAODFileName.Data()) ;
3159  printf("Centrality: Class %s, Option %d, Bin [%d,%d] \n", fCentralityClass.Data(),fCentralityOpt,fCentralityBin[0], fCentralityBin[1]) ;
3160 
3161  printf(" \n") ;
3162 }
3163 
3164 //__________________________________________
3170 //__________________________________________
3172 {
3173  // For LHC11a
3174  // Count number of cells with energy larger than 0.1 in SM3, cut on this number
3175  if ( fRemoveLEDEvents == 1 )
3176  {
3177  Int_t ncellsSM3 = 0;
3178  for(Int_t icell = 0; icell < fInputEvent->GetEMCALCells()->GetNumberOfCells(); icell++)
3179  {
3180  Int_t absID = fInputEvent->GetEMCALCells()->GetCellNumber(icell);
3181  Int_t sm = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
3182  if(fInputEvent->GetEMCALCells()->GetAmplitude(icell) > 0.1 && sm==3) ncellsSM3++;
3183  }
3184 
3185  Int_t ncellcut = 21;
3186  if(GetFiredTriggerClasses().Contains("EMC")) ncellcut = 35;
3187 
3188  if(ncellsSM3 >= ncellcut)
3189  {
3190  AliDebug(1,Form("Reject event with ncells in SM3 %d, cut %d, trig %s",
3191  ncellsSM3,ncellcut,GetFiredTriggerClasses().Data()));
3192  return kTRUE;
3193  }
3194  }
3195  // For testing
3196  // Count number of cells with energy larger than 0.1 any SM, cut on this number
3197  else if ( fRemoveLEDEvents > 1 )
3198  {
3199  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};
3200 
3201  for(Int_t icell = 0; icell < fInputEvent->GetEMCALCells()->GetNumberOfCells(); icell++)
3202  {
3203  Int_t absID = fInputEvent->GetEMCALCells()->GetCellNumber(icell);
3204  Int_t sm = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
3205  if(fInputEvent->GetEMCALCells()->GetAmplitude(icell) > 0.1) ncellsSM[sm]++;
3206  }
3207 
3208  Int_t ncellcut = 21;
3209  if(GetFiredTriggerClasses().Contains("EMC")) ncellcut = 35;
3210 
3211  for(Int_t ism = 0; ism < GetCaloUtils()->GetEMCALGeometry()->GetNumberOfSuperModules(); ism++)
3212  {
3213  if(ncellsSM[ism] >= ncellcut)
3214  {
3215  AliDebug(1,Form("Reject event with ncells in SM%d %d, cut %d, trig %s",
3216  ism,ncellsSM[ism],ncellcut,GetFiredTriggerClasses().Data()));
3217 
3218  return kTRUE;
3219  }
3220  }
3221  }
3222 
3223  return kFALSE;
3224 }
3225 
3226 //_________________________________________________________
3230 //_________________________________________________________
3232 {
3233  if(label < 0) return ;
3234 
3235  AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fInputEvent) ;
3236  if(!evt) return ;
3237 
3238  TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
3239  if(!arr) return ;
3240 
3241  if(label < arr->GetEntriesFast())
3242  {
3243  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
3244  if(!particle) return ;
3245 
3246  if(label == particle->Label()) return ; // label already OK
3247  //else printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
3248  }
3249  //else printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
3250 
3251  // loop on the particles list and check if there is one with the same label
3252  for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
3253  {
3254  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
3255  if(!particle) continue ;
3256 
3257  if(label == particle->Label())
3258  {
3259  label = ind;
3260  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
3261  return;
3262  }
3263  }
3264 
3265  label = -1;
3266 
3267  //printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label not found set to -1 \n");
3268 }
3269 
3270 //___________________________________
3272 //___________________________________
3274 {
3275  if(fCTSTracks) fCTSTracks -> Clear();
3276  if(fEMCALClusters) fEMCALClusters -> Clear("C");
3277  if(fPHOSClusters) fPHOSClusters -> Clear("C");
3278 
3279  fV0ADC[0] = 0; fV0ADC[1] = 0;
3280  fV0Mul[0] = 0; fV0Mul[1] = 0;
3281 
3282  if(fNonStandardJets) fNonStandardJets -> Clear("C");
3283  fBackgroundJets->Reset();
3284 }
3285 
3286 //___________________________________________
3290 //___________________________________________
3292 {
3293  fEventTrigMinBias = kFALSE;
3294  fEventTrigCentral = kFALSE;
3295  fEventTrigSemiCentral = kFALSE;
3296  fEventTrigEMCALL0 = kFALSE;
3297  fEventTrigEMCALL1Gamma1 = kFALSE;
3298  fEventTrigEMCALL1Gamma2 = kFALSE;
3299  fEventTrigEMCALL1Jet1 = kFALSE;
3300  fEventTrigEMCALL1Jet2 = kFALSE;
3301 
3302  AliDebug(1,Form("Select trigger mask bit %d - Trigger Event %s - Select <%s>",
3304 
3305  if(fEventTriggerMask <=0 )// in case no mask set
3306  {
3307  // EMC triggered event? Which type?
3308  if( GetFiredTriggerClasses().Contains("-B-") || GetFiredTriggerClasses().Contains("-S-") || GetFiredTriggerClasses().Contains("-I-") )
3309  {
3310  if ( GetFiredTriggerClasses().Contains("EGA" ) ||
3311  GetFiredTriggerClasses().Contains("EG1" ) )
3312  {
3313  fEventTrigEMCALL1Gamma1 = kTRUE;
3314  if( GetFiredTriggerClasses().Contains("EG1" ) && !fFiredTriggerClassName.Contains("EG1") ) fEventTrigEMCALL1Gamma1 = kFALSE;
3315  }
3316  else if( GetFiredTriggerClasses().Contains("EG2" ) )
3317  {
3318  fEventTrigEMCALL1Gamma2 = kTRUE;
3319  if( !fFiredTriggerClassName.Contains("EG2") ) fEventTrigEMCALL1Gamma2 = kFALSE;
3320  }
3321  else if( GetFiredTriggerClasses().Contains("EJE" ) ||
3322  GetFiredTriggerClasses().Contains("EJ1" ) )
3323  {
3324  fEventTrigEMCALL1Jet1 = kTRUE;
3325  if( GetFiredTriggerClasses().Contains("EJ1" ) && !fFiredTriggerClassName.Contains("EJ1") )
3326  fEventTrigEMCALL1Jet1 = kFALSE;
3327  }
3328  else if( GetFiredTriggerClasses().Contains("EJ2" ) )
3329  {
3330  fEventTrigEMCALL1Jet2 = kTRUE;
3331  if( !fFiredTriggerClassName.Contains("EJ2") ) fEventTrigEMCALL1Jet2 = kFALSE;
3332  }
3333  else if( GetFiredTriggerClasses().Contains("CEMC") &&
3334  !GetFiredTriggerClasses().Contains("EGA" ) &&
3335  !GetFiredTriggerClasses().Contains("EJE" ) &&
3336  !GetFiredTriggerClasses().Contains("EG1" ) &&
3337  !GetFiredTriggerClasses().Contains("EJ1" ) &&
3338  !GetFiredTriggerClasses().Contains("EG2" ) &&
3339  !GetFiredTriggerClasses().Contains("EJ2" ) )
3340  {
3341  fEventTrigEMCALL0 = kTRUE;
3342  }
3343 
3344  //Min bias event trigger?
3345  if (GetFiredTriggerClasses().Contains("CCENT_R2-B-NOPF-ALLNOTRD"))
3346  {
3347  fEventTrigCentral = kTRUE;
3348  }
3349  else if(GetFiredTriggerClasses().Contains("CSEMI_R1-B-NOPF-ALLNOTRD"))
3350  {
3351  fEventTrigSemiCentral = kTRUE;
3352  }
3353  else if((GetFiredTriggerClasses().Contains("CINT") || GetFiredTriggerClasses().Contains("CPBI2_B1") ) &&
3354  GetFiredTriggerClasses().Contains("-NOPF-ALLNOTRD") )
3355  {
3356  fEventTrigMinBias = kTRUE;
3357  }
3358  }
3359  }
3360  else
3361  {
3362  // EMC L1 Gamma
3363  if ( fEventTriggerMask & AliVEvent::kEMCEGA )
3364  {
3365  //printf("EGA trigger bit\n");
3366  if (GetFiredTriggerClasses().Contains("EG"))
3367  {
3368  if (GetFiredTriggerClasses().Contains("EGA")) fEventTrigEMCALL1Gamma1 = kTRUE;
3369  else
3370  {
3371  if(GetFiredTriggerClasses().Contains("EG1")) fEventTrigEMCALL1Gamma1 = kTRUE;
3372  if(GetFiredTriggerClasses().Contains("EG2")) fEventTrigEMCALL1Gamma2 = kTRUE;
3373  }
3374  }
3375  }
3376  // EMC L1 Jet
3377  else if( fEventTriggerMask & AliVEvent::kEMCEJE )
3378  {
3379  //printf("EGA trigger bit\n");
3380  if (GetFiredTriggerClasses().Contains("EJ"))
3381  {
3382  if (GetFiredTriggerClasses().Contains("EJE")) fEventTrigEMCALL1Jet1 = kTRUE;
3383  else
3384  {
3385  if(GetFiredTriggerClasses().Contains("EJ1")) fEventTrigEMCALL1Jet1 = kTRUE;
3386  if(GetFiredTriggerClasses().Contains("EJ2")) fEventTrigEMCALL1Jet2 = kTRUE;
3387  }
3388  }
3389  }
3390  // EMC L0
3391  else if((fEventTriggerMask & AliVEvent::kEMC7) ||
3392  (fEventTriggerMask & AliVEvent::kEMC1) )
3393  {
3394  //printf("L0 trigger bit\n");
3395  fEventTrigEMCALL0 = kTRUE;
3396  }
3397  // Min Bias Pb-Pb
3398  else if( fEventTriggerMask & AliVEvent::kCentral )
3399  {
3400  //printf("MB semi central trigger bit\n");
3401  fEventTrigSemiCentral = kTRUE;
3402  }
3403  // Min Bias Pb-Pb
3404  else if( fEventTriggerMask & AliVEvent::kSemiCentral )
3405  {
3406  //printf("MB central trigger bit\n");
3407  fEventTrigCentral = kTRUE;
3408  }
3409  // Min Bias pp, PbPb, pPb
3410  else if((fEventTriggerMask & AliVEvent::kMB ) ||
3411  (fEventTriggerMask & AliVEvent::kINT7) ||
3412  (fEventTriggerMask & AliVEvent::kINT8) ||
3413  (fEventTriggerMask & AliVEvent::kAnyINT) )
3414  {
3415  //printf("MB trigger bit\n");
3416  fEventTrigMinBias = kTRUE;
3417  }
3418  }
3419 
3420  AliDebug(1,Form("Event bits: \n \t MB %d, Cen %d, Sem %d, L0 %d, L1G1 %d, L1G2 %d, L1J1 %d, L1J2 %d",
3424 
3425  // L1 trigger bit
3426  if( fBitEGA == 0 && fBitEJE == 0 )
3427  {
3428  // Init the trigger bit once, correct depending on AliESD(AOD)CaloTrigger header version
3429 
3430  // Simpler way to do it ...
3431 // if( fInputEvent->GetRunNumber() < 172000 )
3432 // reader->SetEventTriggerL1Bit(4,5); // current LHC11 data
3433 // else
3434 // reader->SetEventTriggerL1Bit(6,8); // LHC12-13 data
3435 
3436  // Old values
3437  fBitEGA = 4;
3438  fBitEJE = 5;
3439 
3440  TFile* file = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
3441 
3442  const TList *clist = file->GetStreamerInfoCache();
3443 
3444  if(clist)
3445  {
3446  TStreamerInfo *cinfo = (TStreamerInfo*)clist->FindObject("AliESDCaloTrigger");
3447  Int_t verid = 5; // newer ESD header version
3448  if(!cinfo)
3449  {
3450  cinfo = (TStreamerInfo*)clist->FindObject("AliAODCaloTrigger");
3451  verid = 3; // newer AOD header version
3452  }
3453 
3454  if(cinfo)
3455  {
3456  Int_t classversionid = cinfo->GetClassVersion();
3457  //printf("********* Header class version %d *********** \n",classversionid);
3458 
3459  if (classversionid >= verid)
3460  {
3461  fBitEGA = 6;
3462  fBitEJE = 8;
3463  }
3464  } else AliInfo("AliCaloTrackReader()::SetEventTriggerBit() - Streamer info for trigger class not available, bit not changed");
3465  } else AliInfo("AliCaloTrackReader::SetEventTriggerBit() - Streamer list not available!, bit not changed");
3466 
3467  } // set once the EJE, EGA trigger bit
3468 }
3469 
3470 //____________________________________________________________
3473 //____________________________________________________________
3474 void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
3475 {
3476  fInputEvent = input;
3477  fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ;
3478  if (fMixedEvent)
3479  fNMixedEvent = fMixedEvent->GetNumberOfEvents() ;
3480 
3481  //Delete previous vertex
3482  if(fVertex)
3483  {
3484  for (Int_t i = 0; i < fNMixedEvent; i++)
3485  {
3486  delete [] fVertex[i] ;
3487  }
3488  delete [] fVertex ;
3489  }
3490 
3491  fVertex = new Double_t*[fNMixedEvent] ;
3492  for (Int_t i = 0; i < fNMixedEvent; i++)
3493  {
3494  fVertex[i] = new Double_t[3] ;
3495  fVertex[i][0] = 0.0 ;
3496  fVertex[i][1] = 0.0 ;
3497  fVertex[i][2] = 0.0 ;
3498  }
3499 }
3500 
3501 
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