Fix Coverity
[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / 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 /* $Id: AliAnalysisTaskCaloFilter.cxx $ */
17
18 //////////////////////////////////////////////////////////
19 // Filter the ESDCaloClusters and ESDCaloCells of EMCAL,
20 // PHOS or both, creating the corresponing AODCaloClusters
21 // and AODCaloCells.
22 // Keep also the AODHeader information and the vertex.
23 // Needed for calorimeter calibration.
24 // Copy of AliAnalysisTaskESDfilter.
25 // Author: Gustavo Conesa Balbastre (INFN - Frascati)
26 //////////////////////////////////////////////////////////
27
28 //Root
29 #include "TGeoManager.h"
30 #include "TFile.h"
31 #include "TNtuple.h"
32 #include "TROOT.h"
33 #include "TInterpreter.h"
34
35 //STEER
36 #include "AliESDEvent.h"
37 #include "AliAODEvent.h"
38 #include "AliLog.h"
39 #include "AliVCluster.h"
40 #include "AliVCaloCells.h"
41 #include "AliVEventHandler.h"
42 #include "AliAnalysisManager.h"
43 #include "AliInputEventHandler.h"
44 #include "AliESDtrackCuts.h"
45 #include "AliTriggerAnalysis.h"
46
47 //EMCAL
48 #include "AliEMCALRecoUtils.h"
49 #include "AliEMCALGeometry.h"
50
51 #include "AliAnalysisTaskCaloFilter.h"
52
53 ClassImp(AliAnalysisTaskCaloFilter)
54   
55 ////////////////////////////////////////////////////////////////////////
56
57 AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter():
58   AliAnalysisTaskSE("CaloFilterTask"), //fCuts(0x0),
59   fCaloFilter(0), fCorrect(kFALSE), 
60   fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEARV1"), 
61   fEMCALRecoUtils(new AliEMCALRecoUtils),
62   fESDtrackCuts(0), fTriggerAnalysis (new AliTriggerAnalysis), fTrackMultEtaCut(0.8),
63   fLoadEMCALMatrices(kFALSE), //fLoadPHOSMatrices(kFALSE),
64   fGeoMatrixSet(kFALSE), fEventNtuple(0),fConfigName(""),fFillAODFile(kTRUE)
65 {
66   // Default constructor
67   fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
68   for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
69   //for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix[i]  = 0 ;
70
71   DefineOutput(1, TNtuple::Class());
72   
73 }
74
75 //__________________________________________________
76 AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter(const char* name):
77   AliAnalysisTaskSE(name), //fCuts(0x0),
78   fCaloFilter(0), fCorrect(kFALSE),
79   fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEARV1"), 
80   fEMCALRecoUtils(new AliEMCALRecoUtils),
81   fESDtrackCuts(0), fTriggerAnalysis (new AliTriggerAnalysis), fTrackMultEtaCut(0.8),
82   fLoadEMCALMatrices(kFALSE), //fLoadPHOSMatrices(kFALSE),
83   fGeoMatrixSet(kFALSE), fEventNtuple(0),fConfigName(""),fFillAODFile(kTRUE)
84 {
85   // Constructor
86   fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
87   for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
88   //for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix[i]  = 0 ;
89
90   DefineOutput(1, TNtuple::Class());
91
92 }
93
94 //__________________________________________________
95 AliAnalysisTaskCaloFilter::~AliAnalysisTaskCaloFilter()
96 {
97   //Destructor.
98         
99   if(fEMCALGeo)        delete fEMCALGeo;        
100   if(fEMCALRecoUtils)  delete fEMCALRecoUtils;
101   if(fESDtrackCuts)    delete fESDtrackCuts;
102   if(fTriggerAnalysis) delete fTriggerAnalysis;
103
104   if(fEventNtuple)     delete fEventNtuple;
105   
106 }
107
108 //-------------------------------------------------------------------
109 void AliAnalysisTaskCaloFilter::Init()
110 {
111   //Init analysis with configuration macro if available
112   
113   if(gROOT->LoadMacro(fConfigName) >=0){
114     printf("Configure analysis with %s\n",fConfigName.Data());
115     AliAnalysisTaskCaloFilter *filter = (AliAnalysisTaskCaloFilter*)gInterpreter->ProcessLine("ConfigCaloFilter()");
116     fEMCALGeoName      = filter->fEMCALGeoName; 
117     fLoadEMCALMatrices = filter->fLoadEMCALMatrices;
118     fFillAODFile       = filter->fFillAODFile;
119     fEMCALRecoUtils    = filter->fEMCALRecoUtils; 
120     fConfigName        = filter->fConfigName;
121     fCaloFilter        = filter->fCaloFilter;
122     fCorrect           = filter->fCorrect;
123     fTrackMultEtaCut   = filter->fTrackMultEtaCut;
124     fESDtrackCuts      = filter->fESDtrackCuts;
125     for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = filter->fEMCALMatrix[i] ;
126   }
127
128
129 //__________________________________________________
130 void AliAnalysisTaskCaloFilter::UserCreateOutputObjects()
131 {
132   // Init geometry 
133         
134   fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;   
135   
136   OpenFile(1);
137   
138   fEventNtuple = new TNtuple("EventSelection","EventSelection", "bPileup:bGoodVertex:bV0AND:trackMult");
139  
140   PostData(1, fEventNtuple);
141
142 }  
143
144 //__________________________________________________
145 void AliAnalysisTaskCaloFilter::UserExec(Option_t */*option*/)
146 {
147   // Execute analysis for current event
148   // Copy input ESD or AOD header, vertex, CaloClusters and CaloCells to output AOD
149
150   if (fDebug > 0)  
151     printf("CaloFilter: Analysing event # %d\n", (Int_t)Entry());
152   
153   AliVEvent* event = InputEvent();
154   if(!event) {
155     printf("AliAnalysisTaskCaloFilter::UserExec - This event does not contain Input?");
156     return;
157   }
158
159   //Magic line to write events to file
160   AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
161    
162   // cast event, depending on input we will have to use one or the other type of event
163   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);  
164   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
165   
166   //-------------------------------------------
167   //Event selection parameters
168   //-------------------------------------------
169   //Is it a pileup event?
170   Bool_t bPileup = event->IsPileupFromSPD(3, 0.8, 3., 2., 5.); //Default values, if not it does not compile
171   //Bool_t bPileup = event->IsPileupFromSPD(); 
172   //if(bPileup) return kFALSE;
173   Int_t  trackMult    = 0;
174   Bool_t bV0AND       = kFALSE;
175   Bool_t bGoodVertex  = kFALSE;
176   if(esdevent){
177     //Get track multiplicity
178     Int_t nTracks   = InputEvent()->GetNumberOfTracks() ;
179     for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
180       AliVTrack * track = (AliVTrack*)InputEvent()->GetTrack(itrack) ; // retrieve track from esd
181       if(!fESDtrackCuts->AcceptTrack((AliESDtrack*)track)) continue;
182       //Count the tracks in eta < 0.8
183       if(TMath::Abs(track->Eta())< fTrackMultEtaCut) trackMult++;
184     }  
185     //V0AND?   
186     bV0AND = fTriggerAnalysis->IsOfflineTriggerFired(esdevent, AliTriggerAnalysis::kV0AND);
187     //if(!bV0AND) return kFALSE;
188     //Well reconstructed vertex
189     bGoodVertex = CheckForPrimaryVertex();
190     //if(!bGoodVertex) return kFALSE;
191     
192   }//ESDs
193   
194   if(fDebug > 0)
195     printf("AliAnalysisTaskCaloFilter::UserExec() - PileUp %d, Good Vertex %d, v0AND %d, Track Mult in |eta| < %2.1f = %d\n",
196            bPileup,bGoodVertex,bV0AND, fTrackMultEtaCut, trackMult);
197   
198   //Put bools with event selection parameters in a TNtuple
199   //Int_t eventSelection[] = {bPileup,bGoodVertex,bV0AND,trackMult};
200   fEventNtuple->Fill(bPileup,bGoodVertex,bV0AND,trackMult);
201   
202   //--------------------------------------------------------------------
203   //Set in AOD General Event Parameters, vertex, runnumber, trigger, etc 
204   //-------------------------------------------------------------------
205   
206   // set arrays and pointers
207   Float_t  posF[3]  ;
208   Double_t pos[3]   ;
209   Double_t covVtx[6];
210   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
211       
212   AliAODHeader* header = AODEvent()->GetHeader();
213   
214   header->SetRunNumber(event->GetRunNumber());
215   
216   if(esdevent){
217     TTree* tree = fInputHandler->GetTree();
218     if (tree) {
219       TFile* file = tree->GetCurrentFile();
220       if (file) header->SetESDFileName(file->GetName());
221     }
222   }
223   else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
224   
225   header->SetBunchCrossNumber(event->GetBunchCrossNumber());
226   header->SetOrbitNumber(event->GetOrbitNumber());
227   header->SetPeriodNumber(event->GetPeriodNumber());
228   header->SetEventType(event->GetEventType());
229   
230   //Centrality
231   if(event->GetCentrality()){
232     header->SetCentrality(new AliCentrality(*(event->GetCentrality())));
233   }
234   else{
235     header->SetCentrality(0);
236   }
237   
238   //Trigger  
239   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
240   if      (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
241   else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
242   header->SetTriggerMask(event->GetTriggerMask()); 
243   header->SetTriggerCluster(event->GetTriggerCluster());
244   if(esdevent){
245     header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());    
246     header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());    
247     header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());    
248   }
249   else if (aodevent){
250     header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());    
251     header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());    
252     header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());    
253   }
254
255   header->SetMagneticField(event->GetMagneticField());
256   //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.); 
257
258   header->SetZDCN1Energy(event->GetZDCN1Energy());
259   header->SetZDCP1Energy(event->GetZDCP1Energy());
260   header->SetZDCN2Energy(event->GetZDCN2Energy());
261   header->SetZDCP2Energy(event->GetZDCP2Energy());
262   header->SetZDCEMEnergy(event->GetZDCEMEnergy(0),event->GetZDCEMEnergy(1));
263   
264   Float_t diamxy[2]={event->GetDiamondX(),event->GetDiamondY()};
265   Float_t diamcov[3];
266   event->GetDiamondCovXY(diamcov);
267   header->SetDiamond(diamxy,diamcov);
268   if      (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
269   else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
270   
271   //
272   //
273   Int_t nVertices = 1 ;/* = prim. vtx*/;
274   Int_t nCaloClus = event->GetNumberOfCaloClusters();
275   
276   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
277   
278   // Access to the AOD container of vertices
279   TClonesArray &vertices = *(AODEvent()->GetVertices());
280   Int_t jVertices=0;
281   
282   // Add primary vertex. The primary tracks will be defined
283   // after the loops on the composite objects (V0, cascades, kinks)
284   event->GetPrimaryVertex()->GetXYZ(pos);
285   Float_t chi = 0;
286   if      (esdevent){
287     esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
288     chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
289   }
290   else if (aodevent){
291     aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
292     chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
293   }
294   
295   AliAODVertex * primary = new(vertices[jVertices++])
296     AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
297   primary->SetName(event->GetPrimaryVertex()->GetName());
298   primary->SetTitle(event->GetPrimaryVertex()->GetTitle());
299   
300   // Access to the AOD container of clusters
301   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
302   Int_t jClusters=0;
303   
304   //-------------------------------------------
305   //Do Corrections in EMCAL 
306   //-------------------------------------------
307   //If EMCAL, and requested, correct energy, position ...
308   //Need to do this in a separate loop before filling the ESDs because of the track matching recalculations
309   if(fCorrect && (fCaloFilter==kEMCAL || fCaloFilter==kBoth) ) {
310     
311     if(!fGeoMatrixSet){
312       if(fLoadEMCALMatrices){
313         printf("AliAnalysisTaskCaloFilter::UserExec() - Load user defined EMCAL geometry matrices\n");
314         for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
315           if(fEMCALMatrix[mod]){
316             if(DebugLevel() > 1) 
317               fEMCALMatrix[mod]->Print();
318             fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;  
319           }
320           fGeoMatrixSet=kTRUE;
321         }//SM loop
322       }//Load matrices
323       else if(!gGeoManager){
324         printf("AliAnalysisTaskCaloFilter::UserExec() - Get geo matrices from data\n");
325         //Still not implemented in AOD, just a workaround to be able to work at least with ESDs 
326         if(!strcmp(event->GetName(),"AliAODEvent")) {
327           if(DebugLevel() > 1) 
328             printf("AliAnalysisTaskCaloFilter Use ideal geometry, values geometry matrix not kept in AODs.\n");
329         }//AOD
330         else {  
331           if(DebugLevel() > 1) printf("AliAnalysisTaskCaloFilter Load Misaligned matrices. \n");
332           AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
333           if(!esd) {
334             printf("AliAnalysisTaskCaloFilter::UserExec() - This event does not contain ESDs?");
335             return;
336           }
337           for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
338             //if(DebugLevel() > 1) 
339             esd->GetEMCALMatrix(mod)->Print();
340             if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
341           } 
342           fGeoMatrixSet=kTRUE;
343         }//ESD
344       }//Load matrices from Data 
345       
346       //Recover time dependent corrections, put then in recalibration histograms. Do it once
347       fEMCALRecoUtils->SetRunDependentCorrections(InputEvent()->GetRunNumber());
348
349     }//first event
350     
351     
352     //Cluster Loop
353     for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
354       
355       AliVCluster * cluster = event->GetCaloCluster(iClust);
356       if(cluster->IsPHOS()) continue ;
357
358       Float_t position[]={0,0,0};
359       if(DebugLevel() > 2)
360         printf("Check cluster %d for bad channels and close to border\n",cluster->GetID());
361       if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;        
362       //      if(!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, cluster, event->GetEMCALCells())) {
363       //        printf("Finally reject\n");
364       //        continue;
365       //      }
366       if(DebugLevel() > 2)
367       { 
368         printf("Filter, before  : i %d, E %f, dispersion %f, m02 %f, m20 %f, distToBad %f\n",iClust,cluster->E(),
369                cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20(), cluster->GetDistanceToBadChannel());
370         cluster->GetPosition(position);
371         printf("Filter, before  : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
372       }
373       
374       //Recalculate distance to bad channels, if new list of bad channels provided
375       fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, event->GetEMCALCells(), cluster);
376
377       if(fEMCALRecoUtils->IsRecalibrationOn())  {
378         fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, cluster, event->GetEMCALCells());
379         fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, event->GetEMCALCells(),cluster);
380         fEMCALRecoUtils->RecalculateClusterPID(cluster);
381       }
382             
383       fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, event->GetEMCALCells(),cluster);
384       
385       if(DebugLevel() > 2)
386       { 
387         printf("Filter, after   : i %d, E %f, dispersion %f, m02 %f, m20 %f, distToBad %f\n",cluster->GetID(),cluster->E(),
388                cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20(), cluster->GetDistanceToBadChannel());
389         cluster->GetPosition(position);
390         printf("Filter, after   : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
391       }    
392       
393       cluster->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(cluster));
394       
395     }
396     //Recalculate track-matching
397     fEMCALRecoUtils->FindMatches(event,0,fEMCALGeo);
398     
399   } // corrections in EMCAL
400   
401   //-------------------------------------------
402   // Now loop on clusters to fill AODs
403   //-------------------------------------------
404   for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
405     
406     AliVCluster * cluster = event->GetCaloCluster(iClust);
407     
408     //Check which calorimeter information we want to keep.
409     
410     if(fCaloFilter!=kBoth){
411       if     (fCaloFilter==kPHOS  && cluster->IsEMCAL()) continue ;
412       else if(fCaloFilter==kEMCAL && cluster->IsPHOS())  continue ;
413     }  
414     
415     //Temporary trick, FIXME
416     Float_t dR = cluster->GetTrackDx();
417     Float_t dZ = cluster->GetTrackDz();
418     if(DebugLevel() > 2)
419       printf("Original residuals : dZ %f, dR %f\n ",dZ, dR);
420     //--------------------------------------------------------------
421     //If EMCAL and corrections done, get the new matching parameters, do not copy noisy clusters
422     if(cluster->IsEMCAL() && fCorrect && esdevent){
423       if(DebugLevel() > 2)
424         printf("Check cluster %d for bad channels and close to border\n",cluster->GetID());
425       if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;        
426       
427       //      if(!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, cluster, event->GetEMCALCells())) {
428       //        printf("Finally reject\n");
429       //        continue;
430       //      }
431
432       if(fEMCALRecoUtils->IsExoticCluster(cluster, InputEvent()->GetEMCALCells(),InputEvent()->GetBunchCrossNumber())) continue;        
433
434       fEMCALRecoUtils->GetMatchedResiduals(cluster->GetID(),dR,dZ);
435       cluster->SetTrackDistance(dR,dZ);
436     }
437     
438     if(DebugLevel() > 2){
439       if(cluster->IsEMCAL()) printf("EMCAL Track-Cluster Residuals : dZ %f, dR %f\n ",dZ, dR);
440       if(cluster->IsPHOS())  printf("PHOS  Track-Cluster Residuals : dZ %f, dR %f\n ",dZ, dR);
441     }
442     
443     //--------------------------------------------------------------
444
445     //Now fill AODs
446     
447     Int_t id       = cluster->GetID();
448     Float_t energy = cluster->E();
449     cluster->GetPosition(posF);
450     
451     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) 
452       AliAODCaloCluster(id,
453                         cluster->GetNLabels(),
454                         cluster->GetLabels(),
455                         energy,
456                         posF,
457                         NULL,
458                         cluster->GetType());
459     
460     caloCluster->SetChi2(dZ);//Temporary trick, FIXME
461     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
462                                 cluster->GetDispersion(),
463                                 cluster->GetM20(), cluster->GetM02(),
464                                 dR,  //Temporary trick, FIXME
465                                 cluster->GetNExMax(),cluster->GetTOF()) ;
466     
467     caloCluster->SetPIDFromESD(cluster->GetPID());
468     caloCluster->SetNCells(cluster->GetNCells());
469     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
470     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
471     
472     if(DebugLevel() > 2)
473     { 
474       printf("Filter, aod     : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",caloCluster->GetID(),caloCluster->E(),
475              caloCluster->GetDispersion(),caloCluster->GetM02(),caloCluster->GetM20());
476       caloCluster->GetPosition(posF);
477       printf("Filter, aod     : i %d, x %f, y %f, z %f\n",caloCluster->GetID(), posF[0], posF[1], posF[2]);
478     }    
479     
480     //Matched tracks, just to know if there was any match, the track pointer is useless.
481     //Temporary trick, FIXME
482     if(esdevent){
483       if(TMath::Abs(dR) < 990 && TMath::Abs(dZ) < 990) { //Default value in PHOS 999, in EMCAL 1024, why?
484         if(DebugLevel() > 2) 
485           printf("*** Cluster Track-Matched *** dR %f, dZ %f\n",caloCluster->GetEmcCpvDistance(),caloCluster->Chi2());
486         caloCluster->AddTrackMatched(new AliAODTrack); 
487       }// fill the array with one entry to signal a possible match
488       //TArrayI* matchedT =     ((AliESDCaloCluster*)cluster)->GetTracksMatched();
489       //if (InputEvent()->GetNumberOfTracks() > 0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {        
490       //  for (Int_t im = 0; im < matchedT->GetSize(); im++) {
491       //    Int_t iESDtrack = matchedT->At(im);;
492       //    if ((AliVTrack*)InputEvent()->GetTrack(iESDtrack) != 0) {
493       //      caloCluster->AddTrackMatched((AliVTrack*)InputEvent()->GetTrack(iESDtrack));
494       //    }
495       //  }
496       //}// There is at least a match with a track
497     }
498   } 
499   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters      
500   // end of loop on calo clusters
501   
502   // fill EMCAL cell info
503   if ((fCaloFilter==kBoth ||  fCaloFilter==kEMCAL) && event->GetEMCALCells()) { // protection against missing ESD information
504     AliVCaloCells &eventEMcells = *(event->GetEMCALCells());
505     Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
506     
507     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
508     aodEMcells.CreateContainer(nEMcell);
509     aodEMcells.SetType(AliVCaloCells::kEMCALCell);
510     Double_t calibFactor = 1.;   
511     for (Int_t iCell = 0; iCell < nEMcell; iCell++) { 
512       Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; 
513       fEMCALGeo->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta); 
514       fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);       
515       
516       if(fCorrect && fEMCALRecoUtils->IsRecalibrationOn()){ 
517         calibFactor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
518       }
519       
520       if(!fEMCALRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
521         aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
522         //printf("GOOD channel\n");
523       }
524       else {
525         aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
526         //printf("BAD channel\n");
527       }
528     }
529     aodEMcells.Sort();
530   }
531   
532   // fill PHOS cell info
533   if ((fCaloFilter==kBoth ||  fCaloFilter==kPHOS) && event->GetPHOSCells()) { // protection against missing ESD information
534     AliVCaloCells &eventPHcells = *(event->GetPHOSCells());
535     Int_t nPHcell = eventPHcells.GetNumberOfCells() ;
536     
537     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
538     aodPHcells.CreateContainer(nPHcell);
539     aodPHcells.SetType(AliVCaloCells::kPHOSCell);
540     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {      
541       aodPHcells.SetCell(iCell,eventPHcells.GetCellNumber(iCell),eventPHcells.GetAmplitude(iCell));
542     }
543     aodPHcells.Sort();
544   }
545   
546    PostData(1, fEventNtuple);
547   
548   //return;
549 }
550
551 //____________________________________________________________________________
552 Bool_t AliAnalysisTaskCaloFilter::CheckForPrimaryVertex(){
553   //Check if the vertex was well reconstructed, copy from V0Reader of conversion group
554   //It only works for ESDs
555   
556   AliESDEvent * event = dynamic_cast<AliESDEvent*> (InputEvent());
557   if(!event) return kFALSE;
558   
559   if(event->GetPrimaryVertexTracks()->GetNContributors() > 0) {
560     return kTRUE;
561   }
562   
563   if(event->GetPrimaryVertexTracks()->GetNContributors() < 1) {
564     // SPD vertex
565     if(event->GetPrimaryVertexSPD()->GetNContributors() > 0) {
566       //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
567       return kTRUE;
568       
569     }
570     if(event->GetPrimaryVertexSPD()->GetNContributors() < 1) {
571       //      cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
572       return kFALSE;
573     }
574   }
575   return kFALSE;
576
577 }
578
579
580 //_____________________________________________________
581 void AliAnalysisTaskCaloFilter::PrintInfo(){
582
583   //Print settings
584
585   printf("TASK: AnalysisCaloFilter \n");
586   printf("\t Not only filter, correct Clusters? %d\n",fCorrect);
587   printf("\t Calorimeter Filtering Option     ? %d\n",fCaloFilter);
588   //printf("\t Use handmade geo matrices?   EMCAL %d, PHOS %d\n",fLoadEMCALMatrices, fLoadPHOSMatrices);
589   printf("\t Use handmade geo matrices?   EMCAL %d, PHOS 0\n",fLoadEMCALMatrices);
590 }
591
592 //_____________________________________________________
593 //void AliAnalysisTaskCaloFilter::LocalInit()
594 //{
595 //      // Local Initialization
596 //      
597 //      // Create cuts/param objects and publish to slot
598 //      const Int_t buffersize = 255;
599 //      char onePar[buffersize] ;
600 //      fCuts = new TList();
601 //  
602 //      snprintf(onePar,buffersize, "Calorimeter Filtering Option %d", fCaloFilter) ;
603 //      fCuts->Add(new TObjString(onePar));
604 //      snprintf(onePar,buffersize, "Not only filter but correct? %d cells;", fCorrect) ;
605 //      fCuts->Add(new TObjString(onePar));
606 //  
607 //      // Post Data
608 //      PostData(2, fCuts);
609 //      
610 //}
611
612
613 //__________________________________________________
614 void AliAnalysisTaskCaloFilter::Terminate(Option_t */*option*/)
615 {
616   // Terminate analysis
617   //
618     if (fDebug > 1) printf("AnalysisCaloFilter: Terminate() \n");
619 }
620