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