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