]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGGA/CaloTasks/AliAnalysisTaskCaloFilter.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGGA / CaloTasks / AliAnalysisTaskCaloFilter.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//////////////////////////////////////////////////////////
17// Filter the ESDCaloClusters and ESDCaloCells of EMCAL,
18// PHOS or both, creating the corresponing AODCaloClusters
19// and AODCaloCells.
20// Keep also the AODHeader information and the vertex.
21// Keep tracks, v0s, VZERO if requested
22// Select events containing a cluster or track avobe a given threshold
23// Copy of AliAnalysisTaskESDfilter.
24// Author: Gustavo Conesa Balbastre (INFN - Frascati)
25//////////////////////////////////////////////////////////
26
27//Root
28#include "TGeoManager.h"
29#include "TFile.h"
30#include "TROOT.h"
31#include "TInterpreter.h"
32
33//STEER
34#include "AliESDEvent.h"
35#include "AliAODEvent.h"
36#include "AliLog.h"
37#include "AliVCluster.h"
38#include "AliVCaloCells.h"
39#include "AliVEventHandler.h"
40#include "AliAODHandler.h"
41#include "AliAnalysisManager.h"
42#include "AliInputEventHandler.h"
43
44//EMCAL
45#include "AliEMCALRecoUtils.h"
46#include "AliEMCALGeometry.h"
47
48#include "AliAnalysisTaskCaloFilter.h"
49
50ClassImp(AliAnalysisTaskCaloFilter)
51
52////////////////////////////////////////////////////////////////////////
53
54AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter():
55AliAnalysisTaskSE("CaloFilterTask"),
56fCaloFilter(0), fEventSelection(),
57fAcceptAllMBEvent(kFALSE),fMBTriggerMask(AliVEvent::kMB),
58fCorrect(kFALSE),
59fEMCALGeo(0x0), fEMCALGeoName("EMCAL_COMPLETE12SMV1"),
60fEMCALRecoUtils(new AliEMCALRecoUtils),
61fLoadEMCALMatrices(kFALSE), //fLoadPHOSMatrices(kFALSE),
62fGeoMatrixSet(kFALSE),
63fConfigName(""), fFillAODFile(kTRUE),
64fFillMCParticles(kFALSE),
65fFillTracks(kFALSE), fFillHybridTracks(kFALSE),
66fFillAllVertices(kFALSE), fFillv0s(kFALSE),
67fFillVZERO(kFALSE),
68fEMCALEnergyCut(0.), fEMCALNcellsCut (0),
69fPHOSEnergyCut(0.), fPHOSNcellsCut (0),
70fTrackPtCut(-1),
71fVzCut(100.), fCheckEventVertex(kTRUE),
72fEvent(0x0),
73fESDEvent(0x0), fAODEvent(0x0)
74{
75 // Default constructor
76
77 fEventSelection[0] = kFALSE;
78 fEventSelection[1] = kFALSE;
79 fEventSelection[2] = kFALSE;
80
81 for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
82 //for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix[i] = 0 ;
83
84}
85
86//__________________________________________________
87AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter(const char* name):
88AliAnalysisTaskSE(name),
89fCaloFilter(0), fEventSelection(),
90fAcceptAllMBEvent(kFALSE),fMBTriggerMask(AliVEvent::kMB),
91fCorrect(kFALSE),
92fEMCALGeo(0x0), fEMCALGeoName("EMCAL_COMPLETE12SMV1"),
93fEMCALRecoUtils(new AliEMCALRecoUtils),
94fLoadEMCALMatrices(kFALSE), //fLoadPHOSMatrices(kFALSE),
95fGeoMatrixSet(kFALSE),
96fConfigName(""), fFillAODFile(kTRUE),
97fFillMCParticles(kFALSE),
98fFillTracks(kFALSE), fFillHybridTracks(kFALSE),
99fFillAllVertices(kFALSE), fFillv0s(kFALSE),
100fFillVZERO(kFALSE),
101fEMCALEnergyCut(0.), fEMCALNcellsCut(0),
102fPHOSEnergyCut(0.), fPHOSNcellsCut(0),
103fTrackPtCut(-1),
104fVzCut(100.), fCheckEventVertex(kTRUE),
105fEvent(0x0),
106fESDEvent(0x0), fAODEvent(0x0)
107{
108 // Constructor
109
110 fEventSelection[0] = kFALSE;
111 fEventSelection[1] = kFALSE;
112 fEventSelection[2] = kFALSE;
113
114 for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
115 //for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix[i] = 0 ;
116
117}
118
119//__________________________________________________
120AliAnalysisTaskCaloFilter::~AliAnalysisTaskCaloFilter()
121{
122 //Destructor.
123
124 if(fEMCALGeo) delete fEMCALGeo;
125 if(fEMCALRecoUtils) delete fEMCALRecoUtils;
126
127}
128
129//_____________________________________________
130Bool_t AliAnalysisTaskCaloFilter::AcceptEvent()
131{
132 // Define conditions to accept the event to be filtered
133
134 if(!AcceptEventVertex()) return kFALSE;
135
136 Bool_t eventSel = kFALSE;
137
138 Bool_t isMB = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fMBTriggerMask);
139
140 if ( isMB && fAcceptAllMBEvent ) eventSel = kTRUE; // accept any MB event
141
142 else if( fEventSelection[0] && AcceptEventEMCAL() ) eventSel = kTRUE; // accept event depending on EMCAL activity
143
144 else if( fEventSelection[1] && AcceptEventPHOS () ) eventSel = kTRUE; // accept event depending on PHOS activity
145
146 else if( fEventSelection[2] && AcceptEventTrack() ) eventSel = kTRUE; // accept event depending on Track activity
147
148 return eventSel ;
149
150}
151
152//__________________________________________________
153Bool_t AliAnalysisTaskCaloFilter::AcceptEventEMCAL()
154{
155 // Accept event given there is a EMCAL cluster with enough energy, and not noisy, exotic
156
157 if(fCaloFilter==kPHOS) return kTRUE; // accept
158
159 if(fEMCALEnergyCut <= 0) return kTRUE; // accept
160
161 Int_t nCluster = InputEvent() -> GetNumberOfCaloClusters();
162 AliVCaloCells * caloCell = InputEvent() -> GetEMCALCells();
163 Int_t bc = InputEvent() -> GetBunchCrossNumber();
164
165 for(Int_t icalo = 0; icalo < nCluster; icalo++)
166 {
167 AliVCluster *clus = (AliVCluster*) (InputEvent()->GetCaloCluster(icalo));
168
169 if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fEMCALNcellsCut ) && ( clus->E() > fEMCALEnergyCut ) &&
170 fEMCALRecoUtils->IsGoodCluster(clus,fEMCALGeo,caloCell,bc))
171 {
172
173 AliDebug(1,Form("Accept : E %2.2f > %2.2f, nCells %d > %d",
174 clus->E(), fEMCALEnergyCut, clus->GetNCells(), fEMCALNcellsCut));
175
176 return kTRUE;
177 }
178
179 }// loop
180
181 AliDebug(1,"Reject");
182
183 //printf("Fired %s\n",((AliESDEvent*)InputEvent())->GetFiredTriggerClasses().Data());
184
185 return kFALSE;
186
187}
188
189//_________________________________________________
190Bool_t AliAnalysisTaskCaloFilter::AcceptEventPHOS()
191{
192 // Accept event given there is a PHOS cluster with enough energy and not noisy/exotic
193
194 if(fCaloFilter==kEMCAL) return kTRUE; // accept
195
196 if(fPHOSEnergyCut <= 0) return kTRUE; // accept
197
198 Int_t nCluster = InputEvent() -> GetNumberOfCaloClusters();
199
200 for(Int_t icalo = 0; icalo < nCluster; icalo++)
201 {
202 AliVCluster *clus = (AliVCluster*) (InputEvent()->GetCaloCluster(icalo));
203
204 if( ( clus->IsPHOS() ) && ( clus->GetNCells() > fPHOSNcellsCut ) && ( clus->E() > fPHOSEnergyCut ))
205 {
206
207 AliDebug(1,Form("Accept : E %2.2f > %2.2f, nCells %d > %d",
208 clus->E(), fPHOSEnergyCut, clus->GetNCells(), fPHOSNcellsCut));
209
210 return kTRUE;
211 }
212
213 }// loop
214
215 AliDebug(1,"Reject");
216
217 return kFALSE;
218
219}
220
221//__________________________________________________
222Bool_t AliAnalysisTaskCaloFilter::AcceptEventTrack()
223{
224 // Accept event if there is a track avobe a certain pT
225
226 if(fTrackPtCut <= 0) return kTRUE; // accept
227
228 Double_t pTrack[3] = {0,0,0};
229
230 for (Int_t nTrack = 0; nTrack < fEvent->GetNumberOfTracks(); ++nTrack)
231 {
232 AliVTrack *track = (AliVTrack*) fEvent->GetTrack(nTrack);
233
234 // Select only hybrid tracks?
235 if(fAODEvent && fFillHybridTracks && !((AliAODTrack*)track)->IsHybridGlobalConstrainedGlobal()) continue;
236
237 track->GetPxPyPz(pTrack) ;
238
239 TLorentzVector momentum(pTrack[0],pTrack[1],pTrack[2],0);
240
241 if(momentum.Pt() > fTrackPtCut)
242 {
243 AliDebug(1,Form("Accept : pT %2.2f > %2.2f",momentum.Pt(), fTrackPtCut));
244
245 return kTRUE;
246 }
247
248 }
249
250 AliDebug(1,"Reject");
251
252 return kFALSE;
253
254}
255
256//___________________________________________________
257Bool_t AliAnalysisTaskCaloFilter::AcceptEventVertex()
258{
259 // Accept event with good vertex
260
261 Double_t v[3];
262 InputEvent()->GetPrimaryVertex()->GetXYZ(v) ;
263
264 if(TMath::Abs(v[2]) > fVzCut)
265 {
266 AliDebug(1,Form("Vz Reject : vz %2.2f > %2.2f",v[2],fVzCut));
267
268 return kFALSE ;
269 }
270
271 return CheckForPrimaryVertex();
272}
273
274//_______________________________________________________
275Bool_t AliAnalysisTaskCaloFilter::CheckForPrimaryVertex()
276{
277 //Check if the vertex was well reconstructed, copy from v0Reader of conversion group
278 //It only works for ESDs
279
280 if(!fCheckEventVertex) return kTRUE;
281
282 // AODs
283 if(!fESDEvent)
284 {
285 // Check that the vertex is not (0,0,0)
286 Double_t v[3];
287 InputEvent()->GetPrimaryVertex()->GetXYZ(v) ;
288
289 if(TMath::Abs(v[2]) < 1e-6 &&
290 TMath::Abs(v[1]) < 1e-6 &&
291 TMath::Abs(v[0]) < 1e-6 )
292 {
293 AliDebug(1,"Reject v(0,0,0)");
294
295 return kFALSE ;
296 }
297
298 return kTRUE;
299 }
300
301 // ESDs
302 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors() > 0)
303 {
304 return kTRUE;
305 }
306
307 if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors() < 1)
308 {
309 // SPD vertex
310 if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors() > 0)
311 {
312 //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
313 return kTRUE;
314
315 }
316 if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors() < 1)
317 {
318 // cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
319 AliDebug(1,"Reject, GetPrimaryVertexSPD()->GetNContributors() < 1");
320
321 return kFALSE;
322 }
323 }
324
325 AliDebug(1,"Reject, GetPrimaryVertexTracks()->GetNContributors() > 1");
326
327 return kFALSE;
328
329}
330
331//__________________________________________________
332void AliAnalysisTaskCaloFilter::CorrectionsInEMCAL()
333{
334 //If EMCAL, and requested, correct energy, position ...
335
336 //Need to do this in a separate loop before filling the ESDs because of the track matching recalculations
337
338 if(fCorrect && (fCaloFilter==kEMCAL || fCaloFilter==kBoth) )
339 {
340 if(!fGeoMatrixSet)
341 {
342 if(fLoadEMCALMatrices)
343 {
344 AliInfo("Load user defined EMCAL geometry matrices");
345 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
346 {
347 if(fEMCALMatrix[mod])
348 {
349 if(DebugLevel() > 1)
350 fEMCALMatrix[mod]->Print();
351 fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
352 }
353 fGeoMatrixSet=kTRUE;
354 }//SM loop
355 }//Load matrices
356 else if(!gGeoManager)
357 {
358 AliInfo("Get geo matrices from data");
359 //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
360 if(!strcmp(InputEvent()->GetName(),"AliAODEvent"))
361 {
362 AliDebug(1,"Use ideal geometry, values geometry matrix not kept in AODs");
363 }//AOD
364 else
365 {
366 AliDebug(1,"Load Misaligned matrices");
367 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
368 if(!esd)
369 {
370 AliInfo("This event does not contain ESDs?");
371 return;
372 }
373 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
374 {
375 if(DebugLevel() > 1)
376 esd->GetEMCALMatrix(mod)->Print();
377 if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
378 }
379 fGeoMatrixSet=kTRUE;
380 }//ESD
381 }//Load matrices from Data
382
383 }//first event
384
385 //Cluster Loop
386 Int_t nCaloClus = InputEvent()->GetNumberOfCaloClusters();
387
388 for (Int_t iClust=0; iClust<nCaloClus; ++iClust)
389 {
390 AliVCluster * cluster = InputEvent()->GetCaloCluster(iClust);
391
392 if(cluster->IsPHOS()) continue ;
393
394 Float_t position[]={0,0,0};
395
396 AliDebug(1,Form("Check cluster %d for bad channels and close to border",cluster->GetID()));
397
398 if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
399
400 AliDebug(2,Form("Filter, before : i %d, E %f, dispersion %f, m02 %f, m20 %f, distToBad %f",iClust,cluster->E(),
401 cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20(), cluster->GetDistanceToBadChannel()));
402 cluster->GetPosition(position);
403 AliDebug(2,Form("Filter, before : i %d, x %f, y %f, z %f",cluster->GetID(), position[0], position[1], position[2]));
404
405 //Recalculate distance to bad channels, if new list of bad channels provided
406 fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, InputEvent()->GetEMCALCells(), cluster);
407
408 if(fEMCALRecoUtils->IsRecalibrationOn())
409 {
410 fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, cluster, InputEvent()->GetEMCALCells());
411 fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, InputEvent()->GetEMCALCells(),cluster);
412 fEMCALRecoUtils->RecalculateClusterPID(cluster);
413 }
414
415 fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, InputEvent()->GetEMCALCells(),cluster);
416
417 AliDebug(2,Form("Filter, after : i %d, E %f, dispersion %f, m02 %f, m20 %f, distToBad %f",cluster->GetID(),cluster->E(),
418 cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20(), cluster->GetDistanceToBadChannel()));
419 cluster->GetPosition(position);
420 AliDebug(1,Form("Filter, after : i %d, x %f, y %f, z %f",cluster->GetID(), position[0], position[1], position[2]));
421
422 cluster->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(cluster));
423
424 }
425
426 //Recalculate track-matching
427 fEMCALRecoUtils->FindMatches(InputEvent(),0,fEMCALGeo);
428
429 } // corrections in EMCAL
430}
431
432//________________________________________________
433void AliAnalysisTaskCaloFilter::FillAODCaloCells()
434{
435 // Fill EMCAL/PHOS cell info
436
437 // EMCAL
438 if ((fCaloFilter==kBoth || fCaloFilter==kEMCAL) && fEvent->GetEMCALCells())
439 { // protection against missing ESD information
440 AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
441 Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
442
443 AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
444 aodEMcells.CreateContainer(nEMcell);
445 aodEMcells.SetType(AliVCaloCells::kEMCALCell);
446 Double_t calibFactor = 1.;
447 for (Int_t iCell = 0; iCell < nEMcell; iCell++)
448 {
449 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
450 fEMCALGeo->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
451 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
452
453 if(fCorrect && fEMCALRecoUtils->IsRecalibrationOn())
454 {
455 calibFactor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
456 }
457
458 if(!fEMCALRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
459 { //Channel is not declared as bad
460 aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
461 eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
462 //printf("GOOD channel\n");
463 }
464 else
465 {
466 aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
467 //printf("BAD channel\n");
468 }
469 }
470 aodEMcells.Sort();
471 }
472
473 // PHOS
474 if ((fCaloFilter==kBoth || fCaloFilter==kPHOS) && fEvent->GetPHOSCells())
475 { // protection against missing ESD information
476 AliVCaloCells &eventPHcells = *(fEvent->GetPHOSCells());
477 Int_t nPHcell = eventPHcells.GetNumberOfCells() ;
478
479 AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
480 aodPHcells.CreateContainer(nPHcell);
481 aodPHcells.SetType(AliVCaloCells::kPHOSCell);
482
483 for (Int_t iCell = 0; iCell < nPHcell; iCell++)
484 {
485 aodPHcells.SetCell(iCell,eventPHcells.GetCellNumber(iCell),eventPHcells.GetAmplitude(iCell),
486 eventPHcells.GetTime(iCell),eventPHcells.GetMCLabel(iCell),eventPHcells.GetEFraction(iCell));
487 }
488
489 aodPHcells.Sort();
490 }
491}
492
493
494//___________________________________________________
495void AliAnalysisTaskCaloFilter::FillAODCaloClusters()
496{
497 // Fill the AOD with caloclusters
498
499 // Access to the AOD container of clusters
500
501 TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
502 Int_t jClusters=0;
503 Float_t posF[3] ;
504
505 Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
506 for (Int_t iClust=0; iClust<nCaloClus; ++iClust)
507 {
508
509 AliVCluster * cluster = fEvent->GetCaloCluster(iClust);
510
511 //Check which calorimeter information we want to keep.
512
513 if(fCaloFilter!=kBoth)
514 {
515 if (fCaloFilter==kPHOS && cluster->IsEMCAL()) continue ;
516 else if(fCaloFilter==kEMCAL && cluster->IsPHOS()) continue ;
517 }
518
519 // Get original residuals, in case of previous recalculation, reset them
520 Float_t dR = cluster->GetTrackDx();
521 Float_t dZ = cluster->GetTrackDz();
522
523 AliDebug(2,Form("Original residuals : dZ %f, dR %f",dZ, dR));
524
525 //--------------------------------------------------------------
526 //If EMCAL and corrections done, get the new matching parameters, do not copy noisy clusters
527 if(cluster->IsEMCAL() && fCorrect)
528 {
529 AliDebug(2,Form("Check cluster %d for bad channels and close to border",cluster->GetID()));
530
531 if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
532
533 if(fEMCALRecoUtils->IsExoticCluster(cluster, InputEvent()->GetEMCALCells(),InputEvent()->GetBunchCrossNumber())) continue;
534
535 fEMCALRecoUtils->GetMatchedResiduals(cluster->GetID(),dR,dZ);
536 cluster->SetTrackDistance(dR,dZ);
537 }
538
539 AliDebug(2,Form("EMCAL? %d, PHOS? %d Track-Cluster Residuals : dZ %f, dR %f",cluster->IsEMCAL(), cluster->IsPHOS(),dZ, dR));
540
541 //--------------------------------------------------------------
542
543 //Now fill AODs
544
545 Int_t id = cluster->GetID();
546 Float_t energy = cluster->E();
547 cluster->GetPosition(posF);
548
549 AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
550 AliAODCaloCluster(id,
551 cluster->GetNLabels(),
552 cluster->GetLabels(),
553 energy,
554 posF,
555 NULL,
556 cluster->GetType());
557
558 caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
559 cluster->GetDispersion(),
560 cluster->GetM20(), cluster->GetM02(),
561 -1,
562 cluster->GetNExMax(),cluster->GetTOF()) ;
563
564 caloCluster->SetPIDFromESD(cluster->GetPID());
565 caloCluster->SetNCells(cluster->GetNCells());
566 caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
567 caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
568 caloCluster->SetTrackDistance(dR, dZ);
569
570 AliDebug(2,Form("Filter, aod : i %d, E %f, dispersion %f, m02 %f, m20 %f",caloCluster->GetID(),caloCluster->E(),
571 caloCluster->GetDispersion(),caloCluster->GetM02(),caloCluster->GetM20()));
572 caloCluster->GetPosition(posF);
573 AliDebug(2,Form("Filter, aod : i %d, x %f, y %f, z %f",caloCluster->GetID(), posF[0], posF[1], posF[2]));
574
575 //Matched tracks, just to know if there was any match, the track pointer is useless if tracks not stored
576 if(TMath::Abs(dR) < 990 && TMath::Abs(dZ) < 990)
577 { //Default value in PHOS 999, in EMCAL 1024, why?
578 caloCluster->AddTrackMatched(new AliAODTrack);
579 }
580 // TO DO, in case Tracks available, think how to put the matched track in AOD
581 }
582
583 caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots
584
585}
586
587//__________________________________________________
588void AliAnalysisTaskCaloFilter::FillAODCaloTrigger()
589{
590 // AOD CaloTrigger copy
591
592 if( !AODEvent() || !fAODEvent ) return;
593
594 AliAODCaloTrigger* triggerEM = AODEvent()->GetCaloTrigger("EMCAL");
595 AliAODCaloTrigger* triggerPH = AODEvent()->GetCaloTrigger("PHOS");
596
597 // Copy from AODs
598
599 AliAODCaloTrigger* inTriggerEM = fAODEvent ->GetCaloTrigger("EMCAL");
600 AliAODCaloTrigger* inTriggerPH = fAODEvent ->GetCaloTrigger("PHOS");
601
602 if(inTriggerPH && (fCaloFilter==kBoth || fCaloFilter==kPHOS)) *triggerPH = *inTriggerPH;
603
604 if(inTriggerEM && (fCaloFilter==kBoth || fCaloFilter==kEMCAL)) *triggerEM = *inTriggerEM;
605}
606
607//______________________________________________
608void AliAnalysisTaskCaloFilter::FillAODHeader()
609{
610 // AOD header copy
611
612 AliAODHeader* header = dynamic_cast<AliAODHeader*>(AODEvent()->GetHeader());
613 if(!header)
614 {
615 AliFatal("Not a standard AOD");
616 return; // not needed but coverity complains
617 }
618
619 // Copy from AODs
620 if(fAODEvent)
621 {
622 *header = *((AliAODHeader*)fAODEvent->GetHeader());
623 return;
624 }
625
626 if(!fESDEvent) return;
627
628 // Copy from ESDs
629
630 header->SetRunNumber(fEvent->GetRunNumber());
631
632 TTree* tree = fInputHandler->GetTree();
633 if (tree)
634 {
635 TFile* file = tree->GetCurrentFile();
636 if (file) header->SetESDFileName(file->GetName());
637 }
638
639 header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
640 header->SetOrbitNumber(fEvent->GetOrbitNumber());
641 header->SetPeriodNumber(fEvent->GetPeriodNumber());
642 header->SetEventType(fEvent->GetEventType());
643
644 //Centrality
645 if(fEvent->GetCentrality())
646 {
647 header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
648 }
649 else
650 {
651 header->SetCentrality(0);
652 }
653
654 //Trigger
655 header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
656 header->SetFiredTriggerClasses(fESDEvent->GetFiredTriggerClasses());
657 header->SetTriggerMask(fEvent->GetTriggerMask());
658 header->SetTriggerCluster(fEvent->GetTriggerCluster());
659 header->SetL0TriggerInputs(fESDEvent->GetHeader()->GetL0TriggerInputs());
660 header->SetL1TriggerInputs(fESDEvent->GetHeader()->GetL1TriggerInputs());
661 header->SetL2TriggerInputs(fESDEvent->GetHeader()->GetL2TriggerInputs());
662
663 header->SetMagneticField(fEvent->GetMagneticField());
664 header->SetMuonMagFieldScale(fESDEvent->GetCurrentDip()/6000.);
665
666 header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
667 header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
668 header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
669 header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
670 header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
671
672 Float_t diamxy[2]={(Float_t)fEvent->GetDiamondX(),(Float_t)fEvent->GetDiamondY()};
673 Float_t diamcov[3];
674 fEvent->GetDiamondCovXY(diamcov);
675 header->SetDiamond(diamxy,diamcov);
676 header->SetDiamondZ(fESDEvent->GetDiamondZ(),fESDEvent->GetSigma2DiamondZ());
677
678}
679
680
681//__________________________________________________
682void AliAnalysisTaskCaloFilter::FillAODMCParticles()
683{
684 // Copy MC particles
685
686 if(!fFillMCParticles) return;
687
688 TClonesArray* inMCParticles = (TClonesArray*) (fAODEvent ->FindListObject("mcparticles"));
689 TClonesArray* ouMCParticles = (TClonesArray*) ( AODEvent()->FindListObject("mcparticles"));
690
691 if( inMCParticles && ouMCParticles ) new (ouMCParticles) TClonesArray(*inMCParticles);
692
693}
694
695//_____________________________________________
696void AliAnalysisTaskCaloFilter::FillAODTracks()
697{
698 // AOD track copy
699
700 if(!fFillTracks) return;
701
702 AliAODTrack* aodTrack(0x0);
703
704 Double_t pos[3] = { 0. };
705 Double_t covTr[21]= { 0. };
706 //Double_t pid[10] = { 0. };
707 Double_t p[3] = { 0. };
708
709 // Copy from AODs
710 if(fAODEvent)
711 {
712 //TClonesArray* inTracks = fAODEvent ->GetTracks();
713 TClonesArray* ouTracks = AODEvent()->GetTracks();
714 //new (ouTracks) TClonesArray(*inTracks);
715
716 //printf("N tracks %d\n",fAODEvent->GetNumberOfTracks());
717 Int_t nCopyTrack = 0;
718 for (Int_t nTrack = 0; nTrack < fAODEvent->GetNumberOfTracks(); ++nTrack)
719 {
720 AliAODTrack *track = dynamic_cast<AliAODTrack*>(fAODEvent->GetTrack(nTrack));
721 if(!track) AliFatal("Not a standard AOD");
722
723 // Select only hybrid tracks?
724 if(fFillHybridTracks && !track->IsHybridGlobalConstrainedGlobal()) continue;
725
726 // Remove PID object to save space
727 //track->SetDetPID(0x0);
728
729 //new((*ouTracks)[nCopyTrack++]) AliAODTrack(*track);
730
731 track->GetPxPyPz(p);
732 Bool_t isDCA = track->GetPosition(pos);
733 track->GetCovMatrix(covTr);
734 //track->GetPID(pid);
735
736 AliAODVertex* primVertex = (AliAODVertex*) AODEvent()->GetVertices()->At(0); // primary vertex, copied previously!!!
737
738 aodTrack = new((*ouTracks)[nCopyTrack++]) AliAODTrack(
739 track->GetID(),
740 track->GetLabel(),
741 p,
742 kTRUE,
743 pos,
744 isDCA,
745 covTr,
746 track->Charge(),
747 track->GetITSClusterMap(),
748 // pid,
749 primVertex,
750 track->GetUsedForVtxFit(),
751 track->GetUsedForPrimVtxFit(),
752 (AliAODTrack::AODTrk_t) track->GetType(),
753 track->GetFilterMap());
754
755
756 aodTrack->SetPIDForTracking(track->GetPIDForTracking());
757 aodTrack->SetIsHybridGlobalConstrainedGlobal(track->IsHybridGlobalConstrainedGlobal());
758 aodTrack->SetIsHybridTPCConstrainedGlobal (track->IsHybridTPCConstrainedGlobal());
759 aodTrack->SetIsGlobalConstrained (track->IsGlobalConstrained());
760 aodTrack->SetIsTPCConstrained (track->IsTPCConstrained());
761
762 aodTrack->SetTPCFitMap (track->GetTPCFitMap());
763 aodTrack->SetTPCClusterMap(track->GetTPCClusterMap());
764 aodTrack->SetTPCSharedMap (track->GetTPCSharedMap());
765
766 aodTrack->SetChi2MatchTrigger(track->GetChi2MatchTrigger());
767
768 // set the DCA values to the AOD track
769
770 aodTrack->SetPxPyPzAtDCA(track->PxAtDCA(),track->PyAtDCA(),track->PzAtDCA());
771 aodTrack->SetXYAtDCA (track->XAtDCA() ,track->YAtDCA());
772
773 aodTrack->SetFlags (track->GetFlags());
774 aodTrack->SetTPCPointsF (track->GetTPCNclsF());
775
776 // Calo
777
778 if(track->IsEMCAL()) aodTrack->SetEMCALcluster(track->GetEMCALcluster());
779 if(track->IsPHOS()) aodTrack->SetPHOScluster (track->GetPHOScluster());
780 aodTrack->SetTrackPhiEtaPtOnEMCal( track->GetTrackPhiOnEMCal(), track->GetTrackPhiOnEMCal(), track->GetTrackPtOnEMCal() );
781
782 }
783
784 //printf("Final N tracks %d\n",nCopyTrack);
785
786 return;
787 }
788
789}
790
791//_________________________________________
792void AliAnalysisTaskCaloFilter::FillAODv0s()
793{
794 // Copy v0s (use if you know what you do, use quite a lot of memory)
795
796 if(!fFillv0s) return;
797
798 // Copy from AODs
799 if(fAODEvent)
800 {
801 TClonesArray* inv0 = fAODEvent ->GetV0s();
802 TClonesArray* ouv0 = AODEvent()->GetV0s();
803
804 //new (ouv0s) TClonesArray(*inv0s);
805
806 Int_t allv0s = inv0->GetEntriesFast();
807
808 for (Int_t nv0s = 0; nv0s < allv0s; ++nv0s)
809 {
810 AliAODv0 *v0 = (AliAODv0*)inv0->At(nv0s);
811
812 new((*ouv0)[nv0s]) AliAODv0(*v0);
813 }
814
815 return;
816 }
817
818}
819
820//____________________________________________
821void AliAnalysisTaskCaloFilter::FillAODVZERO()
822{
823 // Copy VZERO
824
825 if(!fFillVZERO) return;
826
827 AliAODVZERO* vzeroData = AODEvent()->GetVZEROData();
828
829 if(fESDEvent) *vzeroData = *(fESDEvent->GetVZEROData());
830 else *vzeroData = *(fAODEvent->GetVZEROData());
831
832}
833
834//_______________________________________________
835void AliAnalysisTaskCaloFilter::FillAODVertices()
836{
837 // Copy vertices
838
839 // set arrays and pointers
840 Double_t pos[3] ;
841 Double_t covVtx[6];
842 for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
843
844 // Copy from AODs
845 if(fAODEvent)
846 {
847 TClonesArray* inVertices = fAODEvent ->GetVertices();
848 TClonesArray* ouVertices = AODEvent()->GetVertices();
849
850 //new (ouVertices) TClonesArray(*inVertices);
851
852 //Keep only the first 3 vertices if not requested
853 Int_t allVertices = inVertices->GetEntriesFast();
854
855 //printf("n Vertices %d\n",allVertices);
856
857 if(!fFillAllVertices)
858 {
859 if(allVertices > 3) allVertices = 3;
860 }
861
862 //printf("Final n Vertices %d\n",allVertices);
863
864 for (Int_t nVertices = 0; nVertices < allVertices; ++nVertices)
865 {
866 AliAODVertex *vertex = (AliAODVertex*)inVertices->At(nVertices);
867
868 new((*ouVertices)[nVertices]) AliAODVertex(*vertex);
869 }
870
871 return;
872 }
873
874 if(!fESDEvent) return;
875
876 // Copy from ESDs
877
878 // Access to the AOD container of vertices
879 Int_t jVertices=0;
880 TClonesArray &vertices = *(AODEvent()->GetVertices());
881
882 // Add primary vertex. The primary tracks will be defined
883 // after the loops on the composite objects (v0, cascades, kinks)
884 fEvent ->GetPrimaryVertex()->GetXYZ(pos);
885 fESDEvent->GetPrimaryVertex()->GetCovMatrix(covVtx);
886 Float_t chi = fESDEvent->GetPrimaryVertex()->GetChi2toNDF();
887
888 AliAODVertex * primary = new(vertices[jVertices++])
889 AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
890 primary->SetName(fEvent->GetPrimaryVertex()->GetName());
891 primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
892
893}
894
895//____________________________________
896void AliAnalysisTaskCaloFilter::Init()
897{
898 //Init analysis with configuration macro if available
899
900 if(gROOT->LoadMacro(fConfigName) >=0)
901 {
902 AliInfo(Form("Configure analysis with %s",fConfigName.Data()));
903
904 AliAnalysisTaskCaloFilter *filter = (AliAnalysisTaskCaloFilter*)gInterpreter->ProcessLine("ConfigCaloFilter()");
905
906 fEMCALGeoName = filter->fEMCALGeoName;
907 fLoadEMCALMatrices = filter->fLoadEMCALMatrices;
908 fFillAODFile = filter->fFillAODFile;
909 fFillTracks = filter->fFillTracks;
910 fFillHybridTracks = filter->fFillHybridTracks;
911 fFillv0s = filter->fFillv0s;
912 fFillVZERO = filter->fFillVZERO;
913 fFillAllVertices = filter->fFillAllVertices;
914 fEMCALRecoUtils = filter->fEMCALRecoUtils;
915 fConfigName = filter->fConfigName;
916 fCaloFilter = filter->fCaloFilter;
917 fEventSelection[0] = filter->fEventSelection[0];
918 fEventSelection[1] = filter->fEventSelection[1];
919 fEventSelection[2] = filter->fEventSelection[2];
920 fAcceptAllMBEvent = filter->fAcceptAllMBEvent;
921 fMBTriggerMask = filter->fMBTriggerMask;
922 fCorrect = filter->fCorrect;
923 fEMCALEnergyCut = filter->fEMCALEnergyCut;
924 fEMCALNcellsCut = filter->fEMCALNcellsCut;
925 fPHOSEnergyCut = filter->fPHOSEnergyCut;
926 fPHOSNcellsCut = filter->fPHOSNcellsCut;
927 fTrackPtCut = filter->fTrackPtCut;
928 fVzCut = filter->fVzCut;
929 fCheckEventVertex = filter->fCheckEventVertex;
930
931 for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = filter->fEMCALMatrix[i] ;
932 }
933}
934
935//_________________________________________
936void AliAnalysisTaskCaloFilter::PrintInfo()
937{
938 //Print settings
939
940 printf("AnalysisCaloFilter::PrintInfo() \n");
941
942 printf("\t Not only filter, correct Clusters? %d\n",fCorrect);
943 printf("\t Calorimeter Filtering Option ? %d\n",fCaloFilter);
944
945 //printf("\t Use handmade geo matrices? EMCAL %d, PHOS %d\n",fLoadEMCALMatrices, fLoadPHOSMatrices);
946 printf("\t Use handmade geo matrices? EMCAL %d, PHOS 0\n",fLoadEMCALMatrices);
947
948 printf("\t Fill: AOD file? %d Tracks? %d; all Vertex? %d; v0s? %d; VZERO ? %d\n",
949 fFillAODFile,fFillTracks,fFillAllVertices, fFillv0s, fFillVZERO);
950
951 printf("\t Event Selection based : EMCAL? %d, PHOS? %d Tracks? %d - Accept all MB with mask %d? %d\n",
952 fEventSelection[0],fEventSelection[1],fEventSelection[2],fMBTriggerMask, fAcceptAllMBEvent);
953
954 printf("\t \t EMCAL E > %2.2f, EMCAL nCells >= %d, PHOS E > %2.2f, PHOS nCells >= %d, Track pT > %2.2f, |vz| < %2.2f\n",
955 fEMCALEnergyCut,fEMCALNcellsCut,fPHOSEnergyCut,fPHOSNcellsCut, fTrackPtCut,fVzCut);
956}
957
958//_______________________________________________________
959void AliAnalysisTaskCaloFilter::UserCreateOutputObjects()
960{
961 // Init geometry
962
963 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
964
965 if(fFillMCParticles)
966 {
967 TClonesArray * aodMCParticles = new TClonesArray("AliAODMCParticle",500);
968 aodMCParticles->SetName("mcparticles");
969 ((AliAODHandler*)AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler())->AddBranch("TClonesArray", &aodMCParticles);
970 }
971
972}
973
974//____________________________________________________________
975void AliAnalysisTaskCaloFilter::UserExec(Option_t */*option*/)
976{
977 // Execute analysis for current event
978 // Copy input ESD or AOD header, vertex, CaloClusters and CaloCells to output AOD
979
980 AliDebug(1,Form("Analysing event # %d", (Int_t)Entry()));
981
982 fEvent = InputEvent();
983 fAODEvent = dynamic_cast<AliAODEvent*> (fEvent);
984 fESDEvent = dynamic_cast<AliESDEvent*> (fEvent);
985
986 if(!fEvent)
987 {
988 AliInfo("This event does not contain Input?");
989 return;
990 }
991
992 // printf("Start processing : %s\n",fAODEvent->GetFiredTriggerClasses().Data());
993
994 // Select the event
995
996 if(!AcceptEvent()) return ;
997
998 //Magic line to write events to file
999
1000 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
1001
1002 // Reset output AOD
1003
1004 Int_t nVertices = 0;
1005 if(fFillv0s) nVertices = fEvent->GetNumberOfV0s();
1006 Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
1007 Int_t nTracks = fEvent->GetNumberOfTracks();
1008
1009 AODEvent()->ResetStd(nTracks, nVertices, 0, 0, 0, nCaloClus, 0, 0);
1010
1011 // Copy
1012
1013 FillAODHeader();
1014
1015 //
1016 FillAODv0s();
1017
1018 //
1019 FillAODVertices(); // Do it before the track filtering to have the reference to the vertex
1020
1021 //
1022 FillAODVZERO();
1023
1024 //
1025 FillAODTracks();
1026
1027 //
1028 CorrectionsInEMCAL();
1029
1030 //
1031 FillAODCaloClusters();
1032
1033 //
1034 FillAODCaloCells();
1035
1036 //
1037 FillAODCaloTrigger();
1038
1039 //
1040 FillAODMCParticles();
1041
1042 //printf("Filtered event, end processing : %s\n",fAODEvent->GetFiredTriggerClasses().Data());
1043
1044}
1045