PWG4/CaloCalib/AliAnalysisTaskPHOSPi0CalibSelection.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 <cstdlib>
  25
  26 // Root
  27 #include "TLorentzVector.h"
  28 #include "TVector3.h"
  29 #include "TRefArray.h"
  30 #include "TList.h"
  31
  32 // AliRoot
  33 #include "AliAnalysisTaskPHOSPi0CalibSelection.h"
  34 #include "AliAODEvent.h"
  35 #include "AliESDEvent.h"
  36 #include "AliPHOSAodCluster.h"
  37 #include "AliPHOSPIDv1.h"
  38
  39 ClassImp(AliAnalysisTaskPHOSPi0CalibSelection)
  40
  41 AliAnalysisTaskPHOSPi0CalibSelection::AliAnalysisTaskPHOSPi0CalibSelection() :
  42 AliAnalysisTaskSE(),fOutputContainer(0x0),fPhosGeo(0x0),fCalibData(0x0),fHmgg(0x0),
  43   fEmin(0.), fLogWeight(4.5)
  44 {
  45   //Default constructor.
  46
  47   for(Int_t iMod=0; iMod<5; iMod++) {
  48     for(Int_t iX=0; iX<64; iX++) {
  49       for(Int_t iZ=0; iZ<56; iZ++) {
  50         fHmpi0[iMod][iX][iZ]=0;
  51       }
  52     }
  53   }
  54
  55 }
  56
  57 AliAnalysisTaskPHOSPi0CalibSelection::AliAnalysisTaskPHOSPi0CalibSelection(const char* name) :
  58   AliAnalysisTaskSE(name),fOutputContainer(0x0),fPhosGeo(0x0),fCalibData(0x0),fHmgg(0x0),
  59   fEmin(0.), fLogWeight(4.5)
  60 {
  61   //Named constructor which should be used.
  62
  63   DefineOutput(1,TList::Class());
  64
  65   for(Int_t iMod=0; iMod<5; iMod++) {
  66     for(Int_t iX=0; iX<64; iX++) {
  67       for(Int_t iZ=0; iZ<56; iZ++) {
  68         fHmpi0[iMod][iX][iZ]=0;
  69       }
  70     }
  71   }
  72
  73 }
  74
  75 AliAnalysisTaskPHOSPi0CalibSelection::~AliAnalysisTaskPHOSPi0CalibSelection()
  76 {
  77   //Destructor.
  78
  79   if(fOutputContainer){
  80     fOutputContainer->Delete() ;
  81     delete fOutputContainer ;
  82   }
  83
  84   if(fCalibData) delete fCalibData;
  85   if(fPhosGeo)   delete fPhosGeo;
  86
  87 }
  88
  89 void AliAnalysisTaskPHOSPi0CalibSelection::UserCreateOutputObjects()
  90 {
  91   //Create output container
  92   fOutputContainer = new TList();
  93
  94   char hname[128], htitl[128];
  95
  96   for(Int_t iMod=0; iMod<5; iMod++) {
  97     for(Int_t iX=0; iX<64; iX++) {
  98       for(Int_t iZ=0; iZ<56; iZ++) {
  99         sprintf(hname,"%d_%d_%d",iMod,iX,iZ);
 100         sprintf(htitl,"Two-gamma inv. mass for mod %d, cell (%d,%d)",iMod,iX,iZ);
 101         fHmpi0[iMod][iX][iZ] = new TH1F(hname,htitl,100,0.,300.);
 102         fOutputContainer->Add(fHmpi0[iMod][iX][iZ]);
 103       }
 104     }
 105   }
 106
 107   fHmgg = new TH1F("hmgg","2-cluster invariant mass",100,0.,300.);
 108   fOutputContainer->Add(fHmgg);
 109
 110 }
 111
 112 void AliAnalysisTaskPHOSPi0CalibSelection::UserExec(Option_t* /* option */)
 113 {
 114   //Analysis per event.
 115
 116   AliAODEvent* aod = 0x0;
 117   if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) aod = dynamic_cast<AliAODEvent*>(InputEvent());
 118   else  if(!strcmp(InputEvent()->GetName(),"AliESDEvent")) aod = AODEvent();
 119   else {
 120           printf("AliAnalysisTaskPHOSPi0CalibSelection: Unknown event type, STOP!\n");
 121           abort();
 122   }
 123
 124   Double_t v[] = {aod->GetVertex(0)->GetX(),aod->GetVertex(0)->GetY(),aod->GetVertex(0)->GetZ()}; //to check!!
 125   //aod->GetVertex()->GetXYZ(v) ;
 126   TVector3 vtx(v); //Check
 127
 128   printf("Vertex: (%.3f,%.3f,%.3f)\n",vtx.X(),vtx.Y(),vtx.Z());
 129
 130   Int_t runNum = aod->GetRunNumber();
 131   printf("Run number: %d\n",runNum);
 132
 133   //Get the matrix with geometry information
 134   //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
 135   if(!strcmp(InputEvent()->GetName(),"AliAODEvent"))
 136                 printf("Use ideal geometry, values geometry matrix not kept in AODs");
 137   else{
 138           AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
 139           for(Int_t mod=0; mod<5; mod++)
 140                   fPhosGeo->SetMisalMatrix(esd->GetPHOSMatrix(mod),mod) ;
 141   }
 142
 143   printf("Will use fLogWeight %.3f .\n",fLogWeight);
 144
 145   AliPHOSPIDv1 pid;
 146
 147   Int_t mod1 = -1;
 148   TLorentzVector p1;
 149   TLorentzVector p2;
 150   TLorentzVector p12;
 151   Int_t relid[4] ;
 152   Int_t maxId;
 153
 154   TRefArray * caloClustersArr  = new TRefArray();
 155   aod->GetPHOSClusters(caloClustersArr);
 156
 157   const Int_t kNumberOfPhosClusters   = caloClustersArr->GetEntries() ;
 158   printf("CaloClusters: %d\n", kNumberOfPhosClusters);
 159
 160   // PHOS cells
 161   AliAODCaloCells *phsCells = aod->GetPHOSCells();
 162
 163   // loop over PHOS clusters
 164   for(Int_t iClu=0; iClu<kNumberOfPhosClusters; iClu++) {
 165
 166     AliAODCaloCluster *c1 = (AliAODCaloCluster *) caloClustersArr->At(iClu);
 167     if(!c1->IsPHOSCluster()) continue; // EMCAL cluster!
 168
 169     Float_t e1i = c1->E();   // cluster energy before correction
 170     if(e1i<fEmin) continue;
 171
 172     AliPHOSAodCluster clu1(*c1);
 173     clu1.Recalibrate(fCalibData, phsCells);
 174     clu1.EvalAll(fLogWeight,vtx);
 175     clu1.EnergyCorrection(&pid) ;
 176
 177     clu1.GetMomentum(p1,v);
 178
 179     MaxEnergyCellPos(phsCells,&clu1,maxId);
 180     fPhosGeo->AbsToRelNumbering(maxId, relid);
 181
 182     mod1 = relid[0]-1;        // module
 183     Int_t iX = relid[2]-1;    // cluster X-coord
 184     Int_t iZ = relid[3]-1;    // cluster Z-coord
 185         Float_t e1ii = clu1.E(); // cluster energy after correction
 186
 187     for (Int_t jClu=iClu; jClu<kNumberOfPhosClusters; jClu++) {
 188       AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
 189       if(!c2->IsPHOSCluster())   continue; // EMCAL cluster!
 190       if(c2->IsEqual(c1)) continue;
 191
 192       Float_t e2i = c2->E();
 193       if(e2i<fEmin) continue;
 194
 195       AliPHOSAodCluster clu2(*c2);
 196       clu2.Recalibrate(fCalibData, phsCells);
 197       clu2.EvalAll(fLogWeight,vtx);
 198       clu2.EnergyCorrection(&pid) ;
 199
 200       clu2.GetMomentum(p2,v);
 201           Float_t e2ii = clu2.E();
 202
 203       p12 = p1+p2;
 204
 205       fHmgg->Fill(p12.M()*1000);
 206       fHmpi0[mod1][iX][iZ]->Fill(p12.M()*1000);
 207
 208       printf("Mass in (mod%d,%d,%d): %.3f GeV  E1_i=%f E1_ii=%f  E2_i=%f E2_ii=%f\n",
 209              mod1,iX,iZ,p12.M(),e1i,e1ii,e2i,e2ii);
 210     }
 211
 212   } // end of loop over PHOS clusters
 213
 214   delete caloClustersArr;
 215   PostData(1,fOutputContainer);
 216 }
 217
 218 void AliAnalysisTaskPHOSPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells *cells, AliAODCaloCluster* clu, Int_t& maxId)
 219 {
 220   //For a given CaloCluster calculates the absId of the cell
 221   //with maximum energy deposit.
 222
 223   Double_t eMax = -111;
 224
 225   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 226     Int_t cellAbsId = clu->GetCellAbsId(iDig);
 227     Double_t eCell = cells->GetCellAmplitude(cellAbsId)*clu->GetCellAmplitudeFraction(iDig);
 228     if(eCell>eMax)  {
 229       eMax = eCell;
 230       maxId = cellAbsId;
 231     }
 232   }
 233
 234 }
 235
 236 void AliAnalysisTaskPHOSPi0CalibSelection::SetCalibCorrections(AliPHOSCalibData* cdata)
 237 {
 238   //Set new correction factors (~1) to calibration coefficients, delete previous.
 239
 240    if(fCalibData) delete fCalibData;
 241    fCalibData = cdata;
 242
 243 }
 244