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