AliPhysics  ec7afe5 (ec7afe5)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliCaloTrackReader.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // --- ROOT system ---
17 #include <TFile.h>
18 #include <TGeoManager.h>
19 #include <TStreamerInfo.h>
20 
21 // ---- ANALYSIS system ----
22 #include "AliMCEvent.h"
23 #include "AliAODMCHeader.h"
24 #include "AliGenPythiaEventHeader.h"
25 #include "AliGenCocktailEventHeader.h"
26 #include "AliGenHijingEventHeader.h"
27 #include "AliESDEvent.h"
28 #include "AliAODEvent.h"
29 #include "AliVTrack.h"
30 #include "AliVParticle.h"
31 #include "AliMixedEvent.h"
32 //#include "AliTriggerAnalysis.h"
33 #include "AliESDVZERO.h"
34 #include "AliVCaloCells.h"
35 #include "AliAnalysisManager.h"
36 #include "AliInputEventHandler.h"
37 #include "AliAODMCParticle.h"
38 #include "AliStack.h"
39 #include "AliLog.h"
40 #include "AliMultSelection.h"
41 
42 // ---- Detectors ----
43 #include "AliPHOSGeoUtils.h"
44 #include "AliEMCALGeometry.h"
45 #include "AliEMCALRecoUtils.h"
46 
47 // ---- CaloTrackCorr ---
48 #include "AliCalorimeterUtils.h"
49 #include "AliCaloTrackReader.h"
50 
51 // ---- Jets ----
52 #include "AliAODJet.h"
53 #include "AliAODJetEventBackground.h"
54 
58 
59 //________________________________________
61 //________________________________________
63 TObject(), fEventNumber(-1), //fCurrentFileName(""),
64 fDataType(0), fDebug(0),
65 fFiducialCut(0x0), fCheckFidCut(kFALSE),
66 fComparePtHardAndJetPt(0), fPtHardAndJetPtFactor(0),
67 fComparePtHardAndClusterPt(0),fPtHardAndClusterPtFactor(0),
68 fCTSPtMin(0), fEMCALPtMin(0), fPHOSPtMin(0),
69 fCTSPtMax(0), fEMCALPtMax(0), fPHOSPtMax(0),
70 fEMCALBadChMinDist(0), fPHOSBadChMinDist (0),
71 fEMCALNCellsCut(0), fPHOSNCellsCut(0),
72 fUseEMCALTimeCut(1), fUseParamTimeCut(0),
73 fUseTrackTimeCut(0), fAccessTrackTOF(0),
74 fEMCALTimeCutMin(-10000), fEMCALTimeCutMax(10000),
75 fEMCALParamTimeCutMin(), fEMCALParamTimeCutMax(),
76 fTrackTimeCutMin(-10000), fTrackTimeCutMax(10000),
77 fUseTrackDCACut(0),
78 fAODBranchList(0x0),
79 fCTSTracks(0x0), fEMCALClusters(0x0),
80 fDCALClusters(0x0), fPHOSClusters(0x0),
81 fEMCALCells(0x0), fPHOSCells(0x0),
82 fInputEvent(0x0), fOutputEvent(0x0),fMC(0x0),
83 fFillCTS(0), fFillEMCAL(0),
84 fFillDCAL(0), fFillPHOS(0),
85 fFillEMCALCells(0), fFillPHOSCells(0),
86 fRecalculateClusters(kFALSE),fCorrectELinearity(kTRUE),
87 fSelectEmbeddedClusters(kFALSE),
88 fSmearShowerShape(0), fSmearShowerShapeWidth(0), fRandom(),
89 fSmearingFunction(0), fSmearNLMMin(0), fSmearNLMMax(0),
90 fTrackStatus(0), fSelectSPDHitTracks(0),
91 fTrackMultNPtCut(0), fTrackMultEtaCut(0.9),
92 fReadStack(kFALSE), fReadAODMCParticles(kFALSE),
93 fDeltaAODFileName(""), fFiredTriggerClassName(""),
94 
95 fEventTriggerMask(0), fMixEventTriggerMask(0), fEventTriggerAtSE(0),
96 fEventTrigMinBias(0), fEventTrigCentral(0),
97 fEventTrigSemiCentral(0), fEventTrigEMCALL0(0),
98 fEventTrigEMCALL1Gamma1(0), fEventTrigEMCALL1Gamma2(0),
99 fEventTrigEMCALL1Jet1(0), fEventTrigEMCALL1Jet2(0),
100 fBitEGA(0), fBitEJE(0),
101 
102 fEventType(-1),
103 fTaskName(""), fCaloUtils(0x0),
104 fWeightUtils(0x0), fEventWeight(1),
105 fMixedEvent(NULL), fNMixedEvent(0), fVertex(NULL),
106 fListMixedTracksEvents(), fListMixedCaloEvents(),
107 fLastMixedTracksEvent(-1), fLastMixedCaloEvent(-1),
108 fWriteOutputDeltaAOD(kFALSE),
109 fEMCALClustersListName(""), fZvtxCut(0.),
110 fAcceptFastCluster(kFALSE), fRemoveLEDEvents(kFALSE),
111 //Trigger rejection
112 fRemoveBadTriggerEvents(0), fTriggerPatchClusterMatch(1),
113 fTriggerPatchTimeWindow(), fTriggerL0EventThreshold(0),
114 fTriggerL1EventThreshold(0), fTriggerL1EventThresholdFix(0),
115 fTriggerClusterBC(0), fTriggerClusterIndex(0), fTriggerClusterId(0),
116 fIsExoticEvent(0), fIsBadCellEvent(0), fIsBadMaxCellEvent(0),
117 fIsTriggerMatch(0), fIsTriggerMatchOpenCut(),
118 fTriggerClusterTimeRecal(kTRUE), fRemoveUnMatchedTriggers(kTRUE),
119 fDoPileUpEventRejection(kFALSE), fDoV0ANDEventSelection(kFALSE),
120 fDoVertexBCEventSelection(kFALSE),
121 fDoRejectNoTrackEvents(kFALSE),
122 fUseEventsWithPrimaryVertex(kFALSE),
123 //fTriggerAnalysis (0x0),
124 fTimeStampEventSelect(0),
125 fTimeStampEventFracMin(0), fTimeStampEventFracMax(0),
126 fTimeStampRunMin(0), fTimeStampRunMax(0),
127 fNPileUpClusters(-1), fNNonPileUpClusters(-1), fNPileUpClustersCut(3),
128 fVertexBC(-200), fRecalculateVertexBC(0),
129 fUseAliCentrality(0), fCentralityClass(""), fCentralityOpt(0),
130 fEventPlaneMethod(""),
131 fFillInputNonStandardJetBranch(kFALSE),
132 fNonStandardJets(new TClonesArray("AliAODJet",100)), fInputNonStandardJetBranchName("jets"),
133 fFillInputBackgroundJetBranch(kFALSE),
134 fBackgroundJets(0x0),fInputBackgroundJetBranchName("jets"),
135 fAcceptEventsWithBit(0), fRejectEventsWithBit(0), fRejectEMCalTriggerEventsWith2Tresholds(0),
136 fMomentum(), fOutputContainer(0x0), fEnergyHistogramNbins(0),
137 fhNEventsAfterCut(0), fNMCGenerToAccept(0), fMCGenerEventHeaderToAccept("")
138 {
139  for(Int_t i = 0; i < 8; i++) fhEMCALClusterCutsE [i]= 0x0 ;
140  for(Int_t i = 0; i < 7; i++) fhPHOSClusterCutsE [i]= 0x0 ;
141  for(Int_t i = 0; i < 6; i++) fhCTSTrackCutsPt [i]= 0x0 ;
142  for(Int_t j = 0; j < 5; j++) { fMCGenerToAccept [j] = ""; fMCGenerIndexToAccept[j] = -1; }
143 
144  InitParameters();
145 }
146 
147 //_______________________________________
149 //_______________________________________
151 {
152  DeletePointers();
153 }
154 
155 //_______________________________________
158 //_______________________________________
160 {
161  delete fFiducialCut ;
162 
163  if(fAODBranchList)
164  {
165  fAODBranchList->Delete();
166  delete fAODBranchList ;
167  }
168 
169  if(fCTSTracks)
170  {
171  if(fDataType!=kMC)fCTSTracks->Clear() ;
172  else fCTSTracks->Delete() ;
173  delete fCTSTracks ;
174  }
175 
176  if(fEMCALClusters)
177  {
178  if(fDataType!=kMC)fEMCALClusters->Clear("C") ;
179  else fEMCALClusters->Delete() ;
180  delete fEMCALClusters ;
181  }
182 
183  if(fDCALClusters)
184  {
185  if(fDataType!=kMC)fDCALClusters->Clear("C") ;
186  else fDCALClusters->Delete() ;
187  delete fDCALClusters ;
188  }
189 
190  if(fPHOSClusters)
191  {
192  if(fDataType!=kMC)fPHOSClusters->Clear("C") ;
193  else fPHOSClusters->Delete() ;
194  delete fPHOSClusters ;
195  }
196 
197  if(fVertex)
198  {
199  for (Int_t i = 0; i < fNMixedEvent; i++)
200  {
201  delete [] fVertex[i] ;
202 
203  }
204  delete [] fVertex ;
205  }
206 
207  //delete fTriggerAnalysis;
208 
209  if(fNonStandardJets)
210  {
211  if(fDataType!=kMC) fNonStandardJets->Clear("C") ;
212  else fNonStandardJets->Delete() ;
213  delete fNonStandardJets ;
214  }
215  delete fBackgroundJets ;
216 
217  fRejectEventsWithBit.Reset();
218  fAcceptEventsWithBit.Reset();
219 
220  if ( fWeightUtils ) delete fWeightUtils ;
221 
222  // Pointers not owned, done by the analysis frame
223  // if(fInputEvent) delete fInputEvent ;
224  // if(fOutputEvent) delete fOutputEvent ;
225  // if(fMC) delete fMC ;
226  // Pointer not owned, deleted by maker
227  // if (fCaloUtils) delete fCaloUtils ;
228 }
229 
230 //____________________________________________________________
234 //____________________________________________________________
236 {
237  Float_t cut = fTrackDCACut[0]+fTrackDCACut[1]/TMath::Power(pt,fTrackDCACut[2]);
238 
239  if(TMath::Abs(dca) < cut)
240  return kTRUE;
241  else
242  return kFALSE;
243 }
244 
245 //_____________________________________________________
248 //_____________________________________________________
250 {
251  Int_t nAccept = fAcceptEventsWithBit.GetSize();
252 
253  //printf("N accept %d\n", nAccept);
254 
255  if( nAccept <= 0 )
256  return kTRUE ; // accept the event
257 
258  UInt_t trigFired = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
259 
260  for(Int_t ibit = 0; ibit < nAccept; ibit++)
261  {
262  Bool_t accept = (trigFired & fAcceptEventsWithBit.At(ibit));
263 
264  //printf("accept %d, ibit %d, bit %d \n",accept, ibit,fAcceptEventsWithBit.At(ibit));
265  if(accept) return kTRUE ; // accept the event
266  }
267 
268  return kFALSE ; // reject the event
269 }
270 
271 //_____________________________________________________
276 //_____________________________________________________
278 {
279  if( !fMC || fNMCGenerToAccept <= 0 ) return kTRUE;
280 
281  TString genName;
282  Int_t genIndex;
283  genIndex = GetCocktailGeneratorAndIndex(mcLabel, genName);
284  //fMC->GetCocktailGenerator(mcLabel,genName);
285 
286  Bool_t generOK = kFALSE;
287  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
288  {
289  if ( fMCGenerToAccept[ig].Contains(genName) ) generOK = kTRUE;
290 
291  if ( generOK && fMCGenerIndexToAccept[ig] >= 0 && fMCGenerToAccept[ig] != genIndex) generOK = kFALSE;
292  }
293 
294  if ( !generOK ) AliDebug(1, Form("skip label %d, gen %s",mcLabel,genName.Data()) );
295 
296  return generOK;
297 }
298 
299 //_____________________________________________________________________
307 //_____________________________________________________________________
309 {
310  //method that gives the generator for a given particle with label index (or that of the corresponding primary)
311  AliVParticle* mcpart0 = (AliVParticle*) GetMC()->GetTrack(index);
312  Int_t genIndex = -1;
313 
314  if(!mcpart0)
315  {
316  printf("AliMCEvent-BREAK: No valid AliMCParticle at label %i\n",index);
317  return -1;
318  }
319 
320  nameGen = GetGeneratorNameAndIndex(index,genIndex);
321 
322  if(nameGen.Contains("nococktailheader") ) return -1;
323 
324  Int_t lab=index;
325 
326  while(nameGen.IsWhitespace())
327  {
328  AliVParticle* mcpart = (AliVParticle*) GetMC()->GetTrack(lab);
329 
330  if(!mcpart)
331  {
332  printf("AliMCEvent-BREAK: No valid AliMCParticle at label %i\n",lab);
333  break;
334  }
335 
336  Int_t mother=0;
337  mother = mcpart->GetMother();
338 
339  if(mother<0)
340  {
341  printf("AliMCEvent - BREAK: Reached primary particle without valid mother\n");
342  break;
343  }
344 
345  AliVParticle* mcmom = (AliVParticle*) GetMC()->GetTrack(mother);
346  if(!mcmom)
347  {
348  printf("AliMCEvent-BREAK: No valid AliMCParticle mother at label %i\n",mother);
349  break;
350  }
351 
352  lab=mother;
353 
354  nameGen = GetGeneratorNameAndIndex(mother,genIndex);
355  }
356 
357  return genIndex;
358 }
359 //_____________________________________________________________________
367 //_____________________________________________________________________
369 {
370  Int_t nsumpart = GetMC()->GetNumberOfPrimaries();
371 
372  genIndex = -1;
373 
374  TList* lh = GetMC()->GetCocktailList();
375  if(!lh)
376  {
377  TString noheader="nococktailheader";
378  return noheader;
379  }
380 
381  Int_t nh = lh->GetEntries();
382 
383  for (Int_t i = nh-1; i >= 0; i--)
384  {
385  AliGenEventHeader* gh = (AliGenEventHeader*)lh->At(i);
386 
387  TString genname = gh->GetName();
388 
389  Int_t npart=gh->NProduced();
390 
391  if (i == 0) npart = nsumpart;
392 
393  if(index < nsumpart && index >= (nsumpart-npart))
394  {
395  genIndex = i ;
396  return genname;
397  }
398 
399  nsumpart-=npart;
400  }
401 
402  TString empty="";
403  return empty;
404 }
405 
406 
407 //_____________________________________________________
410 //_____________________________________________________
412 {
413  Int_t nReject = fRejectEventsWithBit.GetSize();
414 
415  //printf("N reject %d\n", nReject);
416 
417  if( nReject <= 0 )
418  return kTRUE ; // accept the event
419 
420  UInt_t trigFired = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
421 
422  for(Int_t ibit = 0; ibit < nReject; ibit++)
423  {
424  Bool_t reject = (trigFired & fRejectEventsWithBit.At(ibit));
425 
426  //printf("reject %d, ibit %d, bit %d \n",reject, ibit,fRejectEventsWithBit.At(ibit));
427  if(reject) return kFALSE ; // reject the event
428  }
429 
430  return kTRUE ; // accept the event
431 }
432 
433 //_____________________________________________
437 //_____________________________________________
439 {
440  //-----------------------------------------------------------
441  // Reject events depending on the trigger name
442  //-----------------------------------------------------------
443 
444  AliDebug(1,Form("FiredTriggerClass <%s>, selected class <%s>, compare name %d",
447 
448  if ( fFiredTriggerClassName != "" )
449  {
451  return kFALSE;
452  else
453  AliDebug(1,"Accepted triggered event");
454  }
455 
456  fhNEventsAfterCut->Fill(1.5);
457 
458  //-----------------------------------------------------------
459  // Reject events depending on the event species type
460  //-----------------------------------------------------------
461 
462  // Event types:
463  // kStartOfRun = 1, // START_OF_RUN
464  // kEndOfRun = 2, // END_OF_RUN
465  // kStartOfRunFiles = 3, // START_OF_RUN_FILES
466  // kEndOfRunFiles = 4, // END_OF_RUN_FILES
467  // kStartOfBurst = 5, // START_OF_BURST
468  // kEndOfBurst = 6, // END_OF_BURST
469  // kPhysicsEvent = 7, // PHYSICS_EVENT
470  // kCalibrationEvent = 8, // CALIBRATION_EVENT
471  // kFormatError = 9, // EVENT_FORMAT_ERROR
472  // kStartOfData = 10, // START_OF_DATA
473  // kEndOfData = 11, // END_OF_DATA
474  // kSystemSoftwareTriggerEvent = 12, // SYSTEM_SOFTWARE_TRIGGER_EVENT
475  // kDetectorSoftwareTriggerEvent = 13 // DETECTOR_SOFTWARE_TRIGGER_EVENT
476 
477  Int_t eventType = 0;
478  if(fInputEvent->GetHeader())
479  eventType = ((AliVHeader*)fInputEvent->GetHeader())->GetEventType();
480 
481  AliDebug(3,Form("Event type %d",eventType));
482 
483  // Select only Physics events in data, eventType = 7,
484  // usually done by physics selection
485  // MC not set, eventType =0, I think, so do not apply a value to fEventType by default
486  // LED events have eventType = 8, implemented a selection cut in the past
487  // but not really useful for EMCal data since LED are not reconstructed (unless wrongly assumed as physics)
488 
489  if ( fEventType >= 0 && eventType != fEventType ) return kFALSE ;
490 
491  AliDebug(1,"Pass event species selection");
492 
493  fhNEventsAfterCut->Fill(2.5);
494 
495  //-----------------------------------------------------------------
496  // In case of mixing analysis, select here the trigger of the event
497  //-----------------------------------------------------------------
498 
499  UInt_t isTrigger = kFALSE;
500  UInt_t isMB = kFALSE;
501 
502  if(!fEventTriggerAtSE)
503  {
504  // In case of mixing analysis, accept MB events, not only Trigger
505  // Track and cluster arrays filled for MB in order to create the pool in the corresponding analysis
506  // via de method in the base class FillMixedEventPool()
507 
508  AliAnalysisManager *manager = AliAnalysisManager::GetAnalysisManager();
509  AliInputEventHandler *inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
510 
511  if(!inputHandler) return kFALSE ; // to content coverity
512 
513  isTrigger = inputHandler->IsEventSelected() & fEventTriggerMask;
514  isMB = inputHandler->IsEventSelected() & fMixEventTriggerMask;
515 
516  if(!isTrigger && !isMB) return kFALSE;
517 
518  //printf("Selected triggered event : %s\n",GetFiredTriggerClasses().Data());
519  AliDebug(0,"Pass uninteresting triggered events rejection in case of mixing analysis");
520 
521  fhNEventsAfterCut->Fill(3.5);
522  }
523 
524 
525  //-------------------------------------------------------------------------------------
526  // Reject or accept events depending on the trigger bit
527  //-------------------------------------------------------------------------------------
528 
531 
532  //printf("AliCaloTrackReader::FillInputEvent() - Accept event? %d, Reject event %d? \n",okA,okR);
533 
534  if(!okA || !okR) return kFALSE;
535 
536  AliDebug(1,"Pass event bit rejection");
537 
538  fhNEventsAfterCut->Fill(4.5);
539 
540  //----------------------------------------------------------------------
541  // Do not count events that were likely triggered by an exotic cluster
542  // or out BC cluster
543  //----------------------------------------------------------------------
544 
545  // Set a bit with the event kind, MB, L0, L1 ...
547 
548  // In case of Mixing, avoid checking the triggers in the min bias events
549  if(!fEventTriggerAtSE && (isMB && !isTrigger)) return kTRUE;
550 
552  {
554  {
555  // Reject triggered events when there is coincidence on both EMCal trigger thresholds,
556  // but the requested trigger is the low trigger threshold
557  if(IsEventEMCALL1Jet1 () && IsEventEMCALL1Jet2 () && fFiredTriggerClassName.Contains("EJ2")) return kFALSE;
558  if(IsEventEMCALL1Gamma1() && IsEventEMCALL1Gamma2() && fFiredTriggerClassName.Contains("EG2")) return kFALSE;
559  }
560 
561  //Get Patches that triggered
563 
564  MatchTriggerCluster(patches);
565 
566  patches.Reset();
567 
568  // If requested, remove badly triggeed events, but only when the EMCal trigger bit is set
570  {
571  AliDebug(1,Form("ACCEPT triggered event? \n exotic? %d - bad cell %d - bad Max cell %d - BC %d - Matched %d\n",
573  if (fIsExoticEvent) return kFALSE;
574  else if(fIsBadCellEvent) return kFALSE;
575  else if(fRemoveUnMatchedTriggers && !fIsTriggerMatch) return kFALSE ;
576  else if(fTriggerClusterBC != 0) return kFALSE;
577  AliDebug(1,Form("\t *** YES for %s",GetFiredTriggerClasses().Data()));
578  }
579 
580  AliDebug(1,"Pass EMCal triggered event rejection \n");
581 
582  fhNEventsAfterCut->Fill(5.5);
583  }
584 
585  //-------------------------------------------------------------------------------------
586  // Select events only fired by a certain trigger configuration if it is provided
587  //-------------------------------------------------------------------------------------
588 
589  if (GetFiredTriggerClasses().Contains("FAST") && !GetFiredTriggerClasses().Contains("ALL") && !fAcceptFastCluster)
590  {
591  AliDebug(1,Form("Do not count events from fast cluster, trigger name %s\n",fFiredTriggerClassName.Data()));
592  return kFALSE;
593 
594  fhNEventsAfterCut->Fill(6.5);
595  }
596 
597  //-------------------------------------------------------------------------------------
598  // Reject event if large clusters with large energy
599  // Use only for LHC11a data for the moment, and if input is clusterizer V1 or V1+unfolding
600  // If clusterzer NxN or V2 it does not help
601  //-------------------------------------------------------------------------------------
602 
603  //Int_t run = fInputEvent->GetRunNumber();
604 
605  if ( fRemoveLEDEvents )
606  {
608 
609  if(reject) return kFALSE;
610 
611  AliDebug(1,"Pass LED event rejection");
612 
613  fhNEventsAfterCut->Fill(7.5);
614  } // Remove LED events
615 
616  // All selection criteria passed, accept the event
617  return kTRUE ;
618 }
619 
620 //________________________________________________
625 //________________________________________________
627 {
628  //printf("AliCaloTrackReader::ComparePtHardAndJetPt() - GenHeaderName : %s\n",GetGenEventHeader()->ClassName());
629 
630  if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
631  {
632  TParticle * jet = 0;
633  AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
634  Int_t nTriggerJets = pygeh->NTriggerJets();
635  Float_t ptHard = pygeh->GetPtHard();
636 
637  AliDebug(1,Form("Njets: %d, pT Hard %f",nTriggerJets, ptHard));
638 
639  Float_t tmpjet[]={0,0,0,0};
640  for(Int_t ijet = 0; ijet< nTriggerJets; ijet++)
641  {
642  pygeh->TriggerJet(ijet, tmpjet);
643  jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
644 
645  AliDebug(1,Form("jet %d; pycell jet pT %f",ijet, jet->Pt()));
646 
647  //Compare jet pT and pt Hard
648  if(jet->Pt() > fPtHardAndJetPtFactor * ptHard)
649  {
650  AliInfo(Form("Reject jet event with : pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n",
651  ptHard, jet->Pt(), fPtHardAndJetPtFactor));
652  return kFALSE;
653  }
654  }
655 
656  if(jet) delete jet;
657  }
658 
659  return kTRUE ;
660 }
661 
662 //____________________________________________________
667 //____________________________________________________
669 {
670  if(!strcmp(GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
671  {
672  AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetGenEventHeader();
673  Float_t ptHard = pygeh->GetPtHard();
674 
675  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
676  for (Int_t iclus = 0; iclus < nclusters; iclus++)
677  {
678  AliVCluster * clus = fInputEvent->GetCaloCluster(iclus) ;
679  Float_t ecluster = clus->E();
680 
681  if(ecluster > fPtHardAndClusterPtFactor*ptHard)
682  {
683  AliInfo(Form("Reject : ecluster %2.2f, calo %d, factor %2.2f, ptHard %f",ecluster,clus->GetType(),fPtHardAndClusterPtFactor,ptHard));
684 
685  return kFALSE;
686  }
687  }
688 
689  }
690 
691  return kTRUE ;
692 }
693 
694 //___________________________________________________
697 //___________________________________________________
699 {
700 
701  fhNEventsAfterCut = new TH1I("hNEventsAfterCut", "Number of analyzed events", 19, 0, 19) ;
702  //fhNEventsAfterCut->SetXTitle("Selection");
703  fhNEventsAfterCut->SetYTitle("# events");
704  fhNEventsAfterCut->GetXaxis()->SetBinLabel(1 ,"1=Input");
705  fhNEventsAfterCut->GetXaxis()->SetBinLabel(2 ,"2=Trigger string");
706  fhNEventsAfterCut->GetXaxis()->SetBinLabel(3 ,"3=Event Type");
707  fhNEventsAfterCut->GetXaxis()->SetBinLabel(4 ,"4=Mixing Event");
708  fhNEventsAfterCut->GetXaxis()->SetBinLabel(5 ,"5=Trigger Bit");
709  fhNEventsAfterCut->GetXaxis()->SetBinLabel(6 ,"6=Good EMC Trigger");
710  fhNEventsAfterCut->GetXaxis()->SetBinLabel(7 ,"7=!Fast Cluster");
711  fhNEventsAfterCut->GetXaxis()->SetBinLabel(8 ,"8=!LED");
712  fhNEventsAfterCut->GetXaxis()->SetBinLabel(9 ,"9=Time stamp");
713  fhNEventsAfterCut->GetXaxis()->SetBinLabel(10,"10=Z vertex");
714  fhNEventsAfterCut->GetXaxis()->SetBinLabel(11,"11=Primary vertex");
715  fhNEventsAfterCut->GetXaxis()->SetBinLabel(12,"12=Pile-up");
716  fhNEventsAfterCut->GetXaxis()->SetBinLabel(13,"13=V0AND");
717  fhNEventsAfterCut->GetXaxis()->SetBinLabel(14,"14=Centrality");
718  fhNEventsAfterCut->GetXaxis()->SetBinLabel(15,"15=GenHeader");
719  fhNEventsAfterCut->GetXaxis()->SetBinLabel(16,"16=PtHard-Jet");
720  fhNEventsAfterCut->GetXaxis()->SetBinLabel(17,"17=PtHard-Cluster");
721  fhNEventsAfterCut->GetXaxis()->SetBinLabel(18,"18=Track multi.");
722  fhNEventsAfterCut->GetXaxis()->SetBinLabel(19,"19=TOF BC");
724 
725 
726  if(fFillEMCAL)
727  {
728  for(Int_t i = 0; i < 8; i++)
729  {
730  TString names[] = {"NoCut", "Corrected", "GoodCluster", "NonLinearity",
731  "EnergyAndFidutial", "Time", "NCells", "BadDist"};
732 
733  fhEMCALClusterCutsE[i] = new TH1F(Form("hEMCALReaderClusterCuts_%d_%s",i,names[i].Data()),
734  Form("EMCal %d, %s",i,names[i].Data()),
736  fhEMCALClusterCutsE[i]->SetYTitle("# clusters");
737  fhEMCALClusterCutsE[i]->SetXTitle("#it{E} (GeV)");
739  }
740  }
741 
742  if(fFillPHOS)
743  {
744  for(Int_t i = 0; i < 7; i++)
745  {
746  TString names[] = {"NoCut", "ExcludeCPV", "BorderCut", "FiducialCut", "EnergyCut", "NCells", "BadDist"};
747 
748  fhPHOSClusterCutsE[i] = new TH1F(Form("hPHOSReaderClusterCuts_%d_%s",i,names[i].Data()),
749  Form("PHOS Cut %d, %s",i,names[i].Data()),
751  fhPHOSClusterCutsE[i]->SetYTitle("# clusters");
752  fhPHOSClusterCutsE[i]->SetXTitle("#it{E} (GeV)");
754  }
755  }
756 
757  if(fFillCTS)
758  {
759  for(Int_t i = 0; i < 6; i++)
760  {
761  TString names[] = {"NoCut", "Status", "ESD_AOD", "TOF", "DCA","PtAcceptanceMult"};
762 
763  fhCTSTrackCutsPt[i] = new TH1F(Form("hCTSReaderClusterCuts_%d_%s",i,names[i].Data()),
764  Form("CTS Cut %d, %s",i,names[i].Data()),
766  fhCTSTrackCutsPt[i]->SetYTitle("# tracks");
767  fhCTSTrackCutsPt[i]->SetXTitle("#it{p}_{T} (GeV)");
769  }
770  }
771 
772  return fOutputContainer ;
773 }
774 
775 //_____________________________________________________
777 //_____________________________________________________
779 {
780  TString parList ; //this will be list of parameters used for this analysis.
781 
782  const Int_t buffersize = 255;
783  char onePar[buffersize] ;
784 
785  snprintf(onePar,buffersize,"--- Reader ---:") ;
786  parList+=onePar ;
787  snprintf(onePar,buffersize,"Data type: %d; zvertex cut %2.2f; EMC cluster name: <%s> ",fDataType, fZvtxCut, fEMCALClustersListName.Data()) ;
788  parList+=onePar ;
789  snprintf(onePar,buffersize,"Use detector: EMC %d, DCA %d, PHOS %d, CTS %d, EMCcells %d, PHOScells %d ; ",
791  snprintf(onePar,buffersize,"E-pT window: EMC (%2.1f,%2.1f), PHOS (%2.1f,%2.1f), CTS (%2.1f,%2.1f); ",
793  parList+=onePar ;
794  snprintf(onePar,buffersize,"Dist to bad channel: EMC > %2.1f, PHOS > %2.1f; ",fEMCALBadChMinDist,fPHOSBadChMinDist) ;
795  parList+=onePar ;
796  snprintf(onePar,buffersize,"N cells: EMC > %d, PHOS > %d; ",fEMCALNCellsCut,fPHOSNCellsCut) ;
797  parList+=onePar ;
798  snprintf(onePar,buffersize,"EMC time cut single window (%2.2f,%2.2f); ",fEMCALTimeCutMin,fEMCALTimeCutMax) ;
799  parList+=onePar ;
800  snprintf(onePar,buffersize,"Check: calo fid cut %d; ",fCheckFidCut) ;
801  parList+=onePar ;
802  snprintf(onePar,buffersize,"Track: status %d, SPD hit %d; ",(Int_t) fTrackStatus, fSelectSPDHitTracks) ;
803  parList+=onePar ;
804  snprintf(onePar,buffersize,"multip. eta cut %1.1f; npt cuts %d;",fTrackMultEtaCut, fTrackMultNPtCut) ;
805  parList+=onePar ;
806 
807  if(fUseTrackDCACut)
808  {
809  snprintf(onePar,buffersize,"DCA cut ON, param (%2.4f,%2.4f,%2.4f); ",fTrackDCACut[0],fTrackDCACut[1],fTrackDCACut[2]) ;
810  parList+=onePar ;
811  }
812 
813  snprintf(onePar,buffersize,"Recalculate Clusters = %d, E linearity = %d; ",fRecalculateClusters, fCorrectELinearity) ;
814  parList+=onePar ;
815 
816  snprintf(onePar,buffersize,"SE trigger sel. %d, not? trigger Mask? %d, MB Trigger Mask for mixed %d; ",
818  parList+=onePar ;
819 
820  snprintf(onePar,buffersize,"Select fired trigger %s; Remove Bad trigger event %d, unmatched %d; Accept fastcluster %d, Reject LED %d ",
822  parList+=onePar ;
823 
825  {
826  snprintf(onePar,buffersize,"Accept only labels from: ");
827  parList+=onePar ;
828  for(Int_t i = 0; i< fNMCGenerToAccept; i++)
829  parList+=(fMCGenerToAccept[i]+" ") ;
830 
831  snprintf(onePar,buffersize,"; ");
832  parList+=onePar ;
833  }
834 
836  {
837  snprintf(onePar,buffersize,"EMC M02 smear ON, function %d, param %2.4f, %d<=NLM<=%d; ",
839  parList+=onePar ;
840  }
841 
843  {
844  snprintf(onePar,buffersize,"jet pt / pt hard < %2.1f; ",fPtHardAndJetPtFactor);
845  parList+=onePar ;
846  }
847 
849  {
850  snprintf(onePar,buffersize,"cluster pt / pt hard < %2.2f",fPtHardAndClusterPtFactor);
851  parList+=onePar ;
852  }
853 
854  snprintf(onePar,buffersize,"Centrality: Class %s, Option %d, Bin [%d,%d]; New centrality %d; Event plane method %s; ",
856  parList+=onePar ;
857 
858  return new TObjString(parList) ;
859 }
860 
861 
862 //____________________________________________
864 //____________________________________________
866 {
867  if(fMC)
868  return fMC->Stack();
869  else
870  {
871  AliDebug(1,"Stack is not available");
872  return 0x0 ;
873  }
874 }
875 
876 //______________________________________________
878 //______________________________________________
880 {
881  if(fMC)
882  {
883  return fMC->Header();
884  }
885  else
886  {
887  AliInfo("Header is not available");
888  return 0x0 ;
889  }
890 }
891 
892 //______________________________________________________________
894 //______________________________________________________________
895 AliGenEventHeader* AliCaloTrackReader::GetGenEventHeader() const
896 {
897  if (ReadStack() && fMC)
898  {
899  AliGenEventHeader * eventHeader = fMC->GenEventHeader();
900 
901  if(fMCGenerEventHeaderToAccept=="") return eventHeader ;
902 
903  AliGenCocktailEventHeader *cocktail = dynamic_cast<AliGenCocktailEventHeader *>(eventHeader);
904 
905  if(!cocktail) return 0x0 ;
906 
907  TList *genHeaders = cocktail->GetHeaders();
908 
909  Int_t nGenerators = genHeaders->GetEntries();
910 
911  for(Int_t igen = 0; igen < nGenerators; igen++)
912  {
913  AliGenEventHeader * eventHeader2 = (AliGenEventHeader*)genHeaders->At(igen) ;
914  TString name = eventHeader2->GetName();
915  //printf("ESD Event header %d %s\n",igen,name.Data());
916 
917  if(name.Contains(fMCGenerEventHeaderToAccept,TString::kIgnoreCase))
918  return eventHeader2 ;
919  }
920 
921  return 0x0;
922 
923  }
924  else if(ReadAODMCParticles() && GetAODMCHeader())
925  {
926  Int_t nGenerators = GetAODMCHeader()->GetNCocktailHeaders();
927 
928  if( nGenerators <= 0) return 0x0;
929 
930  if(fMCGenerEventHeaderToAccept=="") return GetAODMCHeader()->GetCocktailHeader(0);
931 
932  for(Int_t igen = 0; igen < nGenerators; igen++)
933  {
934  AliGenEventHeader * eventHeader = GetAODMCHeader()->GetCocktailHeader(igen) ;
935  TString name = eventHeader->GetName();
936  //printf("AOD Event header %d %s\n",igen,name.Data());
937 
938  if(name.Contains(fMCGenerEventHeaderToAccept,TString::kIgnoreCase))
939  return eventHeader ;
940  }
941 
942  return 0x0;
943 
944  }
945  else
946  {
947  return 0x0;
948  }
949 }
950 
951 //_________________________________________________________
954 //_________________________________________________________
956 {
957  AliInfo("Input are not AODs");
958 
959  return NULL ;
960 }
961 
962 //________________________________________________________
965 //________________________________________________________
966 AliAODMCHeader* AliCaloTrackReader::GetAODMCHeader() const
967 {
968  AliInfo("Input are not AODs");
969 
970  return NULL ;
971 }
972 
973 //___________________________________________________________
981 //___________________________________________________________
982 Int_t AliCaloTrackReader::GetVertexBC(const AliVVertex * vtx)
983 {
984  Int_t vertexBC=vtx->GetBC();
985 
986  if(!fRecalculateVertexBC) return vertexBC;
987 
988  // Value not available, recalculate it.
989 
990  Double_t bz = fInputEvent->GetMagneticField();
991  Bool_t bc0 = kFALSE;
992  Int_t ntr = GetCTSTracks()->GetEntriesFast();
993  //printf("N Tracks %d\n",ntr);
994 
995  for(Int_t i = 0 ; i < ntr ; i++)
996  {
997  AliVTrack * track = (AliVTrack*) (GetCTSTracks()->At(i));
998 
999  //Check if has TOF info, if not skip
1000  ULong_t status = track->GetStatus();
1001  Bool_t okTOF = (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ;
1002  vertexBC = track->GetTOFBunchCrossing(bz);
1003  Float_t pt = track->Pt();
1004 
1005  if(!okTOF) continue;
1006 
1007  // Get DCA x, y
1008  Double_t dca[2] = {1e6,1e6};
1009  Double_t covar[3] = {1e6,1e6,1e6};
1010  track->PropagateToDCA(vtx,bz,100.,dca,covar);
1011 
1012  if(AcceptDCA(pt,dca[0]))
1013  {
1014  if (vertexBC !=0 && fVertexBC != AliVTrack::kTOFBCNA) return vertexBC;
1015  else if(vertexBC == 0) bc0 = kTRUE;
1016  }
1017  }
1018 
1019  if( bc0 ) vertexBC = 0 ;
1020  else vertexBC = AliVTrack::kTOFBCNA ;
1021 
1022  return vertexBC;
1023 }
1024 
1025 //_____________________________
1032 //_____________________________
1034 {
1035  return track->GetID();
1036 }
1037 
1038 //_____________________________
1041 //_____________________________
1043 {
1045  {
1046  AliInfo("Cannot access stack and mcparticles at the same time, change them");
1047  fReadStack = kFALSE;
1048  fReadAODMCParticles = kFALSE;
1049  }
1050 
1051  // Activate debug level in AliAnaWeights
1052  if( fWeightUtils->GetDebug() >= 0 )
1053  (AliAnalysisManager::GetAnalysisManager())->AddClassDebug(fWeightUtils->ClassName(), fWeightUtils->GetDebug());
1054 }
1055 
1056 //_______________________________________
1058 //_______________________________________
1060 {
1061  fDataType = kESD ;
1062  fCTSPtMin = 0.1 ;
1063  fEMCALPtMin = 0.1 ;
1064  fPHOSPtMin = 0.1 ;
1065  fCTSPtMax = 1000. ;
1066  fEMCALPtMax = 1000. ;
1067  fPHOSPtMax = 1000. ;
1068 
1069  fEMCALBadChMinDist = 0; // open, 2; // standard
1070  fPHOSBadChMinDist = 0; // open, 2; // standard
1071 
1072  fEMCALNCellsCut = 0; // open, 1; // standard
1073  fPHOSNCellsCut = 0; // open, 2; // standard
1074 
1075  //Track DCA cuts
1076  // dca_xy cut = 0.0105+0.0350/TMath::Power(pt,1.1);
1077  fTrackDCACut[0] = 0.0105;
1078  fTrackDCACut[1] = 0.0350;
1079  fTrackDCACut[2] = 1.1;
1080 
1081  //Do not filter the detectors input by default.
1082  fFillEMCAL = kFALSE;
1083  fFillDCAL = kFALSE;
1084  fFillPHOS = kFALSE;
1085  fFillCTS = kFALSE;
1086  fFillEMCALCells = kFALSE;
1087  fFillPHOSCells = kFALSE;
1088 
1089  fReadStack = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
1090  fReadAODMCParticles = kFALSE; // Check in the constructor of the other readers if it was set or in the configuration file
1091  fDeltaAODFileName = "deltaAODPartCorr.root";
1093  fEventTriggerMask = AliVEvent::kAny;
1094  fMixEventTriggerMask = AliVEvent::kAnyINT;
1095  fEventTriggerAtSE = kTRUE; // Use only events that pass event selection at SE base class
1096 
1097  fAcceptFastCluster = kTRUE;
1098  fEventType = -1;
1099 
1100  //We want tracks fitted in the detectors:
1101  //fTrackStatus=AliESDtrack::kTPCrefit;
1102  //fTrackStatus|=AliESDtrack::kITSrefit;
1103  fTrackStatus = 0;
1104 
1105  fV0ADC[0] = 0; fV0ADC[1] = 0;
1106  fV0Mul[0] = 0; fV0Mul[1] = 0;
1107 
1108  fZvtxCut = 10.;
1109 
1110  fNMixedEvent = 1;
1111 
1112  fPtHardAndJetPtFactor = 7.;
1114 
1115  //Centrality
1116  fUseAliCentrality = kFALSE;
1117  fCentralityClass = "V0M";
1118  fCentralityOpt = 100;
1119  fCentralityBin[0] = fCentralityBin[1]=-1;
1120 
1121  fEventPlaneMethod = "V0";
1122 
1123  // Allocate memory (not sure this is the right place)
1124  fCTSTracks = new TObjArray();
1125  fEMCALClusters = new TObjArray();
1126  fDCALClusters = new TObjArray();
1127  fPHOSClusters = new TObjArray();
1128  //fTriggerAnalysis = new AliTriggerAnalysis;
1129  fAODBranchList = new TList ;
1130  fOutputContainer = new TList ;
1131 
1132  fEnergyHistogramNbins = 200;
1133  fEnergyHistogramLimit[0] = 0 ;
1134  fEnergyHistogramLimit[1] = 100;
1135 
1136  fPileUpParamSPD[0] = 3 ; fPileUpParamSPD[1] = 0.8 ;
1137  fPileUpParamSPD[2] = 3.0 ; fPileUpParamSPD[3] = 2.0 ; fPileUpParamSPD[4] = 5.0;
1138 
1139  // Parametrized time cut (LHC11d)
1142 
1143  // Parametrized time cut (LHC11c)
1144  //fEMCALParamTimeCutMin[0] =-5; fEMCALParamTimeCutMin[1] =-1 ; fEMCALParamTimeCutMin[2] = 1.87; fEMCALParamTimeCutMin[3] = 0.4;
1145  //fEMCALParamTimeCutMax[0] = 3.5; fEMCALParamTimeCutMax[1] = 50; fEMCALParamTimeCutMax[2] = 0.15; fEMCALParamTimeCutMax[3] = 1.6;
1146 
1147  fTimeStampRunMin = -1;
1148  fTimeStampRunMax = 1e12;
1151 
1152  for(Int_t i = 0; i < 19; i++)
1153  {
1154  fEMCalBCEvent [i] = 0;
1155  fEMCalBCEventCut[i] = 0;
1156  fTrackBCEvent [i] = 0;
1157  fTrackBCEventCut[i] = 0;
1158  }
1159 
1160  // Trigger match-rejection
1161  fTriggerPatchTimeWindow[0] = 8;
1162  fTriggerPatchTimeWindow[1] = 9;
1163 
1164  fTriggerClusterBC = -10000 ;
1167  fTriggerClusterIndex = -1;
1168  fTriggerClusterId = -1;
1169 
1170  //Jets
1173  if(!fNonStandardJets) fNonStandardJets = new TClonesArray("AliAODJet",100);
1176  if(!fBackgroundJets) fBackgroundJets = new AliAODJetEventBackground();
1177 
1178  fSmearShowerShapeWidth = 0.005;
1179  fSmearNLMMin = 1;
1180  fSmearNLMMax = 1;
1181 
1182  fWeightUtils = new AliAnaWeights() ;
1183  fEventWeight = 1 ;
1184 
1185  fTrackMultNPtCut = 8;
1186  fTrackMultPtCut[0] = 0.15; fTrackMultPtCut[1] = 0.5; fTrackMultPtCut[2] = 1.0;
1187  fTrackMultPtCut[3] = 2.0 ; fTrackMultPtCut[4] = 4.0; fTrackMultPtCut[5] = 6.0;
1188  fTrackMultPtCut[6] = 8.0 ; fTrackMultPtCut[7] = 10.;
1189  fTrackMultPtCut[8] = 15.0; fTrackMultPtCut[9] = 20.;
1190 }
1191 
1192 //__________________________________________________________________________
1195 //__________________________________________________________________________
1197 {
1198  // Parametrized cut depending on E
1199  if(fUseParamTimeCut)
1200  {
1203  //printf("tof %f, minCut %f, maxCut %f\n",tof,minCut,maxCut);
1204  if( tof < minCut || tof > maxCut ) return kFALSE ;
1205  }
1206 
1207  //In any case, the time should to be larger than the fixed window ...
1208  if( tof < fEMCALTimeCutMin || tof > fEMCALTimeCutMax ) return kFALSE ;
1209 
1210  return kTRUE ;
1211 }
1212 
1213 //________________________________________________
1216 //________________________________________________
1218 {
1219  return fInputEvent->IsPileupFromSPD((Int_t) fPileUpParamSPD[0] , fPileUpParamSPD[1] ,
1220  fPileUpParamSPD[2] , fPileUpParamSPD[3] , fPileUpParamSPD[4] );
1221  //printf("Param : %d, %2.2f, %2.2f, %2.2f, %2.2f\n",(Int_t) fPileUpParamSPD[0], fPileUpParamSPD[1], fPileUpParamSPD[2], fPileUpParamSPD[3], fPileUpParamSPD[4]);
1222 }
1223 
1224 //__________________________________________________
1226 //__________________________________________________
1228 {
1229  if(fNPileUpClusters > fNPileUpClustersCut) return kTRUE ;
1230  else return kFALSE;
1231 }
1232 
1233 //________________________________________________________
1235 //________________________________________________________
1237 {
1238  if( IsPileUpFromSPD() && IsPileUpFromEMCal()) return kTRUE ;
1239  else return kFALSE;
1240 }
1241 
1242 //_______________________________________________________
1244 //_______________________________________________________
1246 {
1247  if( IsPileUpFromSPD() || IsPileUpFromEMCal()) return kTRUE ;
1248  else return kFALSE;
1249 }
1250 
1251 //___________________________________________________________
1253 //___________________________________________________________
1255 {
1256  if( IsPileUpFromSPD() && !IsPileUpFromEMCal()) return kTRUE ;
1257  else return kFALSE;
1258 }
1259 
1260 //___________________________________________________________
1262 //___________________________________________________________
1264 {
1265  if( !IsPileUpFromSPD() && IsPileUpFromEMCal()) return kTRUE ;
1266  else return kFALSE;
1267 }
1268 
1269 //______________________________________________________________
1271 //______________________________________________________________
1273 {
1274  if( !IsPileUpFromSPD() && !IsPileUpFromEMCal()) return kTRUE ;
1275  else return kFALSE;
1276 }
1277 
1278 //___________________________________________________________________________________
1286 //___________________________________________________________________________________
1287 Bool_t AliCaloTrackReader::FillInputEvent(Int_t iEntry, const char * /*curFileName*/)
1288 {
1289  fEventNumber = iEntry;
1290  fTriggerClusterIndex = -1;
1291  fTriggerClusterId = -1;
1292  fIsTriggerMatch = kFALSE;
1293  fTriggerClusterBC = -10000;
1294  fIsExoticEvent = kFALSE;
1295  fIsBadCellEvent = kFALSE;
1296  fIsBadMaxCellEvent = kFALSE;
1297 
1298  fIsTriggerMatchOpenCut[0] = kFALSE ;
1299  fIsTriggerMatchOpenCut[1] = kFALSE ;
1300  fIsTriggerMatchOpenCut[2] = kFALSE ;
1301 
1302  //fCurrentFileName = TString(currentFileName);
1303  if(!fInputEvent)
1304  {
1305  AliInfo("Input event not available, skip event analysis");
1306  return kFALSE;
1307  }
1308 
1309  fhNEventsAfterCut->Fill(0.5);
1310 
1311  //-----------------------------------------------
1312  // Select the event depending on the trigger type
1313  // and other event characteristics
1314  // like the goodness of the EMCal trigger
1315  //-----------------------------------------------
1316 
1317  Bool_t accept = CheckEventTriggers();
1318  if(!accept) return kFALSE;
1319 
1320  AliDebug(1,"Pass Event trigger selection");
1321 
1322 
1323  //------------------------------------------------------
1324  // Event rejection depending on time stamp
1325  //------------------------------------------------------
1326 
1328  {
1329  AliESDEvent* esd = dynamic_cast<AliESDEvent*> (fInputEvent);
1330  if(esd)
1331  {
1332  Int_t timeStamp = esd->GetTimeStamp();
1333  Float_t timeStampFrac = 1.*(timeStamp-fTimeStampRunMin) / (fTimeStampRunMax-fTimeStampRunMin);
1334 
1335  //printf("stamp0 %d, max0 %d, frac %f\n", timeStamp-fTimeStampRunMin,fTimeStampRunMax-fTimeStampRunMin, timeStampFrac);
1336 
1337  if(timeStampFrac < fTimeStampEventFracMin || timeStampFrac > fTimeStampEventFracMax) return kFALSE;
1338  }
1339 
1340  AliDebug(1,"Pass Time Stamp rejection");
1341  fhNEventsAfterCut->Fill(8.5);
1342  }
1343 
1344  //------------------------------------------------------
1345  // Event rejection depending on vertex, pileup, v0and
1346  //------------------------------------------------------
1347 
1348  FillVertexArray();
1349 
1350  //Reject events with Z vertex too large, only for SE analysis, if not, cut on the analysis code
1351  if(!GetMixedEvent() && TMath::Abs(fVertex[0][2]) > fZvtxCut) return kFALSE;
1352 
1353  fhNEventsAfterCut->Fill(9.5);
1354 
1356  {
1357  if( !CheckForPrimaryVertex() ) return kFALSE; // algorithm in ESD/AOD Readers
1358  if( TMath::Abs(fVertex[0][0] ) < 1.e-6 &&
1359  TMath::Abs(fVertex[0][1] ) < 1.e-6 &&
1360  TMath::Abs(fVertex[0][2] ) < 1.e-6 ) return kFALSE;
1361  }
1362 
1363  AliDebug(1,"Pass primary vertex rejection");
1364 
1365  fhNEventsAfterCut->Fill(10.5);
1366 
1367  //printf("Reader : IsPileUp %d, Multi %d\n",IsPileUpFromSPD(),fInputEvent->IsPileupFromSPDInMultBins());
1368 
1370  {
1371  // Do not analyze events with pileup
1372  Bool_t bPileup = IsPileUpFromSPD();
1373  //IsPileupFromSPDInMultBins() // method to try
1374  //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]);
1375  if(bPileup) return kFALSE;
1376 
1377  AliDebug(1,"Pass Pile-Up event rejection");
1378 
1379  fhNEventsAfterCut->Fill(11.5);
1380  }
1381 
1383  {
1384  AliVVZERO* v0 = fInputEvent->GetVZEROData();
1385 
1386  Bool_t bV0AND = ((v0->GetV0ADecision()==1) && (v0->GetV0CDecision()==1));
1387  //bV0AND = fTriggerAnalysis->IsOfflineTriggerFired((AliESDEvent*)fInputEvent, AliTriggerAnalysis::kV0AND);
1388  //printf("V0AND event? %d\n",bV0AND);
1389 
1390  if(!bV0AND)
1391  {
1392  AliDebug(1,"Reject event by V0AND");
1393  return kFALSE;
1394  }
1395 
1396  AliDebug(1,"Pass V0AND event rejection");
1397 
1398  fhNEventsAfterCut->Fill(12.5);
1399  }
1400 
1401  //------------------------------------------------------
1402  // Check if there is a centrality value, PbPb analysis,
1403  // and if a centrality bin selection is requested
1404  //------------------------------------------------------
1405 
1406  //If we need a centrality bin, we select only those events in the corresponding bin.
1407  Int_t cen = -1;
1408  if ( fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
1409  {
1410  cen = GetEventCentrality();
1411 
1412  if(cen > fCentralityBin[1] || cen <= fCentralityBin[0]) return kFALSE; //reject events out of bin.
1413 
1414  AliDebug(1,"Pass centrality rejection");
1415 
1416  fhNEventsAfterCut->Fill(13.5);
1417  }
1418 
1419  //----------------------------------------------------------------
1420  // Reject the event if the event header name is not
1421  // the one requested among the possible generators.
1422  // Needed in case of cocktail MC generation with multiple options.
1423  //----------------------------------------------------------------
1425  return kFALSE;
1426 
1427  fhNEventsAfterCut->Fill(14.5);
1428 
1429  //---------------------------------------------------------------------------
1430  // In case of analysis of events with jets, skip those with jet pt > 5 pt hard
1431  // To be used on for MC data in pT hard bins
1432  //---------------------------------------------------------------------------
1433 
1435  {
1436  if(!ComparePtHardAndJetPt()) return kFALSE ;
1437  AliDebug(1,"Pass Pt Hard - Jet rejection");
1438  fhNEventsAfterCut->Fill(15.5);
1439  }
1440 
1442  {
1443  if(!ComparePtHardAndClusterPt()) return kFALSE ;
1444  AliDebug(1,"Pass Pt Hard - Cluster rejection");
1445  fhNEventsAfterCut->Fill(16.5);
1446  }
1447 
1448  //------------------------------------------------------------------
1449  // Recover the weight assigned to the event, if provided
1450  // right now only for pT-hard bins and centrality depedent weights
1451  //------------------------------------------------------------------
1453  {
1455 
1456  fWeightUtils->SetPythiaEventHeader(((AliGenPythiaEventHeader*)GetGenEventHeader()));
1457 
1459  }
1460 
1461  //-------------------------------------------------------
1462  // Get the main vertex BC, in case not available
1463  // it is calculated in FillCTS checking the BC of tracks
1464  //------------------------------------------------------
1465  fVertexBC = fInputEvent->GetPrimaryVertex()->GetBC();
1466 
1467  //-----------------------------------------------
1468  // Fill the arrays with cluster/tracks/cells data
1469  //-----------------------------------------------
1470 
1471  if(fFillCTS)
1472  {
1473  FillInputCTS();
1474  //Accept events with at least one track
1475  if(fTrackMult[0] == 0 && fDoRejectNoTrackEvents) return kFALSE ;
1476 
1477  fhNEventsAfterCut->Fill(17.5);
1478 
1479  AliDebug(1,"Pass rejection of null track events");
1480  }
1481 
1483  {
1484  if(fVertexBC != 0 && fVertexBC != AliVTrack::kTOFBCNA) return kFALSE ;
1485 
1486  AliDebug(1,"Pass rejection of events with vertex at BC!=0");
1487 
1488  fhNEventsAfterCut->Fill(18.5);
1489  }
1490 
1491  if(fFillEMCALCells)
1493 
1494  if(fFillPHOSCells)
1496 
1497  if(fFillEMCAL || fFillDCAL)
1498  FillInputEMCAL();
1499 
1500  if(fFillPHOS)
1501  FillInputPHOS();
1502 
1503  FillInputVZERO();
1504 
1505  //one specified jet branch
1510 
1511  AliDebug(1,"Event accepted for analysis");
1512 
1513  return kTRUE ;
1514 }
1515 
1516 //__________________________________________________
1519 //__________________________________________________
1521 {
1522  if(fUseAliCentrality)
1523  {
1524  if ( !GetCentrality() ) return -1;
1525 
1526  AliDebug(1,Form("Cent. Percentile: V0M %2.2f, CL0 %2.2f, CL1 %2.2f; selected class %s",
1527  GetCentrality()->GetCentralityPercentile("V0M"),
1528  GetCentrality()->GetCentralityPercentile("CL0"),
1529  GetCentrality()->GetCentralityPercentile("CL1"),
1530  fCentralityClass.Data()));
1531 
1532  if (fCentralityOpt == 100) return (Int_t) GetCentrality()->GetCentralityPercentile(fCentralityClass); // 100 bins max
1533  else if(fCentralityOpt == 10) return GetCentrality()->GetCentralityClass10(fCentralityClass);// 10 bins max
1534  else if(fCentralityOpt == 20) return GetCentrality()->GetCentralityClass5(fCentralityClass); // 20 bins max
1535  else
1536  {
1537  AliInfo(Form("Unknown centrality option %d, use 10, 20 or 100\n",fCentralityOpt));
1538  return -1;
1539  }
1540  }
1541  else
1542  {
1543  if ( !GetMultSelCen() ) return -1;
1544 
1545  AliDebug(1,Form("Mult. Percentile: V0M %2.2f, CL0 %2.2f, CL1 %2.2f; selected class %s",
1546  GetMultSelCen()->GetMultiplicityPercentile("V0M",1),
1547  GetMultSelCen()->GetMultiplicityPercentile("CL0",1),
1548  GetMultSelCen()->GetMultiplicityPercentile("CL1",1),
1549  fCentralityClass.Data()));
1550 
1551  return GetMultSelCen()->GetMultiplicityPercentile(fCentralityClass, kTRUE); // returns centrality only for events used in calibration
1552 
1553  // equivalent to
1554  //GetMultSelCen()->GetMultiplicityPercentile("V0M", kFALSE); // returns centrality for any event
1555  //Int_t qual = GetMultSelCen()->GetEvSelCode(); if (qual ! = 0) cent = qual;
1556  }
1557 }
1558 
1559 //_____________________________________________________
1562 //_____________________________________________________
1564 {
1565  if( !GetEventPlane() ) return -1000;
1566 
1567  Float_t ep = GetEventPlane()->GetEventplane(GetEventPlaneMethod(), GetInputEvent());
1568 
1569  if(GetEventPlaneMethod()=="Q" && (ep < 0 || ep > TMath::Pi()))
1570  {
1571  AliDebug(1,Form("Bad EP for <Q> method : %f\n",ep));
1572  return -1000;
1573  }
1574  else if(GetEventPlaneMethod().Contains("V0") )
1575  {
1576  if((ep > TMath::Pi()/2 || ep < -TMath::Pi()/2))
1577  {
1578  AliDebug(1,Form("Bad EP for <%s> method : %f\n",GetEventPlaneMethod().Data(), ep));
1579  return -1000;
1580  }
1581 
1582  ep+=TMath::Pi()/2; // put same range as for <Q> method
1583  }
1584 
1585  AliDebug(3,Form("Event plane angle %f",ep));
1586 
1587 // if(fDebug > 0 )
1588 // {
1589 // if (ep > TMath::Pi()) printf("AliCaloTrackReader::GetEventPlaneAngle() - Too large angle = %f\n",ep);
1590 // else if(ep < 0 ) printf("AliCaloTrackReader::GetEventPlaneAngle() - Negative angle = %f\n" ,ep);
1591 // }
1592 
1593  return ep;
1594 }
1595 
1596 //__________________________________________________________
1598 //__________________________________________________________
1600 {
1601  vertex[0] = fVertex[0][0];
1602  vertex[1] = fVertex[0][1];
1603  vertex[2] = fVertex[0][2];
1604 }
1605 
1606 //__________________________________________________________________________
1608 //__________________________________________________________________________
1609 void AliCaloTrackReader::GetVertex(Double_t vertex[3], Int_t evtIndex) const
1610 {
1611  vertex[0] = fVertex[evtIndex][0];
1612  vertex[1] = fVertex[evtIndex][1];
1613  vertex[2] = fVertex[evtIndex][2];
1614 }
1615 
1616 //________________________________________
1619 //________________________________________
1621 {
1622  // Delete previous vertex
1623  if(fVertex)
1624  {
1625  for (Int_t i = 0; i < fNMixedEvent; i++)
1626  {
1627  delete [] fVertex[i] ;
1628  }
1629  delete [] fVertex ;
1630  }
1631 
1632  fVertex = new Double_t*[fNMixedEvent] ;
1633  for (Int_t i = 0; i < fNMixedEvent; i++)
1634  {
1635  fVertex[i] = new Double_t[3] ;
1636  fVertex[i][0] = 0.0 ;
1637  fVertex[i][1] = 0.0 ;
1638  fVertex[i][2] = 0.0 ;
1639  }
1640 
1641  if (!fMixedEvent)
1642  { // Single event analysis
1643  if(fDataType!=kMC)
1644  {
1645 
1646  if(fInputEvent->GetPrimaryVertex())
1647  {
1648  fInputEvent->GetPrimaryVertex()->GetXYZ(fVertex[0]);
1649  }
1650  else
1651  {
1652  AliWarning("NULL primary vertex");
1653  fVertex[0][0]=0.; fVertex[0][1]=0.; fVertex[0][2]=0.;
1654  }//Primary vertex pointer do not exist
1655 
1656  } else
1657  {// MC read event
1658  fVertex[0][0]=0.; fVertex[0][1]=0.; fVertex[0][2]=0.;
1659  }
1660 
1661  AliDebug(1,Form("Single Event Vertex : %f,%f,%f\n",fVertex[0][0],fVertex[0][1],fVertex[0][2]));
1662 
1663  } else
1664  { // MultiEvent analysis
1665  for (Int_t iev = 0; iev < fNMixedEvent; iev++)
1666  {
1667  if (fMixedEvent->GetVertexOfEvent(iev))
1668  fMixedEvent->GetVertexOfEvent(iev)->GetXYZ(fVertex[iev]);
1669  else
1670  AliWarning("No vertex found");
1671 
1672  AliDebug(1,Form("Multi Event %d Vertex : %f,%f,%f",iev,fVertex[iev][0],fVertex[iev][1],fVertex[iev][2]));
1673  }
1674  }
1675 }
1676 
1677 //_____________________________________
1683 //_____________________________________
1685 {
1686  AliDebug(1,"Begin");
1687 
1688  Double_t pTrack[3] = {0,0,0};
1689 
1690  Int_t nTracks = fInputEvent->GetNumberOfTracks() ;
1691  Int_t nstatus = 0;
1692  Double_t bz = GetInputEvent()->GetMagneticField();
1693 
1694  for(Int_t i = 0; i < 19; i++)
1695  {
1696  fTrackBCEvent [i] = 0;
1697  fTrackBCEventCut[i] = 0;
1698  }
1699 
1700  for(Int_t iptCut = 0; iptCut < fTrackMultNPtCut; iptCut++ )
1701  {
1702  fTrackMult [iptCut] = 0;
1703  fTrackSumPt[iptCut] = 0;
1704  }
1705 
1706  Bool_t bc0 = kFALSE;
1707  if(fRecalculateVertexBC) fVertexBC = AliVTrack::kTOFBCNA;
1708 
1709  for (Int_t itrack = 0; itrack < nTracks; itrack++)
1710  {
1711  AliVTrack * track = (AliVTrack*)fInputEvent->GetTrack(itrack) ; // retrieve track from esd
1712 
1713  if ( !AcceptParticleMCLabel( TMath::Abs(track->GetLabel()) ) ) continue ;
1714 
1715  fhCTSTrackCutsPt[0]->Fill(track->Pt());
1716 
1717  //Select tracks under certain conditions, TPCrefit, ITSrefit ... check the set bits
1718  ULong_t status = track->GetStatus();
1719 
1720  if (fTrackStatus && !((status & fTrackStatus) == fTrackStatus))
1721  continue ;
1722 
1723  fhCTSTrackCutsPt[1]->Fill(track->Pt());
1724 
1725  nstatus++;
1726 
1727  //-------------------------
1728  // Select the tracks depending on cuts of AOD or ESD
1729  if(!SelectTrack(track, pTrack)) continue ;
1730 
1731  fhCTSTrackCutsPt[2]->Fill(track->Pt());
1732 
1733  //-------------------------
1734  // TOF cuts
1735  Bool_t okTOF = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
1736  Double_t tof = -1000;
1737  Int_t trackBC = -1000 ;
1738 
1739  if(fAccessTrackTOF)
1740  {
1741  if(okTOF)
1742  {
1743  trackBC = track->GetTOFBunchCrossing(bz);
1744  SetTrackEventBC(trackBC+9);
1745 
1746  tof = track->GetTOFsignal()*1e-3;
1747 
1748  //After selecting tracks with small DCA, pointing to vertex, set vertex BC depeding on tracks BC
1750  {
1751  if (trackBC != 0 && trackBC != AliVTrack::kTOFBCNA) fVertexBC = trackBC;
1752  else if(trackBC == 0) bc0 = kTRUE;
1753  }
1754 
1755  //In any case, the time should to be larger than the fixed window ...
1756  if( fUseTrackTimeCut && (trackBC !=0 || tof < fTrackTimeCutMin || tof > fTrackTimeCutMax) )
1757  {
1758  //printf("Remove track time %f and bc = %d\n",tof,trackBC);
1759  continue ;
1760  }
1761  //else printf("Accept track time %f and bc = %d\n",tof,trackBC);
1762  }
1763  }
1764 
1765  fhCTSTrackCutsPt[3]->Fill(track->Pt());
1766 
1767  //---------------------
1768  // DCA cuts
1769  //
1770  fMomentum.SetPxPyPzE(pTrack[0],pTrack[1],pTrack[2],0);
1771 
1772  if(fUseTrackDCACut)
1773  {
1774  Float_t dcaTPC =-999;
1775  //In case of AODs, TPC tracks cannot be propagated back to primary vertex,
1776  if( fDataType == kAOD ) dcaTPC = ((AliAODTrack*) track)->DCA();
1777 
1778  //normal way to get the dca, cut on dca_xy
1779  if(dcaTPC==-999)
1780  {
1781  Double_t dca[2] = {1e6,1e6};
1782  Double_t covar[3] = {1e6,1e6,1e6};
1783  Bool_t okDCA = track->PropagateToDCA(fInputEvent->GetPrimaryVertex(),bz,100.,dca,covar);
1784  if( okDCA) okDCA = AcceptDCA(fMomentum.Pt(),dca[0]);
1785  if(!okDCA)
1786  {
1787  //printf("AliCaloTrackReader::FillInputCTS() - Reject track pt %2.2f, dca_xy %2.4f\n",fMomentum.Pt(),dca[0]);
1788  continue ;
1789  }
1790  }
1791  }// DCA cuts
1792 
1793  fhCTSTrackCutsPt[4]->Fill(track->Pt());
1794 
1795  //-------------------------
1796  // Kinematic/acceptance cuts
1797  //
1798  // Count the tracks in eta < 0.9 and different pT cuts
1799  Float_t ptTrack = fMomentum.Pt();
1800  if(TMath::Abs(track->Eta())< fTrackMultEtaCut)
1801  {
1802  for(Int_t iptCut = 0; iptCut < fTrackMultNPtCut; iptCut++ )
1803  {
1804  if(ptTrack > fTrackMultPtCut[iptCut])
1805  {
1806  fTrackMult [iptCut]++;
1807  fTrackSumPt[iptCut]+=ptTrack;
1808  }
1809  }
1810  }
1811 
1812  if(fCTSPtMin > ptTrack || fCTSPtMax < ptTrack) continue ;
1813 
1814  // Check effect of cuts on track BC
1815  if(fAccessTrackTOF && okTOF) SetTrackEventBCcut(trackBC+9);
1816 
1817  if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS)) continue;
1818 
1819  fhCTSTrackCutsPt[5]->Fill(track->Pt());
1820 
1821  // ------------------------------
1822  // Add selected tracks to array
1823  AliDebug(2,Form("Selected tracks pt %3.2f, phi %3.2f deg, eta %3.2f",
1824  fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
1825 
1826  fCTSTracks->Add(track);
1827 
1828  // TODO, check if remove
1829  if (fMixedEvent) track->SetID(itrack);
1830 
1831  }// track loop
1832 
1833  if( fRecalculateVertexBC && (fVertexBC == 0 || fVertexBC == AliVTrack::kTOFBCNA))
1834  {
1835  if( bc0 ) fVertexBC = 0 ;
1836  else fVertexBC = AliVTrack::kTOFBCNA ;
1837  }
1838 
1839  AliDebug(1,Form("AOD entries %d, input tracks %d, pass status %d, multipliticy %d", fCTSTracks->GetEntriesFast(), nTracks, nstatus, fTrackMult[0]));//fCTSTracksNormalInputEntries);
1840 }
1841 
1842 //_______________________________________________________________________________
1860 //_______________________________________________________________________________
1862 {
1863  // Accept clusters with the proper label
1864  if ( clus->GetLabel() >= 0 ) // -1 corresponds to noisy MC
1865  {
1866  if ( !AcceptParticleMCLabel(clus->GetLabel()) ) return ;
1867  }
1868 
1869  // TODO, not sure if needed anymore
1870  Int_t vindex = 0 ;
1871  if (fMixedEvent)
1872  vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
1873 
1874  clus->GetMomentum(fMomentum, fVertex[vindex]);
1875 
1876  // No correction/cut applied yet
1877  fhEMCALClusterCutsE[0]->Fill(clus->E());
1878 
1879  //if( (fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10 )
1880  AliDebug(2,Form("Input cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
1881  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
1882 
1883  //---------------------------
1884  // Embedding case
1885  // TO BE REVISED
1887  {
1888  if(clus->GetNLabels()==0 || clus->GetLabel() < 0) return;
1889  //else printf("Embedded cluster, %d, n label %d label %d \n",iclus,clus->GetNLabels(),clus->GetLabel());
1890  }
1891 
1892  //--------------------------------------
1893  // Apply some corrections in the cluster
1894  //
1896  {
1897  //Recalibrate the cluster energy
1899  {
1901 
1902  clus->SetE(energy);
1903  //printf("Recalibrated Energy %f\n",clus->E());
1904 
1907 
1908  } // recalculate E
1909 
1910  //Recalculate distance to bad channels, if new list of bad channels provided
1912 
1913  //Recalculate cluster position
1914  if(GetCaloUtils()->IsRecalculationOfClusterPositionOn())
1915  {
1917  //clus->GetPosition(pos);
1918  //printf("After Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1919  }
1920 
1921  // Recalculate TOF
1922  if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn())
1923  {
1924  Double_t tof = clus->GetTOF();
1925  Float_t frac =-1;
1926  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
1927 
1928  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
1929 
1930  //additional L1 phase shift
1931  if(GetCaloUtils()->GetEMCALRecoUtils()->IsL1PhaseInTimeRecalibrationOn())
1932  {
1933  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absIdMax),
1934  fInputEvent->GetBunchCrossNumber(), tof);
1935  }
1936 
1937  clus->SetTOF(tof);
1938 
1939  }// Time recalibration
1940  }
1941 
1942  // Check effect of corrections
1943  fhEMCALClusterCutsE[1]->Fill(clus->E());
1944 
1945  //-----------------------------------------------------------------
1946  // Reject clusters with bad channels, close to borders and exotic
1947  //
1948  Bool_t goodCluster = GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,
1949  GetCaloUtils()->GetEMCALGeometry(),
1950  GetEMCALCells(),fInputEvent->GetBunchCrossNumber());
1951 
1952  if(!goodCluster)
1953  {
1954  //if( (fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10 )
1955  AliDebug(2,Form("Bad cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
1956  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
1957 
1958  return;
1959  }
1960 
1961  // Check effect of bad cluster removal
1962  fhEMCALClusterCutsE[2]->Fill(clus->E());
1963 
1964  //Float_t pos[3];
1965  //clus->GetPosition(pos);
1966  //printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
1967 
1968  //--------------------------------------
1969  // Correct non linearity or smear energy
1970  //
1972  {
1974 
1975  //if( (fDebug > 5 && fMomentum.E() > 0.1) || fDebug > 10 )
1976  AliDebug(5,Form("Correct Non Lin: Old E %3.2f, New E %3.2f",
1977  fMomentum.E(),clus->E()));
1978 
1979  // In case of MC analysis, to match resolution/calibration in real data
1980  // Not needed anymore, just leave for MC studies on systematics
1981  if( GetCaloUtils()->GetEMCALRecoUtils()->IsClusterEnergySmeared() )
1982  {
1983  Float_t rdmEnergy = GetCaloUtils()->GetEMCALRecoUtils()->SmearClusterEnergy(clus);
1984 
1985  //if( (fDebug > 5 && fMomentum.E() > 0.1) || fDebug > 10 )
1986  AliDebug(5,Form("Smear energy: Old E %3.2f, New E %3.2f",clus->E(),rdmEnergy));
1987 
1988  clus->SetE(rdmEnergy);
1989  }
1990  }
1991 
1992  clus->GetMomentum(fMomentum, fVertex[vindex]);
1993 
1994  // Check effect linearity correction, energy smearing
1995  fhEMCALClusterCutsE[3]->Fill(clus->E());
1996 
1997  // Check the event BC depending on EMCal clustr before final cuts
1998  Double_t tof = clus->GetTOF()*1e9;
1999 
2000  Int_t bc = TMath::Nint(tof/50) + 9;
2001  //printf("tof %2.2f, bc+5=%d\n",tof,bc);
2002 
2003  SetEMCalEventBC(bc);
2004 
2005  //--------------------------------------
2006  // Apply some kinematical/acceptance cuts
2007  //
2008  if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E()) return ;
2009 
2010  // Select cluster fiducial region
2011  //
2012  Bool_t bEMCAL = kFALSE;
2013  Bool_t bDCAL = kFALSE;
2014  if(fCheckFidCut)
2015  {
2016  if(fFillEMCAL && fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) bEMCAL = kTRUE ;
2017  if(fFillDCAL && fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kDCAL )) bDCAL = kTRUE ;
2018  }
2019  else
2020  {
2021  bEMCAL = kTRUE;
2022  }
2023 
2024  //---------------------------------------------------------------------
2025  // Mask all cells in collumns facing ALICE thick material if requested
2026  //
2028  {
2029  Int_t absId = -1;
2030  Int_t iSupMod = -1;
2031  Int_t iphi = -1;
2032  Int_t ieta = -1;
2033  Bool_t shared = kFALSE;
2034  GetCaloUtils()->GetEMCALRecoUtils()->GetMaxEnergyCell(GetCaloUtils()->GetEMCALGeometry(), GetEMCALCells(),clus,absId,iSupMod,ieta,iphi,shared);
2035 
2036  AliDebug(2,Form("Masked collumn: cluster E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2037  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2038 
2039 
2040  if(GetCaloUtils()->MaskFrameCluster(iSupMod, ieta)) return;
2041  }
2042 
2043  // Check effect of energy and fiducial cuts
2044  fhEMCALClusterCutsE[4]->Fill(clus->E());
2045 
2046 
2047  //------------------------------------------
2048  // Apply time cut, count EMCal BC before cut
2049  //
2050  SetEMCalEventBCcut(bc);
2051 
2052  if(!IsInTimeWindow(tof,clus->E()))
2053  {
2054  fNPileUpClusters++ ;
2055  if(fUseEMCALTimeCut)
2056  {
2057  AliDebug(2,Form("Out of time window E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f, time %e",
2058  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta(),tof));
2059 
2060  return ;
2061  }
2062  }
2063  else
2065 
2066  // Check effect of time cut
2067  fhEMCALClusterCutsE[5]->Fill(clus->E());
2068 
2069 
2070  //----------------------------------------------------
2071  // Apply N cells cut
2072  //
2073  if(clus->GetNCells() <= fEMCALNCellsCut && fDataType != AliCaloTrackReader::kMC) return ;
2074 
2075  // Check effect of n cells cut
2076  fhEMCALClusterCutsE[6]->Fill(clus->E());
2077 
2078  //----------------------------------------------------
2079  // Apply distance to bad channel cut
2080  //
2081  Double_t distBad = clus->GetDistanceToBadChannel() ; //Distance to bad channel
2082 
2083  if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2084 
2085  if(distBad < fEMCALBadChMinDist) return ;
2086 
2087  // Check effect distance to bad channel cut
2088  fhEMCALClusterCutsE[7]->Fill(clus->E());
2089 
2090  //----------------------------------------------------
2091  // Smear the SS to try to match data and simulations,
2092  // do it only for simulations.
2093  //
2094  Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(clus, GetEMCALCells());
2095  // Int_t nMaxima = clus->GetNExMax(); // For Run2
2096  if( fSmearShowerShape && clus->GetNCells() > 2 &&
2097  nMaxima >= fSmearNLMMin && nMaxima <= fSmearNLMMax )
2098  {
2099  AliDebug(2,Form("Smear shower shape - Original: %2.4f\n", clus->GetM02()));
2101  {
2102  clus->SetM02( clus->GetM02() + fRandom.Landau(0, fSmearShowerShapeWidth) );
2103  }
2105  {
2106  if(iclus%3 == 0 && clus->GetM02() > 0.1) clus->SetM02( clus->GetM02() + fRandom.Landau(0.05, fSmearShowerShapeWidth) ); //fSmearShowerShapeWidth = 0.035
2107  }
2108  else if (fSmearingFunction == kNoSmearing)
2109  {
2110  clus->SetM02( clus->GetM02() );
2111  }
2112  //clus->SetM02( fRandom.Landau(clus->GetM02(), fSmearShowerShapeWidth) );
2113  AliDebug(2,Form("Width %2.4f Smeared : %2.4f\n", fSmearShowerShapeWidth,clus->GetM02()));
2114  }
2115 
2116  //--------------------------------------------------------
2117  // Fill the corresponding array with the selected clusters
2118  // Usually just filling EMCal array with upper or lower clusters is enough,
2119  // but maybe we want to do EMCal-DCal correlations.
2120 
2121  //if((fDebug > 2 && fMomentum.E() > 0.1) || fDebug > 10)
2122  AliDebug(2,Form("Selected clusters (EMCAL%d, DCAL%d), E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2123  bEMCAL,bDCAL,fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2124 
2125 
2126  if (bEMCAL) fEMCALClusters->Add(clus);
2127  else if(bDCAL ) fDCALClusters ->Add(clus);
2128 
2129  // TODO, not sure if needed anymore
2130  if (fMixedEvent)
2131  clus->SetID(iclus) ;
2132 }
2133 
2134 //_______________________________________
2141 //_______________________________________
2143 {
2144  AliDebug(1,"Begin");
2145 
2146  // First recalibrate cells, time or energy
2147  // if(GetCaloUtils()->IsRecalibrationOn())
2148  // GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCells(GetCaloUtils()->GetEMCALGeometry(),
2149  // GetEMCALCells(),
2150  // fInputEvent->GetBunchCrossNumber());
2151 
2152  fNPileUpClusters = 0; // Init counter
2153  fNNonPileUpClusters = 0; // Init counter
2154  for(Int_t i = 0; i < 19; i++)
2155  {
2156  fEMCalBCEvent [i] = 0;
2157  fEMCalBCEventCut[i] = 0;
2158  }
2159 
2160  //Loop to select clusters in fiducial cut and fill container with aodClusters
2161  if(fEMCALClustersListName=="")
2162  {
2163  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2164  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2165  {
2166  AliVCluster * clus = 0;
2167  if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
2168  {
2169  if (clus->IsEMCAL())
2170  {
2171  FillInputEMCALAlgorithm(clus, iclus);
2172  }//EMCAL cluster
2173  }// cluster exists
2174  }// cluster loop
2175 
2176  //Recalculate track matching
2178 
2179  }//Get the clusters from the input event
2180  else
2181  {
2182  TClonesArray * clusterList = 0x0;
2183 
2184  if (fInputEvent->FindListObject(fEMCALClustersListName))
2185  {
2186  clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
2187  }
2188  else if(fOutputEvent)
2189  {
2190  clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
2191  }
2192 
2193  if(!clusterList)
2194  {
2195  AliWarning(Form("Wrong name of list with clusters? <%s>",fEMCALClustersListName.Data()));
2196  return;
2197  }
2198 
2199  Int_t nclusters = clusterList->GetEntriesFast();
2200  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2201  {
2202  AliVCluster * clus = dynamic_cast<AliVCluster*> (clusterList->At(iclus));
2203  //printf("E %f\n",clus->E());
2204  if (clus) FillInputEMCALAlgorithm(clus, iclus);
2205  else AliWarning("Null cluster in list!");
2206  }// cluster loop
2207 
2208  // Recalculate the pile-up time, in case long time clusters removed during clusterization
2209  //printf("Input event INIT : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
2210 
2211  fNPileUpClusters = 0; // Init counter
2212  fNNonPileUpClusters = 0; // Init counter
2213  for(Int_t i = 0; i < 19; i++)
2214  {
2215  fEMCalBCEvent [i] = 0;
2216  fEMCalBCEventCut[i] = 0;
2217  }
2218 
2219  for (Int_t iclus = 0; iclus < fInputEvent->GetNumberOfCaloClusters(); iclus++)
2220  {
2221  AliVCluster * clus = 0;
2222 
2223  if ( (clus = fInputEvent->GetCaloCluster(iclus)) )
2224  {
2225  if (clus->IsEMCAL())
2226  {
2227 
2228  Float_t frac =-1;
2229  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fEMCALCells, clus,frac);
2230  Double_t tof = clus->GetTOF();
2231  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2232  //additional L1 phase shift
2233  if(GetCaloUtils()->GetEMCALRecoUtils()->IsL1PhaseInTimeRecalibrationOn())
2234  {
2235  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absIdMax), fInputEvent->GetBunchCrossNumber(), tof);
2236  }
2237 
2238  tof*=1e9;
2239 
2240  //printf("Input event cluster : AbsIdMax %d, E %2.2f, time %2.2f \n", absIdMax,clus->E(),tof);
2241 
2242  //Reject clusters with bad channels, close to borders and exotic;
2243  if(!GetCaloUtils()->GetEMCALRecoUtils()->IsGoodCluster(clus,GetCaloUtils()->GetEMCALGeometry(),GetEMCALCells(),fInputEvent->GetBunchCrossNumber())) continue;
2244 
2245  Int_t bc = TMath::Nint(tof/50) + 9;
2246  SetEMCalEventBC(bc);
2247 
2248  if(fEMCALPtMin > clus->E() || fEMCALPtMax < clus->E()) continue ;
2249 
2250  clus->GetMomentum(fMomentum, fVertex[0]);
2251 
2252  if(fCheckFidCut && !fFiducialCut->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) return ;
2253 
2254  SetEMCalEventBCcut(bc);
2255 
2256  if(!IsInTimeWindow(tof,clus->E()))
2257  fNPileUpClusters++ ;
2258  else
2260 
2261  }
2262  }
2263  }
2264 
2265  //printf("Input event : Pile-up clusters %d, NO pile-up %d\n",fNPileUpClusters,fNNonPileUpClusters);
2266 
2267  // Recalculate track matching, not necessary if already done in the reclusterization task.
2268  // in case it was not done ...
2270 
2271  }
2272 
2273  AliDebug(1,Form("AOD entries %d, n pile-up clusters %d, n non pile-up %d", fEMCALClusters->GetEntriesFast(),fNPileUpClusters,fNNonPileUpClusters));
2274 }
2275 
2276 //_______________________________________
2278 //_______________________________________
2280 {
2281  AliDebug(1,"Begin");
2282 
2283  // Loop to select clusters in fiducial cut and fill container with aodClusters
2284  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2285  TString genName;
2286  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2287  {
2288  AliVCluster * clus = fInputEvent->GetCaloCluster(iclus) ;
2289  if ( !clus ) continue ;
2290 
2291  if ( !clus->IsPHOS() ) continue ;
2292 
2293  if(clus->GetLabel() >=0 ) // -1 corresponds to noisy MC
2294  {
2295  if ( !AcceptParticleMCLabel(clus->GetLabel()) ) continue ;
2296  }
2297 
2298  fhPHOSClusterCutsE[0]->Fill(clus->E());
2299 
2300  // Skip CPV input
2301  if( clus->GetType() == AliVCluster::kPHOSCharged ) continue ;
2302 
2303  fhPHOSClusterCutsE[1]->Fill(clus->E());
2304 
2305  //---------------------------------------------
2306  // Old feature, try to rely on PHOS tender
2307  //
2309  {
2310  // Recalibrate the cluster energy
2312  {
2313  Float_t energy = GetCaloUtils()->RecalibrateClusterEnergy(clus, (AliAODCaloCells*)GetPHOSCells());
2314  clus->SetE(energy);
2315  }
2316  }
2317 
2318  //----------------------------------------------------------------------------------
2319  // Check if the cluster contains any bad channel and if close to calorimeter borders
2320  //
2321  // Old feature, try to rely on PHOS tender
2322  if( GetCaloUtils()->ClusterContainsBadChannel(kPHOS,clus->GetCellsAbsId(), clus->GetNCells()))
2323  continue;
2324 
2325  if(!GetCaloUtils()->CheckCellFiducialRegion(clus, fInputEvent->GetPHOSCells()))
2326  continue;
2327 
2328  // TODO, add exotic cut???
2329 
2330  fhPHOSClusterCutsE[2]->Fill(clus->E());
2331 
2332  // TODO Dead code? remove?
2333  Int_t vindex = 0 ;
2334  if (fMixedEvent)
2335  vindex = fMixedEvent->EventIndexForCaloCluster(iclus);
2336 
2337  clus->GetMomentum(fMomentum, fVertex[vindex]);
2338 
2339  //----------------------------------------------------------------------------------
2340  // Remove clusters close to borders
2341  //
2343  continue ;
2344 
2345  fhPHOSClusterCutsE[3]->Fill(clus->E());
2346 
2347  //----------------------------------------------------------------------------------
2348  // Remove clusters with too low energy
2349  //
2350  if (fPHOSPtMin > fMomentum.E() || fPHOSPtMax < fMomentum.E() )
2351  continue ;
2352 
2353  fhPHOSClusterCutsE[4]->Fill(clus->E());
2354 
2355  //----------------------------------------------------
2356  // Apply N cells cut
2357  //
2358  if(clus->GetNCells() <= fPHOSNCellsCut && fDataType != AliCaloTrackReader::kMC) return ;
2359 
2360  // Check effect of n cells cut
2361  fhPHOSClusterCutsE[5]->Fill(clus->E());
2362 
2363  //----------------------------------------------------
2364  // Apply distance to bad channel cut
2365  //
2366  Double_t distBad = clus->GetDistanceToBadChannel() ; //Distance to bad channel
2367 
2368  if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2369 
2370  if(distBad < fPHOSBadChMinDist) return ;
2371 
2372  // Check effect distance to bad channel cut
2373  fhPHOSClusterCutsE[6]->Fill(clus->E());
2374 
2375  // TODO, add time cut
2376 
2377  //----------------------------------------------------------------------------------
2378  // Add selected clusters to array
2379  //
2380  //if(fDebug > 2 && fMomentum.E() > 0.1)
2381  AliDebug(2,Form("Selected clusters E %3.2f, pt %3.2f, phi %3.2f deg, eta %3.2f",
2382  fMomentum.E(),fMomentum.Pt(),RadToDeg(GetPhi(fMomentum.Phi())),fMomentum.Eta()));
2383 
2384  fPHOSClusters->Add(clus);
2385 
2386  // TODO Dead code? remove?
2387  if (fMixedEvent)
2388  clus->SetID(iclus) ;
2389 
2390  } // esd/aod cluster loop
2391 
2392  AliDebug(1,Form("AOD entries %d",fPHOSClusters->GetEntriesFast())) ;
2393 }
2394 
2395 //____________________________________________
2397 //____________________________________________
2399 {
2400  fEMCALCells = fInputEvent->GetEMCALCells();
2401 }
2402 
2403 //___________________________________________
2405 //___________________________________________
2407 {
2408  fPHOSCells = fInputEvent->GetPHOSCells();
2409 }
2410 
2411 //_______________________________________
2414 //_______________________________________
2416 {
2417  AliVVZERO* v0 = fInputEvent->GetVZEROData();
2418  //printf("Init V0: ADC (%d,%d), Multiplicity (%d,%d) \n",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]);
2419 
2420  if (v0)
2421  {
2422  AliESDVZERO* esdV0 = dynamic_cast<AliESDVZERO*> (v0);
2423  for (Int_t i = 0; i < 32; i++)
2424  {
2425  if(esdV0)
2426  {//Only available in ESDs
2427  fV0ADC[0] += (Int_t)esdV0->GetAdcV0C(i);
2428  fV0ADC[1] += (Int_t)esdV0->GetAdcV0A(i);
2429  }
2430 
2431  fV0Mul[0] += (Int_t)v0->GetMultiplicityV0C(i);
2432  fV0Mul[1] += (Int_t)v0->GetMultiplicityV0A(i);
2433  }
2434 
2435  AliDebug(1,Form("ADC (%d,%d), Multiplicity (%d,%d)",fV0ADC[0],fV0ADC[1],fV0Mul[0],fV0Mul[1]));
2436  }
2437  else
2438  {
2439  AliDebug(1,"Cannot retrieve V0 ESD! Run w/ null V0 charges");
2440  }
2441 }
2442 
2443 //_________________________________________________
2447 //_________________________________________________
2449 {
2450  AliDebug(2,"Begin");
2451 
2452  //
2453  //check if branch name is given
2454  if(!fInputNonStandardJetBranchName.Length())
2455  {
2456  fInputEvent->Print();
2457  AliFatal("No non-standard jet branch name specified. Specify among existing ones.");
2458  }
2459 
2460  fNonStandardJets = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(fInputNonStandardJetBranchName.Data()));
2461 
2462  if(!fNonStandardJets)
2463  {
2464  //check if jet branch exist; exit if not
2465  fInputEvent->Print();
2466 
2467  AliFatal(Form("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fInputNonStandardJetBranchName.Data()));
2468  }
2469  else
2470  {
2471  AliDebug(1,Form("AOD input jets %d", fNonStandardJets->GetEntriesFast()));
2472  }
2473 }
2474 
2475 //_________________________________________________
2479 //_________________________________________________
2481 {
2482  AliDebug(1,"Begin");
2483  //
2484  //check if branch name is given
2485  if(!fInputBackgroundJetBranchName.Length())
2486  {
2487  fInputEvent->Print();
2488 
2489  AliFatal("No background jet branch name specified. Specify among existing ones.");
2490  }
2491 
2492  fBackgroundJets = (AliAODJetEventBackground*)(fInputEvent->FindListObject(fInputBackgroundJetBranchName.Data()));
2493 
2494  if(!fBackgroundJets)
2495  {
2496  //check if jet branch exist; exit if not
2497  fInputEvent->Print();
2498 
2499  AliFatal(Form("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fInputBackgroundJetBranchName.Data()));
2500  }
2501  else
2502  {
2503  AliDebug(1,"FillInputBackgroundJets");
2504  fBackgroundJets->Print("");
2505  }
2506 }
2507 
2508 //____________________________________________________________________
2514 //____________________________________________________________________
2516 {
2517  // init some variables
2518  Int_t trigtimes[30], globCol, globRow,ntimes, i;
2519  Int_t absId = -1; //[100];
2520  Int_t nPatch = 0;
2521 
2522  TArrayI patches(0);
2523 
2524  // get object pointer
2525  AliVCaloTrigger *caloTrigger = GetInputEvent()->GetCaloTrigger( "EMCAL" );
2526 
2527  if(!caloTrigger)
2528  {
2529  AliError("Trigger patches input (AliVCaloTrigger) not available in data!");
2530  return patches;
2531  }
2532 
2533  //printf("CaloTrigger Entries %d\n",caloTrigger->GetEntries() );
2534 
2535  // class is not empty
2536  if( caloTrigger->GetEntries() > 0 )
2537  {
2538  // must reset before usage, or the class will fail
2539  caloTrigger->Reset();
2540 
2541  // go throuth the trigger channels
2542  while( caloTrigger->Next() )
2543  {
2544  // get position in global 2x2 tower coordinates
2545  caloTrigger->GetPosition( globCol, globRow );
2546 
2547  //L0
2548  if(IsEventEMCALL0())
2549  {
2550  // get dimension of time arrays
2551  caloTrigger->GetNL0Times( ntimes );
2552 
2553  // no L0s in this channel
2554  // presence of the channel in the iterator still does not guarantee that L0 was produced!!
2555  if( ntimes < 1 )
2556  continue;
2557 
2558  // get timing array
2559  caloTrigger->GetL0Times( trigtimes );
2560  //printf("Get L0 patch : n times %d - trigger time window %d - %d\n",ntimes, tmin,tmax);
2561 
2562  // go through the array
2563  for( i = 0; i < ntimes; i++ )
2564  {
2565  // check if in cut - 8,9 shall be accepted in 2011
2566  if( trigtimes[i] >= tmin && trigtimes[i] <= tmax )
2567  {
2568  //printf("Accepted trigger time %d \n",trigtimes[i]);
2569  //if(nTrig > 99) continue;
2570  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(globCol,globRow, absId);
2571  //printf("pass the time cut globCol %d, globRow %d absId %d\n",globCol,globRow, absId);
2572  patches.Set(nPatch+1);
2573  patches.AddAt(absId,nPatch++);
2574  }
2575  } // trigger time array
2576  }//L0
2577  else if(IsEventEMCALL1()) // L1
2578  {
2579  Int_t bit = 0;
2580  caloTrigger->GetTriggerBits(bit);
2581 
2582  Int_t sum = 0;
2583  caloTrigger->GetL1TimeSum(sum);
2584  //fBitEGA-=2;
2585  Bool_t isEGA1 = ((bit >> fBitEGA ) & 0x1) && IsEventEMCALL1Gamma1() ;
2586  Bool_t isEGA2 = ((bit >> (fBitEGA+1)) & 0x1) && IsEventEMCALL1Gamma2() ;
2587  Bool_t isEJE1 = ((bit >> fBitEJE ) & 0x1) && IsEventEMCALL1Jet1 () ;
2588  Bool_t isEJE2 = ((bit >> (fBitEJE+1)) & 0x1) && IsEventEMCALL1Jet2 () ;
2589 
2590  //if((bit>> fBitEGA )&0x1) printf("Trig Bit %d - bit %d - EG1 %d - EG2 %d\n",fBitEGA ,bit,IsEventEMCALL1Gamma1(),IsEventEMCALL1Gamma2());
2591  //if((bit>>(fBitEGA+1))&0x1) printf("Trig Bit %d - bit %d - EG1 %d - EG2 %d\n",fBitEGA+1,bit,IsEventEMCALL1Gamma1(),IsEventEMCALL1Gamma2());
2592 
2593  if(!isEGA1 && !isEJE1 && !isEGA2 && !isEJE2) continue;
2594 
2595  Int_t patchsize = -1;
2596  if (isEGA1 || isEGA2) patchsize = 2;
2597  else if (isEJE1 || isEJE2) patchsize = 16;
2598 
2599  //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",
2600  // bit,sum,patchsize,isEGA1,isEGA2,isEJE1,isEJE2,fBitEGA,fBitEJE,IsEventEMCALL1Gamma(),IsEventEMCALL1Jet());
2601 
2602 
2603  // add 2x2 (EGA) or 16x16 (EJE) patches
2604  for(Int_t irow=0; irow < patchsize; irow++)
2605  {
2606  for(Int_t icol=0; icol < patchsize; icol++)
2607  {
2608  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(globCol+icol,globRow+irow, absId);
2609  //printf("pass the time cut globCol %d, globRow %d absId %d\n",globCol,globRow, absId);
2610  patches.Set(nPatch+1);
2611  patches.AddAt(absId,nPatch++);
2612  }
2613  }
2614 
2615  } // L1
2616 
2617  } // trigger iterator
2618  } // go through triggers
2619 
2620  if(patches.GetSize()<=0) AliInfo(Form("No patch found! for triggers: %s and selected <%s>",
2622  //else printf(">>>>> N patches %d, test %d,first %d, last %d\n",patches.GetSize(), nPatch, patches.At(0), patches.At(patches.GetSize()-1));
2623 
2624  return patches;
2625 }
2626 
2627 //____________________________________________________________
2631 //____________________________________________________________
2633 {
2634  // Init info from previous event
2635  fTriggerClusterIndex = -1;
2636  fTriggerClusterId = -1;
2637  fTriggerClusterBC = -10000;
2638  fIsExoticEvent = kFALSE;
2639  fIsBadCellEvent = kFALSE;
2640  fIsBadMaxCellEvent = kFALSE;
2641 
2642  fIsTriggerMatch = kFALSE;
2643  fIsTriggerMatchOpenCut[0] = kFALSE;
2644  fIsTriggerMatchOpenCut[1] = kFALSE;
2645  fIsTriggerMatchOpenCut[2] = kFALSE;
2646 
2647  // Do only analysis for triggered events
2648  if(!IsEventEMCALL1() && !IsEventEMCALL0())
2649  {
2650  fTriggerClusterBC = 0;
2651  return;
2652  }
2653 
2654  //printf("***** Try to match trigger to cluster %d **** L0 %d, L1 %d\n",fTriggerPatchClusterMatch,IsEventEMCALL0(),IsEventEMCALL1());
2655 
2656  //Recover the list of clusters
2657  TClonesArray * clusterList = 0;
2658  if (fInputEvent->FindListObject(fEMCALClustersListName))
2659  {
2660  clusterList = dynamic_cast<TClonesArray*> (fInputEvent->FindListObject(fEMCALClustersListName));
2661  }
2662  else if(fOutputEvent)
2663  {
2664  clusterList = dynamic_cast<TClonesArray*> (fOutputEvent->FindListObject(fEMCALClustersListName));
2665  }
2666 
2667  // Get number of clusters and of trigger patches
2668  Int_t nclusters = fInputEvent->GetNumberOfCaloClusters();
2669  if(clusterList)
2670  nclusters = clusterList->GetEntriesFast();
2671 
2672  Int_t nPatch = patches.GetSize();
2673  Float_t exoDiffTime = GetCaloUtils()->GetEMCALRecoUtils()->GetExoticCellDiffTimeCut();
2674 
2675  //Init some variables used in the cluster loop
2676  Float_t tofPatchMax = 100000;
2677  Float_t ePatchMax =-1;
2678 
2679  Float_t tofMax = 100000;
2680  Float_t eMax =-1;
2681 
2682  Int_t clusMax =-1;
2683  Int_t idclusMax =-1;
2684  Bool_t badClMax = kFALSE;
2685  Bool_t badCeMax = kFALSE;
2686  Bool_t exoMax = kFALSE;
2687 // Int_t absIdMaxTrig= -1;
2688  Int_t absIdMaxMax = -1;
2689 
2690  Int_t nOfHighECl = 0 ;
2691 
2692  //
2693  // Check what is the trigger threshold
2694  // set minimu energym of candidate for trigger cluster
2695  //
2697 
2698  Float_t triggerThreshold = fTriggerL1EventThreshold;
2699  if(IsEventEMCALL0()) triggerThreshold = fTriggerL0EventThreshold;
2700  Float_t minE = triggerThreshold / 2.;
2701 
2702  // This method is not really suitable for JET trigger
2703  // but in case, reduce the energy cut since we do not trigger on high energy particle
2704  if(IsEventEMCALL1Jet() || minE < 1) minE = 1;
2705 
2706  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));
2707 
2708  //
2709  // Loop on the clusters, check if there is any that falls into one of the patches
2710  //
2711  for (Int_t iclus = 0; iclus < nclusters; iclus++)
2712  {
2713  AliVCluster * clus = 0;
2714  if(clusterList) clus = (AliVCluster*) clusterList->At(iclus);
2715  else clus = fInputEvent->GetCaloCluster(iclus);
2716 
2717  if ( !clus ) continue ;
2718 
2719  if ( !clus->IsEMCAL() ) continue ;
2720 
2721  //Skip clusters with too low energy to be triggering
2722  if ( clus->E() < minE ) continue ;
2723 
2724  Float_t frac = -1;
2725  Int_t absIdMax = GetCaloUtils()->GetMaxEnergyCell(fInputEvent->GetEMCALCells(), clus,frac);
2726 
2727  Bool_t badCluster = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),
2728  clus->GetCellsAbsId(),clus->GetNCells());
2729  UShort_t cellMax[] = {(UShort_t) absIdMax};
2730  Bool_t badCell = GetCaloUtils()->GetEMCALRecoUtils()->ClusterContainsBadChannel(GetCaloUtils()->GetEMCALGeometry(),cellMax,1);
2731 
2732  // if cell is bad, it can happen that time calibration is not available,
2733  // when calculating if it is exotic, this can make it to be exotic by default
2734  // open it temporarily for this cluster
2735  if(badCell)
2736  GetCaloUtils()->GetEMCALRecoUtils()->SetExoticCellDiffTimeCut(10000000);
2737 
2738  Bool_t exotic = GetCaloUtils()->GetEMCALRecoUtils()->IsExoticCluster(clus, fInputEvent->GetEMCALCells());
2739 
2740  // Set back the time cut on exotics
2741  if(badCell)
2742  GetCaloUtils()->GetEMCALRecoUtils()->SetExoticCellDiffTimeCut(exoDiffTime);
2743 
2744  // Energy threshold for exotic Ecross typically at 4 GeV,
2745  // for lower energy, check that there are more than 1 cell in the cluster
2746  if(!exotic && clus->GetNCells() < 2) exotic = kTRUE;
2747 
2748  Float_t energy = clus->E();
2749  Int_t idclus = clus->GetID();
2750 
2751  Double_t tof = clus->GetTOF();
2752  if(GetCaloUtils()->GetEMCALRecoUtils()->IsTimeRecalibrationOn() && fTriggerClusterTimeRecal){
2753  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTime(absIdMax,fInputEvent->GetBunchCrossNumber(),tof);
2754  //additional L1 phase shift
2755  if(GetCaloUtils()->GetEMCALRecoUtils()->IsL1PhaseInTimeRecalibrationOn()){
2756  GetCaloUtils()->GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absIdMax), fInputEvent->GetBunchCrossNumber(), tof);
2757  }
2758  }
2759  tof *=1.e9;
2760 
2761  //printf("cluster %d, ID %d, E %2.2f, tof %2.2f, AbsId max %d, exotic %d, bad Cluster %d, bad Cell %d\n",
2762  // iclus,idclus, energy,tof,absIdMax, exotic, badCluster,badCell);
2763 
2764  // Find the highest energy cluster, avobe trigger threshold
2765  // in the event in case no match to trigger is found
2766  if( energy > eMax )
2767  {
2768  tofMax = tof;
2769  eMax = energy;
2770  badClMax = badCluster;
2771  badCeMax = badCell;
2772  exoMax = exotic;
2773  clusMax = iclus;
2774  idclusMax = idclus;
2775  absIdMaxMax = absIdMax;
2776  }
2777 
2778  // count the good clusters in the event avobe the trigger threshold
2779  // to check the exotic events
2780  if(!badCluster && !exotic)
2781  nOfHighECl++;
2782 
2783  // Find match to trigger
2784  if(fTriggerPatchClusterMatch && nPatch>0)
2785  {
2786  for(Int_t iabsId =0; iabsId < nPatch; iabsId++)
2787  {
2788  Int_t absIDCell[4];
2789  GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patches.At(iabsId), absIDCell);
2790  //if(tof > 75 ) printf("E %2.2f TOF %2.2f Trigger patch %d, cells : %d, %d, %d, %d\n",
2791  // clus->E(),tof,patches.At(iabsId), absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
2792 
2793  for(Int_t ipatch = 0; ipatch < 4; ipatch++)
2794  {
2795  if(absIdMax == absIDCell[ipatch])
2796  {
2797  //printf("*** Patches : absId %d, E %2.1f, tof %f \n",absIdMax,clus->E(), tof);
2798  if(energy > ePatchMax)
2799  {
2800  tofPatchMax = tof;
2801  ePatchMax = energy;
2802  fIsBadCellEvent = badCluster;
2803  fIsBadMaxCellEvent = badCell;
2804  fIsExoticEvent = exotic;
2805  fTriggerClusterIndex = iclus;
2806  fTriggerClusterId = idclus;
2807  fIsTriggerMatch = kTRUE;
2808 // absIdMaxTrig = absIdMax;
2809  }
2810  }
2811  }// cell patch loop
2812  }// trigger patch loop
2813  } // Do trigger patch matching
2814 
2815  }// Cluster loop
2816 
2817  // If there was no match, assign as trigger
2818  // the highest energy cluster in the event
2819  if(!fIsTriggerMatch)
2820  {
2821  tofPatchMax = tofMax;
2822  ePatchMax = eMax;
2823  fIsBadCellEvent = badClMax;
2824  fIsBadMaxCellEvent = badCeMax;
2825  fIsExoticEvent = exoMax;
2826  fTriggerClusterIndex = clusMax;
2827  fTriggerClusterId = idclusMax;
2828  }
2829 
2830  Double_t tofPatchMaxUS = TMath::Abs(tofPatchMax);
2831 
2832  if (tofPatchMaxUS < 28 ) fTriggerClusterBC = 0 ;
2833  else if(tofPatchMaxUS < 75 ) fTriggerClusterBC = 1 ;
2834  else if(tofPatchMaxUS < 125) fTriggerClusterBC = 2 ;
2835  else if(tofPatchMaxUS < 175) fTriggerClusterBC = 3 ;
2836  else if(tofPatchMaxUS < 225) fTriggerClusterBC = 4 ;
2837  else if(tofPatchMaxUS < 275) fTriggerClusterBC = 5 ;
2838  else
2839  {
2840  //printf("AliCaloTrackReader::MatchTriggerCluster() - Large BC - tof %2.3f - Index %d\n",tofPatchMaxUS,fTriggerClusterIndex);
2842  else
2843  {
2844  fTriggerClusterIndex = -2;
2845  fTriggerClusterId = -2;
2846  }
2847  }
2848 
2849  if(tofPatchMax < 0) fTriggerClusterBC*=-1;
2850 
2851 
2852  // 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",
2853  // fTriggerClusterIndex, fTriggerClusterId,ePatchMax, tofPatchMax,
2854  // fTriggerClusterBC, fIsBadCellEvent,fIsBadMaxCellEvent,fIsExoticEvent, fIsTriggerMatch, nOfHighECl,absIdMaxMax);
2855  //
2856  // 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",
2857  // clusMax, idclusMax, eMax,tofMax, badClMax, badCeMax,exoMax);
2858 
2859  //Redo matching but open cuts
2860  if(!fIsTriggerMatch && fTriggerClusterId >= 0)
2861  {
2862  // Open time patch time
2863  TArrayI patchOpen = GetTriggerPatches(7,10);
2864 
2865  Int_t patchAbsIdOpenTime = -1;
2866  for(Int_t iabsId =0; iabsId < patchOpen.GetSize(); iabsId++)
2867  {
2868  Int_t absIDCell[4];
2869  patchAbsIdOpenTime = patchOpen.At(iabsId);
2870  GetCaloUtils()->GetEMCALGeometry()->GetCellIndexFromFastORIndex(patchAbsIdOpenTime, absIDCell);
2871  //if(tof > 75 ) printf("E %2.2f TOF %2.2f Trigger patch %d, cells : %d, %d, %d, %d\n",
2872  // clus->E(),tof,patches.At(iabsId), absIDCell[0],absIDCell[1],absIDCell[2],absIDCell[3]);
2873 
2874  for(Int_t ipatch = 0; ipatch < 4; ipatch++)
2875  {
2876  if(absIdMaxMax == absIDCell[ipatch])
2877  {
2878  fIsTriggerMatchOpenCut[0] = kTRUE;
2879  break;
2880  }
2881  }// cell patch loop
2882  }// trigger patch loop
2883 
2884  // Check neighbour patches
2885  Int_t patchAbsId = -1;
2886  Int_t globalCol = -1;
2887  Int_t globalRow = -1;
2888  GetCaloUtils()->GetEMCALGeometry()->GetFastORIndexFromCellIndex(absIdMaxMax, patchAbsId);
2889  GetCaloUtils()->GetEMCALGeometry()->GetPositionInEMCALFromAbsFastORIndex(patchAbsId,globalCol,globalRow);
2890 
2891  // Check patches with strict time cut
2892  Int_t patchAbsIdNeigh = -1;
2893  for(Int_t icol = globalCol-1; icol <= globalCol+1; icol++)
2894  {
2895  if(icol < 0 || icol > 47) continue;
2896 
2897  for(Int_t irow = globalRow; irow <= globalRow+1; irow++)
2898  {
2899  if(irow < 0 || irow > 63) continue;
2900 
2901  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(icol, irow, patchAbsIdNeigh);
2902 
2903  if ( patchAbsIdNeigh < 0 ) continue;
2904 
2905  for(Int_t iabsId =0; iabsId < patches.GetSize(); iabsId++)
2906  {
2907  if(patchAbsIdNeigh == patches.At(iabsId))
2908  {
2909  fIsTriggerMatchOpenCut[1] = kTRUE;
2910  break;
2911  }
2912  }// trigger patch loop
2913 
2914  }// row
2915  }// col
2916 
2917  // Check patches with open time cut
2918  Int_t patchAbsIdNeighOpenTime = -1;
2919  for(Int_t icol = globalCol-1; icol <= globalCol+1; icol++)
2920  {
2921  if(icol < 0 || icol > 47) continue;
2922 
2923  for(Int_t irow = globalRow; irow <= globalRow+1; irow++)
2924  {
2925  if(irow < 0 || irow > 63) continue;
2926 
2927  GetCaloUtils()->GetEMCALGeometry()->GetAbsFastORIndexFromPositionInEMCAL(icol, irow, patchAbsIdNeighOpenTime);
2928 
2929  if ( patchAbsIdNeighOpenTime < 0 ) continue;
2930 
2931  for(Int_t iabsId =0; iabsId < patchOpen.GetSize(); iabsId++)
2932  {
2933  if(patchAbsIdNeighOpenTime == patchOpen.At(iabsId))
2934  {
2935  fIsTriggerMatchOpenCut[2] = kTRUE;
2936  break;
2937  }
2938  }// trigger patch loop
2939 
2940  }// row
2941  }// col
2942 
2943  // printf("No match, new match: Open time %d-%d, open Neigh %d-%d, both open %d-%d\n",fIsTriggerMatchOpenCut[0],patchAbsIdOpenTime,
2944  // fIsTriggerMatchOpenCut[1],patchAbsIdNeigh,
2945  // fIsTriggerMatchOpenCut[2],patchAbsIdNeighOpenTime);
2946 
2947  patchOpen.Reset();
2948 
2949  }// No trigger match found
2950  //printf("Trigger BC %d, Id %d, Index %d\n",fTriggerClusterBC,fTriggerClusterId,fTriggerClusterIndex);
2951 }
2952 
2953 //_________________________________________________________
2957 //_________________________________________________________
2959 {
2961  {
2962  // get object pointer
2963  AliVCaloTrigger *caloTrigger = GetInputEvent()->GetCaloTrigger( "EMCAL" );
2964 
2965  if ( fBitEGA == 6 )
2966  {
2967  if (IsEventEMCALL1Gamma1()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(1);
2968  else if(IsEventEMCALL1Gamma2()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(3);
2969  else if(IsEventEMCALL1Jet1 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(0);
2970  else if(IsEventEMCALL1Jet2 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(2);
2971 
2972  // printf("L1 trigger Threshold Jet1 %f, Gamma1 %f, Jet2 %f, Gamma2 %f\n",
2973  // 0.07874*caloTrigger->GetL1Threshold(0),
2974  // 0.07874*caloTrigger->GetL1Threshold(1),
2975  // 0.07874*caloTrigger->GetL1Threshold(2),
2976  // 0.07874*caloTrigger->GetL1Threshold(3));
2977  }
2978  else
2979  {
2980  // Old AOD data format, in such case, order of thresholds different!!!
2981  if (IsEventEMCALL1Gamma1()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(0);
2982  else if(IsEventEMCALL1Gamma2()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(2);
2983  else if(IsEventEMCALL1Jet1 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(1);
2984  else if(IsEventEMCALL1Jet2 ()) fTriggerL1EventThreshold = 0.07874*caloTrigger->GetL1Threshold(3);
2985 
2986  // printf("L1 trigger Threshold Jet1 %f, Gamma1 %f, Jet2 %f, Gamma2 %f\n",
2987  // 0.07874*caloTrigger->GetL1Threshold(1),
2988  // 0.07874*caloTrigger->GetL1Threshold(0),
2989  // 0.07874*caloTrigger->GetL1Threshold(3),
2990  // 0.07874*caloTrigger->GetL1Threshold(2));
2991  }
2992  }
2993 
2994  // Set L0 threshold, if not set by user
2996  {
2997  // Revise for periods > LHC11d
2998  Int_t runNumber = fInputEvent->GetRunNumber();
2999  if (runNumber < 146861) fTriggerL0EventThreshold = 3. ;
3000  else if(runNumber < 154000) fTriggerL0EventThreshold = 4. ;
3001  else if(runNumber < 165000) fTriggerL0EventThreshold = 5.5;
3002  else fTriggerL0EventThreshold = 2 ;
3003  }
3004 }
3005 
3006 //________________________________________________________
3008 //________________________________________________________
3009 void AliCaloTrackReader::Print(const Option_t * opt) const
3010 {
3011  if(! opt)
3012  return;
3013 
3014  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
3015  printf("Task name : %s\n", fTaskName.Data()) ;
3016  printf("Data type : %d\n", fDataType) ;
3017  printf("CTS Min pT : %2.1f GeV/c\n", fCTSPtMin) ;
3018  printf("EMCAL Min pT : %2.1f GeV/c\n", fEMCALPtMin) ;
3019  printf("PHOS Min pT : %2.1f GeV/c\n", fPHOSPtMin) ;
3020  printf("CTS Max pT : %2.1f GeV/c\n", fCTSPtMax) ;
3021  printf("EMCAL Max pT : %2.1f GeV/c\n", fEMCALPtMax) ;
3022  printf("PHOS Max pT : %2.1f GeV/c\n", fPHOSPtMax) ;
3023  printf("EMCAL Bad Dist > %2.1f \n" , fEMCALBadChMinDist) ;
3024  printf("PHOS Bad Dist > %2.1f \n" , fPHOSBadChMinDist) ;
3025  printf("EMCAL N cells > %d \n" , fEMCALNCellsCut) ;
3026  printf("PHOS N cells > %d \n" , fPHOSNCellsCut) ;
3027  printf("EMCAL Time Cut: %3.1f < TOF < %3.1f\n", fEMCALTimeCutMin, fEMCALTimeCutMax);
3028  printf("Use CTS = %d\n", fFillCTS) ;
3029  printf("Use EMCAL = %d\n", fFillEMCAL) ;
3030  printf("Use DCAL = %d\n", fFillDCAL) ;
3031  printf("Use PHOS = %d\n", fFillPHOS) ;
3032  printf("Use EMCAL Cells = %d\n", fFillEMCALCells) ;
3033  printf("Use PHOS Cells = %d\n", fFillPHOSCells) ;
3034  printf("Track status = %d\n", (Int_t) fTrackStatus) ;
3035 
3036  printf("Track Mult Eta Cut = %2.2f\n", fTrackMultEtaCut) ;
3037 
3038  printf("Track Mult Pt Cuts:") ;
3039  for(Int_t icut = 0; icut < fTrackMultNPtCut; icut++) printf(" %2.2f GeV;",fTrackMultPtCut[icut]);
3040  printf(" \n") ;
3041 
3042  printf("Write delta AOD = %d\n", fWriteOutputDeltaAOD) ;
3043  printf("Recalculate Clusters = %d, E linearity = %d\n", fRecalculateClusters, fCorrectELinearity) ;
3044 
3045  printf("Use Triggers selected in SE base class %d; If not what Trigger Mask? %d; MB Trigger Mask for mixed %d \n",
3047 
3049  printf("Compare jet pt and pt hard to accept event, factor = %2.2f",fPtHardAndJetPtFactor);
3050 
3052  printf("Compare cluster pt and pt hard to accept event, factor = %2.2f",fPtHardAndClusterPtFactor);
3053 
3054  printf("Read Kine from, stack? %d, AOD ? %d \n", fReadStack, fReadAODMCParticles) ;
3055  printf("Delta AOD File Name = %s\n", fDeltaAODFileName.Data()) ;
3056  printf("Centrality: Class %s, Option %d, Bin [%d,%d] \n", fCentralityClass.Data(),fCentralityOpt,fCentralityBin[0], fCentralityBin[1]) ;
3057 
3058  printf(" \n") ;
3059 }
3060 
3061 //__________________________________________
3065 //__________________________________________
3067 {
3068  // Count number of cells with energy larger than 0.1 in SM3, cut on this number
3069  Int_t ncellsSM3 = 0;
3070  for(Int_t icell = 0; icell < fInputEvent->GetEMCALCells()->GetNumberOfCells(); icell++)
3071  {
3072  Int_t absID = fInputEvent->GetEMCALCells()->GetCellNumber(icell);
3073  Int_t sm = GetCaloUtils()->GetEMCALGeometry()->GetSuperModuleNumber(absID);
3074  if(fInputEvent->GetEMCALCells()->GetAmplitude(icell) > 0.1 && sm==3) ncellsSM3++;
3075  }
3076 
3077  Int_t ncellcut = 21;
3078  if(GetFiredTriggerClasses().Contains("EMC")) ncellcut = 35;
3079 
3080  if(ncellsSM3 >= ncellcut)
3081  {
3082  AliDebug(1,Form("Reject event with ncells in SM3 %d, cut %d, trig %s",
3083  ncellsSM3,ncellcut,GetFiredTriggerClasses().Data()));
3084  return kTRUE;
3085  }
3086 
3087  return kFALSE;
3088 }
3089 
3090 //_________________________________________________________
3094 //_________________________________________________________
3096 {
3097  if(label < 0) return ;
3098 
3099  AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fInputEvent) ;
3100  if(!evt) return ;
3101 
3102  TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
3103  if(!arr) return ;
3104 
3105  if(label < arr->GetEntriesFast())
3106  {
3107  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
3108  if(!particle) return ;
3109 
3110  if(label == particle->Label()) return ; // label already OK
3111  //else printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
3112  }
3113  //else printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
3114 
3115  // loop on the particles list and check if there is one with the same label
3116  for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
3117  {
3118  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
3119  if(!particle) continue ;
3120 
3121  if(label == particle->Label())
3122  {
3123  label = ind;
3124  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
3125  return;
3126  }
3127  }
3128 
3129  label = -1;
3130 
3131  //printf("AliCaloTrackReader::RemapMCLabelForAODs() - Label not found set to -1 \n");
3132 }
3133 
3134 //___________________________________
3136 //___________________________________
3138 {
3139  if(fCTSTracks) fCTSTracks -> Clear();
3140  if(fEMCALClusters) fEMCALClusters -> Clear("C");
3141  if(fPHOSClusters) fPHOSClusters -> Clear("C");
3142 
3143  fV0ADC[0] = 0; fV0ADC[1] = 0;
3144  fV0Mul[0] = 0; fV0Mul[1] = 0;
3145 
3146  if(fNonStandardJets) fNonStandardJets -> Clear("C");
3147  fBackgroundJets->Reset();
3148 }
3149 
3150 //___________________________________________
3154 //___________________________________________
3156 {
3157  fEventTrigMinBias = kFALSE;
3158  fEventTrigCentral = kFALSE;
3159  fEventTrigSemiCentral = kFALSE;
3160  fEventTrigEMCALL0 = kFALSE;
3161  fEventTrigEMCALL1Gamma1 = kFALSE;
3162  fEventTrigEMCALL1Gamma2 = kFALSE;
3163  fEventTrigEMCALL1Jet1 = kFALSE;
3164  fEventTrigEMCALL1Jet2 = kFALSE;
3165 
3166  AliDebug(1,Form("Select trigger mask bit %d - Trigger Event %s - Select <%s>",
3168 
3169  if(fEventTriggerMask <=0 )// in case no mask set
3170  {
3171  // EMC triggered event? Which type?
3172  if( GetFiredTriggerClasses().Contains("-B-") || GetFiredTriggerClasses().Contains("-S-") || GetFiredTriggerClasses().Contains("-I-") )
3173  {
3174  if ( GetFiredTriggerClasses().Contains("EGA" ) ||
3175  GetFiredTriggerClasses().Contains("EG1" ) )
3176  {
3177  fEventTrigEMCALL1Gamma1 = kTRUE;
3178  if( GetFiredTriggerClasses().Contains("EG1" ) && !fFiredTriggerClassName.Contains("EG1") ) fEventTrigEMCALL1Gamma1 = kFALSE;
3179  }
3180  else if( GetFiredTriggerClasses().Contains("EG2" ) )
3181  {
3182  fEventTrigEMCALL1Gamma2 = kTRUE;
3183  if( !fFiredTriggerClassName.Contains("EG2") ) fEventTrigEMCALL1Gamma2 = kFALSE;
3184  }
3185  else if( GetFiredTriggerClasses().Contains("EJE" ) ||
3186  GetFiredTriggerClasses().Contains("EJ1" ) )
3187  {
3188  fEventTrigEMCALL1Jet1 = kTRUE;
3189  if( GetFiredTriggerClasses().Contains("EJ1" ) && !fFiredTriggerClassName.Contains("EJ1") )
3190  fEventTrigEMCALL1Jet1 = kFALSE;
3191  }
3192  else if( GetFiredTriggerClasses().Contains("EJ2" ) )
3193  {
3194  fEventTrigEMCALL1Jet2 = kTRUE;
3195  if( !fFiredTriggerClassName.Contains("EJ2") ) fEventTrigEMCALL1Jet2 = kFALSE;
3196  }
3197  else if( GetFiredTriggerClasses().Contains("CEMC") &&
3198  !GetFiredTriggerClasses().Contains("EGA" ) &&
3199  !GetFiredTriggerClasses().Contains("EJE" ) &&
3200  !GetFiredTriggerClasses().Contains("EG1" ) &&
3201  !GetFiredTriggerClasses().Contains("EJ1" ) &&
3202  !GetFiredTriggerClasses().Contains("EG2" ) &&
3203  !GetFiredTriggerClasses().Contains("EJ2" ) )
3204  {
3205  fEventTrigEMCALL0 = kTRUE;
3206  }
3207 
3208  //Min bias event trigger?
3209  if (GetFiredTriggerClasses().Contains("CCENT_R2-B-NOPF-ALLNOTRD"))
3210  {
3211  fEventTrigCentral = kTRUE;
3212  }
3213  else if(GetFiredTriggerClasses().Contains("CSEMI_R1-B-NOPF-ALLNOTRD"))
3214  {
3215  fEventTrigSemiCentral = kTRUE;
3216  }
3217  else if((GetFiredTriggerClasses().Contains("CINT") || GetFiredTriggerClasses().Contains("CPBI2_B1") ) &&
3218  GetFiredTriggerClasses().Contains("-NOPF-ALLNOTRD") )
3219  {
3220  fEventTrigMinBias = kTRUE;
3221  }
3222  }
3223  }
3224  else
3225  {
3226  // EMC L1 Gamma
3227  if ( fEventTriggerMask & AliVEvent::kEMCEGA )
3228  {
3229  //printf("EGA trigger bit\n");
3230  if (GetFiredTriggerClasses().Contains("EG"))
3231  {
3232  if (GetFiredTriggerClasses().Contains("EGA")) fEventTrigEMCALL1Gamma1 = kTRUE;
3233  else
3234  {
3235  if(GetFiredTriggerClasses().Contains("EG1")) fEventTrigEMCALL1Gamma1 = kTRUE;
3236  if(GetFiredTriggerClasses().Contains("EG2")) fEventTrigEMCALL1Gamma2 = kTRUE;
3237  }
3238  }
3239  }
3240  // EMC L1 Jet
3241  else if( fEventTriggerMask & AliVEvent::kEMCEJE )
3242  {
3243  //printf("EGA trigger bit\n");
3244  if (GetFiredTriggerClasses().Contains("EJ"))
3245  {
3246  if (GetFiredTriggerClasses().Contains("EJE")) fEventTrigEMCALL1Jet1 = kTRUE;
3247  else
3248  {
3249  if(GetFiredTriggerClasses().Contains("EJ1")) fEventTrigEMCALL1Jet1 = kTRUE;
3250  if(GetFiredTriggerClasses().Contains("EJ2")) fEventTrigEMCALL1Jet2 = kTRUE;
3251  }
3252  }
3253  }
3254  // EMC L0
3255  else if((fEventTriggerMask & AliVEvent::kEMC7) ||
3256  (fEventTriggerMask & AliVEvent::kEMC1) )
3257  {
3258  //printf("L0 trigger bit\n");
3259  fEventTrigEMCALL0 = kTRUE;
3260  }
3261  // Min Bias Pb-Pb
3262  else if( fEventTriggerMask & AliVEvent::kCentral )
3263  {
3264  //printf("MB semi central trigger bit\n");
3265  fEventTrigSemiCentral = kTRUE;
3266  }
3267  // Min Bias Pb-Pb
3268  else if( fEventTriggerMask & AliVEvent::kSemiCentral )
3269  {
3270  //printf("MB central trigger bit\n");
3271  fEventTrigCentral = kTRUE;
3272  }
3273  // Min Bias pp, PbPb, pPb
3274  else if((fEventTriggerMask & AliVEvent::kMB ) ||
3275  (fEventTriggerMask & AliVEvent::kINT7) ||
3276  (fEventTriggerMask & AliVEvent::kINT8) ||
3277  (fEventTriggerMask & AliVEvent::kAnyINT) )
3278  {
3279  //printf("MB trigger bit\n");
3280  fEventTrigMinBias = kTRUE;
3281  }
3282  }
3283 
3284  AliDebug(1,Form("Event bits: \n \t MB %d, Cen %d, Sem %d, L0 %d, L1G1 %d, L1G2 %d, L1J1 %d, L1J2 %d",
3288 
3289  // L1 trigger bit
3290  if( fBitEGA == 0 && fBitEJE == 0 )
3291  {
3292  // Init the trigger bit once, correct depending on AliESD(AOD)CaloTrigger header version
3293 
3294  // Simpler way to do it ...
3295 // if( fInputEvent->GetRunNumber() < 172000 )
3296 // reader->SetEventTriggerL1Bit(4,5); // current LHC11 data
3297 // else
3298 // reader->SetEventTriggerL1Bit(6,8); // LHC12-13 data
3299 
3300  // Old values
3301  fBitEGA = 4;
3302  fBitEJE = 5;
3303 
3304  TFile* file = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
3305 
3306  const TList *clist = file->GetStreamerInfoCache();
3307 
3308  if(clist)
3309  {
3310  TStreamerInfo *cinfo = (TStreamerInfo*)clist->FindObject("AliESDCaloTrigger");
3311  Int_t verid = 5; // newer ESD header version
3312  if(!cinfo)
3313  {
3314  cinfo = (TStreamerInfo*)clist->FindObject("AliAODCaloTrigger");
3315  verid = 3; // newer AOD header version
3316  }
3317 
3318  if(cinfo)
3319  {
3320  Int_t classversionid = cinfo->GetClassVersion();
3321  //printf("********* Header class version %d *********** \n",classversionid);
3322 
3323  if (classversionid >= verid)
3324  {
3325  fBitEGA = 6;
3326  fBitEJE = 8;
3327  }
3328  } else AliInfo("AliCaloTrackReader()::SetEventTriggerBit() - Streamer info for trigger class not available, bit not changed");
3329  } else AliInfo("AliCaloTrackReader::SetEventTriggerBit() - Streamer list not available!, bit not changed");
3330 
3331  } // set once the EJE, EGA trigger bit
3332 }
3333 
3334 //____________________________________________________________
3337 //____________________________________________________________
3338 void AliCaloTrackReader::SetInputEvent(AliVEvent* const input)
3339 {
3340  fInputEvent = input;
3341  fMixedEvent = dynamic_cast<AliMixedEvent*>(GetInputEvent()) ;
3342  if (fMixedEvent)
3343  fNMixedEvent = fMixedEvent->GetNumberOfEvents() ;
3344 
3345  //Delete previous vertex
3346  if(fVertex)
3347  {
3348  for (Int_t i = 0; i < fNMixedEvent; i++)
3349  {
3350  delete [] fVertex[i] ;
3351  }
3352  delete [] fVertex ;
3353  }
3354 
3355  fVertex = new Double_t*[fNMixedEvent] ;
3356  for (Int_t i = 0; i < fNMixedEvent; i++)
3357  {
3358  fVertex[i] = new Double_t[3] ;
3359  fVertex[i][0] = 0.0 ;
3360  fVertex[i][1] = 0.0 ;
3361  fVertex[i][2] = 0.0 ;
3362  }
3363 }
3364 
3365 
Bool_t IsPileUpFromSPD() const
Double_t fEventWeight
Weight assigned to the event when filling histograms.
Bool_t fUseTrackDCACut
Do DCA selection.
TArrayI GetTriggerPatches(Int_t tmin, Int_t tmax)
virtual void FillInputVZERO()
virtual AliMCEvent * GetMC() const
Int_t fV0ADC[2]
Integrated V0 signal.
Bool_t 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.
double Double_t
Definition: External.C:58
virtual void FillInputNonStandardJets()
Double_t fEMCALTimeCutMax
Remove clusters/cells with time larger than this value, in ns.
Int_t fBitEJE
Trigger bit on VCaloTrigger for EJE.
TLorentzVector fMomentum
! Temporal TLorentzVector container, avoid declaration of TLorentzVectors per event.
Double_t fTimeStampRunMin
Minimum value of time stamp in run.
Bool_t fDoPileUpEventRejection
Select pile-up events by SPD.
void SetCentrality(Float_t cen)
Definition: AliAnaWeights.h:61
TObjArray * fPHOSClusters
Temporal array with PHOS CaloClusters.
TString GetGeneratorNameAndIndex(Int_t index, Int_t &genIndex) const
Double_t fTrackDCACut[3]
Remove tracks with DCA larger than cut, parameters of function stored here.
Bool_t fUseEventsWithPrimaryVertex
Select events with primary vertex.
virtual AliVCaloCells * GetPHOSCells() const
Bool_t fIsBadCellEvent
Bad cell triggered event flag, any cell in cluster is bad.
Bool_t IsCorrectionOfClusterEnergyOn() const
Bool_t fEventTrigEMCALL1Gamma1
Event is L1-Gamma, threshold 1 on its name, it should correspond kEMCEGA.
virtual AliHeader * GetHeader() const
Bool_t IsEventEMCALL1Jet1() const
Bool_t IsEventEMCALL1Gamma1() const
AliEMCALRecoUtils * GetEMCALRecoUtils() const
Bool_t ReadAODMCParticles() const
Int_t fSmearNLMMin
Do smearing for clusters with at least this value.
TString fEventPlaneMethod
Name of event plane method, by default "Q".
Int_t fTrackBCEventCut[19]
Fill one entry per event if there is a track in a given BC, depend on track pT, acceptance cut...
Bool_t fDoVertexBCEventSelection
Select events with vertex on BC=0 or -100.
AliMixedEvent * fMixedEvent
! Mixed event object. This class is not the owner.
virtual AliVEvent * GetInputEvent() const
Bool_t IsPileUpFromSPDAndNotEMCal() const
Check if event is from pile-up determined by SPD and not by EMCal.
Bool_t IsEventEMCALL1() const
Float_t fPtHardAndClusterPtFactor
Factor between ptHard and cluster pT to reject/accept event.
Bool_t fAcceptFastCluster
Accept events from fast cluster, exclude these events for LHC11a.
energy
Definition: HFPtSpectrum.C:44
Bool_t fSmearShowerShape
Smear shower shape (use in MC).
AliVEvent * fInputEvent
! pointer to esd or aod input.
Calculate the weight to the event to be applied when filling histograms.
Definition: AliAnaWeights.h:34
virtual void SetInputEvent(AliVEvent *input)
Double_t fTrackTimeCutMax
Remove tracks with time larger than this value, in ns.
Int_t fTriggerClusterBC
Event triggered by a cluster in BC -5 0 to 5.
virtual AliMultSelection * GetMultSelCen() const
TString fEMCALClustersListName
Alternative list of clusters produced elsewhere and not from InputEvent.
Float_t fPHOSBadChMinDist
Minimal distance to bad channel to accept cluster in PHOS, cm.
Bool_t fWriteOutputDeltaAOD
Write the created delta AOD objects into file.
Int_t fV0Mul[2]
Integrated V0 Multiplicity.
void RecalculateClusterDistanceToBadChannel(AliVCaloCells *cells, AliVCluster *clu)
Float_t fEMCALPtMin
pT Threshold on emcal clusters.
Bool_t fTriggerClusterTimeRecal
In case cluster already calibrated, do not try to recalibrate even if recalib on in AliEMCALRecoUtils...
Double_t ** fVertex
! Vertex array 3 dim for each mixed event buffer.
TH1F * fhCTSTrackCutsPt[6]
! Control histogram on the different CTS tracks selection cuts, pT
Int_t GetDebug() const
Definition: AliAnaWeights.h:48
UInt_t fEventTriggerMask
Select this triggerered event.
virtual Int_t GetEventCentrality() const
virtual AliGenEventHeader * GetGenEventHeader() const
virtual AliCentrality * GetCentrality() const
AliVCaloCells * fPHOSCells
! Temporal array with PHOS AliVCaloCells.
TList * fOutputContainer
! Output container with cut control histograms.
Bool_t IsEventEMCALL1Gamma2() const
virtual AliMixedEvent * GetMixedEvent() const
Int_t fEnergyHistogramNbins
Binning of the control histograms, min and max window.
Bool_t fFillInputBackgroundJetBranch
Flag to use data from background jets.
void SetTrackEventBC(Int_t bc)
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
Int_t fTrackMult[10]
Track multiplicity, count for different pT cuts.
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.
Float_t fTrackMultPtCut[10]
Track multiplicity and sum pt cuts list.
Bool_t IsPileUpFromSPDOrEMCal() const
Check if event is from pile-up determined by SPD or EMCal.
Bool_t IsPileUpFromEMCalAndNotSPD() const
Check if event is from pile-up determined by EMCal, not by SPD.
Bool_t fEventTrigEMCALL1Gamma2
Event is L1-Gamma, threshold 2 on its name, it should correspond kEMCEGA.
Bool_t fEventTrigEMCALL1Jet1
Event is L1-Gamma, threshold 1 on its name, it should correspond kEMCEGA.
Bool_t fFillEMCALCells
Use data from EMCAL.
Float_t fEnergyHistogramLimit[2]
Binning of the control histograms, number of bins.
Int_t fTriggerClusterId
Id of trigger cluster (cluster->GetID()).
Int_t fPHOSNCellsCut
Accept for the analysis PHOS clusters with more than fNCellsCut cells.
Float_t fCTSPtMax
pT Threshold on charged particles.
Int_t GetCocktailGeneratorAndIndex(Int_t index, TString &nameGen) const
Double_t fEMCALParamTimeCutMax[4]
Remove clusters/cells with time larger than parametrized value, in ns.
virtual void FillInputEMCALCells()
Connects the array with EMCAL cells and the pointer.
Bool_t fFillPHOSCells
Use data from PHOS.
virtual AliVCaloCells * GetEMCALCells() const
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.
int Int_t
Definition: External.C:63
void RecalculateClusterShowerShapeParameters(AliVCaloCells *cells, AliVCluster *clu)
Bool_t fSelectSPDHitTracks
Ensure that track hits SPD layers.
TString fFiredTriggerClassName
Name of trigger event type used to do the analysis.
virtual TClonesArray * GetAODMCParticles() const
Definition: External.C:204
TString GetFiredTriggerClasses() const
Bool_t ReadStack() const
unsigned int UInt_t
Definition: External.C:33
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.
float Float_t
Definition: External.C:68
virtual Bool_t CheckForPrimaryVertex() const
Int_t fEMCalBCEvent[19]
Fill one entry per event if there is a cluster in a given BC.
TString fInputBackgroundJetBranchName
Name of background jet branch.
virtual AliEventplane * GetEventPlane() const
Bool_t fEventTriggerAtSE
Select triggered event at SE base task or here.
Int_t fNMCGenerToAccept
Number of MC generators that should not be included in analysis.
Float_t fEMCALPtMax
pT Threshold on emcal clusters.
Int_t fEventType
Set the event species: 7 physics, 0 MC, 8 LED (not useful now)
virtual Bool_t ComparePtHardAndJetPt()
void RecalculateClusterTrackMatching(AliVEvent *event, TObjArray *clusterArray=0x0, AliMCEvent *mc=0x0)
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.
virtual TObjString * GetListOfParameters()
Save parameters used for analysis in a string.
Float_t fTrackMultEtaCut
Track multiplicity eta cut.
TString fMCGenerToAccept[5]
List with name of generators that should not be included.
Int_t fCentralityBin[2]
Minimum and maximum value of the centrality for the analysis.
virtual void FillInputPHOSCells()
Connects the array with PHOS cells and the pointer.
Base class for event, clusters and tracks filtering and preparation for the analysis.
TObjArray * fDCALClusters
Temporal array with DCAL CaloClusters, not needed in the normal case, use just EMCal array with DCal ...
Bool_t fUseTrackTimeCut
Do time cut selection.
virtual TString GetEventPlaneMethod() const
virtual Int_t GetTrackID(AliVTrack *track)
TArrayI fAcceptEventsWithBit
Accept events if trigger bit is on.
unsigned long ULong_t
Definition: External.C:38
TH1F * fhEMCALClusterCutsE[8]
! Control histogram on the different EMCal cluster selection cuts, E
virtual 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.
Int_t GetNumberOfLocalMaxima(AliVCluster *cluster, AliVCaloCells *cells)
Find the number of local maxima in cluster.
Bool_t IsWeightSettingOn() const
Definition: AliAnaWeights.h:78
Int_t fTrackBCEvent[19]
Fill one entry per event if there is a track in a given BC.
virtual void FillInputEMCALAlgorithm(AliVCluster *clus, Int_t iclus)
Bool_t fComparePtHardAndJetPt
In MonteCarlo, jet events, reject fake events with wrong jet energy.
Int_t fNNonPileUpClusters
Number of clusters with time below 20 ns.
Int_t fNMixedEvent
Number of events in mixed event buffer.
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Int_t fCentralityOpt
Option for the returned value of the centrality, possible options 5, 10, 100.
Float_t fZvtxCut
Cut on vertex position.
void RecalculateClusterPosition(AliVCaloCells *cells, AliVCluster *clu)
AliCaloTrackReader()
Constructor. Initialize parameters.
Bool_t fTriggerPatchClusterMatch
Search for the trigger patch and check if associated cluster was the trigger.
Bool_t fUseParamTimeCut
Use simple or parametrized time cut.
Int_t GetNMaskCellColumns() const
AliFiducialCut * fFiducialCut
Acceptance cuts.
Bool_t fRejectEMCalTriggerEventsWith2Tresholds
Reject events EG2 also triggered by EG1 or EJ2 also triggered by EJ1.
TH1F * fhPHOSClusterCutsE[7]
! Control histogram on the different PHOS cluster selection cuts, E
virtual void FillInputEMCAL()
Bool_t fFillEMCAL
Use data from EMCAL.
Bool_t fFillPHOS
Use data from PHOS.
AliVCaloCells * fEMCALCells
! Temporal array with EMCAL AliVCaloCells.
virtual AliAODMCHeader * GetAODMCHeader() const
Int_t fEMCALNCellsCut
Accept for the analysis EMCAL clusters with more than fNCellsCut cells.
void CorrectClusterEnergy(AliVCluster *cl)
Correct cluster energy non linearity.
Float_t fCTSPtMin
pT Threshold on charged particles.
virtual AliStack * GetStack() const
Int_t fSmearNLMMax
Do smearing for clusters with at maximum this value.
Float_t fTriggerL0EventThreshold
L0 Threshold to look for triggered events, set outside.
Float_t fTimeStampEventFracMax
Maximum value of time stamp fraction event.
virtual void ResetLists()
Reset lists, called in AliAnaCaloTrackCorrMaker.
Int_t fNPileUpClusters
Number of clusters with time avobe 20 ns.
AliAODEvent * fOutputEvent
! pointer to aod output.
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Bool_t fRecalculateClusters
Correct clusters, recalculate them if recalibration parameters is given.
Bool_t fDoRejectNoTrackEvents
Reject events with no selected tracks in event.
Bool_t IsEventEMCALL1Jet2() const
Float_t fEMCALBadChMinDist
Minimal distance to bad channel to accept cluster in EMCal, cell units.
Bool_t 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
void SetPythiaEventHeader(AliGenPythiaEventHeader *py)
Definition: AliAnaWeights.h:83
TObjArray * fCTSTracks
Temporal array with tracks.
Int_t fTriggerPatchTimeWindow[2]
Trigger patch selection window.
unsigned short UShort_t
Definition: External.C:28
Bool_t AcceptParticleMCLabel(Int_t mcLabel) const
AliMCEvent * fMC
! Monte Carlo Event Handler.
TH1I * fhNEventsAfterCut
! Each bin represents number of events resulting after a given selection cut: vertex, trigger, ...
Float_t fPtHardAndJetPtFactor
Factor between ptHard and jet pT to reject/accept event.
TString fMCGenerEventHeaderToAccept
Accept events that contain at least this event header name.
const char Option_t
Definition: External.C:48
TString fDeltaAODFileName
Delta AOD file name.
Int_t GetVertexBC() const
Float_t fTrackSumPt[10]
Track sum pT, count for different pT cuts.
Bool_t IsPileUpFromEMCal() const
Check if event is from pile-up determined by EMCal.
virtual Double_t GetWeight()
Bool_t fFillInputNonStandardJetBranch
Flag to use data from non standard jets.
TClonesArray * fNonStandardJets
! Temporal array with jets.
Int_t fEMCalBCEventCut[19]
Fill one entry per event if there is a cluster in a given BC, depend on cluster E, acceptance cut.
Int_t fTriggerClusterIndex
Index in clusters array of trigger cluster.
Bool_t fTriggerL1EventThresholdFix
L1 Threshold is fix and set outside.
bool Bool_t
Definition: External.C:53
Float_t fSmearShowerShapeWidth
Smear shower shape landau function "width" (use in MC).
Int_t fBitEGA
Trigger bit on VCaloTrigger for EGA.
Int_t fTrackMultNPtCut
Track multiplicty, number of pt cuts.
virtual void FillInputBackgroundJets()
Bool_t fDoV0ANDEventSelection
Select events depending on V0AND.
AliAODJetEventBackground * fBackgroundJets
! Background jets.
Bool_t fCorrectELinearity
Correct cluster linearity, always on.
virtual TList * GetCreateControlHistograms()
Int_t fVertexBC
Vertex BC.
Bool_t fIsExoticEvent
Exotic trigger event flag.
TString fInputNonStandardJetBranchName
Name of non standard jet branch.
TArrayI fRejectEventsWithBit
Reject events if trigger bit is on.
Bool_t fUseEMCALTimeCut
Do time cut selection.
Float_t GetPhi(Float_t phi) const
Shift phi angle in case of negative value 360 degrees. Example TLorenzVector::Phi defined in -pi to p...
Int_t fDataType
Select MC: Kinematics, Data: ESD/AOD, MCData: Both.
virtual Bool_t SelectTrack(AliVTrack *, Double_t *)
Int_t GetMaxEnergyCell(AliVCaloCells *cells, AliVCluster *clu, Float_t &fraction) const
For a given CaloCluster, it gets the absId of the cell with maximum energy deposit.
Bool_t fFillDCAL
Use data from DCAL, not needed in the normal case, use just EMCal array with DCal limits...
void SetTrackEventBCcut(Int_t bc)
Bool_t fIsTriggerMatch
Could match the event to a trigger patch?
Float_t fTriggerL1EventThreshold
L1 Threshold to look for triggered events, set in data.
TString fCentralityClass
Name of selected centrality class.
Int_t fMCGenerIndexToAccept[5]
List with index of generators that should not be included.
ULong_t fTrackStatus
Track selection bit, select tracks refitted in TPC, ITS ...
Bool_t IsPileUpFromSPDAndEMCal() const
Check if event is from pile-up determined by SPD and EMCal.
virtual void FillInputPHOS()
Fill the array with PHOS filtered clusters.
virtual TObjArray * GetCTSTracks() const
Bool_t fUseAliCentrality
Select as centrality estimator AliCentrality (Run1) or AliMultSelection (Run1 and Run2) ...
Int_t fSmearingFunction
Choice of smearing function. 0 no smearing. 1 smearing from Gustavo (Landau center at 0)...
Bool_t fRecalculateVertexBC
Recalculate vertex BC from tracks pointing to vertex.
Float_t RecalibrateClusterEnergy(AliVCluster *cluster, AliVCaloCells *cells)
Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that co...
void SetEMCalEventBC(Int_t bc)
Bool_t reject
Bool_t CheckCellFiducialRegion(AliVCluster *cluster, AliVCaloCells *cells) const
Bool_t IsEventEMCALL0() const