PHOS/macros/pi0Calib/AliAnalysisTaskPi0CalibSelection.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: Boris Polishchuk                                               *
   5  *                                                                        *
   6  * Permission to use, copy, modify and distribute this software and its   *
   7  * documentation strictly for non-commercial purposes is hereby granted   *
   8  * without fee, provided that the above copyright notice appears in all   *
   9  * copies and that both the copyright notice and this permission notice   *
  10  * appear in the supporting documentation. The authors make no claims     *
  11  * about the suitability of this software for any purpose. It is          *
  12  * provided "as is" without express or implied warranty.                  *
  13  **************************************************************************/
  14
  15 //---------------------------------------------------------------------------//
  16 //                                                                           //
  17 // Fill histograms (one per cell) with two-cluster invariant mass            //
  18 // using calibration coefficients of the previous iteration.                 //
  19 // Histogram for a given cell is filled if the most energy of one cluster    //
  20 // is deposited in this cell and the other cluster could be anywhere in PHOS.//
  21 //                                                                           //
  22 //---------------------------------------------------------------------------//
  23
  24 #include "AliAnalysisTaskPi0CalibSelection.h"
  25 #include "TLorentzVector.h"
  26 #include "TVector3.h"
  27 #include "AliESDEvent.h"
  28 #include "AliPHOSPIDv1.h"
  29 #include "AliPHOSEsdCluster.h"
  30 #include "TRefArray.h"
  31 #include "TList.h"
  32
  33 ClassImp(AliAnalysisTaskPi0CalibSelection)
  34
  35 AliAnalysisTaskPi0CalibSelection::AliAnalysisTaskPi0CalibSelection() :
  36 AliAnalysisTaskSE(),fOutputContainer(0x0),fPhosGeo(0x0),fCalibData(0x0),
  37 fHmgg(0x0),fEmin(0.),fLogWeight(4.5)
  38 {
  39   //Default constructor.
  40
  41   for(Int_t iMod=0; iMod<5; iMod++) {
  42     for(Int_t iX=0; iX<64; iX++) {
  43       for(Int_t iZ=0; iZ<56; iZ++) {
  44         fHmpi0[iMod][iX][iZ]=0;
  45       }
  46     }
  47   }
  48
  49 }
  50
  51 AliAnalysisTaskPi0CalibSelection::AliAnalysisTaskPi0CalibSelection(const char* name) :
  52   AliAnalysisTaskSE(name),fOutputContainer(0x0),fPhosGeo(0x0),fCalibData(0x0),
  53   fHmgg(0x0),fEmin(0.),fLogWeight(4.5)
  54 {
  55   //Named constructor which should be used.
  56
  57   DefineOutput(1,TList::Class());
  58
  59   for(Int_t iMod=0; iMod<5; iMod++) {
  60     for(Int_t iX=0; iX<64; iX++) {
  61       for(Int_t iZ=0; iZ<56; iZ++) {
  62         fHmpi0[iMod][iX][iZ]=0;
  63       }
  64     }
  65   }
  66
  67 }
  68
  69 AliAnalysisTaskPi0CalibSelection::~AliAnalysisTaskPi0CalibSelection()
  70 {
  71   //Destructor.
  72
  73   if(fOutputContainer){
  74     fOutputContainer->Delete() ;
  75     delete fOutputContainer ;
  76   }
  77
  78   if(fCalibData) delete fCalibData;
  79   if(fPhosGeo)   delete fPhosGeo;
  80 }
  81
  82 void AliAnalysisTaskPi0CalibSelection::UserCreateOutputObjects()
  83 {
  84   fOutputContainer = new TList();
  85
  86   char hname[128], htitl[128];
  87
  88   for(Int_t iMod=0; iMod<5; iMod++) {
  89     for(Int_t iX=0; iX<64; iX++) {
  90       for(Int_t iZ=0; iZ<56; iZ++) {
  91         snprintf(hname,128,"%d_%d_%d",iMod,iX,iZ);
  92         snprintf(htitl,128,"Two-gamma inv. mass for mod %d, cell (%d,%d)",iMod,iX,iZ);
  93         fHmpi0[iMod][iX][iZ] = new TH1F(hname,htitl,100,0.,300.);
  94         fOutputContainer->Add(fHmpi0[iMod][iX][iZ]);
  95       }
  96     }
  97   }
  98
  99   fHmgg = new TH1F("hmgg","2-cluster invariant mass",100,0.,300.);
 100   fOutputContainer->Add(fHmgg);
 101
 102   fCalibData = new AliPHOSCalibData();
 103   fPhosGeo =  AliPHOSGeometry::GetInstance("IHEP") ;
 104 }
 105
 106 void AliAnalysisTaskPi0CalibSelection::UserExec(Option_t* /* option */)
 107 {
 108   //Analysis per event.
 109
 110   AliVEvent* event = InputEvent();
 111   AliESDEvent* esd = static_cast<AliESDEvent*>(event) ;
 112
 113   Double_t v[3] ; //vertex ;
 114   esd->GetVertex()->GetXYZ(v) ;
 115   TVector3 vtx(v);
 116   printf("Vertex: (%.3f,%.3f,%.3f)\n",vtx.X(),vtx.Y(),vtx.Z());
 117
 118   Int_t runNum = esd->GetRunNumber();
 119   printf("Run number: %d\n",runNum);
 120
 121   for(Int_t mod=0; mod<5; mod++)
 122     fPhosGeo->SetMisalMatrix(esd->GetPHOSMatrix(mod),mod) ;
 123
 124   Float_t logWeight = fLogWeight;
 125   printf("Will use logWeight %.3f .\n",logWeight);
 126
 127   AliPHOSPIDv1 pid;
 128
 129   Int_t mod1 = -1;
 130   TLorentzVector p1;
 131   TLorentzVector p2;
 132   TLorentzVector p12;
 133   Int_t relid[4] ;
 134   Int_t maxId;
 135
 136   TRefArray * caloClustersArr  = new TRefArray();
 137   esd->GetPHOSClusters(caloClustersArr);
 138
 139   const Int_t kNumberOfPhosClusters   = caloClustersArr->GetEntries() ;
 140   printf("CaloClusters: %d\n", kNumberOfPhosClusters);
 141
 142   // PHOS cells
 143   AliESDCaloCells *phsCells = esd->GetPHOSCells();
 144
 145   // loop over PHOS clusters
 146   for(Int_t iClu=0; iClu<kNumberOfPhosClusters; iClu++) {
 147
 148     AliESDCaloCluster *c1 = (AliESDCaloCluster *) caloClustersArr->At(iClu);
 149     if(!c1->IsPHOS()) continue; // EMCAL cluster!
 150
 151     Float_t E1_i = c1->E();   // cluster energy before correction
 152     if(E1_i<fEmin) continue;
 153
 154     AliPHOSEsdCluster clu1(*c1);
 155     clu1.Recalibrate(fCalibData, phsCells);
 156     clu1.EvalAll(logWeight,vtx);
 157     clu1.EnergyCorrection() ;
 158
 159     clu1.GetMomentum(p1,v);
 160
 161     MaxEnergyCellPos(phsCells,&clu1,maxId);
 162     fPhosGeo->AbsToRelNumbering(maxId, relid);
 163
 164     mod1 = relid[0]-1;        // module
 165     Int_t iX = relid[2]-1;    // cluster X-coord
 166     Int_t iZ = relid[3]-1;    // cluster Z-coord
 167     Float_t E1_ii = clu1.E(); // cluster energy after correction
 168
 169     for (Int_t jClu=iClu; jClu<kNumberOfPhosClusters; jClu++) {
 170       AliESDCaloCluster *c2 = (AliESDCaloCluster *) caloClustersArr->At(jClu);
 171       if(!c2->IsPHOS())   continue; // EMCAL cluster!
 172       if(c2->IsEqual(c1)) continue;
 173
 174       Float_t E2_i = c2->E();
 175       if(E2_i<fEmin) continue;
 176
 177       AliPHOSEsdCluster clu2(*c2);
 178       clu2.Recalibrate(fCalibData, phsCells);
 179       clu2.EvalAll(logWeight,vtx);
 180       clu2.EnergyCorrection() ;
 181
 182       clu2.GetMomentum(p2,v);
 183       Float_t E2_ii = clu2.E();
 184
 185       p12 = p1+p2;
 186
 187       fHmgg->Fill(p12.M()*1000);
 188       fHmpi0[mod1][iX][iZ]->Fill(p12.M()*1000);
 189
 190       printf("Mass in (mod%d,%d,%d): %.3f GeV  E1_i=%f E1_ii=%f  E2_i=%f E2_ii=%f\n",
 191              mod1,iX,iZ,p12.M(),E1_i,E1_ii,E2_i,E2_ii);
 192     }
 193
 194   } // end of loop over PHOS clusters
 195
 196   delete caloClustersArr;
 197   PostData(1,fOutputContainer);
 198 }
 199
 200 void AliAnalysisTaskPi0CalibSelection::MaxEnergyCellPos(AliESDCaloCells *cells, AliESDCaloCluster* clu, Int_t& maxId)
 201 {
 202   //For a given CaloCluster calculates the absId of the cell
 203   //with maximum energy deposit.
 204
 205   Double_t eMax = -111;
 206
 207   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 208     Int_t cellAbsId = clu->GetCellAbsId(iDig);
 209     Double_t eCell = cells->GetCellAmplitude(cellAbsId)*clu->GetCellAmplitudeFraction(iDig);
 210     if(eCell>eMax)  {
 211       eMax = eCell;
 212       maxId = cellAbsId;
 213     }
 214   }
 215
 216 }
 217
 218 void AliAnalysisTaskPi0CalibSelection::SetCalibCorrections(AliPHOSCalibData* cdata)
 219 {
 220   //Set new correction factors (~1) to calibration coefficients, delete previous.
 221
 222   if(fCalibData) delete fCalibData;
 223   fCalibData = cdata;
 224
 225 }