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