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