AliPhysics  1811c8f (1811c8f)
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  for (Int_t iCell = 0; iCell < nEMcell; iCell++)
475  {
476  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
477  fEMCALGeo->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
478  fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
479 
481  {
482  calibFactor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
483  }
484 
485  if(!fEMCALRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
486  { //Channel is not declared as bad
487  aodEMcells.SetCell(iCell,
488  eventEMcells.GetCellNumber(iCell),
489  eventEMcells.GetAmplitude(iCell)*calibFactor,
490  eventEMcells.GetTime(iCell),
491  eventEMcells.GetMCLabel(iCell),
492  eventEMcells.GetEFraction(iCell),
493  eventEMcells.GetHighGain(iCell));
494  //printf("GOOD channel\n");
495  }
496  else
497  {
498  aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
499  //printf("BAD channel\n");
500  }
501  }
502  aodEMcells.Sort();
503  }
504 
505  // PHOS
506  if ((fCaloFilter==kBoth || fCaloFilter==kPHOS) && fEvent->GetPHOSCells())
507  { // protection against missing ESD information
508  AliVCaloCells &eventPHcells = *(fEvent->GetPHOSCells());
509  Int_t nPHcell = eventPHcells.GetNumberOfCells() ;
510 
511  AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
512  aodPHcells.CreateContainer(nPHcell);
513  aodPHcells.SetType(AliVCaloCells::kPHOSCell);
514 
515  for (Int_t iCell = 0; iCell < nPHcell; iCell++)
516  {
517  aodPHcells.SetCell(iCell,
518  eventPHcells.GetCellNumber(iCell),
519  eventPHcells.GetAmplitude(iCell),
520  eventPHcells.GetTime(iCell),
521  eventPHcells.GetMCLabel(iCell),
522  eventPHcells.GetEFraction(iCell),
523  eventPHcells.GetHighGain(iCell));
524  }
525 
526  aodPHcells.Sort();
527  }
528 }
529 
530 //___________________________________________________
533 //___________________________________________________
535 {
536  TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
537  Int_t jClusters=0;
538  Float_t posF[3] ;
539 
540  Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
541  for (Int_t iClust=0; iClust<nCaloClus; ++iClust)
542  {
543  AliVCluster * cluster = fEvent->GetCaloCluster(iClust);
544 
545  // Check which calorimeter information we want to keep.
546 
547  if(fCaloFilter!=kBoth)
548  {
549  if (fCaloFilter==kPHOS && cluster->IsEMCAL()) continue ;
550  else if(fCaloFilter==kEMCAL && cluster->IsPHOS()) continue ;
551  }
552 
553  // Get original residuals, in case of previous recalculation, reset them
554  Float_t dR = cluster->GetTrackDx();
555  Float_t dZ = cluster->GetTrackDz();
556 
557  AliDebug(2,Form("Original residuals : dZ %f, dR %f",dZ, dR));
558 
559  //--------------------------------------------------------------
560  // If EMCAL and corrections done, get the new matching parameters, do not copy noisy clusters
561  if(cluster->IsEMCAL() && fCorrect)
562  {
563  AliDebug(2,Form("Check cluster %d for bad channels and close to border",cluster->GetID()));
564 
565  if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
566 
567  if(fEMCALRecoUtils->IsExoticCluster(cluster, InputEvent()->GetEMCALCells(),InputEvent()->GetBunchCrossNumber())) continue;
568 
569  fEMCALRecoUtils->GetMatchedResiduals(cluster->GetID(),dR,dZ);
570  cluster->SetTrackDistance(dR,dZ);
571  }
572 
573  AliDebug(2,Form("EMCAL? %d, PHOS? %d Track-Cluster Residuals : dZ %f, dR %f",cluster->IsEMCAL(), cluster->IsPHOS(),dZ, dR));
574 
575  //--------------------------------------------------------------
576 
577  // Now fill AODs
578 
579  Int_t id = cluster->GetID();
580  Float_t energy = cluster->E();
581  cluster->GetPosition(posF);
582 
583  AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
584  AliAODCaloCluster(id,
585  cluster->GetNLabels(),
586  cluster->GetLabels(),
587  energy,
588  posF,
589  NULL,
590  cluster->GetType());
591 
592  caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
593  cluster->GetDispersion(),
594  cluster->GetM20(), cluster->GetM02(),
595  -1,
596  cluster->GetNExMax(),cluster->GetTOF()) ;
597 
598  caloCluster->SetPIDFromESD(cluster->GetPID());
599  caloCluster->SetNCells(cluster->GetNCells());
600  caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
601  caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
602  caloCluster->SetTrackDistance(dR, dZ);
603 
604  caloCluster->SetCellsMCEdepFractionMap(cluster->GetCellsMCEdepFractionMap());
605  caloCluster->SetClusterMCEdepFraction (cluster->GetClusterMCEdepFraction ());
606 
607  AliDebug(2,Form("Filter, aod : i %d, E %f, dispersion %f, m02 %f, m20 %f",caloCluster->GetID(),caloCluster->E(),
608  caloCluster->GetDispersion(),caloCluster->GetM02(),caloCluster->GetM20()));
609  caloCluster->GetPosition(posF);
610  AliDebug(2,Form("Filter, aod : i %d, x %f, y %f, z %f",caloCluster->GetID(), posF[0], posF[1], posF[2]));
611 
612  //Matched tracks, just to know if there was any match, the track pointer is useless if tracks not stored
613  if(TMath::Abs(dR) < 990 && TMath::Abs(dZ) < 990)
614  { //Default value in PHOS 999, in EMCAL 1024, why?
615  caloCluster->AddTrackMatched(new AliAODTrack);
616  }
617  // TO DO, in case Tracks available, think how to put the matched track in AOD
618  }
619 
620  caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots
621 }
622 
623 //__________________________________________________
625 //__________________________________________________
627 {
628  if( !AODEvent() || !fAODEvent ) return;
629 
630  AliAODCaloTrigger* triggerEM = AODEvent()->GetCaloTrigger("EMCAL");
631  AliAODCaloTrigger* triggerPH = AODEvent()->GetCaloTrigger("PHOS");
632 
633  // Copy from AODs
634 
635  AliAODCaloTrigger* inTriggerEM = fAODEvent ->GetCaloTrigger("EMCAL");
636  AliAODCaloTrigger* inTriggerPH = fAODEvent ->GetCaloTrigger("PHOS");
637 
638  if(inTriggerPH && (fCaloFilter==kBoth || fCaloFilter==kPHOS)) *triggerPH = *inTriggerPH;
639 
640  if(inTriggerEM && (fCaloFilter==kBoth || fCaloFilter==kEMCAL)) *triggerEM = *inTriggerEM;
641 }
642 
643 //______________________________________________
645 //______________________________________________
647 {
648  AliAODHeader* header = dynamic_cast<AliAODHeader*>(AODEvent()->GetHeader());
649  if(!header)
650  {
651  AliFatal("Not a standard AOD");
652  return; // not needed but coverity complains
653  }
654 
655  // Copy from AODs
656  if(fAODEvent)
657  {
658  *header = *((AliAODHeader*)fAODEvent->GetHeader());
659  return;
660  }
661 
662  if(!fESDEvent) return;
663 
664  // Copy from ESDs
665 
666  header->SetRunNumber(fEvent->GetRunNumber());
667 
668  TTree* tree = fInputHandler->GetTree();
669  if (tree)
670  {
671  TFile* file = tree->GetCurrentFile();
672  if (file) header->SetESDFileName(file->GetName());
673  }
674 
675  header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
676  header->SetOrbitNumber(fEvent->GetOrbitNumber());
677  header->SetPeriodNumber(fEvent->GetPeriodNumber());
678  header->SetEventType(fEvent->GetEventType());
679 
680  //Centrality
681  if(fEvent->GetCentrality())
682  {
683  header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
684  }
685  else
686  {
687  header->SetCentrality(0);
688  }
689 
690  //Trigger
691  header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
692  header->SetFiredTriggerClasses(fESDEvent->GetFiredTriggerClasses());
693  header->SetTriggerMask(fEvent->GetTriggerMask());
694  header->SetTriggerCluster(fEvent->GetTriggerCluster());
695  header->SetL0TriggerInputs(fESDEvent->GetHeader()->GetL0TriggerInputs());
696  header->SetL1TriggerInputs(fESDEvent->GetHeader()->GetL1TriggerInputs());
697  header->SetL2TriggerInputs(fESDEvent->GetHeader()->GetL2TriggerInputs());
698 
699  header->SetMagneticField(fEvent->GetMagneticField());
700  header->SetMuonMagFieldScale(fESDEvent->GetCurrentDip()/6000.);
701 
702  header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
703  header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
704  header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
705  header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
706  header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
707 
708  Float_t diamxy[2]={(Float_t)fEvent->GetDiamondX(),(Float_t)fEvent->GetDiamondY()};
709  Float_t diamcov[3];
710  fEvent->GetDiamondCovXY(diamcov);
711  header->SetDiamond(diamxy,diamcov);
712  header->SetDiamondZ(fESDEvent->GetDiamondZ(),fESDEvent->GetSigma2DiamondZ());
713 }
714 
715 //__________________________________________________
717 //__________________________________________________
719 {
720  if(!fFillMCParticles) return;
721 
722  TClonesArray* inMCParticles = (TClonesArray*) (fAODEvent ->FindListObject("mcparticles"));
723  TClonesArray* ouMCParticles = (TClonesArray*) ( AODEvent()->FindListObject("mcparticles"));
724 
725  if( inMCParticles && ouMCParticles ) new (ouMCParticles) TClonesArray(*inMCParticles);
726 }
727 
728 //_____________________________________________
730 //_____________________________________________
732 {
733  if(!fFillTracks) return;
734 
735  AliAODTrack* aodTrack(0x0);
736 
737  Double_t pos[3] = { 0. };
738  Double_t covTr[21]= { 0. };
739  //Double_t pid[10] = { 0. };
740  Double_t p[3] = { 0. };
741 
742  // Copy from AODs
743  if(fAODEvent)
744  {
745  //TClonesArray* inTracks = fAODEvent ->GetTracks();
746  TClonesArray* ouTracks = AODEvent()->GetTracks();
747  //new (ouTracks) TClonesArray(*inTracks);
748 
749  //printf("N tracks %d\n",fAODEvent->GetNumberOfTracks());
750  Int_t nCopyTrack = 0;
751  for (Int_t nTrack = 0; nTrack < fAODEvent->GetNumberOfTracks(); ++nTrack)
752  {
753  AliAODTrack *track = dynamic_cast<AliAODTrack*>(fAODEvent->GetTrack(nTrack));
754  if(!track) AliFatal("Not a standard AOD");
755 
756  // Select only hybrid tracks?
757  if(fFillHybridTracks && !track->IsHybridGlobalConstrainedGlobal()) continue;
758 
759  // Remove PID object to save space
760  //track->SetDetPID(0x0);
761 
762  //new((*ouTracks)[nCopyTrack++]) AliAODTrack(*track);
763 
764  track->GetPxPyPz(p);
765  Bool_t isDCA = track->GetPosition(pos);
766  track->GetCovMatrix(covTr);
767  //track->GetPID(pid);
768 
769  AliAODVertex* primVertex = (AliAODVertex*) AODEvent()->GetVertices()->At(0); // primary vertex, copied previously!!!
770 
771  aodTrack = new((*ouTracks)[nCopyTrack++]) AliAODTrack(
772  track->GetID(),
773  track->GetLabel(),
774  p,
775  kTRUE,
776  pos,
777  isDCA,
778  covTr,
779  track->Charge(),
780  track->GetITSClusterMap(),
781  // pid,
782  primVertex,
783  track->GetUsedForVtxFit(),
784  track->GetUsedForPrimVtxFit(),
785  (AliAODTrack::AODTrk_t) track->GetType(),
786  track->GetFilterMap());
787 
788 
789  aodTrack->SetPIDForTracking(track->GetPIDForTracking());
790  aodTrack->SetIsHybridGlobalConstrainedGlobal(track->IsHybridGlobalConstrainedGlobal());
791  aodTrack->SetIsHybridTPCConstrainedGlobal (track->IsHybridTPCConstrainedGlobal());
792  aodTrack->SetIsGlobalConstrained (track->IsGlobalConstrained());
793  aodTrack->SetIsTPCConstrained (track->IsTPCConstrained());
794 
795  aodTrack->SetTPCFitMap (track->GetTPCFitMap());
796  aodTrack->SetTPCClusterMap(track->GetTPCClusterMap());
797  aodTrack->SetTPCSharedMap (track->GetTPCSharedMap());
798 
799  aodTrack->SetChi2MatchTrigger(track->GetChi2MatchTrigger());
800 
801  // set the DCA values to the AOD track
802 
803  aodTrack->SetPxPyPzAtDCA(track->PxAtDCA(),track->PyAtDCA(),track->PzAtDCA());
804  aodTrack->SetXYAtDCA (track->XAtDCA() ,track->YAtDCA());
805 
806  aodTrack->SetFlags (track->GetFlags());
807  aodTrack->SetTPCPointsF (track->GetTPCNclsF());
808 
809  // Calo
810 
811  if(track->IsEMCAL()) aodTrack->SetEMCALcluster(track->GetEMCALcluster());
812  if(track->IsPHOS()) aodTrack->SetPHOScluster (track->GetPHOScluster());
813  aodTrack->SetTrackPhiEtaPtOnEMCal( track->GetTrackPhiOnEMCal(), track->GetTrackPhiOnEMCal(), track->GetTrackPtOnEMCal() );
814 
815  }
816 
817  //printf("Final N tracks %d\n",nCopyTrack);
818 
819  return;
820  }
821 }
822 
823 //_________________________________________
826 //_________________________________________
828 {
829  if(!fFillv0s) return;
830 
831  // Copy from AODs
832  if(fAODEvent)
833  {
834  TClonesArray* inv0 = fAODEvent ->GetV0s();
835  TClonesArray* ouv0 = AODEvent()->GetV0s();
836 
837  //new (ouv0s) TClonesArray(*inv0s);
838 
839  Int_t allv0s = inv0->GetEntriesFast();
840 
841  for (Int_t nv0s = 0; nv0s < allv0s; ++nv0s)
842  {
843  AliAODv0 *v0 = (AliAODv0*)inv0->At(nv0s);
844 
845  new((*ouv0)[nv0s]) AliAODv0(*v0);
846  }
847 
848  return;
849  }
850 }
851 
852 //____________________________________________
854 //____________________________________________
856 {
857  if(!fFillVZERO) return;
858 
859  AliAODVZERO* vzeroData = AODEvent()->GetVZEROData();
860 
861  if(fESDEvent) *vzeroData = *(fESDEvent->GetVZEROData());
862  else *vzeroData = *(fAODEvent->GetVZEROData());
863 }
864 
865 //_______________________________________________
867 //_______________________________________________
869 {
870  // set arrays and pointers
871  Double_t pos[3] ;
872  Double_t covVtx[6];
873  for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
874 
875  // Copy from AODs
876  if(fAODEvent)
877  {
878  TClonesArray* inVertices = fAODEvent ->GetVertices();
879  TClonesArray* ouVertices = AODEvent()->GetVertices();
880 
881  //new (ouVertices) TClonesArray(*inVertices);
882 
883  //Keep only the first 3 vertices if not requested
884  Int_t allVertices = inVertices->GetEntriesFast();
885 
886  //printf("n Vertices %d\n",allVertices);
887 
888  if(!fFillAllVertices)
889  {
890  if(allVertices > 3) allVertices = 3;
891  }
892 
893  //printf("Final n Vertices %d\n",allVertices);
894 
895  for (Int_t nVertices = 0; nVertices < allVertices; ++nVertices)
896  {
897  AliAODVertex *vertex = (AliAODVertex*)inVertices->At(nVertices);
898 
899  new((*ouVertices)[nVertices]) AliAODVertex(*vertex);
900  }
901 
902  return;
903  }
904 
905  if(!fESDEvent) return;
906 
907  // Copy from ESDs
908 
909  // Access to the AOD container of vertices
910  Int_t jVertices=0;
911  TClonesArray &vertices = *(AODEvent()->GetVertices());
912 
913  // Add primary vertex. The primary tracks will be defined
914  // after the loops on the composite objects (v0, cascades, kinks)
915  fEvent ->GetPrimaryVertex()->GetXYZ(pos);
916  fESDEvent->GetPrimaryVertex()->GetCovMatrix(covVtx);
917  Float_t chi = fESDEvent->GetPrimaryVertex()->GetChi2toNDF();
918 
919  AliAODVertex * primary = new(vertices[jVertices++])
920  AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
921  primary->SetName(fEvent->GetPrimaryVertex()->GetName());
922  primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
923 }
924 
925 //____________________________________
930 //____________________________________
932 {
933  if(gROOT->LoadMacro(fConfigName) >=0)
934  {
935  AliInfo(Form("Configure analysis with %s",fConfigName.Data()));
936 
937  AliAnalysisTaskCaloFilter *filter = (AliAnalysisTaskCaloFilter*)gInterpreter->ProcessLine("ConfigCaloFilter()");
938 
939  fEMCALGeoName = filter->fEMCALGeoName;
941  fFillAODFile = filter->fFillAODFile;
942  fFillTracks = filter->fFillTracks;
944  fFillv0s = filter->fFillv0s;
945  fFillVZERO = filter->fFillVZERO;
947  fEMCALRecoUtils = filter->fEMCALRecoUtils;
948  fConfigName = filter->fConfigName;
949  fCaloFilter = filter->fCaloFilter;
950  fEventSelection[0] = filter->fEventSelection[0];
951  fEventSelection[1] = filter->fEventSelection[1];
952  fEventSelection[2] = filter->fEventSelection[2];
954  fMBTriggerMask = filter->fMBTriggerMask;
955  fCorrect = filter->fCorrect;
958  fPHOSEnergyCut = filter->fPHOSEnergyCut;
959  fPHOSNcellsCut = filter->fPHOSNcellsCut;
960  fTrackPtCut = filter->fTrackPtCut;
961  fVzCut = filter->fVzCut;
963 
964  for(Int_t i = 0; i < 22; i++) fEMCALMatrix[i] = filter->fEMCALMatrix[i] ;
965  }
966 }
967 
968 //_________________________________________
970 //_________________________________________
972 {
973  printf("AnalysisCaloFilter::PrintInfo() \n");
974 
975  printf("\t Not only filter, correct Clusters? %d\n",fCorrect);
976  printf("\t Calorimeter Filtering Option ? %d\n",fCaloFilter);
977 
978  //printf("\t Use handmade geo matrices? EMCAL %d, PHOS %d\n",fLoadEMCALMatrices, fLoadPHOSMatrices);
979  printf("\t Use handmade geo matrices? EMCAL %d, PHOS 0\n",fLoadEMCALMatrices);
980 
981  printf("\t Fill: AOD file? %d Tracks? %d; all Vertex? %d; v0s? %d; VZERO ? %d\n",
983 
984  printf("\t Event Selection based : EMCAL? %d, PHOS? %d Tracks? %d - Accept all MB with mask %d? %d\n",
986 
987  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",
989 }
990 
991 //_______________________________________________________
993 //_______________________________________________________
995 {
996  fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
997 
998  if(fFillMCParticles)
999  {
1000  TClonesArray * aodMCParticles = new TClonesArray("AliAODMCParticle",500);
1001  aodMCParticles->SetName("mcparticles");
1002  ((AliAODHandler*)AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler())->AddBranch("TClonesArray", &aodMCParticles);
1003  }
1004 }
1005 
1006 //____________________________________________________________
1010 //____________________________________________________________
1012 {
1013  AliDebug(1,Form("Analysing event # %d", (Int_t)Entry()));
1014 
1015  fEvent = InputEvent();
1016  fAODEvent = dynamic_cast<AliAODEvent*> (fEvent);
1017  fESDEvent = dynamic_cast<AliESDEvent*> (fEvent);
1018 
1019  if(!fEvent)
1020  {
1021  AliInfo("This event does not contain Input?");
1022  return;
1023  }
1024 
1025  // printf("Start processing : %s\n",fAODEvent->GetFiredTriggerClasses().Data());
1026 
1027  // Select the event
1028 
1029  if(!AcceptEvent()) return ;
1030 
1031  //Magic line to write events to file
1032 
1033  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
1034 
1035  // Reset output AOD
1036 
1037  Int_t nVertices = 0;
1038  if(fFillv0s) nVertices = fEvent->GetNumberOfV0s();
1039  Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
1040  Int_t nTracks = fEvent->GetNumberOfTracks();
1041 
1042  AODEvent()->ResetStd(nTracks, nVertices, 0, 0, 0, nCaloClus, 0, 0);
1043 
1044  // Copy
1045 
1046  FillAODHeader();
1047 
1048  //
1049  FillAODv0s();
1050 
1051  //
1052  FillAODVertices(); // Do it before the track filtering to have the reference to the vertex
1053 
1054  //
1055  FillAODVZERO();
1056 
1057  //
1058  FillAODTracks();
1059 
1060  //
1062 
1063  //
1065 
1066  //
1067  FillAODCaloCells();
1068 
1069  //
1071 
1072  //
1074 
1075  //printf("Filtered event, end processing : %s\n",fAODEvent->GetFiredTriggerClasses().Data());
1076 }
1077 
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:44
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.
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.
Int_t GetEMCALChannelStatus(Int_t iSM, Int_t iCol, Int_t iRow) const
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.