PWGGA/PHOSTasks/UserTasks/AliAnalysisTaskSEPHOSpPbPi0.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2006, 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 //
  18 // AliAnalysisTaskSE for the gamma and pi0 from pPb collision analysis
  19 //
  20 // Author: H-S. Zhu, hongsheng.zhu@cern.ch
  21 //                   hszhu@iopp.ccnu.edu.cn
  22 ///////////////////////////////////////////////////////////////////////////
  23
  24 #include <TF1.h>
  25 #include <TH2I.h>
  26 #include <TList.h>
  27 #include <TMath.h>
  28 #include <TArray.h>
  29 #include <TRandom.h>
  30 #include <TVector3.h>
  31 #include <TGeoManager.h>
  32 #include <TClonesArray.h>
  33
  34 #include "AliInputEventHandler.h"
  35 #include "AliAnalysisManager.h"
  36 #include "AliMCEventHandler.h"
  37 #include "AliMCEvent.h"
  38 #include "AliStack.h"
  39 #include "AliVEvent.h"
  40 #include "AliAODEvent.h"
  41 #include "AliESDEvent.h"
  42 #include "AliAODHeader.h"
  43 #include "AliESDHeader.h"
  44 #include "AliVCluster.h"
  45 #include "AliVCaloCells.h"
  46 #include "AliAODCaloCells.h"
  47 #include "AliESDCaloCells.h"
  48 #include "AliPHOSGeometry.h"
  49 #include "AliPHOSGeoUtils.h"
  50 #include "AliPHOSCalibData.h"
  51 #include "AliPHOSAodCluster.h"
  52 #include "AliPHOSEsdCluster.h"
  53 #include "AliOADBContainer.h"
  54 #include "AliPHOSpPbPi0Header.h"
  55 #include "AliCaloClusterInfo.h"
  56 #include "AliAnalysisTaskSEPHOSpPbPi0.h"
  57
  58 ClassImp(AliAnalysisTaskSEPHOSpPbPi0)
  59
  60 const Float_t AliAnalysisTaskSEPHOSpPbPi0::kLogWeight = 4.5;
  61 Bool_t   AliAnalysisTaskSEPHOSpPbPi0::fgRemovePileup    = kFALSE;
  62 Bool_t   AliAnalysisTaskSEPHOSpPbPi0::fgUseFiducialCut  = kFALSE;
  63 Double_t AliAnalysisTaskSEPHOSpPbPi0::fgDecaliWidth     = 0.055;
  64 Double_t AliAnalysisTaskSEPHOSpPbPi0::fgCuts[5]         = { 0.3, 2., 0.2, 2.5, 7e-8 };
  65
  66 //________________________________________________________________________
  67 AliAnalysisTaskSEPHOSpPbPi0::AliAnalysisTaskSEPHOSpPbPi0():
  68   AliAnalysisTaskSE(), fIsMC(kFALSE), fCentralityBin(10), fBufferSize(10), fRunNumber(-1),
  69   fPHOSGeo(0), fPHOSCalibData(0), fListQA(0), fListRD(0), fListMC(0), fHeader(0), fCaloClArr(0)
  70 {
  71   //
  72   // Default constructor
  73   //
  74   for (Int_t i=0; i<10; i++) { for (Int_t j=0; j<10; j++) fEventList[i][j] = 0; }
  75
  76   // Init bad channel map
  77   for (Int_t i=0; i<5; i++)
  78     fPHOSBadMap[i] = new TH2I(Form("PHOS_BadMap_mod%d", i), Form("PHOS_BadMap_mod%d", i), 64, 0., 64., 56, 0., 56.);
  79 }
  80 //________________________________________________________________________
  81 AliAnalysisTaskSEPHOSpPbPi0::AliAnalysisTaskSEPHOSpPbPi0(const char *name):
  82   AliAnalysisTaskSE(name), fIsMC(kFALSE), fCentralityBin(10), fBufferSize(10), fRunNumber(-1),
  83   fPHOSGeo(0), fPHOSCalibData(0), fListQA(0), fListRD(0), fListMC(0), fHeader(0), fCaloClArr(0)
  84 {
  85   // Constructor
  86   for (Int_t i=0; i<10; i++) { for (Int_t j=0; j<10; j++) fEventList[i][j] = 0; }
  87
  88   // Init bad channel map
  89   for (Int_t i=0; i<5; i++)
  90     fPHOSBadMap[i] = new TH2I(Form("PHOS_BadMap_mod%d", i), Form("PHOS_BadMap_mod%d", i), 64, 0., 64., 56, 0., 56.);
  91
  92   DefineOutput(1, TList::Class());
  93   DefineOutput(2, TList::Class());
  94   DefineOutput(3, TList::Class());
  95 }
  96
  97 //_____________________________________________________________________________
  98 AliAnalysisTaskSEPHOSpPbPi0::~AliAnalysisTaskSEPHOSpPbPi0()
  99 {
 100   //
 101   // Default destructor
 102   //
 103   if (fListQA)      { delete fListQA;      fListQA     = NULL; }
 104   if (fListRD)      { delete fListRD;      fListRD     = NULL; }
 105   if (fListMC)      { delete fListMC;      fListMC     = NULL; }
 106   if (fHeader)      { delete fHeader;      fHeader     = NULL; }
 107   if (fCaloClArr)   { delete fCaloClArr;   fCaloClArr  = NULL; }
 108
 109 }
 110 //________________________________________________________________________
 111 void AliAnalysisTaskSEPHOSpPbPi0::UserCreateOutputObjects()
 112 {
 113   // Create the output container
 114
 115   AliPHOSpPbPi0Header::SetIsMC(fIsMC);
 116   AliPHOSpPbPi0Header::SetUseFiducialCut(fgUseFiducialCut);
 117
 118   if (!fHeader)     fHeader     = new AliPHOSpPbPi0Header();
 119   if (!fCaloClArr)  fCaloClArr  = new TClonesArray("AliCaloClusterInfo", 0);
 120   if (!fListQA)     fListQA     = new TList();
 121   if (!fListRD)     fListRD     = new TList();
 122   if (!fListMC)     fListMC     = new TList();
 123
 124   fHeader->SetNCent(fCentralityBin.GetSize()-1);
 125   fHeader->CreateHistograms(fListQA, fListRD, fListMC);
 126
 127   // Post output data.
 128   PostData(1, fListQA);
 129   PostData(2, fListRD);
 130   if (fIsMC) PostData(3, fListMC);
 131
 132   return;
 133 }
 134
 135 //________________________________________________________________________
 136 void AliAnalysisTaskSEPHOSpPbPi0::UserExec(Option_t *)
 137 {
 138   // Main loop
 139   // Called for each event
 140
 141   if (fIsMC) {
 142     if (MCEvent()) {
 143       if (MCEvent()->GetNumberOfTracks()<=0)
 144            { AliError("MC event not found. Nothing done!"); return; }
 145     } else { AliError("MC event not found. Nothing done!"); return; }
 146   }
 147
 148   AliAODEvent *aod   = 0x0;
 149   AliESDEvent *esd   = 0x0;
 150
 151   if (((TString)fInputEvent->IsA()->GetName())=="AliAODEvent") {
 152     aod = dynamic_cast<AliAODEvent*>(fInputEvent);
 153     if (!aod) { AliError("AOD event not found. Nothing done!");   return; }
 154     if (!fIsMC && (aod->GetHeader()->GetEventType()!=7))          return; // check event type; should be PHYSICS = 7 for data and 0 for MC
 155   } else {
 156     esd = dynamic_cast<AliESDEvent*>(fInputEvent);
 157     if (!esd) { AliError("ESD event not found. Nothing done!");   return; }
 158     if (!fIsMC && (esd->GetHeader()->GetEventType()!=7))          return; // check event type; should be PHYSICS = 7 for data and 0 for MC
 159   }
 160
 161   // Fill Event info
 162   fHeader->SetEventInfo(fInputHandler);
 163   fHeader->FillHistosEvent(fListQA);
 164
 165   // PHOS Geometry and Misalignment initialization at the first time it runs
 166   if(fRunNumber != fInputEvent->GetRunNumber()) {
 167     fRunNumber = fInputEvent->GetRunNumber();
 168     fHeader->SetIspARun(fRunNumber>195344 && fRunNumber<197388);  // flag for pA collisions
 169     PHOSInitialize(esd);
 170   }
 171
 172   // Event Selection
 173   if (!fHeader->IsSelected())                     return;
 174   if (fgRemovePileup && fHeader->IsPileupSPD())   return;
 175
 176   // Fill PHOS cells QA histograms
 177   fHeader->FillHistosCaloCellsQA(fListQA, fInputEvent->GetPHOSCells(), fPHOSGeo);
 178
 179   // Fill PHOS cluster Clones Array
 180   FillCaloClusterInfo(aod, esd);
 181   aod = 0x0;
 182
 183   if (!fCaloClArr->GetEntriesFast()) return;
 184
 185   // vertex bining and centrality bining
 186   Int_t zvtx = (Int_t)((fHeader->Vz() + 10.)/2.);   if (zvtx<0) zvtx = 0;   if (zvtx>9) zvtx = 9;
 187   Int_t cent = TMath::BinarySearch<Float_t>(fCentralityBin.GetSize()-1, fCentralityBin.GetArray(), fHeader->Centrality());
 188   if (!fEventList[zvtx][cent]) fEventList[zvtx][cent] = new TList();
 189   TList *eventList = fEventList[zvtx][cent];
 190
 191   // Fill cluster histograms
 192   fHeader->FillHistosCaloCluster(fListQA, fCaloClArr, cent);
 193
 194   // Fill pi0 histograms
 195   fHeader->FillHistosPi0(fListRD, fCaloClArr, cent);
 196
 197   // Fill mixed pi0 histograms
 198   fHeader->FillHistosMixPi0(fListRD, fCaloClArr, eventList, cent);
 199
 200   // Fill MC info
 201   AliStack *stack = 0x0;
 202   if (fIsMC) {
 203     if (esd) {
 204       if (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) {
 205         if (static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent())
 206           stack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
 207       }
 208       fHeader->FillHistosMC(fListMC, stack, fPHOSGeo, cent);
 209     } else
 210       fHeader->FillHistosMC(fListMC, MCEvent(), fPHOSGeo, cent);
 211   }
 212   esd = 0x0;
 213
 214   // Fill event list for mixing
 215   if (fCaloClArr->GetEntriesFast()>0) {
 216     eventList->AddFirst(fCaloClArr);   fCaloClArr = 0x0;
 217     if (eventList->GetSize()>fBufferSize[cent]) { // Remove redundant events
 218       TClonesArray *tmp = static_cast<TClonesArray*>(eventList->Last());
 219       eventList->RemoveLast();
 220       delete tmp;
 221     }
 222   }
 223
 224   return;
 225 }
 226
 227 //________________________________________________________________________
 228 void AliAnalysisTaskSEPHOSpPbPi0::Terminate(Option_t *)
 229 {
 230   // Terminate analysis
 231
 232   // add the correction matrix
 233
 234   return;
 235 }
 236
 237 //________________________________________________________________________
 238 void AliAnalysisTaskSEPHOSpPbPi0::PHOSInitialize(AliESDEvent* const esd)
 239 {
 240   // Initialize PHOS Geometry ,misalignment and calibration
 241 //TGeoManager::Import("$ALICE_ROOT/test/QA/geometry.root");
 242
 243   fPHOSGeo =  AliPHOSGeometry::GetInstance("IHEP");
 244   AliOADBContainer geomContainer("phosGeo");
 245   geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
 246   TObjArray *matrixes = (TObjArray*)geomContainer.GetObject(fRunNumber,"PHOSRotationMatrixes");
 247   for (Int_t mod=0; mod<5; mod++) {
 248     if (!matrixes->At(mod)) continue;
 249     else fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod);
 250   }
 251
 252   // sets the PHOS Misalignment vertex if ESD
 253   if (esd) {
 254     for (Int_t mod=0; mod<5; mod++) {
 255       const TGeoHMatrix* modMatrix = fInputEvent->GetPHOSMatrix(mod);
 256       if (!modMatrix) continue;
 257       else fPHOSGeo->SetMisalMatrix(modMatrix, mod);
 258     }
 259   }
 260
 261   TRandom random;   random.SetSeed(0); // the seed is set to the current  machine clock
 262   fPHOSCalibData = new AliPHOSCalibData();   fPHOSCalibData->SetName("PHOSCalibData");
 263   for(Int_t iMod=1; iMod<6; iMod++) {
 264     for(Int_t iCol=1; iCol<57; iCol++) {
 265       for(Int_t iRow=1; iRow<65; iRow++) {
 266         Float_t value = fIsMC ? random.Gaus(1., fgDecaliWidth) : 0.;  // ADC channel Emc
 267         fPHOSCalibData->SetADCchannelEmc(iMod, iCol, iRow, value);
 268       }
 269     }
 270   }
 271 }
 272
 273 //________________________________________________________________________
 274 void AliAnalysisTaskSEPHOSpPbPi0::FillCaloClusterInfo(AliAODEvent* const aod, AliESDEvent* const esd)
 275 {
 276   // Fill calo cluster info
 277   if (fCaloClArr) fCaloClArr->Clear();
 278   else fCaloClArr = new TClonesArray("AliCaloClusterInfo", 0);
 279
 280   Int_t    nclsts       = fInputEvent->GetNumberOfCaloClusters();
 281   Int_t    countN       = fCaloClArr->GetEntriesFast();
 282   Int_t    relID[4]     = {0,0,0,0 }; // module = relID[0]; cellX = relID[2]; cellZ = relID[3];
 283   Float_t  position[3]  = {0.,0.,0.};
 284   Double_t vtx[3]       = {0,0,0};   fHeader->GetXYZ(vtx);   TVector3 vtxVector(vtx);
 285   Double_t magfield     = fInputEvent->GetMagneticField();
 286   TF1 *nonLinCorr = new TF1("Non-linear", "0.0241+1.0504*x+0.000249*x*x", 0., 50.); // GCB's non-linear correction function
 287
 288   TClonesArray &caloRef = *fCaloClArr;
 289   TLorentzVector momentum;
 290   AliVCluster        *clust       = 0x0;
 291   AliCaloClusterInfo *caloCluster = 0x0;
 292   for (Int_t iclst=0; iclst<nclsts; iclst++) {  // loop over all clusters
 293     clust = fInputEvent->GetCaloCluster(iclst);
 294     if (!(clust && clust->IsPHOS() && clust->E()>fgCuts[0]))                 { clust=0; continue; }
 295     if (!(clust->GetNCells()>(Int_t)fgCuts[1] && clust->GetM02()>fgCuts[2])) { clust=0; continue; } // To remove exotic clusters
 296     if (!(clust->GetDistanceToBadChannel()>fgCuts[3]))                       { clust=0; continue; }
 297 //  if (!(TMath::Abs(clust->GetTOF()<fgCuts[4])))                            { clust=0; continue; } // remove clusters from pileup or same benches
 298     clust->GetPosition(position); TVector3 global(position); fPHOSGeo->GlobalPos2RelId(global,relID);
 299     if (relID[0] == 2)                                                       { clust=0; continue; } // !remove module 2
 300     if (!IsGoodCaloCluster(relID[0], relID[2], relID[3]))                    { clust=0; continue; } // calo cluster selection
 301
 302     caloCluster = new AliCaloClusterInfo(clust, esd, relID, magfield);
 303     if (fgUseFiducialCut && !caloCluster->IsInFiducialRegion(relID[2], relID[3])) { delete caloCluster; caloCluster=0; clust=0; continue; } // Fiducial cut
 304
 305     if (aod) { // TODO recalibration for AOD
 306       if (fIsMC) {
 307         AliPHOSAodCluster aodClust(*(AliAODCaloCluster*) (clust));
 308         aodClust.Recalibrate(fPHOSCalibData, aod->GetPHOSCells());
 309         aodClust.EvalAll(kLogWeight, vtxVector);        // recalculate all cluster parameters
 310         aodClust.SetE(nonLinCorr->Eval(aodClust.E()));  // non-linear correction
 311         aodClust.GetMomentum(momentum, vtx);
 312       } else clust->GetMomentum(momentum, vtx);
 313       caloCluster->SetLorentzVector(momentum);
 314     } else { // for ESD
 315       AliPHOSEsdCluster esdClust(*(AliESDCaloCluster*) (clust));
 316       esdClust.Recalibrate(fPHOSCalibData, esd->GetPHOSCells());
 317       esdClust.EvalAll(kLogWeight, vtxVector);       // recalculate all cluster parameters
 318       esdClust.SetE(nonLinCorr->Eval(esdClust.E()));  // non-linear correction
 319       esdClust.GetMomentum(momentum, vtx);
 320       caloCluster->SetLorentzVector(momentum);
 321     }
 322
 323     if (fIsMC && !(momentum.E()>fgCuts[0])) { delete caloCluster; caloCluster=0; clust=0; continue; } // check for MC
 324
 325     Double_t disp = caloCluster->TestDisp();
 326     if (disp < 2.5*2.5) caloCluster->SetPIDBit(BIT(1)); // set Disp1 bit
 327     if (disp < 1.5*1.5) caloCluster->SetPIDBit(BIT(3)); // set Disp2 bit
 328
 329     clust = 0;
 330
 331     new(caloRef[countN++]) AliCaloClusterInfo(*caloCluster);
 332     delete caloCluster; caloCluster=0;
 333   }  // end loop of all clusters
 334 }
 335
 336 //________________________________________________________________________
 337 Bool_t AliAnalysisTaskSEPHOSpPbPi0::IsGoodCaloCluster(Int_t iMod, Int_t cellX, Int_t cellZ)
 338 {
 339   // !Check whether this cluster is not in bad channel
 340
 341   if (!(iMod>0 && iMod<5 && fPHOSBadMap[iMod]))          return kTRUE;   //No bad maps for this Module
 342   if (fPHOSBadMap[iMod]->GetBinContent(cellX,cellZ)>0)   return kFALSE;
 343
 344   return kTRUE;
 345 }