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   //Do Corrections in EMCAL
 200   //If EMCAL, and requested, correct energy, position ...
 201   //Need to do this in a separate loop before filling the ESDs because of the track matching recalculations
 202   if(fCorrect && (fCaloFilter==kEMCAL || fCaloFilter==kBoth) ) {
 203     //Cluster Loop
 204     for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
 205
 206       AliVCluster * cluster = event->GetCaloCluster(iClust);
 207       if(cluster->IsPHOS()) continue ;
 208
 209       Float_t position[]={0,0,0};
 210       if(DebugLevel() > 2)
 211         printf("Check cluster %d for bad channels and close to border\n",cluster->GetID());
 212       if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
 213       //      if(!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, cluster, event->GetEMCALCells())) {
 214       //        printf("Finally reject\n");
 215       //        continue;
 216       //      }
 217       if(DebugLevel() > 2)
 218       {
 219         printf("Filter, before  : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClust,cluster->E(),
 220                cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20());
 221         cluster->GetPosition(position);
 222         printf("Filter, before  : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
 223       }
 224
 225       if(fEMCALRecoUtils->IsRecalibrationOn())  {
 226         fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, cluster, event->GetEMCALCells());
 227         fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, event->GetEMCALCells(),cluster);
 228         fEMCALRecoUtils->RecalculateClusterPID(cluster);
 229       }
 230       cluster->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(cluster));
 231
 232       fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, event->GetEMCALCells(),cluster);
 233
 234       if(DebugLevel() > 2)
 235       {
 236         printf("Filter, after   : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",cluster->GetID(),cluster->E(),
 237                cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20());
 238         cluster->GetPosition(position);
 239         printf("Filter, after   : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
 240       }
 241
 242     }
 243     //Recalculate track-matching
 244     fEMCALRecoUtils->FindMatches(event);
 245
 246   } // corrections in EMCAL
 247
 248   //Now loop on clusters to fill AODs
 249   for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
 250
 251     AliVCluster * cluster = event->GetCaloCluster(iClust);
 252
 253     //Check which calorimeter information we want to keep.
 254
 255     if(fCaloFilter!=kBoth){
 256       if     (fCaloFilter==kPHOS  && cluster->IsEMCAL()) continue ;
 257       else if(fCaloFilter==kEMCAL && cluster->IsPHOS())  continue ;
 258     }
 259
 260     //Temporary trick, FIXME
 261     Float_t dR = cluster->GetTrackDx();
 262     Float_t dZ = cluster->GetTrackDz();
 263     if(DebugLevel() > 2)
 264       printf("Original residuals : dZ %f, dR %f\n ",dZ, dR);
 265     //--------------------------------------------------------------
 266     //If EMCAL and corrections done, get the new matching parameters, do not copy noisy clusters
 267     if(cluster->IsEMCAL() && fCorrect){
 268       if(DebugLevel() > 2)
 269         printf("Check cluster %d for bad channels and close to border\n",cluster->GetID());
 270       if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
 271       //      if(!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, cluster, event->GetEMCALCells())) {
 272       //        printf("Finally reject\n");
 273       //        continue;
 274       //      }
 275
 276       fEMCALRecoUtils->GetMatchedResiduals(cluster->GetID(),dR,dZ);
 277       if(DebugLevel() > 2)
 278         printf("Corrected Residuals : dZ %f, dR %f\n ",dZ, dR);
 279
 280     }
 281     //--------------------------------------------------------------
 282
 283     //Now fill AODs
 284
 285     Int_t id       = cluster->GetID();
 286     Float_t energy = cluster->E();
 287     cluster->GetPosition(posF);
 288
 289     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
 290       AliAODCaloCluster(id,
 291                         0,
 292                         0x0,
 293                         energy,
 294                         posF,
 295                         NULL,
 296                         cluster->GetType());
 297
 298     caloCluster->SetChi2(dZ);//Temporary trick, FIXME
 299     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
 300                                 cluster->GetDispersion(),
 301                                 cluster->GetM20(), cluster->GetM02(),
 302                                 dR,  //Temporary trick, FIXME
 303                                 cluster->GetNExMax(),cluster->GetTOF()) ;
 304
 305     caloCluster->SetPIDFromESD(cluster->GetPID());
 306     caloCluster->SetNCells(cluster->GetNCells());
 307     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
 308     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
 309
 310     if(DebugLevel() > 2)
 311     {
 312       printf("Filter, aod     : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",caloCluster->GetID(),caloCluster->E(),
 313              caloCluster->GetDispersion(),caloCluster->GetM02(),caloCluster->GetM20());
 314       caloCluster->GetPosition(posF);
 315       printf("Filter, aod     : i %d, x %f, y %f, z %f\n",caloCluster->GetID(), posF[0], posF[1], posF[2]);
 316     }
 317
 318     //Matched tracks, just to know if there was any match, the track pointer is useless.
 319     //Temporary trick, FIXME
 320     if(bESD){
 321       if(TMath::Abs(dR) < 990 && TMath::Abs(dZ) < 990) { //Default value in PHOS 999, in EMCAL 1024, why?
 322         if(DebugLevel() > 2)
 323           printf("*** Cluster Track-Matched *** dR %f, dZ %f\n",caloCluster->GetEmcCpvDistance(),caloCluster->Chi2());
 324         caloCluster->AddTrackMatched(0x0);
 325       }// fill the array with one entry to signal a possible match
 326       //TArrayI* matchedT =     ((AliESDCaloCluster*)cluster)->GetTracksMatched();
 327       //if (InputEvent()->GetNumberOfTracks() > 0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {
 328       //  for (Int_t im = 0; im < matchedT->GetSize(); im++) {
 329       //    Int_t iESDtrack = matchedT->At(im);;
 330       //    if ((AliVTrack*)InputEvent()->GetTrack(iESDtrack) != 0) {
 331       //      caloCluster->AddTrackMatched((AliVTrack*)InputEvent()->GetTrack(iESDtrack));
 332       //    }
 333       //  }
 334       //}// There is at least a match with a track
 335     }
 336   }
 337   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
 338   // end of loop on calo clusters
 339
 340   // fill EMCAL cell info
 341   if ((fCaloFilter==kBoth ||  fCaloFilter==kEMCAL) && event->GetEMCALCells()) { // protection against missing ESD information
 342     AliVCaloCells &eventEMcells = *(event->GetEMCALCells());
 343     Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
 344
 345     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
 346     aodEMcells.CreateContainer(nEMcell);
 347     aodEMcells.SetType(AliVCaloCells::kEMCALCell);
 348     Double_t calibFactor = 1.;
 349     for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
 350       Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 351       fEMCALGeo->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
 352       fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 353
 354       if(fCorrect && fEMCALRecoUtils->IsRecalibrationOn()){
 355         calibFactor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 356       }
 357
 358       if(!fEMCALRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
 359         aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
 360         //printf("GOOD channel\n");
 361       }
 362       else {
 363         aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
 364         //printf("BAD channel\n");
 365
 366       }
 367     }
 368     aodEMcells.Sort();
 369   }
 370
 371   // fill PHOS cell info
 372   if ((fCaloFilter==kBoth ||  fCaloFilter==kPHOS) && event->GetPHOSCells()) { // protection against missing ESD information
 373     AliVCaloCells &eventPHcells = *(event->GetPHOSCells());
 374     Int_t nPHcell = eventPHcells.GetNumberOfCells() ;
 375
 376     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
 377     aodPHcells.CreateContainer(nPHcell);
 378     aodPHcells.SetType(AliVCaloCells::kPHOSCell);
 379     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
 380       aodPHcells.SetCell(iCell,eventPHcells.GetCellNumber(iCell),eventPHcells.GetAmplitude(iCell));
 381     }
 382     aodPHcells.Sort();
 383   }
 384
 385
 386   return;
 387 }
 388
 389 //_____________________________________________________
 390 void AliAnalysisTaskCaloFilter::PrintInfo(){
 391
 392   //Print settings
 393
 394   printf("TASK: AnalysisCaloFilter \n");
 395   printf("\t Not only filter, correct Clusters? %d\n",fCorrect);
 396   printf("\t Calorimeter Filtering Option     ? %d\n",fCaloFilter);
 397
 398 }
 399
 400 //_____________________________________________________
 401 //void AliAnalysisTaskCaloFilter::LocalInit()
 402 //{
 403 //      // Local Initialization
 404 //
 405 //      // Create cuts/param objects and publish to slot
 406 //      const Int_t buffersize = 255;
 407 //      char onePar[buffersize] ;
 408 //      fCuts = new TList();
 409 //
 410 //      snprintf(onePar,buffersize, "Calorimeter Filtering Option %d", fCaloFilter) ;
 411 //      fCuts->Add(new TObjString(onePar));
 412 //      snprintf(onePar,buffersize, "Not only filter but correct? %d cells;", fCorrect) ;
 413 //      fCuts->Add(new TObjString(onePar));
 414 //
 415 //      // Post Data
 416 //      PostData(2, fCuts);
 417 //
 418 //}
 419
 420
 421 //__________________________________________________
 422 void AliAnalysisTaskCaloFilter::Terminate(Option_t */*option*/)
 423 {
 424   // Terminate analysis
 425   //
 426     if (fDebug > 1) printf("AnalysisCaloFilter: Terminate() \n");
 427 }
 428