PWG4/CaloCalib/AliAnalysisTaskPHOSPi0CalibSelection.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Boris Polishchuk                                               *
   5  * Adapted to AOD reading by Gustavo Conesa  *
   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 //                                                                           //
  18 // Fill histograms (one per cell) with two-cluster invariant mass            //
  19 // using calibration coefficients of the previous iteration.                 //
  20 // Histogram for a given cell is filled if the most energy of one cluster    //
  21 // is deposited in this cell and the other cluster could be anywhere in PHOS.//
  22 //                                                                           //
  23 //---------------------------------------------------------------------------//
  24
  25 #include <cstdlib>
  26
  27 // Root
  28 #include "TLorentzVector.h"
  29 #include "TVector3.h"
  30 #include "TRefArray.h"
  31 #include "TList.h"
  32
  33 // AliRoot
  34 #include "AliAnalysisTaskPHOSPi0CalibSelection.h"
  35 #include "AliAODEvent.h"
  36 #include "AliESDEvent.h"
  37 #include "AliESDCaloCluster.h"
  38 #include "AliESDCaloCells.h"
  39 #include "AliPHOSAodCluster.h"
  40 #include "AliPHOSPIDv1.h"
  41
  42 ClassImp(AliAnalysisTaskPHOSPi0CalibSelection)
  43
  44 //__________________________________________________
  45 AliAnalysisTaskPHOSPi0CalibSelection::AliAnalysisTaskPHOSPi0CalibSelection() :
  46 AliAnalysisTaskSE(),fOutputContainer(0x0),fPhosGeo(0x0),fCalibData(0x0),fHmgg(0x0),
  47 fEmin(0.), fLogWeight(4.5), fCopyAOD(kFALSE)
  48 {
  49   //Default constructor.
  50
  51   for(Int_t iMod=0; iMod<5; iMod++) {
  52     for(Int_t iX=0; iX<64; iX++) {
  53       for(Int_t iZ=0; iZ<56; iZ++) {
  54         fHmpi0[iMod][iX][iZ]=0;
  55       }
  56     }
  57   }
  58
  59 }
  60
  61 //__________________________________________________
  62 AliAnalysisTaskPHOSPi0CalibSelection::AliAnalysisTaskPHOSPi0CalibSelection(const char* name) :
  63   AliAnalysisTaskSE(name),fOutputContainer(0x0),fPhosGeo(0x0),fCalibData(0x0),fHmgg(0x0),
  64   fEmin(0.), fLogWeight(4.5), fCopyAOD(kFALSE)
  65 {
  66   //Named constructor which should be used.
  67
  68   DefineOutput(1,TList::Class());
  69
  70   for(Int_t iMod=0; iMod<5; iMod++) {
  71     for(Int_t iX=0; iX<64; iX++) {
  72       for(Int_t iZ=0; iZ<56; iZ++) {
  73         fHmpi0[iMod][iX][iZ]=0;
  74       }
  75     }
  76   }
  77
  78 }
  79
  80 //__________________________________________________
  81 AliAnalysisTaskPHOSPi0CalibSelection::~AliAnalysisTaskPHOSPi0CalibSelection()
  82 {
  83   //Destructor.
  84
  85   if(fOutputContainer){
  86     fOutputContainer->Delete() ;
  87     delete fOutputContainer ;
  88   }
  89
  90   if(fCalibData) delete fCalibData;
  91   if(fPhosGeo)   delete fPhosGeo;
  92
  93 }
  94
  95 //__________________________________________________
  96 void AliAnalysisTaskPHOSPi0CalibSelection::CreateAODFromAOD()
  97 {
  98   // Copy AOD header, vertex, CaloClusters and CaloCells to output AOD
  99
 100   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
 101
 102   // set arrays and pointers
 103   Float_t posF[3];
 104   Double_t pos[3];
 105
 106   Double_t covVtx[6];
 107
 108   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 109
 110   // Update the header
 111   AliAODHeader*headerin = aod->GetHeader();
 112   AliAODHeader* header = AODEvent()->GetHeader();
 113   header->SetRunNumber(headerin->GetRunNumber());
 114   header->SetBunchCrossNumber(headerin->GetBunchCrossNumber());
 115   header->SetOrbitNumber(headerin->GetOrbitNumber());
 116   header->SetPeriodNumber(headerin->GetPeriodNumber());
 117   header->SetEventType(headerin->GetEventType());
 118   header->SetMuonMagFieldScale(headerin->GetMuonMagFieldScale());
 119   header->SetCentrality(headerin->GetCentrality());
 120
 121   header->SetTriggerMask(headerin->GetTriggerMask());
 122   header->SetTriggerCluster(headerin->GetTriggerCluster());
 123   header->SetMagneticField(headerin->GetMagneticField());
 124   header->SetZDCN1Energy(headerin->GetZDCN1Energy());
 125   header->SetZDCP1Energy(headerin->GetZDCP1Energy());
 126   header->SetZDCN2Energy(headerin->GetZDCN2Energy());
 127   header->SetZDCP2Energy(headerin->GetZDCP2Energy());
 128   header->SetZDCEMEnergy(headerin->GetZDCEMEnergy(0),headerin->GetZDCEMEnergy(1));
 129   Float_t diamxy[2]={aod->GetDiamondX(),aod->GetDiamondY()};
 130   Float_t diamcov[3]; aod->GetDiamondCovXY(diamcov);
 131   header->SetDiamond(diamxy,diamcov);
 132   //
 133   //
 134   Int_t nVertices = 1 ;/* = prim. vtx*/;
 135   Int_t nCaloClus = aod->GetNCaloClusters();
 136
 137   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
 138
 139   // Access to the AOD container of vertices
 140   TClonesArray &vertices = *(AODEvent()->GetVertices());
 141   Int_t jVertices=0;
 142
 143   // Add primary vertex. The primary tracks will be defined
 144   // after the loops on the composite objects (V0, cascades, kinks)
 145   const AliAODVertex *vtx = aod->GetPrimaryVertex();
 146
 147   vtx->GetXYZ(pos); // position
 148   vtx->GetCovMatrix(covVtx); //covariance matrix
 149
 150   AliAODVertex * primary = new(vertices[jVertices++])
 151     AliAODVertex(pos, covVtx, vtx->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
 152   primary->SetName(vtx->GetName());
 153   primary->SetTitle(vtx->GetTitle());
 154
 155   // Access to the AOD container of clusters
 156   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
 157   Int_t jClusters=0;
 158
 159   for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
 160
 161     AliAODCaloCluster * cluster = aod->GetCaloCluster(iClust);
 162
 163     //Check if it is a PHOS cluster
 164     if(!cluster->IsPHOSCluster())  continue ;
 165
 166     Int_t id       = cluster->GetID();
 167     Float_t energy = cluster->E();
 168     cluster->GetPosition(posF);
 169     Char_t ttype   = cluster->GetType();
 170
 171     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
 172       AliAODCaloCluster(id,
 173                         0,
 174                         0x0,
 175                         energy,
 176                         posF,
 177                         NULL,
 178                         ttype);
 179
 180     caloCluster->SetCaloCluster(cluster->GetDistToBadChannel(),
 181                                 cluster->GetDispersion(),
 182                                 cluster->GetM20(), cluster->GetM02(),
 183                                 cluster->GetEmcCpvDistance(),
 184                                 cluster->GetNExMax(),cluster->GetTOF()) ;
 185
 186     caloCluster->SetPIDFromESD(cluster->PID());
 187     caloCluster->SetNCells(cluster->GetNCells());
 188     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
 189     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
 190
 191   }
 192
 193   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
 194   // end of loop on calo clusters
 195
 196   // fill PHOS cell info
 197   if (aod->GetPHOSCells()) { // protection against missing AOD information
 198     AliAODCaloCells &aodinPHcells = *(aod->GetPHOSCells());
 199     Int_t nPHcell = aodinPHcells.GetNumberOfCells() ;
 200
 201     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
 202     aodPHcells.CreateContainer(nPHcell);
 203     aodPHcells.SetType(AliAODCaloCells::kPHOS);
 204     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
 205       aodPHcells.SetCell(iCell,aodinPHcells.GetCellNumber(iCell),aodinPHcells.GetAmplitude(iCell));
 206     }
 207     aodPHcells.Sort();
 208   }
 209
 210   return;
 211 }
 212
 213 //__________________________________________________
 214 void AliAnalysisTaskPHOSPi0CalibSelection::CreateAODFromESD()
 215 {
 216
 217   // Copy Header, Vertex, CaloClusters and CaloCells from ESDs to AODs
 218
 219   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 220
 221   // set arrays and pointers
 222   Float_t posF[3];
 223   Double_t pos[3];
 224
 225   Double_t covVtx[6];
 226
 227   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 228
 229   // Update the header
 230
 231   AliAODHeader* header = AODEvent()->GetHeader();
 232   header->SetRunNumber(esd->GetRunNumber());
 233   header->SetBunchCrossNumber(esd->GetBunchCrossNumber());
 234   header->SetOrbitNumber(esd->GetOrbitNumber());
 235   header->SetPeriodNumber(esd->GetPeriodNumber());
 236   header->SetEventType(esd->GetEventType());
 237   header->SetMuonMagFieldScale(-999.); // FIXME
 238   header->SetCentrality(-999.);        // FIXME
 239
 240
 241   header->SetTriggerMask(esd->GetTriggerMask());
 242   header->SetTriggerCluster(esd->GetTriggerCluster());
 243   header->SetMagneticField(esd->GetMagneticField());
 244   header->SetZDCN1Energy(esd->GetZDCN1Energy());
 245   header->SetZDCP1Energy(esd->GetZDCP1Energy());
 246   header->SetZDCN2Energy(esd->GetZDCN2Energy());
 247   header->SetZDCP2Energy(esd->GetZDCP2Energy());
 248   header->SetZDCEMEnergy(esd->GetZDCEMEnergy(0),esd->GetZDCEMEnergy(1));
 249   Float_t diamxy[2]={esd->GetDiamondX(),esd->GetDiamondY()};
 250   Float_t diamcov[3]; esd->GetDiamondCovXY(diamcov);
 251   header->SetDiamond(diamxy,diamcov);
 252   //
 253   //
 254   Int_t nVertices = 1 ;/* = prim. vtx*/;
 255   Int_t nCaloClus = esd->GetNumberOfCaloClusters();
 256
 257   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
 258
 259   // Access to the AOD container of vertices
 260   TClonesArray &vertices = *(AODEvent()->GetVertices());
 261   Int_t jVertices=0;
 262
 263   // Add primary vertex. The primary tracks will be defined
 264   // after the loops on the composite objects (V0, cascades, kinks)
 265   const AliESDVertex *vtx = esd->GetPrimaryVertex();
 266
 267   vtx->GetXYZ(pos); // position
 268   vtx->GetCovMatrix(covVtx); //covariance matrix
 269
 270   AliAODVertex * primary = new(vertices[jVertices++])
 271     AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
 272   primary->SetName(vtx->GetName());
 273   primary->SetTitle(vtx->GetTitle());
 274
 275   // Access to the AOD container of clusters
 276   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
 277   Int_t jClusters=0;
 278
 279   for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
 280
 281     AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
 282
 283     //Check which calorimeter information we want to keep.
 284     if(!cluster->IsPHOS())  continue ;
 285
 286     Int_t id       = cluster->GetID();
 287     Float_t energy = cluster->E();
 288     cluster->GetPosition(posF);
 289
 290     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++])
 291       AliAODCaloCluster(id,
 292                         0,
 293                         0x0,
 294                         energy,
 295                         posF,
 296                         NULL,
 297                         AliAODCluster::kPHOSNeutral);
 298
 299     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
 300                                 cluster->GetClusterDisp(),
 301                                 cluster->GetM20(), cluster->GetM02(),
 302                                 cluster->GetEmcCpvDistance(),
 303                                 cluster->GetNExMax(),cluster->GetTOF()) ;
 304
 305     caloCluster->SetPIDFromESD(cluster->GetPid());
 306     caloCluster->SetNCells(cluster->GetNCells());
 307     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
 308     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
 309
 310   }
 311   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters
 312   // end of loop on calo clusters
 313
 314   // fill PHOS cell info
 315   if( esd->GetPHOSCells()) { // protection against missing ESD information
 316     AliESDCaloCells &esdPHcells = *(esd->GetPHOSCells());
 317     Int_t nPHcell = esdPHcells.GetNumberOfCells() ;
 318
 319     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
 320     aodPHcells.CreateContainer(nPHcell);
 321     aodPHcells.SetType(AliAODCaloCells::kPHOS);
 322     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {
 323       aodPHcells.SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell));
 324     }
 325     aodPHcells.Sort();
 326   }
 327
 328 }
 329
 330 //__________________________________________________
 331 void AliAnalysisTaskPHOSPi0CalibSelection::UserCreateOutputObjects()
 332 {
 333   //Create output container
 334   fOutputContainer = new TList();
 335
 336   char hname[128], htitl[128];
 337
 338   for(Int_t iMod=0; iMod<5; iMod++) {
 339     for(Int_t iX=0; iX<64; iX++) {
 340       for(Int_t iZ=0; iZ<56; iZ++) {
 341         sprintf(hname,"%d_%d_%d",iMod,iX,iZ);
 342         sprintf(htitl,"Two-gamma inv. mass for mod %d, cell (%d,%d)",iMod,iX,iZ);
 343         fHmpi0[iMod][iX][iZ] = new TH1F(hname,htitl,100,0.,300.);
 344         fOutputContainer->Add(fHmpi0[iMod][iX][iZ]);
 345       }
 346     }
 347   }
 348
 349   fHmgg = new TH1F("hmgg","2-cluster invariant mass",100,0.,300.);
 350   fOutputContainer->Add(fHmgg);
 351
 352   fCalibData = new AliPHOSCalibData();
 353   fPhosGeo =  AliPHOSGeometry::GetInstance("IHEP") ;
 354
 355 }
 356
 357 //__________________________________________________
 358 void AliAnalysisTaskPHOSPi0CalibSelection::UserExec(Option_t* /* option */)
 359 {
 360   //Analysis per event.
 361   if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection <<< Event %d >>>\n",(Int_t)Entry());
 362
 363   AliAODEvent* aod = 0x0;
 364   if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
 365     //Input are ESDs
 366     aod = dynamic_cast<AliAODEvent*>(InputEvent());
 367     // Create new AOD with only CaloClusters and CaloCells
 368     if(fCopyAOD) CreateAODFromAOD();
 369   }
 370   else  if(!strcmp(InputEvent()->GetName(),"AliESDEvent")) {
 371     //Input are ESDs
 372     aod = AODEvent();
 373     // Create AOD with CaloClusters and use it as input.
 374     // If filtering task is already executed, this is not needed.
 375     if(fCopyAOD) CreateAODFromESD();
 376   }
 377   else {
 378     printf("AliAnalysisTaskPHOSPi0CalibSelection: Unknown event type, STOP!\n");
 379     abort();
 380   }
 381
 382   Double_t v[] = {aod->GetVertex(0)->GetX(),aod->GetVertex(0)->GetY(),aod->GetVertex(0)->GetZ()}; //to check!!
 383   //aod->GetVertex()->GetXYZ(v) ;
 384   TVector3 vtx(v); //Check
 385
 386   if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Vertex: (%.3f,%.3f,%.3f)\n",vtx.X(),vtx.Y(),vtx.Z());
 387
 388   Int_t runNum = aod->GetRunNumber();
 389   if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
 390
 391   //Get the matrix with geometry information
 392   //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
 393   if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
 394           if(DebugLevel() > 1)
 395                   printf("AliAnalysisTaskPHOSPi0CalibSelection Use ideal geometry, values geometry matrix not kept in AODs.\n");
 396   }
 397   else{
 398           if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Load Misaligned matrices. \n");
 399           AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
 400           for(Int_t mod=0; mod<5; mod++){
 401                 if(esd->GetPHOSMatrix(mod))
 402                         fPhosGeo->SetMisalMatrix(esd->GetPHOSMatrix(mod),mod) ;
 403                 }
 404   }
 405
 406   if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
 407
 408   AliPHOSPIDv1 pid;
 409
 410   Int_t mod1 = -1;
 411   TLorentzVector p1;
 412   TLorentzVector p2;
 413   TLorentzVector p12;
 414   Int_t relid[4] ;
 415   Int_t maxId;
 416
 417   TRefArray * caloClustersArr  = new TRefArray();
 418   aod->GetPHOSClusters(caloClustersArr);
 419
 420   const Int_t kNumberOfPhosClusters   = caloClustersArr->GetEntries() ;
 421   if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection CaloClusters: %d\n", kNumberOfPhosClusters);
 422
 423   // PHOS cells
 424   AliAODCaloCells *phsCells = aod->GetPHOSCells();
 425
 426   // loop over PHOS clusters
 427   for(Int_t iClu=0; iClu<kNumberOfPhosClusters; iClu++) {
 428
 429     AliAODCaloCluster *c1 = (AliAODCaloCluster *) caloClustersArr->At(iClu);
 430     if(!c1->IsPHOSCluster()) continue; // EMCAL cluster!
 431
 432     Float_t e1i = c1->E();   // cluster energy before correction
 433     if(e1i<fEmin) continue;
 434
 435     AliPHOSAodCluster clu1(*c1);
 436     clu1.Recalibrate(fCalibData, phsCells);
 437     clu1.EvalAll(fLogWeight,vtx);
 438     clu1.EnergyCorrection(&pid) ;
 439
 440     clu1.GetMomentum(p1,v);
 441
 442     MaxEnergyCellPos(phsCells,&clu1,maxId);
 443     fPhosGeo->AbsToRelNumbering(maxId, relid);
 444
 445     mod1 = relid[0]-1;        // module
 446     Int_t iX = relid[2]-1;    // cluster X-coord
 447     Int_t iZ = relid[3]-1;    // cluster Z-coord
 448         Float_t e1ii = clu1.E(); // cluster energy after correction
 449
 450     for (Int_t jClu=iClu; jClu<kNumberOfPhosClusters; jClu++) {
 451       AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
 452       if(!c2->IsPHOSCluster())   continue; // EMCAL cluster!
 453       if(c2->IsEqual(c1)) continue;
 454
 455       Float_t e2i = c2->E();
 456       if(e2i<fEmin) continue;
 457
 458       AliPHOSAodCluster clu2(*c2);
 459       clu2.Recalibrate(fCalibData, phsCells);
 460       clu2.EvalAll(fLogWeight,vtx);
 461       clu2.EnergyCorrection(&pid) ;
 462
 463       clu2.GetMomentum(p2,v);
 464           Float_t e2ii = clu2.E();
 465
 466       p12 = p1+p2;
 467
 468       fHmgg->Fill(p12.M()*1000);
 469       fHmpi0[mod1][iX][iZ]->Fill(p12.M()*1000);
 470
 471           if(DebugLevel() > 1) printf("AliAnalysisTaskPHOSPi0CalibSelection Mass in (mod%d,%d,%d): %.3f GeV  E1_i=%f E1_ii=%f  E2_i=%f E2_ii=%f\n",
 472              mod1,iX,iZ,p12.M(),e1i,e1ii,e2i,e2ii);
 473     }
 474
 475   } // end of loop over PHOS clusters
 476
 477   delete caloClustersArr;
 478   PostData(1,fOutputContainer);
 479 }
 480
 481 //__________________________________________________
 482 void AliAnalysisTaskPHOSPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells *cells, AliAODCaloCluster* clu, Int_t& maxId)
 483 {
 484   //For a given CaloCluster calculates the absId of the cell
 485   //with maximum energy deposit.
 486
 487   Double_t eMax = -111;
 488
 489   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 490     Int_t cellAbsId = clu->GetCellAbsId(iDig);
 491     Double_t eCell = cells->GetCellAmplitude(cellAbsId)*clu->GetCellAmplitudeFraction(iDig);
 492     if(eCell>eMax)  {
 493       eMax = eCell;
 494       maxId = cellAbsId;
 495     }
 496   }
 497
 498 }
 499
 500 //__________________________________________________
 501 void AliAnalysisTaskPHOSPi0CalibSelection::SetCalibCorrections(AliPHOSCalibData* cdata)
 502 {
 503   //Set new correction factors (~1) to calibration coefficients, delete previous.
 504
 505    if(fCalibData) delete fCalibData;
 506    fCalibData = cdata;
 507
 508 }
 509