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