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->SetTimeDependentCorrections(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->IsRejectExoticCluster() && fEMCALRecoUtils->IsExoticCluster(cluster)) 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(0x0);
 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