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