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 #include "AliAnalysisTaskCaloFilter.h"
  29 #include "AliESDEvent.h"
  30 #include "AliAODEvent.h"
  31 #include "AliLog.h"
  32 #include "AliVCluster.h"
  33 #include "AliVCaloCells.h"
  34 #include "AliEMCALRecoUtils.h"
  35 #include "AliEMCALGeometry.h"
  36 #include "AliVEventHandler.h"
  37 #include "AliAnalysisManager.h"
  38 #include "AliInputEventHandler.h"
  39 #include "AliESDtrackCuts.h"
  40
  41 ClassImp(AliAnalysisTaskCaloFilter)
  42
  43 ////////////////////////////////////////////////////////////////////////
  44
  45 AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter():
  46   AliAnalysisTaskSE(), //fCuts(0x0),
  47   fCaloFilter(0), fCorrect(kFALSE),
  48   fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEAR"),
  49   fEMCALRecoUtils(new AliEMCALRecoUtils),
  50   fESDtrackCuts(0), fTrackMultEtaCut(0.9)
  51 {
  52   // Default constructor
  53   fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
  54
  55 }
  56
  57 //__________________________________________________
  58 AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter(const char* name):
  59   AliAnalysisTaskSE(name), //fCuts(0x0),
  60   fCaloFilter(0), fCorrect(kFALSE),
  61   fEMCALGeo(0x0),fEMCALGeoName("EMCAL_FIRSTYEAR"),
  62   fEMCALRecoUtils(new AliEMCALRecoUtils),
  63   fESDtrackCuts(0), fTrackMultEtaCut(0.9)
  64 {
  65   // Constructor
  66   fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
  67
  68 }
  69
  70 //__________________________________________________
  71 AliAnalysisTaskCaloFilter::~AliAnalysisTaskCaloFilter()
  72 {
  73   //Destructor.
  74
  75   if(fEMCALGeo)       delete fEMCALGeo;
  76   if(fEMCALRecoUtils) delete fEMCALRecoUtils;
  77   if(fESDtrackCuts)   delete fESDtrackCuts;
  78 }
  79
  80 //__________________________________________________
  81 void AliAnalysisTaskCaloFilter::UserCreateOutputObjects()
  82 {
  83   // Init geometry
  84
  85   fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
  86
  87 }
  88
  89 //__________________________________________________
  90 void AliAnalysisTaskCaloFilter::UserExec(Option_t */*option*/)
  91 {
  92   // Execute analysis for current event
  93   //
  94
  95   if (fDebug > 0)
  96     printf("CaloFilter: Analysing event # %d\n", (Int_t)Entry());
  97
  98   // Copy input ESD or AOD header, vertex, CaloClusters and CaloCells to output AOD
  99
 100   AliVEvent* event = InputEvent();
 101   if(!event) {
 102     printf("AliAnalysisTaskCaloFilter::CreateAODFromESD() - This event does not contain Input?");
 103     return;
 104   }
 105
 106   //Magic line to write events to file
 107   AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
 108
 109   Bool_t bAOD = kFALSE;
 110   if(!strcmp(event->GetName(),"AliAODEvent")) bAOD=kTRUE;
 111   Bool_t bESD = kFALSE;
 112   if(!strcmp(event->GetName(),"AliESDEvent")) bESD=kTRUE;
 113
 114   //Get track multiplicity
 115   Int_t trackMult = 0;
 116   if(bESD){
 117     Int_t nTracks   = InputEvent()->GetNumberOfTracks() ;
 118     for (Int_t itrack =  0; itrack <  nTracks; itrack++) {////////////// track loop
 119       AliVTrack * track = (AliVTrack*)InputEvent()->GetTrack(itrack) ; // retrieve track from esd
 120       if(!fESDtrackCuts->AcceptTrack((AliESDtrack*)track)) continue;
 121       //Count the tracks in eta < 0.9
 122       if(TMath::Abs(track->Eta())< fTrackMultEtaCut) trackMult++;
 123     }
 124   }
 125
 126   // set arrays and pointers
 127   Float_t posF[3];
 128   Double_t pos[3];
 129
 130   Double_t covVtx[6];
 131
 132   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 133
 134   AliAODHeader* header = AODEvent()->GetHeader();
 135
 136   header->SetRunNumber(event->GetRunNumber());
 137   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
 138   if(bESD)
 139     header->SetFiredTriggerClasses(((AliESDEvent*)event)->GetFiredTriggerClasses());
 140   header->SetTriggerMask(event->GetTriggerMask());
 141   header->SetTriggerCluster(event->GetTriggerCluster());
 142
 143   header->SetBunchCrossNumber(event->GetBunchCrossNumber());
 144   header->SetOrbitNumber(event->GetOrbitNumber());
 145   header->SetPeriodNumber(event->GetPeriodNumber());
 146   header->SetEventType(event->GetEventType());
 147   header->SetMuonMagFieldScale(-999.); // FIXME
 148   //printf("Track Multiplicity for eta < %f: %d \n",fTrackMultEtaCut,trackMult);
 149   header->SetCentrality((Double_t)trackMult);        // FIXME
 150   //printf("Centrality %f\n",header->GetCentrality());
 151
 152   header->SetTriggerMask(event->GetTriggerMask());
 153   header->SetTriggerCluster(event->GetTriggerCluster());
 154   header->SetMagneticField(event->GetMagneticField());
 155   header->SetZDCN1Energy(event->GetZDCN1Energy());
 156   header->SetZDCP1Energy(event->GetZDCP1Energy());
 157   header->SetZDCN2Energy(event->GetZDCN2Energy());
 158   header->SetZDCP2Energy(event->GetZDCP2Energy());
 159   header->SetZDCEMEnergy(event->GetZDCEMEnergy(0),event->GetZDCEMEnergy(1));
 160   Float_t diamxy[2]={event->GetDiamondX(),event->GetDiamondY()};
 161   Float_t diamcov[3]; event->GetDiamondCovXY(diamcov);
 162   header->SetDiamond(diamxy,diamcov);
 163   if(bESD){
 164     header->SetDiamondZ(((AliESDEvent*)event)->GetDiamondZ(),((AliESDEvent*)event)->GetSigma2DiamondZ());
 165   }
 166   //
 167   //
 168   Int_t nVertices = 1 ;/* = prim. vtx*/;
 169   Int_t nCaloClus = event->GetNumberOfCaloClusters();
 170
 171   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
 172
 173   // Access to the AOD container of vertices
 174   TClonesArray &vertices = *(AODEvent()->GetVertices());
 175   Int_t jVertices=0;
 176
 177   // Add primary vertex. The primary tracks will be defined
 178   // after the loops on the composite objects (V0, cascades, kinks)
 179   event->GetPrimaryVertex()->GetXYZ(pos);
 180   Float_t chi = 0;
 181   if      (bESD){
 182     ((AliESDEvent*)event)->GetPrimaryVertex()->GetCovMatrix(covVtx);
 183     chi = ((AliESDEvent*)event)->GetPrimaryVertex()->GetChi2toNDF();
 184   }
 185   else if (bAOD){
 186     ((AliAODEvent*)event)->GetPrimaryVertex()->GetCovMatrix(covVtx);
 187     chi = ((AliAODEvent*)event)->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
 188   }
 189
 190   AliAODVertex * primary = new(vertices[jVertices++])
 191     AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
 192   primary->SetName(event->GetPrimaryVertex()->GetName());
 193   primary->SetTitle(event->GetPrimaryVertex()->GetTitle());
 194
 195   // Access to the AOD container of clusters
 196   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
 197   Int_t jClusters=0;
 198
 199   for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
 200
 201     AliVCluster * cluster = event->GetCaloCluster(iClust);
 202
 203     //Check which calorimeter information we want to keep.
 204
 205     if(fCaloFilter!=kBoth){
 206       if     (fCaloFilter==kPHOS  && cluster->IsEMCAL()) continue ;
 207       else if(fCaloFilter==kEMCAL && cluster->IsPHOS())  continue ;
 208     }
 209
 210     //--------------------------------------------------------------
 211     //If EMCAL, and requested, correct energy, position ...
 212     if(cluster->IsEMCAL() && fCorrect){
 213       Float_t position[]={0,0,0};
 214       if(DebugLevel() > 2)
 215         printf("Check cluster %d for bad channels and close to border\n",cluster->GetID());
 216       if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
 217 //      if(!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, cluster, event->GetEMCALCells())) {
 218 //        printf("Finally reject\n");
 219 //        continue;
 220 //      }
 221       if(DebugLevel() > 2)
 222       {
 223         printf("Filter, before  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClust,cluster->E(),
 224                cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20());
 225         cluster->GetPosition(position);
 226         printf("Filter, before  : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
 227       }
 228
 229       if(fEMCALRecoUtils->IsRecalibrationOn())  {
 230         fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, cluster, event->GetEMCALCells());
 231         fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, event->GetEMCALCells(),cluster);
 232         fEMCALRecoUtils->RecalculateClusterPID(cluster);
 233
 234       }
 235       cluster->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(cluster));
 236
 237       fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, event->GetEMCALCells(),cluster);
 238
 239       if(DebugLevel() > 2)
 240       {
 241         printf("Filter, after   : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",cluster->GetID(),cluster->E(),
 242                cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20());
 243         cluster->GetPosition(position);
 244         printf("Filter, after   : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
 245       }
 246
 247     }
 248     //--------------------------------------------------------------
 249
 250     //Now fill AODs
 251     Int_t id       = cluster->GetID();
 252     Float_t energy = cluster->E();
 253     cluster->GetPosition(posF);
 254
 255     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
 256       AliAODCaloCluster(id,
 257                         0,
 258                         0x0,
 259                         energy,
 260                         posF,
 261                         NULL,
 262                         cluster->GetType());
 263
 264     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
 265                                 cluster->GetDispersion(),
 266                                 cluster->GetM20(), cluster->GetM02(),
 267                                 cluster->GetEmcCpvDistance(),
 268                                 cluster->GetNExMax(),cluster->GetTOF()) ;
 269
 270     caloCluster->SetPIDFromESD(cluster->GetPID());
 271     caloCluster->SetNCells(cluster->GetNCells());
 272     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
 273     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
 274
 275     if(DebugLevel() > 2)
 276     {
 277       printf("Filter, aod       : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",caloCluster->GetID(),caloCluster->E(),
 278              caloCluster->GetDispersion(),caloCluster->GetM02(),caloCluster->GetM20());
 279       caloCluster->GetPosition(posF);
 280       printf("Filter, aod       : i %d, x %f, y %f, z %f\n",caloCluster->GetID(), posF[0], posF[1], posF[2]);
 281     }
 282
 283     //Matched tracks, just to know if there was any match, the track pointer is useless.
 284     if(bESD){
 285       TArrayI* matchedT =       ((AliESDCaloCluster*)cluster)->GetTracksMatched();
 286       if (InputEvent()->GetNumberOfTracks() > 0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {
 287         for (Int_t im = 0; im < matchedT->GetSize(); im++) {
 288           Int_t iESDtrack = matchedT->At(im);;
 289           if ((AliVTrack*)InputEvent()->GetTrack(iESDtrack) != 0) {
 290             caloCluster->AddTrackMatched((AliVTrack*)InputEvent()->GetTrack(iESDtrack));
 291           }
 292         }
 293       }// There is at least a match with a track
 294     }
 295   }
 296   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
 297   // end of loop on calo clusters
 298
 299   // fill EMCAL cell info
 300   if ((fCaloFilter==kBoth ||  fCaloFilter==kEMCAL) && event->GetEMCALCells()) { // protection against missing ESD information
 301     AliVCaloCells &eventEMcells = *(event->GetEMCALCells());
 302     Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
 303
 304     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
 305     aodEMcells.CreateContainer(nEMcell);
 306     aodEMcells.SetType(AliVCaloCells::kEMCALCell);
 307     Double_t calibFactor = 1.;
 308     for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
 309       Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 310       fEMCALGeo->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
 311       fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 312
 313       if(fCorrect && fEMCALRecoUtils->IsRecalibrationOn()){
 314         calibFactor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 315       }
 316
 317       if(!fEMCALRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
 318         aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
 319         //printf("GOOD channel\n");
 320       }
 321       else {
 322         aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
 323         //printf("BAD channel\n");
 324
 325       }
 326     }
 327     aodEMcells.Sort();
 328   }
 329
 330   // fill PHOS cell info
 331   if ((fCaloFilter==kBoth ||  fCaloFilter==kPHOS) && event->GetPHOSCells()) { // protection against missing ESD information
 332     AliVCaloCells &eventPHcells = *(event->GetPHOSCells());
 333     Int_t nPHcell = eventPHcells.GetNumberOfCells() ;
 334
 335     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
 336     aodPHcells.CreateContainer(nPHcell);
 337     aodPHcells.SetType(AliVCaloCells::kPHOSCell);
 338     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
 339       aodPHcells.SetCell(iCell,eventPHcells.GetCellNumber(iCell),eventPHcells.GetAmplitude(iCell));
 340     }
 341     aodPHcells.Sort();
 342   }
 343
 344
 345   return;
 346 }
 347
 348 //_____________________________________________________
 349 void AliAnalysisTaskCaloFilter::PrintInfo(){
 350
 351   //Print settings
 352
 353   printf("TASK: AnalysisCaloFilter \n");
 354   printf("\t Not only filter, correct Clusters? %d\n",fCorrect);
 355   printf("\t Calorimeter Filtering Option     ? %d\n",fCaloFilter);
 356
 357 }
 358
 359 //_____________________________________________________
 360 //void AliAnalysisTaskCaloFilter::LocalInit()
 361 //{
 362 //      // Local Initialization
 363 //
 364 //      // Create cuts/param objects and publish to slot
 365 //      const Int_t buffersize = 255;
 366 //      char onePar[buffersize] ;
 367 //      fCuts = new TList();
 368 //
 369 //      snprintf(onePar,buffersize, "Calorimeter Filtering Option %d", fCaloFilter) ;
 370 //      fCuts->Add(new TObjString(onePar));
 371 //      snprintf(onePar,buffersize, "Not only filter but correct? %d cells;", fCorrect) ;
 372 //      fCuts->Add(new TObjString(onePar));
 373 //
 374 //      // Post Data
 375 //      PostData(2, fCuts);
 376 //
 377 //}
 378
 379
 380 //__________________________________________________
 381 void AliAnalysisTaskCaloFilter::Terminate(Option_t */*option*/)
 382 {
 383   // Terminate analysis
 384   //
 385     if (fDebug > 1) printf("AnalysisCaloFilter: Terminate() \n");
 386 }
 387