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