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 // Keep tracks, v0s, VZERO if requested
  22 // Select events containing a cluster or track avobe a given threshold
  23 // Copy of AliAnalysisTaskESDfilter.
  24 // Author: Gustavo Conesa Balbastre (INFN - Frascati)
  25 //////////////////////////////////////////////////////////
  26
  27 //Root
  28 #include "TGeoManager.h"
  29 #include "TFile.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 "AliAODHandler.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliInputEventHandler.h"
  43
  44 //EMCAL
  45 #include "AliEMCALRecoUtils.h"
  46 #include "AliEMCALGeometry.h"
  47
  48 #include "AliAnalysisTaskCaloFilter.h"
  49
  50 ClassImp(AliAnalysisTaskCaloFilter)
  51
  52 ////////////////////////////////////////////////////////////////////////
  53
  54 AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter():
  55 AliAnalysisTaskSE("CaloFilterTask"),
  56 fCaloFilter(0),           fEventSelection(),
  57 fAcceptAllMBEvent(kFALSE),fMBTriggerMask(AliVEvent::kMB),
  58 fCorrect(kFALSE),
  59 fEMCALGeo(0x0),           fEMCALGeoName("EMCAL_COMPLETE12SMV1"),
  60 fEMCALRecoUtils(new AliEMCALRecoUtils),
  61 fLoadEMCALMatrices(kFALSE), //fLoadPHOSMatrices(kFALSE),
  62 fGeoMatrixSet(kFALSE),
  63 fConfigName(""),          fFillAODFile(kTRUE),
  64 fFillMCParticles(kFALSE),
  65 fFillTracks(kFALSE),      fFillHybridTracks(kFALSE),
  66 fFillAllVertices(kFALSE), fFillv0s(kFALSE),
  67 fFillVZERO(kFALSE),
  68 fEMCALEnergyCut(0.),      fEMCALNcellsCut (0),
  69 fPHOSEnergyCut(0.),       fPHOSNcellsCut (0),
  70 fTrackPtCut(-1),
  71 fVzCut(100.),             fCheckEventVertex(kTRUE),
  72 fEvent(0x0),
  73 fESDEvent(0x0),           fAODEvent(0x0)
  74 {
  75   // Default constructor
  76
  77   fEventSelection[0] = kFALSE;
  78   fEventSelection[1] = kFALSE;
  79   fEventSelection[2] = kFALSE;
  80
  81   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
  82   //for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix[i]  = 0 ;
  83
  84 }
  85
  86 //__________________________________________________
  87 AliAnalysisTaskCaloFilter::AliAnalysisTaskCaloFilter(const char* name):
  88 AliAnalysisTaskSE(name),
  89 fCaloFilter(0),           fEventSelection(),
  90 fAcceptAllMBEvent(kFALSE),fMBTriggerMask(AliVEvent::kMB),
  91 fCorrect(kFALSE),
  92 fEMCALGeo(0x0),           fEMCALGeoName("EMCAL_COMPLETE12SMV1"),
  93 fEMCALRecoUtils(new AliEMCALRecoUtils),
  94 fLoadEMCALMatrices(kFALSE), //fLoadPHOSMatrices(kFALSE),
  95 fGeoMatrixSet(kFALSE),
  96 fConfigName(""),          fFillAODFile(kTRUE),
  97 fFillMCParticles(kFALSE),
  98 fFillTracks(kFALSE),      fFillHybridTracks(kFALSE),
  99 fFillAllVertices(kFALSE), fFillv0s(kFALSE),
 100 fFillVZERO(kFALSE),
 101 fEMCALEnergyCut(0.),      fEMCALNcellsCut(0),
 102 fPHOSEnergyCut(0.),       fPHOSNcellsCut(0),
 103 fTrackPtCut(-1),
 104 fVzCut(100.),             fCheckEventVertex(kTRUE),
 105 fEvent(0x0),
 106 fESDEvent(0x0),           fAODEvent(0x0)
 107 {
 108   // Constructor
 109
 110   fEventSelection[0] = kFALSE;
 111   fEventSelection[1] = kFALSE;
 112   fEventSelection[2] = kFALSE;
 113
 114   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
 115   //for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix[i]  = 0 ;
 116
 117 }
 118
 119 //__________________________________________________
 120 AliAnalysisTaskCaloFilter::~AliAnalysisTaskCaloFilter()
 121 {
 122   //Destructor.
 123
 124   if(fEMCALGeo)        delete fEMCALGeo;
 125   if(fEMCALRecoUtils)  delete fEMCALRecoUtils;
 126
 127 }
 128
 129 //_____________________________________________
 130 Bool_t AliAnalysisTaskCaloFilter::AcceptEvent()
 131 {
 132   // Define conditions to accept the event to be filtered
 133
 134   if(!AcceptEventVertex()) return kFALSE;
 135
 136   Bool_t eventSel = kFALSE;
 137
 138   Bool_t isMB = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fMBTriggerMask);
 139
 140   if     ( isMB && fAcceptAllMBEvent )                eventSel = kTRUE; // accept any MB event
 141
 142   else if( fEventSelection[0] && AcceptEventEMCAL() ) eventSel = kTRUE; // accept event depending on EMCAL activity
 143
 144   else if( fEventSelection[1] && AcceptEventPHOS () ) eventSel = kTRUE; // accept event depending on PHOS  activity
 145
 146   else if( fEventSelection[2] && AcceptEventTrack() ) eventSel = kTRUE; // accept event depending on Track activity
 147
 148   return eventSel ;
 149
 150 }
 151
 152 //__________________________________________________
 153 Bool_t AliAnalysisTaskCaloFilter::AcceptEventEMCAL()
 154 {
 155   // Accept event given there is a EMCAL cluster with enough energy, and not noisy, exotic
 156
 157   if(fCaloFilter==kPHOS)   return kTRUE; // accept
 158
 159   if(fEMCALEnergyCut <= 0) return kTRUE; // accept
 160
 161   Int_t           nCluster = InputEvent() -> GetNumberOfCaloClusters();
 162   AliVCaloCells * caloCell = InputEvent() -> GetEMCALCells();
 163   Int_t           bc       = InputEvent() -> GetBunchCrossNumber();
 164
 165   for(Int_t icalo = 0; icalo < nCluster; icalo++)
 166   {
 167     AliVCluster *clus = (AliVCluster*) (InputEvent()->GetCaloCluster(icalo));
 168
 169     if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fEMCALNcellsCut ) && ( clus->E() > fEMCALEnergyCut ) &&
 170        fEMCALRecoUtils->IsGoodCluster(clus,fEMCALGeo,caloCell,bc))
 171     {
 172
 173       if (fDebug > 0) printf("AliAnalysisTaskCaloFilter::AcceptEventEMCAL() - Accept :  E %2.2f > %2.2f, nCells %d > %d \n",
 174                              clus->E(), fEMCALEnergyCut, clus->GetNCells(), fEMCALNcellsCut);
 175
 176       return kTRUE;
 177     }
 178
 179   }// loop
 180
 181   if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::AcceptEventEMCAL() - Reject \n");
 182
 183   //printf("Fired %s\n",((AliESDEvent*)InputEvent())->GetFiredTriggerClasses().Data());
 184
 185   return kFALSE;
 186
 187 }
 188
 189 //_________________________________________________
 190 Bool_t AliAnalysisTaskCaloFilter::AcceptEventPHOS()
 191 {
 192   // Accept event given there is a PHOS cluster with enough energy and not noisy/exotic
 193
 194   if(fCaloFilter==kEMCAL) return kTRUE; // accept
 195
 196   if(fPHOSEnergyCut <= 0) return kTRUE; // accept
 197
 198   Int_t nCluster = InputEvent() -> GetNumberOfCaloClusters();
 199
 200   for(Int_t icalo = 0; icalo < nCluster; icalo++)
 201   {
 202     AliVCluster *clus = (AliVCluster*) (InputEvent()->GetCaloCluster(icalo));
 203
 204     if( ( clus->IsPHOS() ) && ( clus->GetNCells() > fPHOSNcellsCut ) && ( clus->E() > fPHOSEnergyCut ))
 205     {
 206
 207       if (fDebug > 0) printf("AliAnalysisTaskCaloFilter::AcceptEventPHOS() - Accept :  E %2.2f > %2.2f, nCells %d > %d \n",
 208                              clus->E(), fPHOSEnergyCut, clus->GetNCells(), fPHOSNcellsCut);
 209
 210       return kTRUE;
 211     }
 212
 213   }// loop
 214
 215   if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::AcceptEventPHOS() - Reject \n");
 216
 217   return kFALSE;
 218
 219 }
 220
 221 //__________________________________________________
 222 Bool_t AliAnalysisTaskCaloFilter::AcceptEventTrack()
 223 {
 224   // Accept event if there is a track avobe a certain pT
 225
 226   if(fTrackPtCut <= 0) return kTRUE; // accept
 227
 228   Double_t pTrack[3] = {0,0,0};
 229
 230   for (Int_t nTrack = 0; nTrack < fEvent->GetNumberOfTracks(); ++nTrack)
 231   {
 232     AliVTrack *track = (AliVTrack*) fEvent->GetTrack(nTrack);
 233
 234     // Select only hybrid tracks?
 235     if(fAODEvent && fFillHybridTracks && !((AliAODTrack*)track)->IsHybridGlobalConstrainedGlobal()) continue;
 236
 237     track->GetPxPyPz(pTrack) ;
 238
 239     TLorentzVector momentum(pTrack[0],pTrack[1],pTrack[2],0);
 240
 241     if(momentum.Pt() > fTrackPtCut)
 242     {
 243       if (fDebug > 0) printf("AliAnalysisTaskCaloFilter::AcceptEventTrack() - Accept :  pT %2.2f > %2.2f \n",
 244                              momentum.Pt(), fTrackPtCut);
 245
 246       return kTRUE;
 247     }
 248
 249   }
 250
 251   if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::AcceptEventTrack() - Reject \n");
 252
 253   return kFALSE;
 254
 255 }
 256
 257 //___________________________________________________
 258 Bool_t AliAnalysisTaskCaloFilter::AcceptEventVertex()
 259 {
 260   // Accept event with good vertex
 261
 262   Double_t v[3];
 263   InputEvent()->GetPrimaryVertex()->GetXYZ(v) ;
 264
 265   if(TMath::Abs(v[2]) > fVzCut)
 266   {
 267     if (fDebug > 0) printf("AliAnalysisTaskCaloFilter::AcceptEventVertex() - Vz Reject : vz %2.2f > %2.2f\n",v[2],fVzCut);
 268
 269     return kFALSE ;
 270   }
 271
 272   return CheckForPrimaryVertex();
 273 }
 274
 275 //_______________________________________________________
 276 Bool_t AliAnalysisTaskCaloFilter::CheckForPrimaryVertex()
 277 {
 278   //Check if the vertex was well reconstructed, copy from v0Reader of conversion group
 279   //It only works for ESDs
 280
 281   if(!fCheckEventVertex) return kTRUE;
 282
 283   // AODs
 284   if(!fESDEvent)
 285   {
 286     // Check that the vertex is not (0,0,0)
 287     Double_t v[3];
 288     InputEvent()->GetPrimaryVertex()->GetXYZ(v) ;
 289
 290     if(TMath::Abs(v[2]) < 1e-6 &&
 291        TMath::Abs(v[1]) < 1e-6 &&
 292        TMath::Abs(v[0]) < 1e-6 )
 293     {
 294       if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::CheckForPrimaryVertex() - Reject v(0,0,0) \n");
 295
 296       return kFALSE ;
 297     }
 298
 299     return kTRUE;
 300   }
 301
 302   // ESDs
 303   if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors() > 0)
 304   {
 305     return kTRUE;
 306   }
 307
 308   if(fESDEvent->GetPrimaryVertexTracks()->GetNContributors() < 1)
 309   {
 310     // SPD vertex
 311     if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors() > 0)
 312     {
 313       //cout<<"spd vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
 314       return kTRUE;
 315
 316     }
 317     if(fESDEvent->GetPrimaryVertexSPD()->GetNContributors() < 1)
 318     {
 319       //      cout<<"bad vertex type::"<< fESDEvent->GetPrimaryVertex()->GetName() << endl;
 320       if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::CheckForPrimaryVertex() - Reject, GetPrimaryVertexSPD()->GetNContributors() < 1 \n");
 321
 322       return kFALSE;
 323     }
 324   }
 325
 326   if (fDebug > 0)  printf("AliAnalysisTaskCaloFilter::CheckForPrimaryVertex() - Reject, GetPrimaryVertexTracks()->GetNContributors() > 1 \n");
 327
 328   return kFALSE;
 329
 330 }
 331
 332 //__________________________________________________
 333 void AliAnalysisTaskCaloFilter::CorrectionsInEMCAL()
 334 {
 335   //If EMCAL, and requested, correct energy, position ...
 336
 337   //Need to do this in a separate loop before filling the ESDs because of the track matching recalculations
 338
 339   if(fCorrect && (fCaloFilter==kEMCAL || fCaloFilter==kBoth) )
 340   {
 341     if(!fGeoMatrixSet)
 342     {
 343       if(fLoadEMCALMatrices)
 344       {
 345         printf("AliAnalysisTaskCaloFilter::UserExec() - Load user defined EMCAL geometry matrices\n");
 346         for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 347           if(fEMCALMatrix[mod]){
 348             if(DebugLevel() > 1)
 349               fEMCALMatrix[mod]->Print();
 350             fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
 351           }
 352           fGeoMatrixSet=kTRUE;
 353         }//SM loop
 354       }//Load matrices
 355       else if(!gGeoManager)
 356       {
 357         printf("AliAnalysisTaskCaloFilter::UserExec() - Get geo matrices from data\n");
 358         //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
 359         if(!strcmp(InputEvent()->GetName(),"AliAODEvent"))
 360         {
 361           if(DebugLevel() > 1)
 362             printf("AliAnalysisTaskCaloFilter Use ideal geometry, values geometry matrix not kept in AODs.\n");
 363         }//AOD
 364         else
 365         {
 366           if(DebugLevel() > 1) printf("AliAnalysisTaskCaloFilter Load Misaligned matrices. \n");
 367           AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
 368           if(!esd)
 369           {
 370             printf("AliAnalysisTaskCaloFilter::UserExec() - This event does not contain ESDs?");
 371             return;
 372           }
 373           for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
 374           {
 375             //if(DebugLevel() > 1)
 376             esd->GetEMCALMatrix(mod)->Print();
 377             if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
 378           }
 379           fGeoMatrixSet=kTRUE;
 380         }//ESD
 381       }//Load matrices from Data
 382
 383     }//first event
 384
 385     //Cluster Loop
 386     Int_t nCaloClus = InputEvent()->GetNumberOfCaloClusters();
 387
 388     for (Int_t iClust=0; iClust<nCaloClus; ++iClust)
 389     {
 390       AliVCluster * cluster = InputEvent()->GetCaloCluster(iClust);
 391
 392       if(cluster->IsPHOS()) continue ;
 393
 394       Float_t position[]={0,0,0};
 395       if(DebugLevel() > 2)
 396         printf("Check cluster %d for bad channels and close to border\n",cluster->GetID());
 397       if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
 398
 399       if(DebugLevel() > 2)
 400       {
 401         printf("Filter, before  : i %d, E %f, dispersion %f, m02 %f, m20 %f, distToBad %f\n",iClust,cluster->E(),
 402                cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20(), cluster->GetDistanceToBadChannel());
 403         cluster->GetPosition(position);
 404         printf("Filter, before  : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
 405       }
 406
 407       //Recalculate distance to bad channels, if new list of bad channels provided
 408       fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, InputEvent()->GetEMCALCells(), cluster);
 409
 410       if(fEMCALRecoUtils->IsRecalibrationOn())
 411       {
 412         fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, cluster, InputEvent()->GetEMCALCells());
 413         fEMCALRecoUtils->RecalculateClusterShowerShapeParameters(fEMCALGeo, InputEvent()->GetEMCALCells(),cluster);
 414         fEMCALRecoUtils->RecalculateClusterPID(cluster);
 415       }
 416
 417       fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, InputEvent()->GetEMCALCells(),cluster);
 418
 419       if(DebugLevel() > 2)
 420       {
 421         printf("Filter, after   : i %d, E %f, dispersion %f, m02 %f, m20 %f, distToBad %f\n",cluster->GetID(),cluster->E(),
 422                cluster->GetDispersion(),cluster->GetM02(),cluster->GetM20(), cluster->GetDistanceToBadChannel());
 423         cluster->GetPosition(position);
 424         printf("Filter, after   : i %d, x %f, y %f, z %f\n",cluster->GetID(), position[0], position[1], position[2]);
 425       }
 426
 427       cluster->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(cluster));
 428
 429     }
 430
 431     //Recalculate track-matching
 432     fEMCALRecoUtils->FindMatches(InputEvent(),0,fEMCALGeo);
 433
 434   } // corrections in EMCAL
 435 }
 436
 437 //________________________________________________
 438 void AliAnalysisTaskCaloFilter::FillAODCaloCells()
 439 {
 440   // Fill EMCAL/PHOS cell info
 441
 442   // EMCAL
 443   if ((fCaloFilter==kBoth ||  fCaloFilter==kEMCAL) && fEvent->GetEMCALCells())
 444   { // protection against missing ESD information
 445     AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
 446     Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
 447
 448     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
 449     aodEMcells.CreateContainer(nEMcell);
 450     aodEMcells.SetType(AliVCaloCells::kEMCALCell);
 451     Double_t calibFactor = 1.;
 452     for (Int_t iCell = 0; iCell < nEMcell; iCell++)
 453     {
 454       Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 455       fEMCALGeo->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
 456       fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 457
 458       if(fCorrect && fEMCALRecoUtils->IsRecalibrationOn())
 459       {
 460         calibFactor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 461       }
 462
 463       if(!fEMCALRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
 464       { //Channel is not declared as bad
 465         aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
 466                            eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
 467         //printf("GOOD channel\n");
 468       }
 469       else
 470       {
 471         aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
 472         //printf("BAD channel\n");
 473       }
 474     }
 475     aodEMcells.Sort();
 476   }
 477
 478   // PHOS
 479   if ((fCaloFilter==kBoth ||  fCaloFilter==kPHOS) && fEvent->GetPHOSCells())
 480   { // protection against missing ESD information
 481     AliVCaloCells &eventPHcells = *(fEvent->GetPHOSCells());
 482     Int_t nPHcell = eventPHcells.GetNumberOfCells() ;
 483
 484     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
 485     aodPHcells.CreateContainer(nPHcell);
 486     aodPHcells.SetType(AliVCaloCells::kPHOSCell);
 487
 488     for (Int_t iCell = 0; iCell < nPHcell; iCell++)
 489     {
 490       aodPHcells.SetCell(iCell,eventPHcells.GetCellNumber(iCell),eventPHcells.GetAmplitude(iCell),
 491                          eventPHcells.GetTime(iCell),eventPHcells.GetMCLabel(iCell),eventPHcells.GetEFraction(iCell));
 492     }
 493
 494     aodPHcells.Sort();
 495   }
 496 }
 497
 498
 499 //___________________________________________________
 500 void AliAnalysisTaskCaloFilter::FillAODCaloClusters()
 501 {
 502   // Fill the AOD with caloclusters
 503
 504   // Access to the AOD container of clusters
 505
 506   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
 507   Int_t jClusters=0;
 508   Float_t  posF[3]  ;
 509
 510   Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
 511   for (Int_t iClust=0; iClust<nCaloClus; ++iClust)
 512   {
 513
 514     AliVCluster * cluster = fEvent->GetCaloCluster(iClust);
 515
 516     //Check which calorimeter information we want to keep.
 517
 518     if(fCaloFilter!=kBoth)
 519     {
 520       if     (fCaloFilter==kPHOS  && cluster->IsEMCAL()) continue ;
 521       else if(fCaloFilter==kEMCAL && cluster->IsPHOS())  continue ;
 522     }
 523
 524     // Get original residuals, in case of previous recalculation, reset them
 525     Float_t dR = cluster->GetTrackDx();
 526     Float_t dZ = cluster->GetTrackDz();
 527
 528     if(DebugLevel() > 2)
 529       printf("Original residuals : dZ %f, dR %f\n ",dZ, dR);
 530
 531     //--------------------------------------------------------------
 532     //If EMCAL and corrections done, get the new matching parameters, do not copy noisy clusters
 533     if(cluster->IsEMCAL() && fCorrect)
 534     {
 535       if(DebugLevel() > 2)
 536         printf("Check cluster %d for bad channels and close to border\n",cluster->GetID());
 537
 538       if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo,cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
 539
 540       if(fEMCALRecoUtils->IsExoticCluster(cluster, InputEvent()->GetEMCALCells(),InputEvent()->GetBunchCrossNumber())) continue;
 541
 542       fEMCALRecoUtils->GetMatchedResiduals(cluster->GetID(),dR,dZ);
 543       cluster->SetTrackDistance(dR,dZ);
 544     }
 545
 546     if(DebugLevel() > 2)
 547     {
 548       if(cluster->IsEMCAL()) printf("EMCAL Track-Cluster Residuals : dZ %f, dR %f\n ",dZ, dR);
 549       if(cluster->IsPHOS())  printf("PHOS  Track-Cluster Residuals : dZ %f, dR %f\n ",dZ, dR);
 550     }
 551
 552     //--------------------------------------------------------------
 553
 554     //Now fill AODs
 555
 556     Int_t id       = cluster->GetID();
 557     Float_t energy = cluster->E();
 558     cluster->GetPosition(posF);
 559
 560     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
 561     AliAODCaloCluster(id,
 562                       cluster->GetNLabels(),
 563                       cluster->GetLabels(),
 564                       energy,
 565                       posF,
 566                       NULL,
 567                       cluster->GetType());
 568
 569     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
 570                                 cluster->GetDispersion(),
 571                                 cluster->GetM20(), cluster->GetM02(),
 572                                 -1,
 573                                 cluster->GetNExMax(),cluster->GetTOF()) ;
 574
 575     caloCluster->SetPIDFromESD(cluster->GetPID());
 576     caloCluster->SetNCells(cluster->GetNCells());
 577     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
 578     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
 579     caloCluster->SetTrackDistance(dR, dZ);
 580
 581     if(DebugLevel() > 2)
 582     {
 583       printf("Filter, aod     : i %d, E %f, dispersion %f, m02 %f, m20 %f\n",caloCluster->GetID(),caloCluster->E(),
 584              caloCluster->GetDispersion(),caloCluster->GetM02(),caloCluster->GetM20());
 585       caloCluster->GetPosition(posF);
 586       printf("Filter, aod     : i %d, x %f, y %f, z %f\n",caloCluster->GetID(), posF[0], posF[1], posF[2]);
 587     }
 588
 589     //Matched tracks, just to know if there was any match, the track pointer is useless if tracks not stored
 590     if(TMath::Abs(dR) < 990 && TMath::Abs(dZ) < 990)
 591     { //Default value in PHOS 999, in EMCAL 1024, why?
 592       caloCluster->AddTrackMatched(new AliAODTrack);
 593     }
 594     // TO DO, in case Tracks available, think how to put the matched track in AOD
 595   }
 596
 597   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots
 598
 599 }
 600
 601 //__________________________________________________
 602 void AliAnalysisTaskCaloFilter::FillAODCaloTrigger()
 603 {
 604   // AOD CaloTrigger copy
 605
 606   if( !AODEvent() || !fAODEvent ) return;
 607
 608   AliAODCaloTrigger* triggerEM   = AODEvent()->GetCaloTrigger("EMCAL");
 609   AliAODCaloTrigger* triggerPH   = AODEvent()->GetCaloTrigger("PHOS");
 610
 611   // Copy from AODs
 612
 613   AliAODCaloTrigger* inTriggerEM = fAODEvent ->GetCaloTrigger("EMCAL");
 614   AliAODCaloTrigger* inTriggerPH = fAODEvent ->GetCaloTrigger("PHOS");
 615
 616   if(inTriggerPH && (fCaloFilter==kBoth || fCaloFilter==kPHOS))  *triggerPH = *inTriggerPH;
 617
 618   if(inTriggerEM && (fCaloFilter==kBoth || fCaloFilter==kEMCAL)) *triggerEM = *inTriggerEM;
 619 }
 620
 621 //______________________________________________
 622 void AliAnalysisTaskCaloFilter::FillAODHeader()
 623 {
 624   // AOD header copy
 625
 626   AliAODHeader* header = dynamic_cast<AliAODHeader*>(AODEvent()->GetHeader());
 627   if(!header) AliFatal("Not a standard AOD");
 628
 629   // Copy from AODs
 630   if(fAODEvent)
 631   {
 632     *header = *((AliAODHeader*)fAODEvent->GetHeader());
 633     return;
 634   }
 635
 636   if(!fESDEvent) return;
 637
 638   // Copy from ESDs
 639
 640   header->SetRunNumber(fEvent->GetRunNumber());
 641
 642   TTree* tree = fInputHandler->GetTree();
 643   if (tree)
 644   {
 645     TFile* file = tree->GetCurrentFile();
 646     if (file) header->SetESDFileName(file->GetName());
 647   }
 648
 649   header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
 650   header->SetOrbitNumber(fEvent->GetOrbitNumber());
 651   header->SetPeriodNumber(fEvent->GetPeriodNumber());
 652   header->SetEventType(fEvent->GetEventType());
 653
 654   //Centrality
 655   if(fEvent->GetCentrality())
 656   {
 657     header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
 658   }
 659   else
 660   {
 661     header->SetCentrality(0);
 662   }
 663
 664   //Trigger
 665   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
 666   header->SetFiredTriggerClasses(fESDEvent->GetFiredTriggerClasses());
 667   header->SetTriggerMask(fEvent->GetTriggerMask());
 668   header->SetTriggerCluster(fEvent->GetTriggerCluster());
 669   header->SetL0TriggerInputs(fESDEvent->GetHeader()->GetL0TriggerInputs());
 670   header->SetL1TriggerInputs(fESDEvent->GetHeader()->GetL1TriggerInputs());
 671   header->SetL2TriggerInputs(fESDEvent->GetHeader()->GetL2TriggerInputs());
 672
 673   header->SetMagneticField(fEvent->GetMagneticField());
 674   header->SetMuonMagFieldScale(fESDEvent->GetCurrentDip()/6000.);
 675
 676   header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 677   header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 678   header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 679   header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 680   header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
 681
 682   Float_t diamxy[2]={(Float_t)fEvent->GetDiamondX(),(Float_t)fEvent->GetDiamondY()};
 683   Float_t diamcov[3];
 684   fEvent->GetDiamondCovXY(diamcov);
 685   header->SetDiamond(diamxy,diamcov);
 686   header->SetDiamondZ(fESDEvent->GetDiamondZ(),fESDEvent->GetSigma2DiamondZ());
 687
 688 }
 689
 690
 691 //__________________________________________________
 692 void AliAnalysisTaskCaloFilter::FillAODMCParticles()
 693 {
 694   // Copy MC particles
 695
 696   if(!fFillMCParticles) return;
 697
 698   TClonesArray* inMCParticles = (TClonesArray*) (fAODEvent  ->FindListObject("mcparticles"));
 699   TClonesArray* ouMCParticles = (TClonesArray*) ( AODEvent()->FindListObject("mcparticles"));
 700
 701   if( inMCParticles &&  ouMCParticles ) new (ouMCParticles) TClonesArray(*inMCParticles);
 702
 703 }
 704
 705 //_____________________________________________
 706 void AliAnalysisTaskCaloFilter::FillAODTracks()
 707 {
 708   // AOD track copy
 709
 710   if(!fFillTracks) return;
 711
 712   AliAODTrack* aodTrack(0x0);
 713
 714   Double_t pos[3]   = { 0. };
 715   Double_t covTr[21]= { 0. };
 716   //Double_t pid[10]  = { 0. };
 717   Double_t p[3]     = { 0. };
 718
 719   // Copy from AODs
 720   if(fAODEvent)
 721   {
 722     //TClonesArray* inTracks = fAODEvent  ->GetTracks();
 723     TClonesArray* ouTracks = AODEvent()->GetTracks();
 724                 //new (ouTracks) TClonesArray(*inTracks);
 725
 726     //printf("N tracks %d\n",fAODEvent->GetNumberOfTracks());
 727     Int_t nCopyTrack = 0;
 728     for (Int_t nTrack = 0; nTrack < fAODEvent->GetNumberOfTracks(); ++nTrack)
 729     {
 730       AliAODTrack *track = dynamic_cast<AliAODTrack*>(fAODEvent->GetTrack(nTrack));
 731       if(!track) AliFatal("Not a standard AOD");
 732
 733       // Select only hybrid tracks?
 734       if(fFillHybridTracks && !track->IsHybridGlobalConstrainedGlobal()) continue;
 735
 736       // Remove PID object to save space
 737       //track->SetDetPID(0x0);
 738
 739       //new((*ouTracks)[nCopyTrack++]) AliAODTrack(*track);
 740
 741       track->GetPxPyPz(p);
 742       Bool_t isDCA  = track->GetPosition(pos);
 743       track->GetCovMatrix(covTr);
 744       //track->GetPID(pid);
 745
 746       AliAODVertex* primVertex = (AliAODVertex*) AODEvent()->GetVertices()->At(0); // primary vertex, copied previously!!!
 747
 748       aodTrack = new((*ouTracks)[nCopyTrack++]) AliAODTrack(
 749                                                             track->GetID(),
 750                                                             track->GetLabel(),
 751                                                             p,
 752                                                             kTRUE,
 753                                                             pos,
 754                                                             isDCA,
 755                                                             covTr,
 756                                                             track->Charge(),
 757                                                             track->GetITSClusterMap(),
 758                                                             // pid,
 759                                                             primVertex,
 760                                                             track->GetUsedForVtxFit(),
 761                                                             track->GetUsedForPrimVtxFit(),
 762                                                             (AliAODTrack::AODTrk_t) track->GetType(),
 763                                                             track->GetFilterMap());
 764
 765
 766       aodTrack->SetPIDForTracking(track->GetPIDForTracking());
 767       aodTrack->SetIsHybridGlobalConstrainedGlobal(track->IsHybridGlobalConstrainedGlobal());
 768       aodTrack->SetIsHybridTPCConstrainedGlobal   (track->IsHybridTPCConstrainedGlobal());
 769       aodTrack->SetIsGlobalConstrained            (track->IsGlobalConstrained());
 770       aodTrack->SetIsTPCConstrained               (track->IsTPCConstrained());
 771
 772       aodTrack->SetTPCFitMap    (track->GetTPCFitMap());
 773       aodTrack->SetTPCClusterMap(track->GetTPCClusterMap());
 774       aodTrack->SetTPCSharedMap (track->GetTPCSharedMap());
 775
 776       aodTrack->SetChi2MatchTrigger(track->GetChi2MatchTrigger());
 777
 778       // set the DCA values to the AOD track
 779
 780       aodTrack->SetPxPyPzAtDCA(track->PxAtDCA(),track->PyAtDCA(),track->PzAtDCA());
 781       aodTrack->SetXYAtDCA    (track->XAtDCA() ,track->YAtDCA());
 782
 783       aodTrack->SetFlags      (track->GetFlags());
 784       aodTrack->SetTPCPointsF (track->GetTPCNclsF());
 785
 786       // Calo
 787
 788       if(track->IsEMCAL()) aodTrack->SetEMCALcluster(track->GetEMCALcluster());
 789       if(track->IsPHOS())  aodTrack->SetPHOScluster (track->GetPHOScluster());
 790       aodTrack->SetTrackPhiEtaPtOnEMCal( track->GetTrackPhiOnEMCal(),  track->GetTrackPhiOnEMCal(),  track->GetTrackPtOnEMCal() );
 791
 792     }
 793
 794     //printf("Final N tracks %d\n",nCopyTrack);
 795
 796     return;
 797   }
 798
 799 }
 800
 801 //_________________________________________
 802 void AliAnalysisTaskCaloFilter::FillAODv0s()
 803 {
 804   // Copy v0s (use if you know what you do, use quite a lot of memory)
 805
 806   if(!fFillv0s) return;
 807
 808   // Copy from AODs
 809   if(fAODEvent)
 810   {
 811     TClonesArray* inv0 = fAODEvent  ->GetV0s();
 812     TClonesArray* ouv0 =  AODEvent()->GetV0s();
 813
 814                 //new (ouv0s) TClonesArray(*inv0s);
 815
 816     Int_t allv0s = inv0->GetEntriesFast();
 817
 818     for (Int_t nv0s = 0; nv0s < allv0s; ++nv0s)
 819     {
 820       AliAODv0 *v0 = (AliAODv0*)inv0->At(nv0s);
 821
 822       new((*ouv0)[nv0s]) AliAODv0(*v0);
 823     }
 824
 825     return;
 826   }
 827
 828 }
 829
 830 //____________________________________________
 831 void AliAnalysisTaskCaloFilter::FillAODVZERO()
 832 {
 833   // Copy VZERO
 834
 835   if(!fFillVZERO) return;
 836
 837   AliAODVZERO* vzeroData = AODEvent()->GetVZEROData();
 838
 839   if(fESDEvent) *vzeroData = *(fESDEvent->GetVZEROData());
 840   else          *vzeroData = *(fAODEvent->GetVZEROData());
 841
 842 }
 843
 844 //_______________________________________________
 845 void AliAnalysisTaskCaloFilter::FillAODVertices()
 846 {
 847   // Copy vertices
 848
 849   // set arrays and pointers
 850   Double_t pos[3]   ;
 851   Double_t covVtx[6];
 852   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 853
 854   // Copy from AODs
 855   if(fAODEvent)
 856   {
 857     TClonesArray* inVertices = fAODEvent  ->GetVertices();
 858     TClonesArray* ouVertices = AODEvent()->GetVertices();
 859
 860                 //new (ouVertices) TClonesArray(*inVertices);
 861
 862     //Keep only the first 3 vertices if not requested
 863     Int_t allVertices = inVertices->GetEntriesFast();
 864
 865     //printf("n Vertices %d\n",allVertices);
 866
 867     if(!fFillAllVertices)
 868     {
 869       if(allVertices >  3) allVertices = 3;
 870     }
 871
 872     //printf("Final n Vertices %d\n",allVertices);
 873
 874     for (Int_t nVertices = 0; nVertices < allVertices; ++nVertices)
 875     {
 876       AliAODVertex *vertex = (AliAODVertex*)inVertices->At(nVertices);
 877
 878       new((*ouVertices)[nVertices]) AliAODVertex(*vertex);
 879     }
 880
 881     return;
 882   }
 883
 884   if(!fESDEvent) return;
 885
 886   // Copy from ESDs
 887
 888   // Access to the AOD container of vertices
 889   Int_t jVertices=0;
 890   TClonesArray &vertices = *(AODEvent()->GetVertices());
 891
 892   // Add primary vertex. The primary tracks will be defined
 893   // after the loops on the composite objects (v0, cascades, kinks)
 894   fEvent   ->GetPrimaryVertex()->GetXYZ(pos);
 895   fESDEvent->GetPrimaryVertex()->GetCovMatrix(covVtx);
 896   Float_t chi = fESDEvent->GetPrimaryVertex()->GetChi2toNDF();
 897
 898   AliAODVertex * primary = new(vertices[jVertices++])
 899   AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
 900   primary->SetName(fEvent->GetPrimaryVertex()->GetName());
 901   primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
 902
 903 }
 904
 905 //____________________________________
 906 void AliAnalysisTaskCaloFilter::Init()
 907 {
 908   //Init analysis with configuration macro if available
 909
 910   if(gROOT->LoadMacro(fConfigName) >=0)
 911   {
 912     printf("Configure analysis with %s\n",fConfigName.Data());
 913
 914     AliAnalysisTaskCaloFilter *filter = (AliAnalysisTaskCaloFilter*)gInterpreter->ProcessLine("ConfigCaloFilter()");
 915
 916     fEMCALGeoName      = filter->fEMCALGeoName;
 917     fLoadEMCALMatrices = filter->fLoadEMCALMatrices;
 918     fFillAODFile       = filter->fFillAODFile;
 919     fFillTracks        = filter->fFillTracks;
 920     fFillHybridTracks  = filter->fFillHybridTracks;
 921     fFillv0s           = filter->fFillv0s;
 922     fFillVZERO         = filter->fFillVZERO;
 923     fFillAllVertices   = filter->fFillAllVertices;
 924     fEMCALRecoUtils    = filter->fEMCALRecoUtils;
 925     fConfigName        = filter->fConfigName;
 926     fCaloFilter        = filter->fCaloFilter;
 927     fEventSelection[0] = filter->fEventSelection[0];
 928     fEventSelection[1] = filter->fEventSelection[1];
 929     fEventSelection[2] = filter->fEventSelection[2];
 930     fAcceptAllMBEvent  = filter->fAcceptAllMBEvent;
 931     fMBTriggerMask     = filter->fMBTriggerMask;
 932     fCorrect           = filter->fCorrect;
 933     fEMCALEnergyCut    = filter->fEMCALEnergyCut;
 934     fEMCALNcellsCut    = filter->fEMCALNcellsCut;
 935     fPHOSEnergyCut     = filter->fPHOSEnergyCut;
 936     fPHOSNcellsCut     = filter->fPHOSNcellsCut;
 937     fTrackPtCut        = filter->fTrackPtCut;
 938     fVzCut             = filter->fVzCut;
 939     fCheckEventVertex  = filter->fCheckEventVertex;
 940
 941     for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = filter->fEMCALMatrix[i] ;
 942   }
 943 }
 944
 945 //_________________________________________
 946 void AliAnalysisTaskCaloFilter::PrintInfo()
 947 {
 948   //Print settings
 949
 950   printf("AnalysisCaloFilter::PrintInfo() \n");
 951
 952   printf("\t Not only filter, correct Clusters? %d\n",fCorrect);
 953   printf("\t Calorimeter Filtering Option     ? %d\n",fCaloFilter);
 954
 955   //printf("\t Use handmade geo matrices?   EMCAL %d, PHOS %d\n",fLoadEMCALMatrices, fLoadPHOSMatrices);
 956   printf("\t Use handmade geo matrices?   EMCAL %d, PHOS 0\n",fLoadEMCALMatrices);
 957
 958   printf("\t Fill: AOD file? %d Tracks? %d; all Vertex? %d; v0s? %d; VZERO ? %d\n",
 959          fFillAODFile,fFillTracks,fFillAllVertices, fFillv0s, fFillVZERO);
 960
 961   printf("\t Event Selection based : EMCAL?  %d, PHOS?  %d Tracks? %d - Accept all MB with mask %d? %d\n",
 962          fEventSelection[0],fEventSelection[1],fEventSelection[2],fMBTriggerMask, fAcceptAllMBEvent);
 963
 964   printf("\t \t EMCAL E > %2.2f, EMCAL nCells >= %d, PHOS E > %2.2f, PHOS nCells >= %d, Track pT > %2.2f, |vz| < %2.2f\n",
 965          fEMCALEnergyCut,fEMCALNcellsCut,fPHOSEnergyCut,fPHOSNcellsCut, fTrackPtCut,fVzCut);
 966 }
 967
 968 //_______________________________________________________
 969 void AliAnalysisTaskCaloFilter::UserCreateOutputObjects()
 970 {
 971   // Init geometry
 972
 973   fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
 974
 975   if(fFillMCParticles)
 976   {
 977     TClonesArray * aodMCParticles = new TClonesArray("AliAODMCParticle",500);
 978                 aodMCParticles->SetName("mcparticles");
 979                 ((AliAODHandler*)AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler())->AddBranch("TClonesArray", &aodMCParticles);
 980   }
 981
 982 }
 983
 984 //____________________________________________________________
 985 void AliAnalysisTaskCaloFilter::UserExec(Option_t */*option*/)
 986 {
 987   // Execute analysis for current event
 988   // Copy input ESD or AOD header, vertex, CaloClusters and CaloCells to output AOD
 989
 990   if (fDebug > 0)
 991     printf("CaloFilter: Analysing event # %d\n", (Int_t)Entry());
 992
 993   fEvent    = InputEvent();
 994   fAODEvent = dynamic_cast<AliAODEvent*> (fEvent);
 995   fESDEvent = dynamic_cast<AliESDEvent*> (fEvent);
 996
 997   if(!fEvent)
 998   {
 999     printf("AliAnalysisTaskCaloFilter::UserExec - This event does not contain Input?");
1000     return;
1001   }
1002
1003   // printf("Start processing : %s\n",fAODEvent->GetFiredTriggerClasses().Data());
1004
1005   // Select the event
1006
1007   if(!AcceptEvent()) return ;
1008
1009   //Magic line to write events to file
1010
1011   AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
1012
1013   // Reset output AOD
1014
1015   Int_t nVertices = 0;
1016   if(fFillv0s) nVertices = fEvent->GetNumberOfV0s();
1017   Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
1018   Int_t nTracks   = fEvent->GetNumberOfTracks();
1019
1020   AODEvent()->ResetStd(nTracks, nVertices, 0, 0, 0, nCaloClus, 0, 0);
1021
1022   // Copy
1023
1024   FillAODHeader();
1025
1026   //
1027   FillAODv0s();
1028
1029   //
1030   FillAODVertices(); // Do it before the track filtering to have the reference to the vertex
1031
1032   //
1033   FillAODVZERO();
1034
1035   //
1036   FillAODTracks();
1037
1038   //
1039   CorrectionsInEMCAL();
1040
1041   //
1042   FillAODCaloClusters();
1043
1044   //
1045   FillAODCaloCells();
1046
1047   //
1048   FillAODCaloTrigger();
1049
1050   //
1051   FillAODMCParticles();
1052
1053   //printf("Filtered event, end processing : %s\n",fAODEvent->GetFiredTriggerClasses().Data());
1054
1055 }
1056