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