AliPhysics  a9863a5 (a9863a5)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliCaloTrackReader.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // --- ROOT system ---
17 #include <TFile.h>
18 #include <TGeoManager.h>
19 #include <TStreamerInfo.h>
20 
21 // ---- ANALYSIS system ----
22 #include "AliMCEvent.h"
23 #include "AliAODMCHeader.h"
24 #include "AliGenPythiaEventHeader.h"
25 #include "AliGenCocktailEventHeader.h"
26 #include "AliGenHijingEventHeader.h"
27 #include "AliESDEvent.h"
28 #include "AliAODEvent.h"
29 #include "AliVTrack.h"
30 #include "AliVParticle.h"
31 #include "AliMixedEvent.h"
32 //#include "AliTriggerAnalysis.h"
33 #include "AliESDVZERO.h"
34 #include "AliVCaloCells.h"
35 #include "AliAnalysisManager.h"
36 #include "AliInputEventHandler.h"
37 #include "AliAODMCParticle.h"
38 #include "AliStack.h"
39 #include "AliLog.h"
40 #include "AliMultSelection.h"
41 
42 // ---- Detectors ----
43 #include "AliPHOSGeoUtils.h"
44 #include "AliEMCALGeometry.h"
45 #include "AliEMCALRecoUtils.h"
46 
47 // ---- CaloTrackCorr ---
48 #include "AliCalorimeterUtils.h"
49 #include "AliCaloTrackReader.h"
50 
51 // ---- Jets ----
52 #include "AliAODJet.h"
53 #include "AliAODJetEventBackground.h"
54 
58 
59 //________________________________________
61 //________________________________________
63 TObject(), fEventNumber(-1), //fCurrentFileName(""),
64 fDataType(0), fDebug(0),
65 fFiducialCut(0x0), fCheckFidCut(kFALSE),
66 fComparePtHardAndJetPt(0), fPtHardAndJetPtFactor(0),
67 fComparePtHardAndClusterPt(0),fPtHardAndClusterPtFactor(0),
68 fCTSPtMin(0), fEMCALPtMin(0), fPHOSPtMin(0),
69 fCTSPtMax(0), fEMCALPtMax(0), fPHOSPtMax(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 fSmearShowerShape(0), fSmearShowerShapeWidth(0), fRandom(),
87 fSmearingFunction(0),
88 fTrackStatus(0), fSelectSPDHitTracks(0),
89 fTrackMult(0), fTrackMultEtaCut(0.9),
90 fReadStack(kFALSE), fReadAODMCParticles(kFALSE),
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),
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(""), fZvtxCut(0.),
108 fAcceptFastCluster(kFALSE), fRemoveLEDEvents(kFALSE),
109 //Trigger rejection
110 fRemoveBadTriggerEvents(0), fTriggerPatchClusterMatch(1),
111 fTriggerPatchTimeWindow(), fTriggerL0EventThreshold(0),
112 fTriggerL1EventThreshold(0), fTriggerL1EventThresholdFix(0),
113 fTriggerClusterBC(0), fTriggerClusterIndex(0), fTriggerClusterId(0),
114 fIsExoticEvent(0), fIsBadCellEvent(0), fIsBadMaxCellEvent(0),
115 fIsTriggerMatch(0), fIsTriggerMatchOpenCut(),
116 fTriggerClusterTimeRecal(kTRUE), fRemoveUnMatchedTriggers(kTRUE),
117 fDoPileUpEventRejection(kFALSE), fDoV0ANDEventSelection(kFALSE),
118 fDoVertexBCEventSelection(kFALSE),
119 fDoRejectNoTrackEvents(kFALSE),
120 fUseEventsWithPrimaryVertex(kFALSE),
121 //fTriggerAnalysis (0x0),
122 fTimeStampEventSelect(0),
123 fTimeStampEventFracMin(0), fTimeStampEventFracMax(0),
124 fTimeStampRunMin(0), fTimeStampRunMax(0),
125 fNPileUpClusters(-1), fNNonPileUpClusters(-1), fNPileUpClustersCut(3),
126 fVertexBC(-200), fRecalculateVertexBC(0),
127 fUseAliCentrality(0), fCentralityClass(""), fCentralityOpt(0),
128 fEventPlaneMethod(""),
129 fAcceptOnlyHIJINGLabels(0), fNMCProducedMin(0), fNMCProducedMax(0),
130 fFillInputNonStandardJetBranch(kFALSE),
131 fNonStandardJets(new TClonesArray("AliAODJet",100)),fInputNonStandardJetBranchName("jets"),
132 fFillInputBackgroundJetBranch(kFALSE),
133 fBackgroundJets(0x0),fInputBackgroundJetBranchName("jets"),
134 fAcceptEventsWithBit(0), fRejectEventsWithBit(0), fRejectEMCalTriggerEventsWith2Tresholds(0),
135 fMomentum()
136 {
137  InitParameters();
138 }
139 
140 //_______________________________________
142 //_______________________________________
144 {
145  DeletePointers();
146 }
147 
148 //_______________________________________
151 //_______________________________________
153 {
154  delete fFiducialCut ;
155 
156  if(fAODBranchList)
157  {
158  fAODBranchList->Delete();
159  delete fAODBranchList ;
160  }
161 
162  if(fCTSTracks)
163  {
164  if(fDataType!=kMC)fCTSTracks->Clear() ;
165  else fCTSTracks->Delete() ;
166  delete fCTSTracks ;
167  }
168 
169  if(fEMCALClusters)
170  {
171  if(fDataType!=kMC)fEMCALClusters->Clear("C") ;
172  else fEMCALClusters->Delete() ;
173  delete fEMCALClusters ;
174  }
175 
176  if(fDCALClusters)
177  {
178  if(fDataType!=kMC)fDCALClusters->Clear("C") ;
179  else fDCALClusters->Delete() ;
180  delete fDCALClusters ;
181  }
182 
183  if(fPHOSClusters)
184  {
185  if(fDataType!=kMC)fPHOSClusters->Clear("C") ;
186  else fPHOSClusters->Delete() ;
187  delete fPHOSClusters ;
188  }
189 
190  if(fVertex)
191  {
192  for (Int_t i = 0; i < fNMixedEvent; i++)
193  {
194  delete [] fVertex[i] ;
195 
196  }
197  delete [] fVertex ;
198  }
199 
200  //delete fTriggerAnalysis;
201 
202  if(fNonStandardJets)
203  {
204  if(fDataType!=kMC) fNonStandardJets->Clear("C") ;
205  else fNonStandardJets->Delete() ;
206  delete fNonStandardJets ;
207  }
208  delete fBackgroundJets ;
209 
210  fRejectEventsWithBit.Reset();
211  fAcceptEventsWithBit.Reset();
212 
213  if ( fWeightUtils ) delete fWeightUtils ;
214 
215  // Pointers not owned, done by the analysis frame
216  // if(fInputEvent) delete fInputEvent ;
217  // if(fOutputEvent) delete fOutputEvent ;
218  // if(fMC) delete fMC ;
219  // Pointer not owned, deleted by maker
220  // if (fCaloUtils) delete fCaloUtils ;
221 }
222 
223 //____________________________________________________________
227 //____________________________________________________________
228 Bool_t AliCaloTrackReader::AcceptDCA(Float_t pt, Float_t dca)
229 {
230  Float_t cut = fTrackDCACut[0]+fTrackDCACut[1]/TMath::Power(pt,fTrackDCACut[2]);
231 
232  if(TMath::Abs(dca) < cut)
233  return kTRUE;
234  else
235  return kFALSE;
236 }
237 
238 //_____________________________________________________
241 //_____________________________________________________
243 {
244  Int_t nAccept = fAcceptEventsWithBit.GetSize();
245 
246  //printf("N accept %d\n", nAccept);
247 
248  if( nAccept <= 0 )
249  return kTRUE ; // accept the event
250 
251  UInt_t trigFired = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
252 
253  for(Int_t ibit = 0; ibit < nAccept; ibit++)
254  {
255  Bool_t accept = (trigFired & fAcceptEventsWithBit.At(ibit));
256 
257  //printf("accept %d, ibit %d, bit %d \n",accept, ibit,fAcceptEventsWithBit.At(ibit));
258  if(accept) return kTRUE ; // accept the event
259  }
260 
261  return kFALSE ; // reject the event
262 }
263 
264 //_____________________________________________________
267 //_____________________________________________________
269 {
270  Int_t nReject = fRejectEventsWithBit.GetSize();
271 
272  //printf("N reject %d\n", nReject);
273 
274  if( nReject <= 0 )
275  return kTRUE ; // accept the event
276 
277  UInt_t trigFired = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
278 
279  for(Int_t ibit = 0; ibit < nReject; ibit++)
280  {
281  Bool_t reject = (trigFired & fRejectEventsWithBit.At(ibit));
282 
283  //printf("reject %d, ibit %d, bit %d \n",reject, ibit,fRejectEventsWithBit.At(ibit));
284  if(reject) return kFALSE ; // reject the event
285  }
286 
287  return kTRUE ; // accept the event
288 }
289 
290 //_____________________________________________
294 //_____________________________________________
296 {
297  //-----------------------------------------------------------
298  // Reject events depending on the trigger name
299  //-----------------------------------------------------------
300 
301  AliDebug(1,Form("FiredTriggerClass <%s>, selected class <%s>, compare name %d",
304 
305  if ( fFiredTriggerClassName != "" )
306  {
308  return kFALSE;
309  else
310  AliDebug(1,"Accepted triggered event");
311  }
312 
313  //-----------------------------------------------------------
314  // Reject events depending on the event species type
315  //-----------------------------------------------------------
316 
317  // Event types:
318  // kStartOfRun = 1, // START_OF_RUN
319  // kEndOfRun = 2, // END_OF_RUN
320  // kStartOfRunFiles = 3, // START_OF_RUN_FILES
321  // kEndOfRunFiles = 4, // END_OF_RUN_FILES
322  // kStartOfBurst = 5, // START_OF_BURST
323  // kEndOfBurst = 6, // END_OF_BURST
324  // kPhysicsEvent = 7, // PHYSICS_EVENT
325  // kCalibrationEvent = 8, // CALIBRATION_EVENT
326  // kFormatError = 9, // EVENT_FORMAT_ERROR
327  // kStartOfData = 10, // START_OF_DATA
328  // kEndOfData = 11, // END_OF_DATA
329  // kSystemSoftwareTriggerEvent = 12, // SYSTEM_SOFTWARE_TRIGGER_EVENT
330  // kDetectorSoftwareTriggerEvent = 13 // DETECTOR_SOFTWARE_TRIGGER_EVENT
331 
332  Int_t eventType = 0;
333  if(fInputEvent->GetHeader())
334  eventType = ((AliVHeader*)fInputEvent->GetHeader())->GetEventType();
335 
336  AliDebug(3,Form("Event type %d",eventType));
337 
338  // Select only Physics events in data, eventType = 7,
339  // usually done by physics selection
340  // MC not set, eventType =0, I think, so do not apply a value to fEventType by default
341  // LED events have eventType = 8, implemented a selection cut in the past
342  // but not really useful for EMCal data since LED are not reconstructed (unless wrongly assumed as physics)
343 
344  if ( fEventType >= 0 && eventType != fEventType ) return kFALSE ;
345 
346  AliDebug(1,"Pass event species selection");
347 
348  //-----------------------------------------------------------------
349  // In case of mixing analysis, select here the trigger of the event
350  //-----------------------------------------------------------------
351 
352  UInt_t isTrigger = kFALSE;
353  UInt_t isMB = kFALSE;
354 
355  if(!fEventTriggerAtSE)
356  {
357  // In case of mixing analysis, accept MB events, not only Trigger
358  // Track and cluster arrays filled for MB in order to create the pool in the corresponding analysis
359  // via de method in the base class FillMixedEventPool()
360 
361  AliAnalysisManager *manager = AliAnalysisManager::GetAnalysisManager();
362  AliInputEventHandler *inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
363 
364  if(!inputHandler) return kFALSE ; // to content coverity
365 
366  isTrigger = inputHandler->IsEventSelected() & fEventTriggerMask;
367  isMB = inputHandler->IsEventSelected() & fMixEventTriggerMask;
368 
369  if(!isTrigger && !isMB) return kFALSE;
370 
371  //printf("Selected triggered event : %s\n",GetFiredTriggerClasses().Data());
372  AliDebug(0,"Pass uninteresting triggered events rejection in case of mixing analysis");
373  }
374 
375  //-------------------------------------------------------------------------------------
376  // Reject or accept events depending on the trigger bit
377  //-------------------------------------------------------------------------------------
378 
379  Bool_t okA = AcceptEventWithTriggerBit();
380  Bool_t okR = RejectEventWithTriggerBit();
381 
382  //printf("AliCaloTrackReader::FillInputEvent() - Accept event? %d, Reject event %d? \n",okA,okR);
383 
384  if(!okA || !okR) return kFALSE;
385 
386  AliDebug(1,"Pass event bit rejection");
387 
388  //----------------------------------------------------------------------
389  // Do not count events that were likely triggered by an exotic cluster
390  // or out BC cluster
391  //----------------------------------------------------------------------
392 
393  // Set a bit with the event kind, MB, L0, L1 ...
395 
396  // In case of Mixing, avoid checking the triggers in the min bias events
397  if(!fEventTriggerAtSE && (isMB && !isTrigger)) return kTRUE;
398 
400  {
402  {
403  // Reject triggered events when there is coincidence on both EMCal trigger thresholds,
404  // but the requested trigger is the low trigger threshold
405  if(IsEventEMCALL1Jet1 () && IsEventEMCALL1Jet2 () && fFiredTriggerClassName.Contains("EJ2")) return kFALSE;
406  if(IsEventEMCALL1Gamma1() && IsEventEMCALL1Gamma2() && fFiredTriggerClassName.Contains("EG2")) return kFALSE;
407  }
408 
409  //Get Patches that triggered
411 
412  MatchTriggerCluster(patches);
413 
414  patches.Reset();
415 
416  // If requested, remove badly triggeed events, but only when the EMCal trigger bit is set
418  {
419  AliDebug(1,Form("ACCEPT triggered event? \n exotic? %d - bad cell %d - bad Max cell %d - BC %d - Matched %d\n",
421  if (fIsExoticEvent) return kFALSE;
422  else if(fIsBadCellEvent) return kFALSE;
423  else if(fRemoveUnMatchedTriggers && !fIsTriggerMatch) return kFALSE ;
424  else if(fTriggerClusterBC != 0) return kFALSE;
425  AliDebug(1,Form("\t *** YES for %s",GetFiredTriggerClasses().Data()));
426  }
427 
428  AliDebug(1,"Pass EMCal triggered event rejection \n");
429  }
430 
431  //-------------------------------------------------------------------------------------
432  // Select events only fired by a certain trigger configuration if it is provided
433  //-------------------------------------------------------------------------------------
434 
435  if (GetFiredTriggerClasses().Contains("FAST") && !GetFiredTriggerClasses().Contains("ALL") && !fAcceptFastCluster)
436  {
437  AliDebug(1,Form("Do not count events from fast cluster, trigger name %s\n",fFiredTriggerClassName.Data()));
438  return kFALSE;
439  }
440 
441  //-------------------------------------------------------------------------------------
442  // Reject event if large clusters with large energy
443  // Use only for LHC11a data for the moment, and if input is clusterizer V1 or V1+unfolding
444  // If clusterzer NxN or V2 it does not help
445  //-------------------------------------------------------------------------------------
446 
447  //Int_t run = fInputEvent->GetRunNumber();
448 
449  if ( fRemoveLEDEvents )
450  {
451  Bool_t reject = RejectLEDEvents();
452 
453  if(reject) return kFALSE;
454 
455  AliDebug(1,"Pass LED event rejection");
456  } // Remove LED events
457 
458  // All selection criteria passed, accept the event
459  return kTRUE ;
460 }
461 
462 //________________________________________________
467 //________________________________________________
469 {
470  //printf("AliCaloTrackReader::ComparePtHardAndJetPt() - GenHeaderName : %s\n",GetGenEventHeader()->ClassName());
471 
472  if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
473  {
474  TParticle * jet = 0;
475  AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
476  Int_t nTriggerJets = pygeh->NTriggerJets();
477  Float_t ptHard = pygeh->GetPtHard();
478 
479  AliDebug(1,Form("Njets: %d, pT Hard %f",nTriggerJets, ptHard));
480 
481  Float_t tmpjet[]={0,0,0,0};
482  for(Int_t ijet = 0; ijet< nTriggerJets; ijet++)
483  {
484  pygeh->TriggerJet(ijet, tmpjet);
485  jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
486 
487  AliDebug(1,Form("jet %d; pycell jet pT %f",ijet, jet->Pt()));
488 
489  //Compare jet pT and pt Hard
490  if(jet->Pt() > fPtHardAndJetPtFactor * ptHard)
491  {
492  AliInfo(Form("Reject jet event with : pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
493  ptHard, jet->Pt(), fPtHardAndJetPtFactor));
494  return kFALSE;
495  }
496  }
497 
498  if(jet) delete jet;
499  }
500 
501  return kTRUE ;
502 }
503 
504 //____________________________________________________
509 //____________________________________________________
511 {
512  if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
513  {
514  AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
515  Float_t ptHard = pygeh->GetPtHard();
516 
517  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
518  for (Int_t iclus = 0; iclus < nclusters; iclus++)
519  {
520  AliVCluster * clus = fInputEvent->GetCaloCluster(iclus) ;
521  Float_t ecluster = clus->E();
522 
523  if(ecluster > fPtHardAndClusterPtFactor*ptHard)
524  {
525  AliInfo(Form("Reject : ecluster %2.2f, calo %d, factor %2.2f, ptHard %f",ecluster,clus->GetType(),fPtHardAndClusterPtFactor,ptHard));
526 
527  return kFALSE;
528  }
529  }
530 
531  }
532 
533  return kTRUE ;
534 }
535 
536 //____________________________________________
538 //____________________________________________
540 {
541  if(fMC)
542  return fMC->Stack();
543  else
544  {
545  AliDebug(1,"Stack is not available");
546  return 0x0 ;
547  }
548 }
549 
550 //______________________________________________
552 //______________________________________________
554 {
555  if(fMC)
556  {
557  return fMC->Header();
558  }
559  else
560  {
561  AliInfo("Header is not available");
562  return 0x0 ;
563  }
564 }
565 
566 //____________________________________________________
570 //____________________________________________________
572 {
573  fNMCProducedMin = 0;
574  fNMCProducedMax = 0;
575 
576  if (ReadStack() && fMC)
577  {
578  AliGenEventHeader * eventHeader = fMC->GenEventHeader();
579 
580  if(!fAcceptOnlyHIJINGLabels) return ;
581 
582  // TODO Check if it works from here ...
583 
584  AliGenCocktailEventHeader *cocktail = dynamic_cast<AliGenCocktailEventHeader *>(eventHeader);
585 
586  if(!cocktail) return ;
587 
588  TList *genHeaders = cocktail->GetHeaders();
589 
590  Int_t nGenerators = genHeaders->GetEntries();
591  //printf("N generators %d \n", nGenerators);
592 
593  for(Int_t igen = 0; igen < nGenerators; igen++)
594  {
595  AliGenEventHeader * eventHeader2 = (AliGenEventHeader*)genHeaders->At(igen) ;
596  TString name = eventHeader2->GetName();
597 
598  //printf("Generator %d: Class Name %s, Name %s, title %s \n",igen, eventHeader2->ClassName(), name.Data(), eventHeader2->GetTitle());
599 
601  fNMCProducedMax+= eventHeader2->NProduced();
602 
603  if(name.Contains("Hijing",TString::kIgnoreCase)) return ;
604  }
605 
606  }
607  else if(ReadAODMCParticles() && GetAODMCHeader())
608  {
609  Int_t nGenerators = GetAODMCHeader()->GetNCocktailHeaders();
610  //printf("AliCaloTrackReader::GetGenEventHeader() - N headers %d\n",nGenerators);
611 
612  if( nGenerators <= 0) return ;
613 
614  if(!fAcceptOnlyHIJINGLabels) return ;
615 
616  for(Int_t igen = 0; igen < nGenerators; igen++)
617  {
618  AliGenEventHeader * eventHeader = GetAODMCHeader()->GetCocktailHeader(igen) ;
619  TString name = eventHeader->GetName();
620 
621  //printf("Generator %d: Class Name %s, Name %s, title %s \n",igen, eventHeader->ClassName(), name.Data(), eventHeader->GetTitle());
622 
624  fNMCProducedMax+= eventHeader->NProduced();
625 
626  if(name.Contains("Hijing",TString::kIgnoreCase)) return ;
627  }
628 
629  }
630 }
631 
632 //______________________________________________________________
635 //______________________________________________________________
636 AliGenEventHeader* AliCaloTrackReader::GetGenEventHeader() const
637 {
638  if (ReadStack() && fMC)
639  {
640  AliGenEventHeader * eventHeader = fMC->GenEventHeader();
641 
642  if(!fAcceptOnlyHIJINGLabels) return eventHeader ;
643 
644  // TODO Check if it works from here ...
645 
646  AliGenCocktailEventHeader *cocktail = dynamic_cast<AliGenCocktailEventHeader *>(eventHeader);
647 
648  if(!cocktail) return 0x0 ;
649 
650  TList *genHeaders = cocktail->GetHeaders();
651 
652  Int_t nGenerators = genHeaders->GetEntries();
653  //printf("N generators %d \n", nGenerators);
654 
655  for(Int_t igen = 0; igen < nGenerators; igen++)
656  {
657  AliGenEventHeader * eventHeader2 = (AliGenEventHeader*)genHeaders->At(igen) ;
658  TString name = eventHeader2->GetName();
659 
660  //printf("Generator %d: Class Name %s, Name %s, title %s \n",igen, eventHeader2->ClassName(), name.Data(), eventHeader2->GetTitle());
661 
662  if(name.Contains("Hijing",TString::kIgnoreCase)) return eventHeader2 ;
663  }
664 
665  return 0x0;
666 
667  }
668  else if(ReadAODMCParticles() && GetAODMCHeader())
669  {
670  Int_t nGenerators = GetAODMCHeader()->GetNCocktailHeaders();
671  //printf("AliCaloTrackReader::GetGenEventHeader() - N headers %d\n",nGenerators);
672 
673  if( nGenerators <= 0) return 0x0;
674 
675  if(!fAcceptOnlyHIJINGLabels) return GetAODMCHeader()->GetCocktailHeader(0);
676 
677  for(Int_t igen = 0; igen < nGenerators; igen++)
678  {
679  AliGenEventHeader * eventHeader = GetAODMCHeader()->GetCocktailHeader(igen) ;
680  TString name = eventHeader->GetName();
681 
682  //printf("Generator %d: Class Name %s, Name %s, title %s \n",igen, eventHeader->ClassName(), name.Data(), eventHeader->GetTitle());
683 
684  if(name.Contains("Hijing",TString::kIgnoreCase)) return eventHeader ;
685  }
686 
687  return 0x0;
688 
689  }
690  else
691  {
692  //printf("AliCaloTrackReader::GetGenEventHeader() - MC header not available! \n");
693  return 0x0;
694  }
695 }
696 
697 //_________________________________________________________
700 //_________________________________________________________
702 {
703  AliInfo("Input are not AODs");
704 
705  return NULL ;
706 }
707 
708 //________________________________________________________
711 //________________________________________________________
712 AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader() const
713 {
714  AliInfo("Input are not AODs");
715 
716  return NULL ;
717 }
718 
719 //___________________________________________________________
727 //___________________________________________________________
728 Int_t AliCaloTrackReader::GetVertexBC(const AliVVertex * vtx)
729 {
730  Int_t vertexBC=vtx->GetBC();
731 
732  if(!fRecalculateVertexBC) return vertexBC;
733 
734  // Value not available, recalculate it.
735 
736  Double_t bz = fInputEvent->GetMagneticField();
737  Bool_t bc0 = kFALSE;
738  Int_t ntr = GetCTSTracks()->GetEntriesFast();
739  //printf("N Tracks %d\n",ntr);
740 
741  for(Int_t i = 0 ; i < ntr ; i++)
742  {
743  AliVTrack * track = (AliVTrack*) (GetCTSTracks()->At(i));
744 
745  //Check if has TOF info, if not skip
746  ULong_t status = track->GetStatus();
747  Bool_t okTOF = (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ;
748  vertexBC = track->GetTOFBunchCrossing(bz);
749  Float_t pt = track->Pt();
750 
751  if(!okTOF) continue;
752 
753  // Get DCA x, y
754  Double_t dca[2] = {1e6,1e6};
755  Double_t covar[3] = {1e6,1e6,1e6};
756  track->PropagateToDCA(vtx,bz,100.,dca,covar);
757 
758  if(AcceptDCA(pt,dca[0]))
759  {
760  if (vertexBC !=0 && fVertexBC != AliVTrack::kTOFBCNA) return vertexBC;
761  else if(vertexBC == 0) bc0 = kTRUE;
762  }
763  }
764 
765  if( bc0 ) vertexBC = 0 ;
766  else vertexBC = AliVTrack::kTOFBCNA ;
767 
768  return vertexBC;
769 }
770 
771 //_____________________________
774 //_____________________________
776 {
778  {
779  AliInfo("Cannot access stack and mcparticles at the same time, change them");
780  fReadStack = kFALSE;
781  fReadAODMCParticles = kFALSE;
782  }
783 
784  // Activate debug level in AliAnaWeights
785  if( fWeightUtils->GetDebug() >= 0 )
786  (AliAnalysisManager::GetAnalysisManager())->AddClassDebug(fWeightUtils->ClassName(), fWeightUtils->GetDebug());
787 }
788 
789 //_______________________________________
791 //_______________________________________
793 {
794  fDataType = kESD ;
795  fCTSPtMin = 0.1 ;
796  fEMCALPtMin = 0.1 ;
797  fPHOSPtMin = 0.1 ;
798  fCTSPtMax = 1000. ;
799  fEMCALPtMax = 1000. ;
800  fPHOSPtMax = 1000. ;
801 
802  //Track DCA cuts
803  // dca_xy cut = 0.0105+0.0350/TMath::Power(pt,1.1);
804  fTrackDCACut[0] = 0.0105;
805  fTrackDCACut[1] = 0.0350;
806  fTrackDCACut[2] = 1.1;
807 
808  //Do not filter the detectors input by default.
809  fFillEMCAL = kFALSE;
810  fFillDCAL = kFALSE;
811  fFillPHOS = kFALSE;
812  fFillCTS = kFALSE;
813  fFillEMCALCells = kFALSE;
814  fFillPHOSCells = kFALSE;
815 
816  fReadStack = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
817  fReadAODMCParticles = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
818  fDeltaAODFileName = "deltaAODPartCorr.root";
820  fEventTriggerMask = AliVEvent::kAny;
821  fMixEventTriggerMask = AliVEvent::kAnyINT;
822  fEventTriggerAtSE = kTRUE; // Use only events that pass event selection at SE base class
823 
824  fAcceptFastCluster = kTRUE;
825  fEventType = -1;
826 
827  //We want tracks fitted in the detectors:
828  //fTrackStatus=AliESDtrack::kTPCrefit;
829  //fTrackStatus|=AliESDtrack::kITSrefit;
830  fTrackStatus = 0;
831 
832  fV0ADC[0] = 0; fV0ADC[1] = 0;
833  fV0Mul[0] = 0; fV0Mul[1] = 0;
834 
835  fZvtxCut = 10.;
836 
837  fNMixedEvent = 1;
838 
841 
842  //Centrality
843  fUseAliCentrality = kFALSE;
844  fCentralityClass = "V0M";
845  fCentralityOpt = 100;
846  fCentralityBin[0] = fCentralityBin[1]=-1;
847 
848  fEventPlaneMethod = "V0";
849 
850  // Allocate memory (not sure this is the right place)
851  fCTSTracks = new TObjArray();
852  fEMCALClusters = new TObjArray();
853  fDCALClusters = new TObjArray();
854  fPHOSClusters = new TObjArray();
855  //fTriggerAnalysis = new AliTriggerAnalysis;
856  fAODBranchList = new TList ;
857 
858  fPileUpParamSPD[0] = 3 ; fPileUpParamSPD[1] = 0.8 ;
859  fPileUpParamSPD[2] = 3.0 ; fPileUpParamSPD[3] = 2.0 ; fPileUpParamSPD[4] = 5.0;
860 
861  // Parametrized time cut (LHC11d)
864 
865  // Parametrized time cut (LHC11c)
866  //fEMCALParamTimeCutMin[0] =-5; fEMCALParamTimeCutMin[1] =-1 ; fEMCALParamTimeCutMin[2] = 1.87; fEMCALParamTimeCutMin[3] = 0.4;
867  //fEMCALParamTimeCutMax[0] = 3.5; fEMCALParamTimeCutMax[1] = 50; fEMCALParamTimeCutMax[2] = 0.15; fEMCALParamTimeCutMax[3] = 1.6;
868 
869  fTimeStampRunMin = -1;
870  fTimeStampRunMax = 1e12;
873 
874  for(Int_t i = 0; i < 19; i++)
875  {
876  fEMCalBCEvent [i] = 0;
877  fEMCalBCEventCut[i] = 0;
878  fTrackBCEvent [i] = 0;
879  fTrackBCEventCut[i] = 0;
880  }
881 
882  // Trigger match-rejection
885 
886  fTriggerClusterBC = -10000 ;
890  fTriggerClusterId = -1;
891 
892  //Jets
895  if(!fNonStandardJets) fNonStandardJets = new TClonesArray("AliAODJet",100);
898  if(!fBackgroundJets) fBackgroundJets = new AliAODJetEventBackground();
899 
900  fSmearShowerShapeWidth = 0.005;
901 
902  fWeightUtils = new AliAnaWeights() ;
903  fEventWeight = 1 ;
904 
905 }
906 
907 //___________________________________________________
912 //___________________________________________________
914 {
915  AliGenHijingEventHeader* hijingHeader = dynamic_cast<AliGenHijingEventHeader *> (GetGenEventHeader());
916 
917  //printf("header %p, label %d\n",hijingHeader,label);
918 
919  if(!hijingHeader || label < 0 ) return kFALSE;
920 
921 
922  //printf("pass a), N produced %d\n",nproduced);
923 
924  if(label >= fNMCProducedMin && label < fNMCProducedMax)
925  {
926  //printf(" accept!, label is smaller than produced, N %d\n",nproduced);
927 
928  return kTRUE;
929  }
930 
931  if(ReadStack())
932  {
933  if(!GetStack()) return kFALSE;
934 
935  Int_t nprimaries = GetStack()->GetNtrack();
936 
937  if(label > nprimaries) return kFALSE;
938 
939  TParticle * mom = GetStack()->Particle(label);
940 
941  Int_t iMom = label;
942  Int_t iParent = mom->GetFirstMother();
943  while(iParent!=-1)
944  {
945  if(iParent >= fNMCProducedMin && iParent < fNMCProducedMax)
946  {
947  //printf("\t accept, mother is %d \n",iParent)
948  return kTRUE;
949  }
950 
951  iMom = iParent;
952  mom = GetStack()->Particle(iMom);
953  iParent = mom->GetFirstMother();
954  }
955 
956  return kFALSE ;
957 
958  } // ESD
959  else
960  {
961  TClonesArray* mcparticles = GetAODMCParticles();
962 
963  if(!mcparticles) return kFALSE;
964 
965  Int_t nprimaries = mcparticles->GetEntriesFast();
966 
967  if(label > nprimaries) return kFALSE;
968 
969  //printf("pass b) N primaries %d \n",nprimaries);
970 
971  if(label >= fNMCProducedMin && label < fNMCProducedMax)
972  {
973  return kTRUE;
974  }
975 
976  // Find grand parent, check if produced in the good range
977  AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
978 
979  Int_t iMom = label;
980  Int_t iParent = mom->GetMother();
981  while(iParent!=-1)
982  {
983  if(iParent >= fNMCProducedMin && iParent < fNMCProducedMax)
984  {
985  //printf("\t accept, mother is %d, with nProduced %d \n",iParent, nproduced);
986  return kTRUE;
987  }
988 
989  iMom = iParent;
990  mom = (AliAODMCParticle *) mcparticles->At(iMom);
991  iParent = mom->GetMother();
992 
993  }
994 
995  //printf("pass c), no match found \n");
996 
997  return kFALSE ;
998 
999  }//AOD
1000 }
1001 
1002 //__________________________________________________________________________
1005 //__________________________________________________________________________
1006 Bool_t AliCaloTrackReader::IsInTimeWindow(Double_t tof, Float_t energy) const
1007 {
1008  // Parametrized cut depending on E
1009  if(fUseParamTimeCut)
1010  {
1011  Float_t minCut= fEMCALParamTimeCutMin[0]+fEMCALParamTimeCutMin[1]*TMath::Exp(-(energy-fEMCALParamTimeCutMin[2])/fEMCALParamTimeCutMin[3]);
1012  Float_t maxCut= fEMCALParamTimeCutMax[0]+fEMCALParamTimeCutMax[1]*TMath::Exp(-(energy-fEMCALParamTimeCutMax[2])/fEMCALParamTimeCutMax[3]);
1013  //printf("tof %f, minCut %f, maxCut %f\n",tof,minCut,maxCut);
1014  if( tof < minCut || tof > maxCut ) return kFALSE ;
1015  }
1016 
1017  //In any case, the time should to be larger than the fixed window ...
1018  if( tof < fEMCALTimeCutMin || tof > fEMCALTimeCutMax ) return kFALSE ;
1019 
1020  return kTRUE ;
1021 }
1022 
1023 //________________________________________________
1026 //________________________________________________
1028 {
1029  return fInputEvent->IsPileupFromSPD((Int_t) fPileUpParamSPD[0] , fPileUpParamSPD[1] ,
1030  fPileUpParamSPD[2] , fPileUpParamSPD[3] , fPileUpParamSPD[4] );
1031  //printf("Param : %d, %2.2f, %2.2f, %2.2f, %2.2f\n",(Int_t) fPileUpParamSPD[0], fPileUpParamSPD[1], fPileUpParamSPD[2], fPileUpParamSPD[3], fPileUpParamSPD[4]);
1032 }
1033 
1034 //__________________________________________________
1036 //__________________________________________________
1038 {
1039  if(fNPileUpClusters > fNPileUpClustersCut) return kTRUE ;
1040  else return kFALSE;
1041 }
1042 
1043 //________________________________________________________
1045 //________________________________________________________
1047 {
1048  if( IsPileUpFromSPD() && IsPileUpFromEMCal()) return kTRUE ;
1049  else return kFALSE;
1050 }
1051 
1052 //_______________________________________________________
1054 //_______________________________________________________
1056 {
1057  if( IsPileUpFromSPD() || IsPileUpFromEMCal()) return kTRUE ;
1058  else return kFALSE;
1059 }
1060 
1061 //___________________________________________________________
1063 //___________________________________________________________
1065 {
1066  if( IsPileUpFromSPD() && !IsPileUpFromEMCal()) return kTRUE ;
1067  else return kFALSE;
1068 }
1069 
1070 //___________________________________________________________
1072 //___________________________________________________________
1074 {
1075  if( !IsPileUpFromSPD() && IsPileUpFromEMCal()) return kTRUE ;
1076  else return kFALSE;
1077 }
1078 
1079 //______________________________________________________________
1081 //______________________________________________________________
1083 {
1084  if( !IsPileUpFromSPD() && !IsPileUpFromEMCal()) return kTRUE ;
1085  else return kFALSE;
1086 }
1087 
1088 //___________________________________________________________________________________
1096 //___________________________________________________________________________________
1097 Bool_t AliCaloTrackReader::FillInputEvent(Int_t iEntry, const char * /*curFileName*/)
1098 {
1099  fEventNumber = iEntry;
1100  fTriggerClusterIndex = -1;
1101  fTriggerClusterId = -1;
1102  fIsTriggerMatch = kFALSE;
1103  fTriggerClusterBC = -10000;
1104  fIsExoticEvent = kFALSE;
1105  fIsBadCellEvent = kFALSE;
1106  fIsBadMaxCellEvent = kFALSE;
1107 
1108  fIsTriggerMatchOpenCut[0] = kFALSE ;
1109  fIsTriggerMatchOpenCut[1] = kFALSE ;
1110  fIsTriggerMatchOpenCut[2] = kFALSE ;
1111 
1112  //fCurrentFileName = TString(currentFileName);
1113  if(!fInputEvent)
1114  {
1115  AliInfo("Input event not available, skip event analysis");
1116  return kFALSE;
1117  }
1118 
1119  //-----------------------------------------------
1120  // Select the event depending on the trigger type
1121  // and other event characteristics
1122  // like the goodness of the EMCal trigger
1123  //-----------------------------------------------
1124 
1125  Bool_t accept = CheckEventTriggers();
1126  if(!accept) return kFALSE;
1127 
1128  AliDebug(1,"Pass Event trigger selection");
1129 
1130  //---------------------------------------------------------------------------
1131  // In case of analysis of events with jets, skip those with jet pt > 5 pt hard
1132  // To be used on for MC data in pT hard bins
1133  //---------------------------------------------------------------------------
1134 
1136  {
1137  if(!ComparePtHardAndJetPt()) return kFALSE ;
1138  AliDebug(1,"Pass Pt Hard - Jet rejection");
1139  }
1140 
1142  {
1143  if(!ComparePtHardAndClusterPt()) return kFALSE ;
1144  AliDebug(1,"Pass Pt Hard - Cluster rejection");
1145  }
1146 
1147  //------------------------------------------------------
1148  // Event rejection depending on time stamp
1149  //------------------------------------------------------
1150 
1152  {
1153  AliESDEvent* esd = dynamic_cast<AliESDEvent*> (fInputEvent);
1154  if(esd)
1155  {
1156  Int_t timeStamp = esd->GetTimeStamp();
1157  Float_t timeStampFrac = 1.*(timeStamp-fTimeStampRunMin) / (fTimeStampRunMax-fTimeStampRunMin);
1158 
1159  //printf("stamp0 %d, max0 %d, frac %f\n", timeStamp-fTimeStampRunMin,fTimeStampRunMax-fTimeStampRunMin, timeStampFrac);
1160 
1161  if(timeStampFrac < fTimeStampEventFracMin || timeStampFrac > fTimeStampEventFracMax) return kFALSE;
1162  }
1163 
1164  AliDebug(1,"Pass Time Stamp rejection");
1165  }
1166 
1167  //------------------------------------------------------
1168  // Event rejection depending on vertex, pileup, v0and
1169  //------------------------------------------------------
1170 
1171  FillVertexArray();
1172 
1173  //Reject events with Z vertex too large, only for SE analysis, if not, cut on the analysis code
1174  if(!GetMixedEvent() && TMath::Abs(fVertex[0][2]) > fZvtxCut) return kFALSE;
1175 
1177  {
1178  if( !CheckForPrimaryVertex() ) return kFALSE; // algorithm in ESD/AOD Readers
1179  if( TMath::Abs(fVertex[0][0] ) < 1.e-6 &&
1180  TMath::Abs(fVertex[0][1] ) < 1.e-6 &&
1181  TMath::Abs(fVertex[0][2] ) < 1.e-6 ) return kFALSE;
1182  }
1183 
1184  AliDebug(1,"Pass primary vertex rejection");
1185 
1186  //printf("Reader : IsPileUp %d, Multi %d\n",IsPileUpFromSPD(),fInputEvent->IsPileupFromSPDInMultBins());
1187 
1189  {
1190  // Do not analyze events with pileup
1191  Bool_t bPileup = IsPileUpFromSPD();
1192  //IsPileupFromSPDInMultBins() // method to try
1193  //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]);
1194  if(bPileup) return kFALSE;
1195 
1196  AliDebug(1,"Pass Pile-Up event rejection");
1197  }
1198 
1200  {
1201  AliVVZERO* v0 = fInputEvent->GetVZEROData();
1202 
1203  Bool_t bV0AND = ((v0->GetV0ADecision()==1) && (v0->GetV0CDecision()==1));
1204  //bV0AND = fTriggerAnalysis->IsOfflineTriggerFired((AliESDEvent*)fInputEvent, AliTriggerAnalysis::kV0AND);
1205  //printf("V0AND event? %d\n",bV0AND);
1206 
1207  if(!bV0AND)
1208  {
1209  AliDebug(1,"Reject event by V0AND");
1210  return kFALSE;
1211  }
1212 
1213  AliDebug(1,"Pass V0AND event rejection");
1214  }
1215 
1216  //------------------------------------------------------
1217  // Check if there is a centrality value, PbPb analysis,
1218  // and if a centrality bin selection is requested
1219  //------------------------------------------------------
1220 
1221  //If we need a centrality bin, we select only those events in the corresponding bin.
1222  Int_t cen = -1;
1223  if ( fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
1224  {
1225  cen = GetEventCentrality();
1226 
1227  if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return kFALSE; //reject events out of bin.
1228 
1229  AliDebug(1,"Pass centrality rejection");
1230  }
1231 
1232  // Recover the weight assigned to the event, if provided
1233  // right now only for pT-hard bins and centrality depedent weights
1235  {
1237 
1239  }
1240  //---------------------------------------------------------------------------
1241  // In case of MC analysis, set the range of interest in the MC particles list
1242  // The particle selection is done inside the FillInputDetector() methods
1243  //---------------------------------------------------------------------------
1245 
1246  //printf("N min %d, N max %d\n",fNMCProducedMin,fNMCProducedMax);
1247 
1248  //-------------------------------------------------------
1249  // Get the main vertex BC, in case not available
1250  // it is calculated in FillCTS checking the BC of tracks
1251  //------------------------------------------------------
1252  fVertexBC = fInputEvent->GetPrimaryVertex()->GetBC();
1253 
1254  //-----------------------------------------------
1255  // Fill the arrays with cluster/tracks/cells data
1256  //-----------------------------------------------
1257 
1258  if(fFillCTS)
1259  {
1260  FillInputCTS();
1261  //Accept events with at least one track
1262  if(fTrackMult == 0 && fDoRejectNoTrackEvents) return kFALSE ;
1263 
1264  AliDebug(1,"Pass rejection of null track events");
1265  }
1266 
1268  {
1269  if(fVertexBC != 0 && fVertexBC != AliVTrack::kTOFBCNA) return kFALSE ;
1270 
1271  AliDebug(1,"Pass rejection of events with vertex at BC!=0");
1272  }
1273 
1274  if(fFillEMCALCells)
1276 
1277  if(fFillPHOSCells)
1279 
1280  if(fFillEMCAL || fFillDCAL)
1281  FillInputEMCAL();
1282 
1283  if(fFillPHOS)
1284  FillInputPHOS();
1285 
1286  FillInputVZERO();
1287 
1288  //one specified jet branch
1293 
1294  AliDebug(1,"Event accepted for analysis");
1295 
1296  return kTRUE ;
1297 }
1298 
1299 //__________________________________________________
1302 //__________________________________________________
1304 {
1305  if(fUseAliCentrality)
1306  {
1307  if ( !GetCentrality() ) return -1;
1308 
1309  if (fCentralityOpt == 100) return (Int_t) GetCentrality()->GetCentralityPercentile(fCentralityClass); // 100 bins max
1310  else if(fCentralityOpt == 10) return GetCentrality()->GetCentralityClass10(fCentralityClass);// 10 bins max
1311  else if(fCentralityOpt == 20) return GetCentrality()->GetCentralityClass5(fCentralityClass); // 20 bins max
1312  else
1313  {
1314  AliInfo(Form("Unknown centrality option %d, use 10, 20 or 100\n",fCentralityOpt));
1315  return -1;
1316  }
1317  }
1318  else
1319  {
1320  if ( !GetMultSelCen() ) return -1;
1321 
1322  return GetMultSelCen()->GetMultiplicityPercentile(fCentralityClass, kTRUE); // returns centrality only for events used in calibration
1323 
1324  // equivalent to
1325  //GetMultSelCen()->GetMultiplicityPercentile("V0M", kFALSE); // returns centrality for any event
1326  //Int_t qual = GetMultSelCen()->GetEvSelCode(); if (qual ! = 0) cent = qual;
1327  }
1328 }
1329 
1330 //_____________________________________________________
1333 //_____________________________________________________
1335 {
1336  if( !GetEventPlane() ) return -1000;
1337 
1338  Float_t ep = GetEventPlane()->GetEventplane(GetEventPlaneMethod(), GetInputEvent());
1339 
1340  if(GetEventPlaneMethod()=="Q" && (ep < 0 || ep > TMath::Pi()))
1341  {
1342  AliDebug(1,Form("Bad EP for <Q> method : %f\n",ep));
1343  return -1000;
1344  }
1345  else if(GetEventPlaneMethod().Contains("V0") )
1346  {
1347  if((ep > TMath::Pi()/2 || ep < -TMath::Pi()/2))
1348  {
1349  AliDebug(1,Form("Bad EP for <%s> method : %f\n",GetEventPlaneMethod().Data(), ep));
1350  return -1000;
1351  }
1352 
1353  ep+=TMath::Pi()/2; // put same range as for <Q> method
1354  }
1355 
1356  AliDebug(3,Form("Event plane angle %f",ep));
1357 
1358 // if(fDebug > 0 )
1359 // {
1360 // if (ep > TMath::Pi()) printf("AliCaloTrackReader::GetEventPlaneAngle() - Too large angle = %f\n",ep);
1361 // else if(ep < 0 ) printf("AliCaloTrackReader::GetEventPlaneAngle() - Negative angle = %f\n" ,ep);
1362 // }
1363 
1364  return ep;
1365 }
1366 
1367 //__________________________________________________________
1369 //__________________________________________________________
1370 void AliCaloTrackReader::GetVertex(Double_t vertex[3]) const
1371 {
1372  vertex[0] = fVertex[0][0];
1373  vertex[1] = fVertex[0][1];
1374  vertex[2] = fVertex[0][2];
1375 }
1376 
1377 //__________________________________________________________________________
1379 //__________________________________________________________________________
1380 void AliCaloTrackReader::GetVertex(Double_t vertex[3], Int_t evtIndex) const
1381 {
1382  vertex[0] = fVertex[evtIndex][0];
1383  vertex[1] = fVertex[evtIndex][1];
1384  vertex[2] = fVertex[evtIndex][2];
1385 }
1386 
1387 //________________________________________
1390 //________________________________________
1392 {
1393  // Delete previous vertex
1394  if(fVertex)
1395  {
1396  for (Int_t i = 0; i < fNMixedEvent; i++)
1397  {
1398  delete [] fVertex[i] ;
1399  }
1400  delete [] fVertex ;
1401  }
1402 
1403  fVertex = new Double_t*[fNMixedEvent] ;
1404  for (Int_t i = 0; i < fNMixedEvent; i++)
1405  {
1406  fVertex[i] = new Double_t[3] ;
1407  fVertex[i][0] = 0.0 ;
1408  fVertex[i][1] = 0.0 ;
1409  fVertex[i][2] = 0.0 ;
1410  }
1411 
1412  if (!fMixedEvent)
1413  { // Single event analysis
1414  if(fDataType!=kMC)
1415  {
1416 
1417  if(fInputEvent->GetPrimaryVertex())
1418  {
1419  fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]);
1420  }
1421  else
1422  {
1423  AliWarning("NULL primary vertex");
1424  fVertex[0][0]=0.; fVertex[0][1]=0.; fVertex[0][2]=0.;
1425  }//Primary vertex pointer do not exist
1426 
1427  } else
1428  {// MC read event
1429  fVertex[0][0]=0.; fVertex[0][1]=0.; fVertex[0][2]=0.;
1430  }
1431 
1432  AliDebug(1,Form("Single Event Vertex : %f,%f,%f\n",fVertex[0][0],fVertex[0][1],fVertex[0][2]));
1433 
1434  } else
1435  { // MultiEvent analysis
1436  for (Int_t iev = 0; iev < fNMixedEvent; iev++)
1437  {
1438  if (fMixedEvent->GetVertexOfEvent(iev))
1439  fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
1440  else
1441  AliWarning("No vertex found");
1442 
1443  AliDebug(1,Form("Multi Event %d Vertex : %f,%f,%f",iev,fVertex[iev][0],fVertex[iev][1],fVertex[iev][2]));
1444  }
1445  }
1446 }
1447 
1448 //_____________________________________
1454 //_____________________________________
1456 {
1457  AliDebug(1,"Begin");
1458 
1459  Double_t pTrack[3] = {0,0,0};
1460 
1461  Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
1462  fTrackMult = 0;
1463  Int_t nstatus = 0;
1464  Double_t bz = GetInputEvent()->GetMagneticField();
1465 
1466  for(Int_t i = 0; i < 19; i++)
1467  {
1468  fTrackBCEvent [i] = 0;
1469  fTrackBCEventCut[i] = 0;
1470  }
1471 
1472  Bool_t bc0 = kFALSE;
1473  if(fRecalculateVertexBC) fVertexBC = AliVTrack::kTOFBCNA;
1474 
1475  for (Int_t itrack = 0; itrack < nTracks; itrack++)
1476  {
1477  AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
1478 
1479  if(fAcceptOnlyHIJINGLabels && !IsHIJINGLabel(track->GetLabel())) continue ;
1480 
1481  //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
1482  ULong_t status = track->GetStatus();
1483 
1484  if (fTrackStatus && !((status & fTrackStatus) == fTrackStatus))
1485  continue ;
1486 
1487  nstatus++;
1488 
1489  // Select the tracks depending on cuts of AOD or ESD
1490  if(!SelectTrack(track, pTrack)) continue ;
1491 
1492  // TOF cuts
1493  Bool_t okTOF = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
1494  Double_t tof = -1000;
1495  Int_t trackBC = -1000 ;
1496 
1497  if(fAccessTrackTOF)
1498  {
1499  if(okTOF)
1500  {
1501  trackBC = track->GetTOFBunchCrossing(bz);
1502  SetTrackEventBC(trackBC+9);
1503 
1504  tof = track->GetTOFsignal()*1e-3;
1505 
1506  //After selecting tracks with small DCA, pointing to vertex, set vertex BC depeding on tracks BC
1508  {
1509  if (trackBC != 0 && trackBC != AliVTrack::kTOFBCNA) fVertexBC = trackBC;
1510  else if(trackBC == 0) bc0 = kTRUE;
1511  }
1512 
1513  //In any case, the time should to be larger than the fixed window ...
1514  if( fUseTrackTimeCut && (trackBC !=0 || tof < fTrackTimeCutMin || tof > fTrackTimeCutMax) )
1515  {
1516  //printf("Remove track time %f and bc = %d\n",tof,trackBC);
1517  continue ;
1518  }
1519  //else printf("Accept track time %f and bc = %d\n",tof,trackBC);
1520  }
1521  }
1522 
1523  fMomentum.SetPxPyPzE(pTrack[0],pTrack[1],pTrack[2],0);
1524 
1525  // DCA cuts
1526 
1527  if(fUseTrackDCACut)
1528  {
1529  Float_t dcaTPC =-999;
1530  //In case of AODs, TPC tracks cannot be propagated back to primary vertex,
1531  if( fDataType == kAOD ) dcaTPC = ((AliAODTrack*) track)->DCA();
1532 
1533  //normal way to get the dca, cut on dca_xy
1534  if(dcaTPC==-999)
1535  {
1536  Double_t dca[2] = {1e6,1e6};
1537  Double_t covar[3] = {1e6,1e6,1e6};
1538  Bool_t okDCA = track->PropagateToDCA(fInputEvent->GetPrimaryVertex(),bz,100.,dca,covar);
1539  if( okDCA) okDCA = AcceptDCA(fMomentum.Pt(),dca[0]);
1540  if(!okDCA)
1541  {
1542  //printf("AliCaloTrackReader::FillInputCTS() - Reject track pt %2.2f, dca_xy %2.4f\n",fMomentum.Pt(),dca[0]);
1543  continue ;
1544  }
1545  }
1546  }// DCA cuts
1547 
1548  //Count the tracks in eta < 0.9
1549  if(TMath::Abs(track->Eta())< fTrackMultEtaCut) fTrackMult++;
1550 
1551  if(fCTSPtMin > fMomentum.Pt() || fCTSPtMax < fMomentum.Pt()) continue ;
1552 
1553  // Check effect of cuts on track BC
1554  if(fAccessTrackTOF && okTOF) SetTrackEventBCcut(trackBC+9);
1555 
1556  if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS)) continue;
1557 
1558  AliDebug(2,Form("Selected tracks pt %3.2f, phi %3.2f deg, eta %3.2f",
1559  fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
1560 
1561  if (fMixedEvent) track->SetID(itrack);
1562 
1563 
1564  fCTSTracks->Add(track);
1565 
1566  }// track loop
1567 
1568  if( fRecalculateVertexBC && (fVertexBC == 0 || fVertexBC == AliVTrack::kTOFBCNA))
1569  {
1570  if( bc0 ) fVertexBC = 0 ;
1571  else fVertexBC = AliVTrack::kTOFBCNA ;
1572  }
1573 
1574  AliDebug(1,Form("AOD entries %d, input tracks %d, pass status %d, multipliticy %d", fCTSTracks->GetEntriesFast(), nTracks, nstatus, fTrackMult));//fCTSTracksNormalInputEntries);
1575 }
1576 
1577 //_______________________________________________________________________________
1584 //_______________________________________________________________________________
1585 void AliCaloTrackReader::FillInputEMCALAlgorithm(AliVCluster * clus, Int_t iclus)
1586 {
1587  // Accept clusters with the proper label, TO DO, use the new more General methods!!!
1588  if(fAcceptOnlyHIJINGLabels && !IsHIJINGLabel( clus->GetLabel() )) return ;
1589 
1590  Int_t vindex = 0 ;
1591  if (fMixedEvent)
1592  vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
1593 
1594  clus->GetMomentum(fMomentum, fVertex[vindex]);
1595 
1596  //if( (fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10 )
1597  AliDebug(2,Form("Input cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
1598  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
1599 
1601  {
1602  //Recalibrate the cluster energy
1603  if(GetCaloUtils()->IsRecalibrationOn())
1604  {
1606 
1607  clus->SetE(energy);
1608  //printf("Recalibrated Energy %f\n",clus->E());
1609 
1612 
1613  } // recalculate E
1614 
1615  //Recalculate distance to bad channels, if new list of bad channels provided
1617 
1618  //Recalculate cluster position
1619  if(GetCaloUtils()->IsRecalculationOfClusterPositionOn())
1620  {
1622  //clus->GetPosition(pos);
1623  //printf("After Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1624  }
1625 
1626  // Recalculate TOF
1627  if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1628  {
1629  Double_t tof = clus->GetTOF();
1630  Float_t frac =-1;
1631  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
1632 
1634  {
1635  tof = fEMCALCells->GetCellTime(absIdMax);
1636  }
1637 
1638  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1639 
1640  //additional L1 phase shift
1641  if(GetCaloUtils()->GetEMCALRecoUtils()->IsL1PhaseInTimeRecalibrationOn()){
1642  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absIdMax), fInputEvent->GetBunchCrossNumber(), tof);
1643  }
1644 
1645  clus->SetTOF(tof);
1646 
1647  }// Time recalibration
1648  }
1649 
1650  //Reject clusters with bad channels, close to borders and exotic;
1651  Bool_t goodCluster = GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,
1652  GetCaloUtils()->GetEMCALGeometry(),
1653  GetEMCALCells(),fInputEvent->GetBunchCrossNumber());
1654 
1655  if(!goodCluster)
1656  {
1657  //if( (fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10 )
1658  AliDebug(2,Form("Bad cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
1659  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
1660 
1661  return;
1662  }
1663 
1664  //Mask all cells in collumns facing ALICE thick material if requested
1665  if(GetCaloUtils()->GetNMaskCellColumns())
1666  {
1667  Int_t absId = -1;
1668  Int_t iSupMod = -1;
1669  Int_t iphi = -1;
1670  Int_t ieta = -1;
1671  Bool_t shared = kFALSE;
1672  GetCaloUtils()->GetEMCALRecoUtils()->GetMaxEnergyCell(GetCaloUtils()->GetEMCALGeometry(), GetEMCALCells(),clus,absId,iSupMod,ieta,iphi,shared);
1673 
1674  AliDebug(2,Form("Masked collumn: cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
1675  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
1676 
1677 
1678  if(GetCaloUtils()->MaskFrameCluster(iSupMod, ieta)) return;
1679  }
1680 
1682  {
1683  if(clus->GetNLabels()==0 || clus->GetLabel() < 0) return;
1684  //else printf("Embedded cluster, %d, n label %d label %d \n",iclus,clus->GetNLabels(),clus->GetLabel());
1685  }
1686 
1687  //Float_t pos[3];
1688  //clus->GetPosition(pos);
1689  //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1690 
1691  //Correct non linearity or smear energy
1693  {
1695 
1696  //if( (fDebug > 5 && fMomentum.E() > 0.1) || fDebug > 10 )
1697  AliDebug(5,Form("Correct Non Lin: Old E %3.2f, New E %3.2f",
1698  fMomentum.E(),clus->E()));
1699 
1700  // In case of MC analysis, to match resolution/calibration in real data
1701  // Not needed anymore, just leave for MC studies on systematics
1702  if( GetCaloUtils()->GetEMCALRecoUtils()->IsClusterEnergySmeared() )
1703  {
1704  Float_t rdmEnergy = GetCaloUtils()->GetEMCALRecoUtils()->SmearClusterEnergy(clus);
1705 
1706  //if( (fDebug > 5 && fMomentum.E() > 0.1) || fDebug > 10 )
1707  AliDebug(5,Form("Smear energy: Old E %3.2f, New E %3.2f",clus->E(),rdmEnergy));
1708 
1709  clus->SetE(rdmEnergy);
1710  }
1711  }
1712 
1713  Double_t tof = clus->GetTOF()*1e9;
1714 
1715  Int_t bc = TMath::Nint(tof/50) + 9;
1716  //printf("tof %2.2f, bc+5=%d\n",tof,bc);
1717 
1718  SetEMCalEventBC(bc);
1719 
1720  if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E())
1721  {
1722  AliDebug(2,Form("Too low energy cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
1723  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
1724 
1725  return ;
1726  }
1727 
1728  clus->GetMomentum(fMomentum, fVertex[vindex]);
1729 
1730  SetEMCalEventBCcut(bc);
1731 
1732  if(!IsInTimeWindow(tof,clus->E()))
1733  {
1734  fNPileUpClusters++ ;
1735  if(fUseEMCALTimeCut)
1736  {
1737  AliDebug(2,Form("Out of time window E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f, time %e",
1738  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta(),tof));
1739 
1740  return ;
1741  }
1742  }
1743  else
1745 
1746  if (fMixedEvent)
1747  clus->SetID(iclus) ;
1748 
1749  // Select cluster fiducial region
1750  Bool_t bEMCAL = kFALSE;
1751  Bool_t bDCAL = kFALSE;
1752  if(fCheckFidCut)
1753  {
1754  if(fFillEMCAL && fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) bEMCAL = kTRUE ;
1755  if(fFillDCAL && fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kDCAL )) bDCAL = kTRUE ;
1756  }
1757  else
1758  {
1759  bEMCAL = kTRUE;
1760  }
1761 
1762  //if((fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10)
1763  AliDebug(2,Form("Selected clusters (EMCAL%d, DCAL%d), E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
1764  bEMCAL,bDCAL,fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
1765 
1766 
1767  // Smear the SS to try to match data and simulations,
1768  // do it only for simulations.
1769  if( fSmearShowerShape && clus->GetNCells() > 2)
1770  {
1771  AliDebug(2,Form("Smear shower shape - Original: %2.4f\n", clus->GetM02()));
1773  {
1774  clus->SetM02( clus->GetM02() + fRandom.Landau(0, fSmearShowerShapeWidth) );
1775  }
1777  {
1778  if(iclus%3 == 0 && clus->GetM02() > 0.1) clus->SetM02( clus->GetM02() + fRandom.Landau(0.05, fSmearShowerShapeWidth) ); //fSmearShowerShapeWidth = 0.035
1779  }
1780  //clus->SetM02( fRandom.Landau(clus->GetM02(), fSmearShowerShapeWidth) );
1781  AliDebug(2,Form("Width %2.4f Smeared : %2.4f\n", fSmearShowerShapeWidth,clus->GetM02()));
1782  }
1783 
1784  // Fill the corresponding array. Usually just filling EMCal array with upper or lower clusters is enough, but maybe we want to do EMCal-DCal correlations.
1785  if (bEMCAL) fEMCALClusters->Add(clus);
1786  else if(bDCAL ) fDCALClusters ->Add(clus);
1787 }
1788 
1789 //_______________________________________
1796 //_______________________________________
1798 {
1799  AliDebug(1,"Begin");
1800 
1801  // First recalibrate cells, time or energy
1802  // if(GetCaloUtils()->IsRecalibrationOn())
1803  // GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCells(GetCaloUtils()->GetEMCALGeometry(),
1804  // GetEMCALCells(),
1805  // fInputEvent->GetBunchCrossNumber());
1806 
1807  fNPileUpClusters = 0; // Init counter
1808  fNNonPileUpClusters = 0; // Init counter
1809  for(Int_t i = 0; i < 19; i++)
1810  {
1811  fEMCalBCEvent [i] = 0;
1812  fEMCalBCEventCut[i] = 0;
1813  }
1814 
1815  //Loop to select clusters in fiducial cut and fill container with aodClusters
1816  if(fEMCALClustersListName=="")
1817  {
1818  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
1819  for (Int_t iclus = 0; iclus < nclusters; iclus++)
1820  {
1821  AliVCluster * clus = 0;
1822  if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1823  {
1824  if (clus->IsEMCAL())
1825  {
1826  FillInputEMCALAlgorithm(clus, iclus);
1827  }//EMCAL cluster
1828  }// cluster exists
1829  }// cluster loop
1830 
1831  //Recalculate track matching
1833 
1834  }//Get the clusters from the input event
1835  else
1836  {
1837  TClonesArray * clusterList = 0x0;
1838 
1839  if (fInputEvent->FindListObject(fEMCALClustersListName))
1840  {
1841  clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
1842  }
1843  else if(fOutputEvent)
1844  {
1845  clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
1846  }
1847 
1848  if(!clusterList)
1849  {
1850  AliWarning(Form("Wrong name of list with clusters? <%s>",fEMCALClustersListName.Data()));
1851  return;
1852  }
1853 
1854  Int_t nclusters = clusterList->GetEntriesFast();
1855  for (Int_t iclus = 0; iclus < nclusters; iclus++)
1856  {
1857  AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
1858  //printf("E %f\n",clus->E());
1859  if (clus) FillInputEMCALAlgorithm(clus, iclus);
1860  else AliWarning("Null cluster in list!");
1861  }// cluster loop
1862 
1863  // Recalculate the pile-up time, in case long time clusters removed during clusterization
1864  //printf("Input event INIT : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
1865 
1866  fNPileUpClusters = 0; // Init counter
1867  fNNonPileUpClusters = 0; // Init counter
1868  for(Int_t i = 0; i < 19; i++)
1869  {
1870  fEMCalBCEvent [i] = 0;
1871  fEMCalBCEventCut[i] = 0;
1872  }
1873 
1874  for (Int_t iclus = 0; iclus < fInputEvent->GetNumberOfCaloClusters(); iclus++)
1875  {
1876  AliVCluster * clus = 0;
1877 
1878  if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
1879  {
1880  if (clus->IsEMCAL())
1881  {
1882 
1883  Float_t frac =-1;
1884  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
1885  Double_t tof = clus->GetTOF();
1886  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1887  //additional L1 phase shift
1888  if(GetCaloUtils()->GetEMCALRecoUtils()->IsL1PhaseInTimeRecalibrationOn()){
1889  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absIdMax), fInputEvent->GetBunchCrossNumber(), tof);
1890  }
1891 
1892  tof*=1e9;
1893 
1894  //printf("Input event cluster : AbsIdMax %d, E %2.2f, time %2.2f \n", absIdMax,clus->E(),tof);
1895 
1896  //Reject clusters with bad channels, close to borders and exotic;
1897  if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber())) continue;
1898 
1899  Int_t bc = TMath::Nint(tof/50) + 9;
1900  SetEMCalEventBC(bc);
1901 
1902  if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E()) continue ;
1903 
1904  clus->GetMomentum(fMomentum, fVertex[0]);
1905 
1906  if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) return ;
1907 
1908  SetEMCalEventBCcut(bc);
1909 
1910  if(!IsInTimeWindow(tof,clus->E()))
1911  fNPileUpClusters++ ;
1912  else
1914 
1915  }
1916  }
1917  }
1918 
1919  //printf("Input event : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
1920 
1921  // Recalculate track matching, not necessary if already done in the reclusterization task.
1922  // in case it was not done ...
1924 
1925  }
1926 
1927  AliDebug(1,Form("AOD entries %d, n pile-up clusters %d, n non pile-up %d", fEMCALClusters->GetEntriesFast(),fNPileUpClusters,fNNonPileUpClusters));
1928 }
1929 
1930 //_______________________________________
1932 //_______________________________________
1934 {
1935  AliDebug(1,"Begin");
1936 
1937  // Loop to select clusters in fiducial cut and fill container with aodClusters
1938  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
1939  for (Int_t iclus = 0; iclus < nclusters; iclus++)
1940  {
1941  AliVCluster * clus = fInputEvent->GetCaloCluster(iclus) ;
1942  if ( !clus ) continue ;
1943 
1944  if ( !clus->IsPHOS() ) continue ;
1945 
1946  // Skip CPV input
1947  if( clus->GetType() == AliVCluster::kPHOSCharged ) continue ;
1948 
1949  if(fAcceptOnlyHIJINGLabels && !IsHIJINGLabel(clus->GetLabel())) continue ;
1950 
1951  // Check if the cluster contains any bad channel and if close to calorimeter borders
1952  if( GetCaloUtils()->ClusterContainsBadChannel(kPHOS,clus->GetCellsAbsId(), clus->GetNCells()))
1953  continue;
1954 
1955  if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells()))
1956  continue;
1957 
1959  {
1960  // Recalibrate the cluster energy
1962  {
1963  Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
1964  clus->SetE(energy);
1965  }
1966  }
1967 
1968  // Dead code? remove?
1969  Int_t vindex = 0 ;
1970  if (fMixedEvent)
1971  vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
1972 
1973  clus->GetMomentum(fMomentum, fVertex[vindex]);
1974 
1976  continue ;
1977 
1978  if (fPHOSPtMin > fMomentum.E() || fPHOSPtMax < fMomentum.E() )
1979  continue ;
1980 
1981  //if(fDebug > 2 && fMomentum.E() > 0.1)
1982  AliDebug(2,Form("Selected clusters E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
1983  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
1984 
1985  // Dead code? remove?
1986  if (fMixedEvent)
1987  clus->SetID(iclus) ;
1988 
1989  fPHOSClusters->Add(clus);
1990 
1991  } // esd cluster loop
1992 
1993  AliDebug(1,Form("AOD entries %d",fPHOSClusters->GetEntriesFast()));
1994 
1995 }
1996 
1997 //____________________________________________
1999 //____________________________________________
2001 {
2002  fEMCALCells = fInputEvent->GetEMCALCells();
2003 }
2004 
2005 //___________________________________________
2007 //___________________________________________
2009 {
2010  fPHOSCells = fInputEvent->GetPHOSCells();
2011 }
2012 
2013 //_______________________________________
2016 //_______________________________________
2018 {
2019  AliVVZERO* v0 = fInputEvent->GetVZEROData();
2020  //printf("Init V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
2021 
2022  if (v0)
2023  {
2024  AliESDVZERO* esdV0 = dynamic_cast<AliESDVZERO*> (v0);
2025  for (Int_t i = 0; i < 32; i++)
2026  {
2027  if(esdV0)
2028  {//Only available in ESDs
2029  fV0ADC[0] += (Int_t)esdV0->GetAdcV0C(i);
2030  fV0ADC[1] += (Int_t)esdV0->GetAdcV0A(i);
2031  }
2032 
2033  fV0Mul[0] += (Int_t)v0->GetMultiplicityV0C(i);
2034  fV0Mul[1] += (Int_t)v0->GetMultiplicityV0A(i);
2035  }
2036 
2037  AliDebug(1,Form("ADC (%d,%d), Multiplicity (%d,%d)",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]));
2038  }
2039  else
2040  {
2041  AliDebug(1,"Cannot retrieve V0 ESD! Run w/ null V0 charges");
2042  }
2043 }
2044 
2045 //_________________________________________________
2049 //_________________________________________________
2051 {
2052  AliDebug(2,"Begin");
2053 
2054  //
2055  //check if branch name is given
2056  if(!fInputNonStandardJetBranchName.Length())
2057  {
2058  fInputEvent->Print();
2059  AliFatal("No non-standard jet branch name specified. Specify among existing ones.");
2060  }
2061 
2062  fNonStandardJets = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fInputNonStandardJetBranchName.Data()));
2063 
2064  if(!fNonStandardJets)
2065  {
2066  //check if jet branch exist; exit if not
2067  fInputEvent->Print();
2068 
2069  AliFatal(Form("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fInputNonStandardJetBranchName.Data()));
2070  }
2071  else
2072  {
2073  AliDebug(1,Form("AOD input jets %d", fNonStandardJets->GetEntriesFast()));
2074  }
2075 }
2076 
2077 //_________________________________________________
2081 //_________________________________________________
2083 {
2084  AliDebug(1,"Begin");
2085  //
2086  //check if branch name is given
2087  if(!fInputBackgroundJetBranchName.Length())
2088  {
2089  fInputEvent->Print();
2090 
2091  AliFatal("No background jet branch name specified. Specify among existing ones.");
2092  }
2093 
2094  fBackgroundJets = (AliAODJetEventBackground*)(fInputEvent->FindListObject(fInputBackgroundJetBranchName.Data()));
2095 
2096  if(!fBackgroundJets)
2097  {
2098  //check if jet branch exist; exit if not
2099  fInputEvent->Print();
2100 
2101  AliFatal(Form("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fInputBackgroundJetBranchName.Data()));
2102  }
2103  else
2104  {
2105  AliDebug(1,"FillInputBackgroundJets");
2106  fBackgroundJets->Print("");
2107  }
2108 }
2109 
2110 //____________________________________________________________________
2116 //____________________________________________________________________
2117 TArrayI AliCaloTrackReader::GetTriggerPatches(Int_t tmin, Int_t tmax )
2118 {
2119  // init some variables
2120  Int_t trigtimes[30], globCol, globRow,ntimes, i;
2121  Int_t absId = -1; //[100];
2122  Int_t nPatch = 0;
2123 
2124  TArrayI patches(0);
2125 
2126  // get object pointer
2127  AliVCaloTrigger *caloTrigger = GetInputEvent()->GetCaloTrigger( "EMCAL" );
2128 
2129  if(!caloTrigger)
2130  {
2131  AliError("Trigger patches input (AliVCaloTrigger) not available in data!");
2132  return patches;
2133  }
2134 
2135  //printf("CaloTrigger Entries %d\n",caloTrigger->GetEntries() );
2136 
2137  // class is not empty
2138  if( caloTrigger->GetEntries() > 0 )
2139  {
2140  // must reset before usage, or the class will fail
2141  caloTrigger->Reset();
2142 
2143  // go throuth the trigger channels
2144  while( caloTrigger->Next() )
2145  {
2146  // get position in global 2x2 tower coordinates
2147  caloTrigger->GetPosition( globCol, globRow );
2148 
2149  //L0
2150  if(IsEventEMCALL0())
2151  {
2152  // get dimension of time arrays
2153  caloTrigger->GetNL0Times( ntimes );
2154 
2155  // no L0s in this channel
2156  // presence of the channel in the iterator still does not guarantee that L0 was produced!!
2157  if( ntimes < 1 )
2158  continue;
2159 
2160  // get timing array
2161  caloTrigger->GetL0Times( trigtimes );
2162  //printf("Get L0 patch : n times %d - trigger time window %d - %d\n",ntimes, tmin,tmax);
2163 
2164  // go through the array
2165  for( i = 0; i < ntimes; i++ )
2166  {
2167  // check if in cut - 8,9 shall be accepted in 2011
2168  if( trigtimes[i] >= tmin && trigtimes[i] <= tmax )
2169  {
2170  //printf("Accepted trigger time %d \n",trigtimes[i]);
2171  //if(nTrig > 99) continue;
2172  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(globCol,globRow, absId);
2173  //printf("pass the time cut globCol %d, globRow %d absId %d\n",globCol,globRow, absId);
2174  patches.Set(nPatch+1);
2175  patches.AddAt(absId,nPatch++);
2176  }
2177  } // trigger time array
2178  }//L0
2179  else if(IsEventEMCALL1()) // L1
2180  {
2181  Int_t bit = 0;
2182  caloTrigger->GetTriggerBits(bit);
2183 
2184  Int_t sum = 0;
2185  caloTrigger->GetL1TimeSum(sum);
2186  //fBitEGA-=2;
2187  Bool_t isEGA1 = ((bit >> fBitEGA ) & 0x1) && IsEventEMCALL1Gamma1() ;
2188  Bool_t isEGA2 = ((bit >> (fBitEGA+1)) & 0x1) && IsEventEMCALL1Gamma2() ;
2189  Bool_t isEJE1 = ((bit >> fBitEJE ) & 0x1) && IsEventEMCALL1Jet1 () ;
2190  Bool_t isEJE2 = ((bit >> (fBitEJE+1)) & 0x1) && IsEventEMCALL1Jet2 () ;
2191 
2192  //if((bit>> fBitEGA )&0x1) printf("Trig Bit %d - bit %d - EG1 %d - EG2 %d\n",fBitEGA ,bit,IsEventEMCALL1Gamma1(),IsEventEMCALL1Gamma2());
2193  //if((bit>>(fBitEGA+1))&0x1) printf("Trig Bit %d - bit %d - EG1 %d - EG2 %d\n",fBitEGA+1,bit,IsEventEMCALL1Gamma1(),IsEventEMCALL1Gamma2());
2194 
2195  if(!isEGA1 && !isEJE1 && !isEGA2 && !isEJE2) continue;
2196 
2197  Int_t patchsize = -1;
2198  if (isEGA1 || isEGA2) patchsize = 2;
2199  else if (isEJE1 || isEJE2) patchsize = 16;
2200 
2201  //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",
2202  // bit,sum,patchsize,isEGA1,isEGA2,isEJE1,isEJE2,fBitEGA,fBitEJE,IsEventEMCALL1Gamma(),IsEventEMCALL1Jet());
2203 
2204 
2205  // add 2x2 (EGA) or 16x16 (EJE) patches
2206  for(Int_t irow=0; irow < patchsize; irow++)
2207  {
2208  for(Int_t icol=0; icol < patchsize; icol++)
2209  {
2210  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(globCol+icol,globRow+irow, absId);
2211  //printf("pass the time cut globCol %d, globRow %d absId %d\n",globCol,globRow, absId);
2212  patches.Set(nPatch+1);
2213  patches.AddAt(absId,nPatch++);
2214  }
2215  }
2216 
2217  } // L1
2218 
2219  } // trigger iterator
2220  } // go through triggers
2221 
2222  if(patches.GetSize()<=0) AliInfo(Form("No patch found! for triggers: %s and selected <%s>",
2224  //else printf(">>>>> N patches %d, test %d,first %d, last %d\n",patches.GetSize(), nPatch, patches.At(0), patches.At(patches.GetSize()-1));
2225 
2226  return patches;
2227 }
2228 
2229 //____________________________________________________________
2233 //____________________________________________________________
2235 {
2236  // Init info from previous event
2237  fTriggerClusterIndex = -1;
2238  fTriggerClusterId = -1;
2239  fTriggerClusterBC = -10000;
2240  fIsExoticEvent = kFALSE;
2241  fIsBadCellEvent = kFALSE;
2242  fIsBadMaxCellEvent = kFALSE;
2243 
2244  fIsTriggerMatch = kFALSE;
2245  fIsTriggerMatchOpenCut[0] = kFALSE;
2246  fIsTriggerMatchOpenCut[1] = kFALSE;
2247  fIsTriggerMatchOpenCut[2] = kFALSE;
2248 
2249  // Do only analysis for triggered events
2250  if(!IsEventEMCALL1() && !IsEventEMCALL0())
2251  {
2252  fTriggerClusterBC = 0;
2253  return;
2254  }
2255 
2256  //printf("***** Try to match trigger to cluster %d **** L0 %d, L1 %d\n",fTriggerPatchClusterMatch,IsEventEMCALL0(),IsEventEMCALL1());
2257 
2258  //Recover the list of clusters
2259  TClonesArray * clusterList = 0;
2260  if (fInputEvent->FindListObject(fEMCALClustersListName))
2261  {
2262  clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
2263  }
2264  else if(fOutputEvent)
2265  {
2266  clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
2267  }
2268 
2269  // Get number of clusters and of trigger patches
2270  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2271  if(clusterList)
2272  nclusters = clusterList->GetEntriesFast();
2273 
2274  Int_t nPatch = patches.GetSize();
2275  Float_t exoDiffTime = GetCaloUtils()->GetEMCALRecoUtils()->GetExoticCellDiffTimeCut();
2276 
2277  //Init some variables used in the cluster loop
2278  Float_t tofPatchMax = 100000;
2279  Float_t ePatchMax =-1;
2280 
2281  Float_t tofMax = 100000;
2282  Float_t eMax =-1;
2283 
2284  Int_t clusMax =-1;
2285  Int_t idclusMax =-1;
2286  Bool_t badClMax = kFALSE;
2287  Bool_t badCeMax = kFALSE;
2288  Bool_t exoMax = kFALSE;
2289 // Int_t absIdMaxTrig= -1;
2290  Int_t absIdMaxMax = -1;
2291 
2292  Int_t nOfHighECl = 0 ;
2293 
2294  //
2295  // Check what is the trigger threshold
2296  // set minimu energym of candidate for trigger cluster
2297  //
2299 
2300  Float_t triggerThreshold = fTriggerL1EventThreshold;
2301  if(IsEventEMCALL0()) triggerThreshold = fTriggerL0EventThreshold;
2302  Float_t minE = triggerThreshold / 2.;
2303 
2304  // This method is not really suitable for JET trigger
2305  // but in case, reduce the energy cut since we do not trigger on high energy particle
2306  if(IsEventEMCALL1Jet() || minE < 1) minE = 1;
2307 
2308  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));
2309 
2310  //
2311  // Loop on the clusters, check if there is any that falls into one of the patches
2312  //
2313  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2314  {
2315  AliVCluster * clus = 0;
2316  if(clusterList) clus = (AliVCluster*) clusterList->At(iclus);
2317  else clus = fInputEvent->GetCaloCluster(iclus);
2318 
2319  if ( !clus ) continue ;
2320 
2321  if ( !clus->IsEMCAL() ) continue ;
2322 
2323  //Skip clusters with too low energy to be triggering
2324  if ( clus->E() < minE ) continue ;
2325 
2326  Float_t frac = -1;
2327  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fInputEvent->GetEMCALCells(), clus,frac);
2328 
2329  Bool_t badCluster = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),
2330  clus->GetCellsAbsId(),clus->GetNCells());
2331  UShort_t cellMax[] = {(UShort_t) absIdMax};
2332  Bool_t badCell = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),cellMax,1);
2333 
2334  // if cell is bad, it can happen that time calibration is not available,
2335  // when calculating if it is exotic, this can make it to be exotic by default
2336  // open it temporarily for this cluster
2337  if(badCell)
2338  GetCaloUtils()->GetEMCALRecoUtils()->SetExoticCellDiffTimeCut(10000000);
2339 
2340  Bool_t exotic = GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCluster(clus, fInputEvent->GetEMCALCells());
2341 
2342  // Set back the time cut on exotics
2343  if(badCell)
2344  GetCaloUtils()->GetEMCALRecoUtils()->SetExoticCellDiffTimeCut(exoDiffTime);
2345 
2346  // Energy threshold for exotic Ecross typically at 4 GeV,
2347  // for lower energy, check that there are more than 1 cell in the cluster
2348  if(!exotic && clus->GetNCells() < 2) exotic = kTRUE;
2349 
2350  Float_t energy = clus->E();
2351  Int_t idclus = clus->GetID();
2352 
2353  Double_t tof = clus->GetTOF();
2354  if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn() && fTriggerClusterTimeRecal){
2355  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2356  //additional L1 phase shift
2357  if(GetCaloUtils()->GetEMCALRecoUtils()->IsL1PhaseInTimeRecalibrationOn()){
2358  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absIdMax), fInputEvent->GetBunchCrossNumber(), tof);
2359  }
2360  }
2361  tof *=1.e9;
2362 
2363  //printf("cluster %d, ID %d, E %2.2f, tof %2.2f, AbsId max %d, exotic %d, bad Cluster %d, bad Cell %d\n",
2364  // iclus,idclus, energy,tof,absIdMax, exotic, badCluster,badCell);
2365 
2366  // Find the highest energy cluster, avobe trigger threshold
2367  // in the event in case no match to trigger is found
2368  if( energy > eMax )
2369  {
2370  tofMax = tof;
2371  eMax = energy;
2372  badClMax = badCluster;
2373  badCeMax = badCell;
2374  exoMax = exotic;
2375  clusMax = iclus;
2376  idclusMax = idclus;
2377  absIdMaxMax = absIdMax;
2378  }
2379 
2380  // count the good clusters in the event avobe the trigger threshold
2381  // to check the exotic events
2382  if(!badCluster && !exotic)
2383  nOfHighECl++;
2384 
2385  // Find match to trigger
2386  if(fTriggerPatchClusterMatch && nPatch>0)
2387  {
2388  for(Int_t iabsId =0; iabsId < nPatch; iabsId++)
2389  {
2390  Int_t absIDCell[4];
2391  GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patches.At(iabsId), absIDCell);
2392  //if(tof > 75 ) printf("E %2.2f TOF %2.2f Trigger patch %d, cells : %d, %d, %d, %d\n",
2393  // clus->E(),tof,patches.At(iabsId), absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
2394 
2395  for(Int_t ipatch = 0; ipatch < 4; ipatch++)
2396  {
2397  if(absIdMax == absIDCell[ipatch])
2398  {
2399  //printf("*** Patches : absId %d, E %2.1f, tof %f \n",absIdMax,clus->E(), tof);
2400  if(energy > ePatchMax)
2401  {
2402  tofPatchMax = tof;
2403  ePatchMax = energy;
2404  fIsBadCellEvent = badCluster;
2405  fIsBadMaxCellEvent = badCell;
2406  fIsExoticEvent = exotic;
2407  fTriggerClusterIndex = iclus;
2408  fTriggerClusterId = idclus;
2409  fIsTriggerMatch = kTRUE;
2410 // absIdMaxTrig = absIdMax;
2411  }
2412  }
2413  }// cell patch loop
2414  }// trigger patch loop
2415  } // Do trigger patch matching
2416 
2417  }// Cluster loop
2418 
2419  // If there was no match, assign as trigger
2420  // the highest energy cluster in the event
2421  if(!fIsTriggerMatch)
2422  {
2423  tofPatchMax = tofMax;
2424  ePatchMax = eMax;
2425  fIsBadCellEvent = badClMax;
2426  fIsBadMaxCellEvent = badCeMax;
2427  fIsExoticEvent = exoMax;
2428  fTriggerClusterIndex = clusMax;
2429  fTriggerClusterId = idclusMax;
2430  }
2431 
2432  Double_t tofPatchMaxUS = TMath::Abs(tofPatchMax);
2433 
2434  if (tofPatchMaxUS < 28 ) fTriggerClusterBC = 0 ;
2435  else if(tofPatchMaxUS < 75 ) fTriggerClusterBC = 1 ;
2436  else if(tofPatchMaxUS < 125) fTriggerClusterBC = 2 ;
2437  else if(tofPatchMaxUS < 175) fTriggerClusterBC = 3 ;
2438  else if(tofPatchMaxUS < 225) fTriggerClusterBC = 4 ;
2439  else if(tofPatchMaxUS < 275) fTriggerClusterBC = 5 ;
2440  else
2441  {
2442  //printf("AliCaloTrackReader::MatchTriggerCluster() - Large BC - tof %2.3f - Index %d\n",tofPatchMaxUS,fTriggerClusterIndex);
2444  else
2445  {
2446  fTriggerClusterIndex = -2;
2447  fTriggerClusterId = -2;
2448  }
2449  }
2450 
2451  if(tofPatchMax < 0) fTriggerClusterBC*=-1;
2452 
2453 
2454  // 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",
2455  // fTriggerClusterIndex, fTriggerClusterId,ePatchMax, tofPatchMax,
2456  // fTriggerClusterBC, fIsBadCellEvent,fIsBadMaxCellEvent,fIsExoticEvent, fIsTriggerMatch, nOfHighECl,absIdMaxMax);
2457  //
2458  // 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",
2459  // clusMax, idclusMax, eMax,tofMax, badClMax, badCeMax,exoMax);
2460 
2461  //Redo matching but open cuts
2462  if(!fIsTriggerMatch && fTriggerClusterId >= 0)
2463  {
2464  // Open time patch time
2465  TArrayI patchOpen = GetTriggerPatches(7,10);
2466 
2467  Int_t patchAbsIdOpenTime = -1;
2468  for(Int_t iabsId =0; iabsId < patchOpen.GetSize(); iabsId++)
2469  {
2470  Int_t absIDCell[4];
2471  patchAbsIdOpenTime = patchOpen.At(iabsId);
2472  GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patchAbsIdOpenTime, absIDCell);
2473  //if(tof > 75 ) printf("E %2.2f TOF %2.2f Trigger patch %d, cells : %d, %d, %d, %d\n",
2474  // clus->E(),tof,patches.At(iabsId), absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
2475 
2476  for(Int_t ipatch = 0; ipatch < 4; ipatch++)
2477  {
2478  if(absIdMaxMax == absIDCell[ipatch])
2479  {
2480  fIsTriggerMatchOpenCut[0] = kTRUE;
2481  break;
2482  }
2483  }// cell patch loop
2484  }// trigger patch loop
2485 
2486  // Check neighbour patches
2487  Int_t patchAbsId = -1;
2488  Int_t globalCol = -1;
2489  Int_t globalRow = -1;
2490  GetCaloUtils()->GetEMCALGeometry()->GetFastORIndexFromCellIndex(absIdMaxMax, patchAbsId);
2491  GetCaloUtils()->GetEMCALGeometry()->GetPositionInEMCALFromAbsFastORIndex(patchAbsId,globalCol,globalRow);
2492 
2493  // Check patches with strict time cut
2494  Int_t patchAbsIdNeigh = -1;
2495  for(Int_t icol = globalCol-1; icol <= globalCol+1; icol++)
2496  {
2497  if(icol < 0 || icol > 47) continue;
2498 
2499  for(Int_t irow = globalRow; irow <= globalRow+1; irow++)
2500  {
2501  if(irow < 0 || irow > 63) continue;
2502 
2503  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(icol, irow, patchAbsIdNeigh);
2504 
2505  if ( patchAbsIdNeigh < 0 ) continue;
2506 
2507  for(Int_t iabsId =0; iabsId < patches.GetSize(); iabsId++)
2508  {
2509  if(patchAbsIdNeigh == patches.At(iabsId))
2510  {
2511  fIsTriggerMatchOpenCut[1] = kTRUE;
2512  break;
2513  }
2514  }// trigger patch loop
2515 
2516  }// row
2517  }// col
2518 
2519  // Check patches with open time cut
2520  Int_t patchAbsIdNeighOpenTime = -1;
2521  for(Int_t icol = globalCol-1; icol <= globalCol+1; icol++)
2522  {
2523  if(icol < 0 || icol > 47) continue;
2524 
2525  for(Int_t irow = globalRow; irow <= globalRow+1; irow++)
2526  {
2527  if(irow < 0 || irow > 63) continue;
2528 
2529  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(icol, irow, patchAbsIdNeighOpenTime);
2530 
2531  if ( patchAbsIdNeighOpenTime < 0 ) continue;
2532 
2533  for(Int_t iabsId =0; iabsId < patchOpen.GetSize(); iabsId++)
2534  {
2535  if(patchAbsIdNeighOpenTime == patchOpen.At(iabsId))
2536  {
2537  fIsTriggerMatchOpenCut[2] = kTRUE;
2538  break;
2539  }
2540  }// trigger patch loop
2541 
2542  }// row
2543  }// col
2544 
2545  // printf("No match, new match: Open time %d-%d, open Neigh %d-%d, both open %d-%d\n",fIsTriggerMatchOpenCut[0],patchAbsIdOpenTime,
2546  // fIsTriggerMatchOpenCut[1],patchAbsIdNeigh,
2547  // fIsTriggerMatchOpenCut[2],patchAbsIdNeighOpenTime);
2548 
2549  patchOpen.Reset();
2550 
2551  }// No trigger match found
2552  //printf("Trigger BC %d, Id %d, Index %d\n",fTriggerClusterBC,fTriggerClusterId,fTriggerClusterIndex);
2553 }
2554 
2555 //_________________________________________________________
2559 //_________________________________________________________
2561 {
2563  {
2564  // get object pointer
2565  AliVCaloTrigger *caloTrigger = GetInputEvent()->GetCaloTrigger( "EMCAL" );
2566 
2567  if ( fBitEGA == 6 )
2568  {
2569  if (IsEventEMCALL1Gamma1()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(1);
2570  else if(IsEventEMCALL1Gamma2()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(3);
2571  else if(IsEventEMCALL1Jet1 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(0);
2572  else if(IsEventEMCALL1Jet2 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(2);
2573 
2574  // printf("L1 trigger Threshold Jet1 %f, Gamma1 %f, Jet2 %f, Gamma2 %f\n",
2575  // 0.07874*caloTrigger->GetL1Threshold(0),
2576  // 0.07874*caloTrigger->GetL1Threshold(1),
2577  // 0.07874*caloTrigger->GetL1Threshold(2),
2578  // 0.07874*caloTrigger->GetL1Threshold(3));
2579  }
2580  else
2581  {
2582  // Old AOD data format, in such case, order of thresholds different!!!
2583  if (IsEventEMCALL1Gamma1()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(0);
2584  else if(IsEventEMCALL1Gamma2()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(2);
2585  else if(IsEventEMCALL1Jet1 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(1);
2586  else if(IsEventEMCALL1Jet2 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(3);
2587 
2588  // printf("L1 trigger Threshold Jet1 %f, Gamma1 %f, Jet2 %f, Gamma2 %f\n",
2589  // 0.07874*caloTrigger->GetL1Threshold(1),
2590  // 0.07874*caloTrigger->GetL1Threshold(0),
2591  // 0.07874*caloTrigger->GetL1Threshold(3),
2592  // 0.07874*caloTrigger->GetL1Threshold(2));
2593  }
2594  }
2595 
2596  // Set L0 threshold, if not set by user
2598  {
2599  // Revise for periods > LHC11d
2600  Int_t runNumber = fInputEvent->GetRunNumber();
2601  if (runNumber < 146861) fTriggerL0EventThreshold = 3. ;
2602  else if(runNumber < 154000) fTriggerL0EventThreshold = 4. ;
2603  else if(runNumber < 165000) fTriggerL0EventThreshold = 5.5;
2604  else fTriggerL0EventThreshold = 2 ;
2605  }
2606 }
2607 
2608 //________________________________________________________
2610 //________________________________________________________
2611 void AliCaloTrackReader::Print(const Option_t * opt) const
2612 {
2613  if(! opt)
2614  return;
2615 
2616  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
2617  printf("Task name : %s\n", fTaskName.Data()) ;
2618  printf("Data type : %d\n", fDataType) ;
2619  printf("CTS Min pT : %2.1f GeV/c\n", fCTSPtMin) ;
2620  printf("EMCAL Min pT : %2.1f GeV/c\n", fEMCALPtMin) ;
2621  printf("PHOS Min pT : %2.1f GeV/c\n", fPHOSPtMin) ;
2622  printf("CTS Max pT : %2.1f GeV/c\n", fCTSPtMax) ;
2623  printf("EMCAL Max pT : %2.1f GeV/c\n", fEMCALPtMax) ;
2624  printf("PHOS Max pT : %2.1f GeV/c\n", fPHOSPtMax) ;
2625  printf("EMCAL Time Cut: %3.1f < TOF < %3.1f\n", fEMCALTimeCutMin, fEMCALTimeCutMax);
2626  printf("Use CTS = %d\n", fFillCTS) ;
2627  printf("Use EMCAL = %d\n", fFillEMCAL) ;
2628  printf("Use DCAL = %d\n", fFillDCAL) ;
2629  printf("Use PHOS = %d\n", fFillPHOS) ;
2630  printf("Use EMCAL Cells = %d\n", fFillEMCALCells) ;
2631  printf("Use PHOS Cells = %d\n", fFillPHOSCells) ;
2632  printf("Track status = %d\n", (Int_t) fTrackStatus) ;
2633  //printf("AODs Track filter mask = %d or hybrid %d (if filter bit comp %d), select : SPD hit %d, primary %d\n",
2634  // (Int_t) fTrackFilterMask, fSelectHybridTracks, (Int_t) fTrackFilterMaskComplementary, fSelectSPDHitTracks,fSelectPrimaryTracks) ;
2635  printf("Track Mult Eta Cut = %d\n", (Int_t) fTrackMultEtaCut) ;
2636  printf("Write delta AOD = %d\n", fWriteOutputDeltaAOD) ;
2637  printf("Recalculate Clusters = %d, E linearity = %d\n", fRecalculateClusters, fCorrectELinearity) ;
2638 
2639  printf("Use Triggers selected in SE base class %d; If not what Trigger Mask? %d; MB Trigger Mask for mixed %d \n",
2641 
2643  printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
2644 
2646  printf("Compare cluster pt and pt hard to accept event, factor = %2.2f",fPtHardAndClusterPtFactor);
2647 
2648  printf("Read Kine from, stack? %d, AOD ? %d \n", fReadStack, fReadAODMCParticles) ;
2649  printf("Delta AOD File Name = %s\n", fDeltaAODFileName.Data()) ;
2650  printf("Centrality: Class %s, Option %d, Bin [%d,%d] \n", fCentralityClass.Data(),fCentralityOpt,fCentralityBin[0], fCentralityBin[1]) ;
2651 
2652  printf(" \n") ;
2653 }
2654 
2655 //__________________________________________
2659 //__________________________________________
2661 {
2662  // Count number of cells with energy larger than 0.1 in SM3, cut on this number
2663  Int_t ncellsSM3 = 0;
2664  for(Int_t icell = 0; icell < fInputEvent->GetEMCALCells()->GetNumberOfCells(); icell++)
2665  {
2666  Int_t absID = fInputEvent->GetEMCALCells()->GetCellNumber(icell);
2667  Int_t sm = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
2668  if(fInputEvent->GetEMCALCells()->GetAmplitude(icell) > 0.1 && sm==3) ncellsSM3++;
2669  }
2670 
2671  Int_t ncellcut = 21;
2672  if(GetFiredTriggerClasses().Contains("EMC")) ncellcut = 35;
2673 
2674  if(ncellsSM3 >= ncellcut)
2675  {
2676  AliDebug(1,Form("Reject event with ncells in SM3 %d, cut %d, trig %s",
2677  ncellsSM3,ncellcut,GetFiredTriggerClasses().Data()));
2678  return kTRUE;
2679  }
2680 
2681  return kFALSE;
2682 }
2683 
2684 //_________________________________________________________
2688 //_________________________________________________________
2690 {
2691  if(label < 0) return ;
2692 
2693  AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fInputEvent) ;
2694  if(!evt) return ;
2695 
2696  TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
2697  if(!arr) return ;
2698 
2699  if(label < arr->GetEntriesFast())
2700  {
2701  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
2702  if(!particle) return ;
2703 
2704  if(label == particle->Label()) return ; // label already OK
2705  //else printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
2706  }
2707  //else printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
2708 
2709  // loop on the particles list and check if there is one with the same label
2710  for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
2711  {
2712  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
2713  if(!particle) continue ;
2714 
2715  if(label == particle->Label())
2716  {
2717  label = ind;
2718  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
2719  return;
2720  }
2721  }
2722 
2723  label = -1;
2724 
2725  //printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label not found set to -1 \n");
2726 }
2727 
2728 //___________________________________
2730 //___________________________________
2732 {
2733  if(fCTSTracks) fCTSTracks -> Clear();
2734  if(fEMCALClusters) fEMCALClusters -> Clear("C");
2735  if(fPHOSClusters) fPHOSClusters -> Clear("C");
2736 
2737  fV0ADC[0] = 0; fV0ADC[1] = 0;
2738  fV0Mul[0] = 0; fV0Mul[1] = 0;
2739 
2740  if(fNonStandardJets) fNonStandardJets -> Clear("C");
2741  fBackgroundJets->Reset();
2742 }
2743 
2744 //___________________________________________
2748 //___________________________________________
2750 {
2751  fEventTrigMinBias = kFALSE;
2752  fEventTrigCentral = kFALSE;
2753  fEventTrigSemiCentral = kFALSE;
2754  fEventTrigEMCALL0 = kFALSE;
2755  fEventTrigEMCALL1Gamma1 = kFALSE;
2756  fEventTrigEMCALL1Gamma2 = kFALSE;
2757  fEventTrigEMCALL1Jet1 = kFALSE;
2758  fEventTrigEMCALL1Jet2 = kFALSE;
2759 
2760  AliDebug(1,Form("Select trigger mask bit %d - Trigger Event %s - Select <%s>",
2762 
2763  if(fEventTriggerMask <=0 )// in case no mask set
2764  {
2765  // EMC triggered event? Which type?
2766  if( GetFiredTriggerClasses().Contains("-B-") || GetFiredTriggerClasses().Contains("-S-") || GetFiredTriggerClasses().Contains("-I-") )
2767  {
2768  if ( GetFiredTriggerClasses().Contains("EGA" ) ||
2769  GetFiredTriggerClasses().Contains("EG1" ) )
2770  {
2771  fEventTrigEMCALL1Gamma1 = kTRUE;
2772  if( GetFiredTriggerClasses().Contains("EG1" ) && !fFiredTriggerClassName.Contains("EG1") ) fEventTrigEMCALL1Gamma1 = kFALSE;
2773  }
2774  else if( GetFiredTriggerClasses().Contains("EG2" ) )
2775  {
2776  fEventTrigEMCALL1Gamma2 = kTRUE;
2777  if( !fFiredTriggerClassName.Contains("EG2") ) fEventTrigEMCALL1Gamma2 = kFALSE;
2778  }
2779  else if( GetFiredTriggerClasses().Contains("EJE" ) ||
2780  GetFiredTriggerClasses().Contains("EJ1" ) )
2781  {
2782  fEventTrigEMCALL1Jet1 = kTRUE;
2783  if( GetFiredTriggerClasses().Contains("EJ1" ) && !fFiredTriggerClassName.Contains("EJ1") )
2784  fEventTrigEMCALL1Jet1 = kFALSE;
2785  }
2786  else if( GetFiredTriggerClasses().Contains("EJ2" ) )
2787  {
2788  fEventTrigEMCALL1Jet2 = kTRUE;
2789  if( !fFiredTriggerClassName.Contains("EJ2") ) fEventTrigEMCALL1Jet2 = kFALSE;
2790  }
2791  else if( GetFiredTriggerClasses().Contains("CEMC") &&
2792  !GetFiredTriggerClasses().Contains("EGA" ) &&
2793  !GetFiredTriggerClasses().Contains("EJE" ) &&
2794  !GetFiredTriggerClasses().Contains("EG1" ) &&
2795  !GetFiredTriggerClasses().Contains("EJ1" ) &&
2796  !GetFiredTriggerClasses().Contains("EG2" ) &&
2797  !GetFiredTriggerClasses().Contains("EJ2" ) )
2798  {
2799  fEventTrigEMCALL0 = kTRUE;
2800  }
2801 
2802  //Min bias event trigger?
2803  if (GetFiredTriggerClasses().Contains("CCENT_R2-B-NOPF-ALLNOTRD"))
2804  {
2805  fEventTrigCentral = kTRUE;
2806  }
2807  else if(GetFiredTriggerClasses().Contains("CSEMI_R1-B-NOPF-ALLNOTRD"))
2808  {
2809  fEventTrigSemiCentral = kTRUE;
2810  }
2811  else if((GetFiredTriggerClasses().Contains("CINT") || GetFiredTriggerClasses().Contains("CPBI2_B1") ) &&
2812  GetFiredTriggerClasses().Contains("-NOPF-ALLNOTRD") )
2813  {
2814  fEventTrigMinBias = kTRUE;
2815  }
2816  }
2817  }
2818  else
2819  {
2820  // EMC L1 Gamma
2821  if ( fEventTriggerMask & AliVEvent::kEMCEGA )
2822  {
2823  //printf("EGA trigger bit\n");
2824  if (GetFiredTriggerClasses().Contains("EG"))
2825  {
2826  if (GetFiredTriggerClasses().Contains("EGA")) fEventTrigEMCALL1Gamma1 = kTRUE;
2827  else
2828  {
2829  if(GetFiredTriggerClasses().Contains("EG1")) fEventTrigEMCALL1Gamma1 = kTRUE;
2830  if(GetFiredTriggerClasses().Contains("EG2")) fEventTrigEMCALL1Gamma2 = kTRUE;
2831  }
2832  }
2833  }
2834  // EMC L1 Jet
2835  else if( fEventTriggerMask & AliVEvent::kEMCEJE )
2836  {
2837  //printf("EGA trigger bit\n");
2838  if (GetFiredTriggerClasses().Contains("EJ"))
2839  {
2840  if (GetFiredTriggerClasses().Contains("EJE")) fEventTrigEMCALL1Jet1 = kTRUE;
2841  else
2842  {
2843  if(GetFiredTriggerClasses().Contains("EJ1")) fEventTrigEMCALL1Jet1 = kTRUE;
2844  if(GetFiredTriggerClasses().Contains("EJ2")) fEventTrigEMCALL1Jet2 = kTRUE;
2845  }
2846  }
2847  }
2848  // EMC L0
2849  else if((fEventTriggerMask & AliVEvent::kEMC7) ||
2850  (fEventTriggerMask & AliVEvent::kEMC1) )
2851  {
2852  //printf("L0 trigger bit\n");
2853  fEventTrigEMCALL0 = kTRUE;
2854  }
2855  // Min Bias Pb-Pb
2856  else if( fEventTriggerMask & AliVEvent::kCentral )
2857  {
2858  //printf("MB semi central trigger bit\n");
2859  fEventTrigSemiCentral = kTRUE;
2860  }
2861  // Min Bias Pb-Pb
2862  else if( fEventTriggerMask & AliVEvent::kSemiCentral )
2863  {
2864  //printf("MB central trigger bit\n");
2865  fEventTrigCentral = kTRUE;
2866  }
2867  // Min Bias pp, PbPb, pPb
2868  else if((fEventTriggerMask & AliVEvent::kMB ) ||
2869  (fEventTriggerMask & AliVEvent::kINT7) ||
2870  (fEventTriggerMask & AliVEvent::kINT8) ||
2871  (fEventTriggerMask & AliVEvent::kAnyINT) )
2872  {
2873  //printf("MB trigger bit\n");
2874  fEventTrigMinBias = kTRUE;
2875  }
2876  }
2877 
2878  AliDebug(1,Form("Event bits: \n \t MB %d, Cen %d, Sem %d, L0 %d, L1G1 %d, L1G2 %d, L1J1 %d, L1J2 %d",
2882 
2883  // L1 trigger bit
2884  if( fBitEGA == 0 && fBitEJE == 0 )
2885  {
2886  // Init the trigger bit once, correct depending on AliESD(AOD)CaloTrigger header version
2887 
2888  // Simpler way to do it ...
2889 // if( fInputEvent->GetRunNumber() < 172000 )
2890 // reader->SetEventTriggerL1Bit(4,5); // current LHC11 data
2891 // else
2892 // reader->SetEventTriggerL1Bit(6,8); // LHC12-13 data
2893 
2894  // Old values
2895  fBitEGA = 4;
2896  fBitEJE = 5;
2897 
2898  TFile* file = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
2899 
2900  const TList *clist = file->GetStreamerInfoCache();
2901 
2902  if(clist)
2903  {
2904  TStreamerInfo *cinfo = (TStreamerInfo*)clist->FindObject("AliESDCaloTrigger");
2905  Int_t verid = 5; // newer ESD header version
2906  if(!cinfo)
2907  {
2908  cinfo = (TStreamerInfo*)clist->FindObject("AliAODCaloTrigger");
2909  verid = 3; // newer AOD header version
2910  }
2911 
2912  if(cinfo)
2913  {
2914  Int_t classversionid = cinfo->GetClassVersion();
2915  //printf("********* Header class version %d *********** \n",classversionid);
2916 
2917  if (classversionid >= verid)
2918  {
2919  fBitEGA = 6;
2920  fBitEJE = 8;
2921  }
2922  } else AliInfo("AliCaloTrackReader()::SetEventTriggerBit() - Streamer info for trigger class not available, bit not changed");
2923  } else AliInfo("AliCaloTrackReader::SetEventTriggerBit() - Streamer list not available!, bit not changed");
2924 
2925  } // set once the EJE, EGA trigger bit
2926 }
2927 
2928 //____________________________________________________________
2931 //____________________________________________________________
2932 void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
2933 {
2934  fInputEvent = input;
2935  fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ;
2936  if (fMixedEvent)
2937  fNMixedEvent = fMixedEvent->GetNumberOfEvents() ;
2938 
2939  //Delete previous vertex
2940  if(fVertex)
2941  {
2942  for (Int_t i = 0; i < fNMixedEvent; i++)
2943  {
2944  delete [] fVertex[i] ;
2945  }
2946  delete [] fVertex ;
2947  }
2948 
2949  fVertex = new Double_t*[fNMixedEvent] ;
2950  for (Int_t i = 0; i < fNMixedEvent; i++)
2951  {
2952  fVertex[i] = new Double_t[3] ;
2953  fVertex[i][0] = 0.0 ;
2954  fVertex[i][1] = 0.0 ;
2955  fVertex[i][2] = 0.0 ;
2956  }
2957 }
2958 
2959 
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)
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
virtual void FillInputVZERO()
Int_t fV0ADC[2]
Integrated V0 signal.
Bool_t fComparePtHardAndClusterPt
In MonteCarlo, jet events, reject events with too large cluster energy.
AliAnaWeights * fWeightUtils
Pointer to AliAnaWeights.
AliCalorimeterUtils * GetCaloUtils() const
Float_t fTimeStampEventFracMin
Minimum value of time stamp fraction event.
Bool_t fReadAODMCParticles
Access kine information from filtered AOD MC particles.
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:50
Bool_t fAcceptOnlyHIJINGLabels
Select clusters or tracks that where generated by HIJING, reject other generators in case of cocktail...
TObjArray * fPHOSClusters
Temporal array with PHOS CaloClusters.
Double_t fTrackDCACut[3]
Remove tracks with DCA larger than cut, parameters of function stored here.
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.
Bool_t ClusterContainsBadChannel(Int_t calo, UShort_t *cellList, Int_t nCells)
virtual AliHeader * GetHeader() const
Bool_t IsEventEMCALL1Jet1() const
Bool_t IsEventEMCALL1Gamma1() const
AliEMCALRecoUtils * GetEMCALRecoUtils() const
Bool_t ReadAODMCParticles() const
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.
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.
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:31
virtual void SetInputEvent(AliVEvent *input)
Double_t fTrackTimeCutMax
Remove tracks with time larger than this value, in ns.
Int_t fTriggerClusterBC
Event triggered by a cluster in BC -5 0 to 5.
virtual AliMultSelection * GetMultSelCen() const
TString fEMCALClustersListName
Alternative list of clusters produced elsewhere and not from InputEvent.
void RecalculateClusterTrackMatching(AliVEvent *event, TObjArray *clusterArray=0x0)
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.
Int_t GetDebug() const
Definition: AliAnaWeights.h:82
UInt_t fEventTriggerMask
Select this triggerered event.
virtual Int_t GetEventCentrality() const
virtual Bool_t IsHIJINGLabel(Int_t label)
virtual AliGenEventHeader * GetGenEventHeader() const
virtual AliCentrality * GetCentrality() const
AliVCaloCells * fPHOSCells
! Temporal array with PHOS AliVCaloCells.
Bool_t IsEventEMCALL1Gamma2() const
virtual AliMixedEvent * GetMixedEvent() const
Bool_t fFillInputBackgroundJetBranch
Flag to use data from background jets.
void SetTrackEventBC(Int_t bc)
virtual Bool_t ComparePtHardAndClusterPt()
TList * fAODBranchList
List with AOD branches created and needed in analysis.
void RemapMCLabelForAODs(Int_t &label)
Bool_t fEventTrigSemiCentral
Event is AliVEvent::kSemiCentral on its name, it should correspond to PbPb.
TRandom3 fRandom
! Random generator.
Bool_t AcceptDCA(Float_t pt, Float_t dca)
Bool_t IsPileUpFromNotSPDAndNotEMCal() const
Check if event not from pile-up determined neither by SPD nor by EMCal.
Bool_t fFillCTS
Use data from CTS.
virtual void InitParameters()
Initialize the parameters with default.
virtual void GetVertex(Double_t v[3]) const
Bool_t IsRecalibrationOn() const
virtual Bool_t FillInputEvent(Int_t iEntry, const char *currentFileName)
Float_t fPHOSPtMin
pT Threshold on phos clusters.
Bool_t fSelectEmbeddedClusters
Use only simulated clusters that come from embedding.
Bool_t fEventTrigMinBias
Event is min bias on its name, it should correspond to AliVEvent::kMB, AliVEvent::kAnyInt.
Float_t fEMCALParamTimeCutMin[4]
Remove clusters/cells with time smaller than parametrized value, in ns.
Int_t fEventNumber
Event number.
void RecalculateClusterPID(AliVCluster *clu)
TString fTaskName
Name of task that executes the analysis.
Bool_t fRemoveBadTriggerEvents
Remove triggered events because trigger was exotic, bad, or out of BC.
Bool_t IsPileUpFromSPDOrEMCal() const
Check if event is from pile-up determined by SPD or EMCal.
Bool_t IsPileUpFromEMCalAndNotSPD() const
Check if event is from pile-up determined by EMCal, not by SPD.
Bool_t fEventTrigEMCALL1Gamma2
Event is L1-Gamma, threshold 2 on its name, it should correspond kEMCEGA.
Bool_t fEventTrigEMCALL1Jet1
Event is L1-Gamma, threshold 1 on its name, it should correspond kEMCEGA.
Bool_t fFillEMCALCells
Use data from EMCAL.
Int_t fTriggerClusterId
Id of trigger cluster (cluster->GetID()).
Float_t fCTSPtMax
pT Threshold on charged particles.
Double_t fEMCALParamTimeCutMax[4]
Remove clusters/cells with time larger than parametrized value, in ns.
virtual void FillInputEMCALCells()
Connects the array with EMCAL cells and the pointer.
Bool_t fFillPHOSCells
Use data from PHOS.
virtual AliVCaloCells * GetEMCALCells() const
Bool_t fRemoveLEDEvents
Remove events where LED was wrongly firing - EMCAL LHC11a.
AliEMCALGeometry * GetEMCALGeometry() const
Bool_t fIsBadMaxCellEvent
Bad cell triggered event flag, only max energy cell is bad.
void RecalculateClusterShowerShapeParameters(AliVCaloCells *cells, AliVCluster *clu)
TString fFiredTriggerClassName
Name of trigger event type used to do the analysis.
virtual TClonesArray * GetAODMCParticles() const
TString GetFiredTriggerClasses() const
Bool_t ReadStack() const
TObjArray * fEMCALClusters
Temporal array with EMCAL CaloClusters.
Bool_t IsEventEMCALL1Jet() const
Float_t fPHOSPtMax
pT Threshold on phos clusters.
Bool_t fAccessTrackTOF
Access the track TOF, in case of problems when accessing GetTOFBunchCrossing.
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.
Float_t fEMCALPtMax
pT Threshold on emcal clusters.
Int_t fEventType
Set the event species: 7 physics, 0 MC, 8 LED (not useful now)
virtual Bool_t ComparePtHardAndJetPt()
Bool_t IsInFiducialCut(Float_t eta, Float_t phi, Int_t det) const
Bool_t fEventTrigCentral
Event is AliVEvent::kCentral on its name, it should correspond to PbPb.
Float_t fTrackMultEtaCut
Track multiplicity eta cut.
Int_t fCentralityBin[2]
Minimum and maximum value of the centrality for the analysis.
virtual void FillInputPHOSCells()
Connects the array with PHOS cells and the pointer.
Base class for event, clusters and tracks filtering and preparation for the analysis.
TObjArray * fDCALClusters
Temporal array with DCAL CaloClusters, not needed in the normal case, use just EMCal array with DCal ...
Bool_t fUseTrackTimeCut
Do time cut selection.
virtual TString GetEventPlaneMethod() const
TArrayI fAcceptEventsWithBit
Accept events if trigger bit is on.
virtual void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
void MatchTriggerCluster(TArrayI patches)
Bool_t fTimeStampEventSelect
Select events within a fraction of data taking time.
Double_t fTimeStampRunMax
Maximum value of time stamp in run.
Double_t fEMCALTimeCutMin
Remove clusters/cells with time smaller than this value, in ns.
Bool_t fEventTrigEMCALL1Jet2
Event is L1-Gamma, threshold 2 on its name, it should correspond kEMCEGA.
void SetEMCalEventBCcut(Int_t bc)
Int_t fNPileUpClustersCut
Cut to select event as pile-up.
Bool_t fCheckFidCut
Do analysis for clusters in defined region.
virtual ~AliCaloTrackReader()
Destructor.
Bool_t fEventTrigEMCALL0
Event is EMCal L0 on its name, it should correspond to AliVEvent::kEMC7, AliVEvent::kEMC1.
Bool_t IsWeightSettingOn() const
Definition: AliAnaWeights.h:80
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.
energy
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.
AliFiducialCut * fFiducialCut
Acceptance cuts.
Bool_t fRejectEMCalTriggerEventsWith2Tresholds
Reject events EG2 also triggered by EG1 or EJ2 also triggered by EJ1.
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
void CorrectClusterEnergy(AliVCluster *cl)
Correct cluster energy non linearity.
Float_t fCTSPtMin
pT Threshold on charged particles.
virtual AliStack * GetStack() const
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.
Int_t fTrackMult
Track multiplicity.
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
Bool_t fReadStack
Access kine information from stack.
Bool_t fIsTriggerMatchOpenCut[3]
Could not match the event to a trigger patch?, retry opening cuts.
UInt_t fMixEventTriggerMask
Select this triggerered event for mixing, tipically kMB or kAnyINT.
Float_t RadToDeg(Float_t rad) const
TFile * file
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
TObjArray * fCTSTracks
Temporal array with tracks.
Int_t fTriggerPatchTimeWindow[2]
Trigger patch selection window.
AliMCEvent * fMC
! Monte Carlo Event Handler.
Float_t fPtHardAndJetPtFactor
Factor between ptHard and jet pT to reject/accept event.
TString fDeltaAODFileName
Delta AOD file name.
Int_t GetVertexBC() const
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.
Float_t fSmearShowerShapeWidth
Smear shower shape landau function "width" (use in MC).
Int_t fBitEGA
Trigger bit on VCaloTrigger for EGA.
virtual void FillInputBackgroundJets()
Bool_t fDoV0ANDEventSelection
Select events depending on V0AND.
AliAODJetEventBackground * fBackgroundJets
! Background jets.
Bool_t fCorrectELinearity
Correct cluster linearity, always on.
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 fNMCProducedMax
In case of cocktail, select particles in the list with label up to this value.
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.
Float_t fSmearingFunction
Choice of smearing function. 0 no smearing. 1 smearing from Gustavo (Landau center at 0)...
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.
Int_t fNMCProducedMin
In case of cocktail, select particles in the list with label from this value.
Bool_t fUseAliCentrality
Select as centrality estimator AliCentrality (Run1) or AliMultSelection (Run1 and Run2) ...
virtual TObjArray * GetCTSTracks() const
Bool_t fRecalculateVertexBC
Recalculate vertex BC from tracks pointing to vertex.
Float_t RecalibrateClusterEnergy(AliVCluster *cluster, AliVCaloCells *cells)
Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that co...
void SetEMCalEventBC(Int_t bc)
Bool_t reject
Bool_t CheckCellFiducialRegion(AliVCluster *cluster, AliVCaloCells *cells) const
Bool_t IsEventEMCALL0() const