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