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