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