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