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     }//first event
 345
 346     //Cluster Loop
 347     for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
 348
 349       AliVCluster * cluster = event->GetCaloCluster(iClust);
 350       if(cluster->IsPHOS()) continue ;
 351
 352       Float_t position[]={0,0,0};
 353       if(DebugLevel() > 2)
 354         printf("Check cluster %d for bad channels and close to border\n",cluster->GetID());
 355       if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
 356       //      if(!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, cluster, event->GetEMCALCells())) {
 357       //        printf("Finally reject\n");
 358       //        continue;
 359       //      }
 360       if(DebugLevel() > 2)
 361       {
 362         printf("Filter, before  : i %d, E %f, dispersion %f, m02 %f, m20 %f, distToBad %f\n",iClust,cluster->E(),
 363                cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20(), cluster->GetDistanceToBadChannel());
 364         cluster->GetPosition(position);
 365         printf("Filter, before  : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
 366       }
 367
 368       //Recalculate distance to bad channels, if new list of bad channels provided
 369       fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, event->GetEMCALCells(), cluster);
 370
 371       if(fEMCALRecoUtils->IsRecalibrationOn())  {
 372         fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, cluster, event->GetEMCALCells());
 373         fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, event->GetEMCALCells(),cluster);
 374         fEMCALRecoUtils->RecalculateClusterPID(cluster);
 375       }
 376
 377       fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, event->GetEMCALCells(),cluster);
 378
 379       if(DebugLevel() > 2)
 380       {
 381         printf("Filter, after   : i %d, E %f, dispersion %f, m02 %f, m20 %f, distToBad %f\n",cluster->GetID(),cluster->E(),
 382                cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20(), cluster->GetDistanceToBadChannel());
 383         cluster->GetPosition(position);
 384         printf("Filter, after   : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
 385       }
 386
 387       cluster->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(cluster));
 388
 389     }
 390     //Recalculate track-matching
 391     fEMCALRecoUtils->FindMatches(event,0,fEMCALGeo);
 392
 393   } // corrections in EMCAL
 394
 395   //-------------------------------------------
 396   // Now loop on clusters to fill AODs
 397   //-------------------------------------------
 398   for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
 399
 400     AliVCluster * cluster = event->GetCaloCluster(iClust);
 401
 402     //Check which calorimeter information we want to keep.
 403
 404     if(fCaloFilter!=kBoth){
 405       if     (fCaloFilter==kPHOS  && cluster->IsEMCAL()) continue ;
 406       else if(fCaloFilter==kEMCAL && cluster->IsPHOS())  continue ;
 407     }
 408
 409     //Temporary trick, FIXME
 410     Float_t dR = cluster->GetTrackDx();
 411     Float_t dZ = cluster->GetTrackDz();
 412     if(DebugLevel() > 2)
 413       printf("Original residuals : dZ %f, dR %f\n ",dZ, dR);
 414     //--------------------------------------------------------------
 415     //If EMCAL and corrections done, get the new matching parameters, do not copy noisy clusters
 416     if(cluster->IsEMCAL() && fCorrect && esdevent){
 417       if(DebugLevel() > 2)
 418         printf("Check cluster %d for bad channels and close to border\n",cluster->GetID());
 419       if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
 420
 421       //      if(!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, cluster, event->GetEMCALCells())) {
 422       //        printf("Finally reject\n");
 423       //        continue;
 424       //      }
 425
 426       if(fEMCALRecoUtils->IsExoticCluster(cluster, InputEvent()->GetEMCALCells(),InputEvent()->GetBunchCrossNumber())) continue;
 427
 428       fEMCALRecoUtils->GetMatchedResiduals(cluster->GetID(),dR,dZ);
 429       cluster->SetTrackDistance(dR,dZ);
 430     }
 431
 432     if(DebugLevel() > 2){
 433       if(cluster->IsEMCAL()) printf("EMCAL Track-Cluster Residuals : dZ %f, dR %f\n ",dZ, dR);
 434       if(cluster->IsPHOS())  printf("PHOS  Track-Cluster Residuals : dZ %f, dR %f\n ",dZ, dR);
 435     }
 436
 437     //--------------------------------------------------------------
 438
 439     //Now fill AODs
 440
 441     Int_t id       = cluster->GetID();
 442     Float_t energy = cluster->E();
 443     cluster->GetPosition(posF);
 444
 445     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
 446       AliAODCaloCluster(id,
 447                         cluster->GetNLabels(),
 448                         cluster->GetLabels(),
 449                         energy,
 450                         posF,
 451                         NULL,
 452                         cluster->GetType());
 453
 454     caloCluster->SetChi2(dZ);//Temporary trick, FIXME
 455     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
 456                                 cluster->GetDispersion(),
 457                                 cluster->GetM20(), cluster->GetM02(),
 458                                 dR,  //Temporary trick, FIXME
 459                                 cluster->GetNExMax(),cluster->GetTOF()) ;
 460
 461     caloCluster->SetPIDFromESD(cluster->GetPID());
 462     caloCluster->SetNCells(cluster->GetNCells());
 463     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
 464     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
 465
 466     if(DebugLevel() > 2)
 467     {
 468       printf("Filter, aod     : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",caloCluster->GetID(),caloCluster->E(),
 469              caloCluster->GetDispersion(),caloCluster->GetM02(),caloCluster->GetM20());
 470       caloCluster->GetPosition(posF);
 471       printf("Filter, aod     : i %d, x %f, y %f, z %f\n",caloCluster->GetID(), posF[0], posF[1], posF[2]);
 472     }
 473
 474     //Matched tracks, just to know if there was any match, the track pointer is useless.
 475     //Temporary trick, FIXME
 476     if(esdevent){
 477       if(TMath::Abs(dR) < 990 && TMath::Abs(dZ) < 990) { //Default value in PHOS 999, in EMCAL 1024, why?
 478         if(DebugLevel() > 2)
 479           printf("*** Cluster Track-Matched *** dR %f, dZ %f\n",caloCluster->GetEmcCpvDistance(),caloCluster->Chi2());
 480         caloCluster->AddTrackMatched(new AliAODTrack);
 481       }// fill the array with one entry to signal a possible match
 482       //TArrayI* matchedT =     ((AliESDCaloCluster*)cluster)->GetTracksMatched();
 483       //if (InputEvent()->GetNumberOfTracks() > 0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {
 484       //  for (Int_t im = 0; im < matchedT->GetSize(); im++) {
 485       //    Int_t iESDtrack = matchedT->At(im);;
 486       //    if ((AliVTrack*)InputEvent()->GetTrack(iESDtrack) != 0) {
 487       //      caloCluster->AddTrackMatched((AliVTrack*)InputEvent()->GetTrack(iESDtrack));
 488       //    }
 489       //  }
 490       //}// There is at least a match with a track
 491     }
 492   }
 493   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
 494   // end of loop on calo clusters
 495
 496   // fill EMCAL cell info
 497   if ((fCaloFilter==kBoth ||  fCaloFilter==kEMCAL) && event->GetEMCALCells()) { // protection against missing ESD information
 498     AliVCaloCells &eventEMcells = *(event->GetEMCALCells());
 499     Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
 500
 501     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
 502     aodEMcells.CreateContainer(nEMcell);
 503     aodEMcells.SetType(AliVCaloCells::kEMCALCell);
 504     Double_t calibFactor = 1.;
 505     for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
 506       Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 507       fEMCALGeo->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
 508       fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 509
 510       if(fCorrect && fEMCALRecoUtils->IsRecalibrationOn()){
 511         calibFactor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 512       }
 513
 514       if(!fEMCALRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
 515         aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
 516                            eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
 517         //printf("GOOD channel\n");
 518       }
 519       else {
 520         aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
 521         //printf("BAD channel\n");
 522       }
 523     }
 524     aodEMcells.Sort();
 525   }
 526
 527   // fill PHOS cell info
 528   if ((fCaloFilter==kBoth ||  fCaloFilter==kPHOS) && event->GetPHOSCells()) { // protection against missing ESD information
 529     AliVCaloCells &eventPHcells = *(event->GetPHOSCells());
 530     Int_t nPHcell = eventPHcells.GetNumberOfCells() ;
 531
 532     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
 533     aodPHcells.CreateContainer(nPHcell);
 534     aodPHcells.SetType(AliVCaloCells::kPHOSCell);
 535     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
 536       aodPHcells.SetCell(iCell,eventPHcells.GetCellNumber(iCell),eventPHcells.GetAmplitude(iCell),
 537                          eventPHcells.GetTime(iCell),eventPHcells.GetMCLabel(iCell),eventPHcells.GetEFraction(iCell));
 538     }
 539     aodPHcells.Sort();
 540   }
 541
 542    PostData(1, fEventNtuple);
 543
 544   //return;
 545 }
 546
 547 //____________________________________________________________________________
 548 Bool_t AliAnalysisTaskCaloFilter::CheckForPrimaryVertex(){
 549   //Check if the vertex was well reconstructed, copy from V0Reader of conversion group
 550   //It only works for ESDs
 551
 552   AliESDEvent * event = dynamic_cast<AliESDEvent*> (InputEvent());
 553   if(!event) return kFALSE;
 554
 555   if(event->GetPrimaryVertexTracks()->GetNContributors() > 0) {
 556     return kTRUE;
 557   }
 558
 559   if(event->GetPrimaryVertexTracks()->GetNContributors() < 1) {
 560     // SPD vertex
 561     if(event->GetPrimaryVertexSPD()->GetNContributors() > 0) {
 562       //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
 563       return kTRUE;
 564
 565     }
 566     if(event->GetPrimaryVertexSPD()->GetNContributors() < 1) {
 567       //      cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
 568       return kFALSE;
 569     }
 570   }
 571   return kFALSE;
 572
 573 }
 574
 575
 576 //_____________________________________________________
 577 void AliAnalysisTaskCaloFilter::PrintInfo(){
 578
 579   //Print settings
 580
 581   printf("TASK: AnalysisCaloFilter \n");
 582   printf("\t Not only filter, correct Clusters? %d\n",fCorrect);
 583   printf("\t Calorimeter Filtering Option     ? %d\n",fCaloFilter);
 584   //printf("\t Use handmade geo matrices?   EMCAL %d, PHOS %d\n",fLoadEMCALMatrices, fLoadPHOSMatrices);
 585   printf("\t Use handmade geo matrices?   EMCAL %d, PHOS 0\n",fLoadEMCALMatrices);
 586 }
 587
 588 //_____________________________________________________
 589 //void AliAnalysisTaskCaloFilter::LocalInit()
 590 //{
 591 //      // Local Initialization
 592 //
 593 //      // Create cuts/param objects and publish to slot
 594 //      const Int_t buffersize = 255;
 595 //      char onePar[buffersize] ;
 596 //      fCuts = new TList();
 597 //
 598 //      snprintf(onePar,buffersize, "Calorimeter Filtering Option %d", fCaloFilter) ;
 599 //      fCuts->Add(new TObjString(onePar));
 600 //      snprintf(onePar,buffersize, "Not only filter but correct? %d cells;", fCorrect) ;
 601 //      fCuts->Add(new TObjString(onePar));
 602 //
 603 //      // Post Data
 604 //      PostData(2, fCuts);
 605 //
 606 //}
 607
 608
 609 //__________________________________________________
 610 void AliAnalysisTaskCaloFilter::Terminate(Option_t */*option*/)
 611 {
 612   // Terminate analysis
 613   //
 614     if (fDebug > 1) printf("AnalysisCaloFilter: Terminate() \n");
 615 }
 616