PWG2/FORWARD/analysis2/AliCentralMultiplicityTask.cxx

   1 //====================================================================
   2 //
   3 // Base class for classes that calculate the multiplicity in the
   4 // central region event-by-event
   5 //
   6 // Inputs:
   7 //   - AliESDEvent
   8 //
   9 // Outputs:
  10 //   - AliAODCentralMult
  11 //
  12 // Histograms
  13 //
  14 // Corrections used
  15 #include "AliCentralMultiplicityTask.h"
  16 #include "AliForwardCorrectionManager.h"
  17 #include "AliForwardUtil.h"
  18 #include "AliLog.h"
  19 #include "AliAODHandler.h"
  20 #include "AliInputEventHandler.h"
  21 #include "AliESDInputHandler.h"
  22 #include "AliAnalysisManager.h"
  23 #include "AliESDEvent.h"
  24 #include "AliMultiplicity.h"
  25 #include <TROOT.h>
  26 #include <TFile.h>
  27 #include <TError.h>
  28 #include <iostream>
  29 #include <iomanip>
  30
  31 //====================================================================
  32 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name)
  33   : AliAnalysisTaskSE(name),
  34     fData(0),
  35     fList(0),
  36     fAODCentral(kFALSE),
  37     fManager()
  38 {
  39
  40   DefineOutput(1, TList::Class());
  41 }
  42 //____________________________________________________________________
  43 void AliCentralMultiplicityTask::UserCreateOutputObjects()
  44 {
  45
  46   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
  47   AliAODHandler*      ah =
  48     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
  49   if (!ah) AliFatal("No AOD output handler set in analysis manager");
  50
  51
  52   TObject* obj = &fAODCentral;
  53   ah->AddBranch("AliAODCentralMult", &obj);
  54
  55   fList = new TList();
  56   PostData(1,fList);
  57
  58 }
  59 //____________________________________________________________________
  60 void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/)
  61 {
  62
  63   AliESDInputHandler* eventHandler =
  64     dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()
  65                                        ->GetInputEventHandler());
  66   if (!eventHandler) {
  67     AliWarning("No inputhandler found for this event!");
  68     return;}
  69
  70   AliESDEvent* esd = eventHandler->GetEvent();
  71
  72   //Selecting only events with |valid vertex| < 10 cm
  73   const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
  74   if (!vertex) return;
  75   if(!(vertex->GetStatus())) return;
  76   if(vertex->GetNContributors() <= 0) return ;
  77   if(vertex->GetZRes() > 0.1 ) return;
  78   Double_t vertexXYZ[3]={0,0,0};
  79   vertex->GetXYZ(vertexXYZ);
  80   if(TMath::Abs(vertexXYZ[2]) > 10) return;
  81
  82   Double_t delta           = 2 ;
  83   Double_t vertexBinDouble = (vertexXYZ[2] + 10) / delta;
  84   //HHD: The vtxbins are 1-10, not 0-9
  85   Int_t    vtxbin          = Int_t(vertexBinDouble + 1) ;
  86
  87   // Make sure AOD is filled
  88   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
  89   AliAODHandler*      ah =
  90     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
  91   if (!ah)
  92     AliFatal("No AOD output handler set in analysis manager");
  93
  94   ah->SetFillAOD(kTRUE);
  95
  96   //Doing analysis
  97   fAODCentral.Clear("");
  98   TH2D *aodHist = &(fAODCentral.GetHistogram());
  99
 100   const AliMultiplicity* spdmult = esd->GetMultiplicity();
 101   //Filling clusters in layer 1 used for tracklets...
 102   for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
 103     aodHist->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
 104
 105   //...and then the unused ones in layer 1
 106   for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++)
 107     aodHist->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),
 108                   spdmult->GetPhiSingle(j));
 109
 110
 111   // Corrections
 112   TH2D* hSecMap     = fManager.GetSecMapCorrection(vtxbin);
 113   TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
 114   if (!hSecMap)     AliFatal("No secondary map!");
 115   if (!hAcceptance) AliFatal("No acceptance!");
 116
 117   aodHist->Divide(hSecMap);
 118
 119   for(Int_t nx = 1; nx <= aodHist->GetNbinsX(); nx++) {
 120     Float_t acccor = hAcceptance->GetBinContent(nx);
 121
 122     Bool_t etabinSeen = kFALSE;
 123     for(Int_t ny = 1; ny <= aodHist->GetNbinsY(); ny++) {
 124       Float_t aodvalue = aodHist->GetBinContent(nx,ny);
 125       Float_t seccor = hSecMap->GetBinContent(nx,ny);
 126       if(seccor > 0.5) etabinSeen = kTRUE;
 127       if(aodvalue < 0.000001) { aodHist->SetBinContent(nx,ny, 0); continue; }
 128
 129       Float_t aodnew   = aodvalue / acccor ;
 130       aodHist->SetBinContent(nx,ny, aodnew);
 131       Float_t aodErr   = aodHist->GetBinError(nx,ny);
 132       Float_t accErr   = hAcceptance->GetBinError(nx);
 133       Float_t error    = aodnew *TMath::Sqrt(TMath::Power(aodErr/aodvalue,2) +
 134                                              TMath::Power(accErr/acccor,2) );
 135       aodHist->SetBinError(nx,ny,error);
 136
 137     }
 138     //Filling underflow bin if we eta bin is in range
 139     if(etabinSeen) aodHist->SetBinContent(nx,0, 1.);
 140   }
 141
 142   PostData(1,fList);
 143 }
 144 //____________________________________________________________________
 145 void AliCentralMultiplicityTask::Terminate(Option_t* /*option*/)
 146 {
 147 }
 148 //____________________________________________________________________
 149 void
 150 AliCentralMultiplicityTask::Print(Option_t* /*option*/) const
 151 {
 152 }
 153 //====================================================================
 154 AliCentralMultiplicityTask::Manager::Manager() :
 155   fAcceptancePath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralAcceptance"),
 156   fSecMapPath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralSecMap"),
 157   fAcceptance(),
 158   fSecmap(),
 159   fAcceptanceName("centralacceptance"),
 160   fSecMapName("centralsecmap")
 161 {
 162
 163
 164 }
 165 //____________________________________________________________________
 166 AliCentralMultiplicityTask::Manager::Manager(const Manager& o)
 167   :fAcceptancePath(o.fAcceptancePath),
 168    fSecMapPath(o.fSecMapPath),
 169    fAcceptance(o.fAcceptance),
 170    fSecmap(o.fSecmap),
 171    fAcceptanceName(o.fAcceptanceName),
 172    fSecMapName(o.fSecMapName)
 173 {}
 174 //____________________________________________________________________
 175 AliCentralMultiplicityTask::Manager&
 176 AliCentralMultiplicityTask::Manager::operator=(const Manager& o)
 177 {
 178   fAcceptancePath = o.fAcceptancePath;
 179   fSecMapPath     = o.fSecMapPath;
 180   fAcceptance     = o.fAcceptance;
 181   fSecmap         = o.fSecmap;
 182   fAcceptanceName = o.fAcceptanceName;
 183   fSecMapName     = o.fSecMapName;
 184   return *this;
 185 }
 186
 187 //____________________________________________________________________
 188 const char*
 189 AliCentralMultiplicityTask::Manager::GetFullFileName(UShort_t what,
 190                                                      UShort_t sys,
 191                                                      UShort_t sNN,
 192                                                      Short_t  field) const
 193 {
 194   return Form("%s/%s",
 195               what == 0 ? GetSecMapPath() : GetAcceptancePath(),
 196               GetFileName(what, sys, sNN, field));
 197 }
 198
 199 //____________________________________________________________________
 200 const char*
 201 AliCentralMultiplicityTask::Manager::GetFileName(UShort_t  what ,
 202                                                  UShort_t  sys,
 203                                                  UShort_t  sNN,
 204                                                  Short_t   field) const
 205 {
 206   // Must be static - otherwise the data may disappear on return from
 207   // this member function
 208   static TString fname = "";
 209
 210   fname = "";
 211   switch(what) {
 212   case 0:  fname.Append(fSecMapName.Data());     break;
 213   case 1:  fname.Append(fAcceptanceName.Data()); break;
 214   default:
 215     ::Error("GetFileName",
 216             "Invalid indentifier %d for central object, must be 0 or 1!", what);
 217     break;
 218   }
 219   fname.Append(Form("_%s_%04dGeV_%c%1dkG.root",
 220                     AliForwardUtil::CollisionSystemString(sys),
 221                     sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
 222
 223   return fname.Data();
 224 }
 225
 226 //____________________________________________________________________
 227 TH2D*
 228 AliCentralMultiplicityTask::Manager::GetSecMapCorrection(UShort_t vtxbin) const
 229 {
 230   if (!fSecmap) {
 231     ::Warning("GetSecMapCorrection","No secondary map defined");
 232     return 0;
 233   }
 234   return fSecmap->GetCorrection(vtxbin);
 235 }
 236 //____________________________________________________________________
 237 TH1D*
 238 AliCentralMultiplicityTask::Manager::GetAcceptanceCorrection(UShort_t vtxbin)
 239   const
 240 {
 241   if (!fAcceptance) {
 242     ::Warning("GetAcceptanceCorrection","No acceptance map defined");
 243     return 0;
 244   }
 245   return fAcceptance->GetCorrection(vtxbin);
 246 }
 247
 248 //____________________________________________________________________
 249 void
 250 AliCentralMultiplicityTask::Manager::Init(UShort_t  sys,
 251                                           UShort_t  sNN,
 252                                           Short_t   field)
 253 {
 254   TFile fsec(GetFullFileName(0,sys,sNN,field));
 255   fSecmap =
 256     dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));
 257   if(!fSecmap) {
 258     ::Error("Init", "no central Secondary Map found!") ;
 259     return;
 260   }
 261   TFile facc(GetFullFileName(1,sys,sNN,field));
 262   fAcceptance =
 263     dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
 264   if(!fAcceptance) {
 265     ::Error("Init", "no central Acceptance found!") ;
 266     return;
 267   }
 268
 269 }
 270 //
 271 // EOF
 272 //