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 "AliFMDEventInspector.h"
  26 #include <TROOT.h>
  27 #include <TFile.h>
  28 #include <TError.h>
  29 #include <iostream>
  30 #include <iomanip>
  31
  32 //====================================================================
  33 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name)
  34   : AliAnalysisTaskSE(name),
  35     fData(0),
  36     fList(0),
  37     fAODCentral(kFALSE),
  38     fManager(),
  39     fUseSecondary(true),
  40     firstEventSeen(false)
  41 {
  42
  43   DefineOutput(1, TList::Class());
  44 }
  45 //____________________________________________________________________
  46 AliCentralMultiplicityTask::AliCentralMultiplicityTask()
  47   : AliAnalysisTaskSE(),
  48     fData(0),
  49     fList(0),
  50     fAODCentral(),
  51     fManager(),
  52     fUseSecondary(true),
  53     firstEventSeen(false)
  54 {
  55 }
  56 //____________________________________________________________________
  57 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const AliCentralMultiplicityTask& o)
  58   : AliAnalysisTaskSE(o),
  59     fData(o.fData),
  60     fList(o.fList),
  61     fAODCentral(o.fAODCentral),
  62     fManager(o.fManager),
  63     fUseSecondary(o.fUseSecondary),
  64     firstEventSeen(o.firstEventSeen)
  65 {
  66 }
  67 //____________________________________________________________________
  68 AliCentralMultiplicityTask&
  69 AliCentralMultiplicityTask::operator=(const AliCentralMultiplicityTask& o)
  70 {
  71   fData          = o.fData;
  72   fList          = o.fList;
  73   fAODCentral    = o.fAODCentral;
  74   fManager       = o.fManager;
  75   fUseSecondary  = o.fUseSecondary;
  76   firstEventSeen = o.firstEventSeen;
  77   return *this;
  78 }
  79 //____________________________________________________________________
  80 void AliCentralMultiplicityTask::UserCreateOutputObjects()
  81 {
  82
  83   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
  84   AliAODHandler*      ah =
  85     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
  86   if (!ah) AliFatal("No AOD output handler set in analysis manager");
  87
  88
  89   TObject* obj = &fAODCentral;
  90   ah->AddBranch("AliAODCentralMult", &obj);
  91
  92   fList = new TList();
  93   PostData(1,fList);
  94
  95 }
  96 //____________________________________________________________________
  97 void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/)
  98 {
  99
 100   AliESDInputHandler* eventHandler =
 101     dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()
 102                                        ->GetInputEventHandler());
 103   if (!eventHandler) {
 104     AliWarning("No inputhandler found for this event!");
 105     return;}
 106
 107   AliESDEvent* esd = eventHandler->GetEvent();
 108
 109   if(!GetManager().IsInit() && !firstEventSeen) {
 110     AliFMDEventInspector inspector;
 111     inspector.ReadRunDetails(esd);
 112     GetManager().Init(inspector.GetCollisionSystem(),
 113                       inspector.GetEnergy(),
 114                       inspector.GetField());
 115
 116     //std::cout<<inspector.GetCollisionSystem()<<"  "<<inspector.GetEnergy()<<"    "<<inspector.GetField()<<std::endl;
 117     AliInfo("Manager of corrections in AliCentralMultiplicityTask init");
 118     firstEventSeen = kTRUE;
 119   }
 120
 121   //Selecting only events with |valid vertex| < 10 cm
 122   const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
 123   if (!vertex) return;
 124   if(!(vertex->GetStatus())) return;
 125   if(vertex->GetNContributors() <= 0) return ;
 126   if(vertex->GetZRes() > 0.1 ) return;
 127   Double_t vertexXYZ[3]={0,0,0};
 128   vertex->GetXYZ(vertexXYZ);
 129   if(TMath::Abs(vertexXYZ[2]) > 10) return;
 130
 131   Double_t delta           = 2 ;
 132   Double_t vertexBinDouble = (vertexXYZ[2] + 10) / delta;
 133   //HHD: The vtxbins are 1-10, not 0-9
 134   Int_t    vtxbin          = Int_t(vertexBinDouble + 1) ;
 135
 136   // Make sure AOD is filled
 137   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 138   AliAODHandler*      ah =
 139     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 140   if (!ah)
 141     AliFatal("No AOD output handler set in analysis manager");
 142
 143   ah->SetFillAOD(kTRUE);
 144
 145   //Doing analysis
 146   fAODCentral.Clear("");
 147   TH2D *aodHist = &(fAODCentral.GetHistogram());
 148
 149   const AliMultiplicity* spdmult = esd->GetMultiplicity();
 150   //Filling clusters in layer 1 used for tracklets...
 151   for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
 152     aodHist->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
 153
 154   //...and then the unused ones in layer 1
 155   for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++)
 156     aodHist->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),
 157                   spdmult->GetPhiSingle(j));
 158
 159
 160   // Corrections
 161   TH2D* hSecMap     = 0;
 162   TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
 163   if (fUseSecondary) hSecMap = fManager.GetSecMapCorrection(vtxbin);
 164   if (fUseSecondary && !hSecMap)  AliFatal("No secondary map!");
 165   if (!hAcceptance)               AliFatal("No acceptance!");
 166
 167   if (hSecMap) aodHist->Divide(hSecMap);
 168
 169   for(Int_t nx = 1; nx <= aodHist->GetNbinsX(); nx++) {
 170     Float_t accCor = hAcceptance->GetBinContent(nx);
 171
 172     Bool_t etabinSeen = kFALSE;
 173     for(Int_t ny = 1; ny <= aodHist->GetNbinsY(); ny++) {
 174       Float_t aodValue = aodHist->GetBinContent(nx,ny);
 175       Float_t secCor   = hSecMap->GetBinContent(nx,ny);
 176       if (secCor > 0.5) etabinSeen = kTRUE;
 177       if (aodValue < 0.000001) { aodHist->SetBinContent(nx,ny, 0); continue; }
 178       if (accCor   < 0.000001) accCor = 1;
 179       Float_t aodNew   = aodValue / accCor ;
 180       aodHist->SetBinContent(nx,ny, aodNew);
 181       Float_t aodErr   = aodHist->GetBinError(nx,ny);
 182       Float_t accErr   = hAcceptance->GetBinError(nx);
 183       Float_t error    = aodNew *TMath::Sqrt(TMath::Power(aodErr/aodValue,2) +
 184                                              TMath::Power(accErr/accCor,2) );
 185       aodHist->SetBinError(nx,ny,error);
 186
 187     }
 188     //Filling underflow bin if we eta bin is in range
 189     if(etabinSeen) aodHist->SetBinContent(nx,0, 1.);
 190   }
 191
 192   PostData(1,fList);
 193 }
 194 //____________________________________________________________________
 195 void AliCentralMultiplicityTask::Terminate(Option_t* /*option*/)
 196 {
 197 }
 198 //____________________________________________________________________
 199 void
 200 AliCentralMultiplicityTask::Print(Option_t* /*option*/) const
 201 {
 202 }
 203 //====================================================================
 204 AliCentralMultiplicityTask::Manager::Manager() :
 205   fAcceptancePath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralAcceptance"),
 206   fSecMapPath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralSecMap"),
 207   fAcceptance(),
 208   fSecmap(),
 209   fAcceptanceName("centralacceptance"),
 210   fSecMapName("centralsecmap"),
 211   fIsInit(kFALSE)
 212 {
 213
 214
 215 }
 216 //____________________________________________________________________
 217 AliCentralMultiplicityTask::Manager::Manager(const Manager& o)
 218   :fAcceptancePath(o.fAcceptancePath),
 219    fSecMapPath(o.fSecMapPath),
 220    fAcceptance(o.fAcceptance),
 221    fSecmap(o.fSecmap),
 222    fAcceptanceName(o.fAcceptanceName),
 223    fSecMapName(o.fSecMapName),
 224    fIsInit(o.fIsInit)
 225 {}
 226 //____________________________________________________________________
 227 AliCentralMultiplicityTask::Manager&
 228 AliCentralMultiplicityTask::Manager::operator=(const Manager& o)
 229 {
 230   fAcceptancePath = o.fAcceptancePath;
 231   fSecMapPath     = o.fSecMapPath;
 232   fAcceptance     = o.fAcceptance;
 233   fSecmap         = o.fSecmap;
 234   fAcceptanceName = o.fAcceptanceName;
 235   fSecMapName     = o.fSecMapName;
 236   fIsInit         = o.fIsInit;
 237   return *this;
 238 }
 239
 240 //____________________________________________________________________
 241 const char*
 242 AliCentralMultiplicityTask::Manager::GetFullFileName(UShort_t what,
 243                                                      UShort_t sys,
 244                                                      UShort_t sNN,
 245                                                      Short_t  field) const
 246 {
 247   return Form("%s/%s",
 248               what == 0 ? GetSecMapPath() : GetAcceptancePath(),
 249               GetFileName(what, sys, sNN, field));
 250 }
 251
 252 //____________________________________________________________________
 253 const char*
 254 AliCentralMultiplicityTask::Manager::GetFileName(UShort_t  what ,
 255                                                  UShort_t  sys,
 256                                                  UShort_t  sNN,
 257                                                  Short_t   field) const
 258 {
 259   // Must be static - otherwise the data may disappear on return from
 260   // this member function
 261   static TString fname = "";
 262
 263   fname = "";
 264   switch(what) {
 265   case 0:  fname.Append(fSecMapName.Data());     break;
 266   case 1:  fname.Append(fAcceptanceName.Data()); break;
 267   default:
 268     ::Error("GetFileName",
 269             "Invalid indentifier %d for central object, must be 0 or 1!", what);
 270     break;
 271   }
 272   fname.Append(Form("_%s_%04dGeV_%c%1dkG.root",
 273                     AliForwardUtil::CollisionSystemString(sys),
 274                     sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
 275
 276   return fname.Data();
 277 }
 278
 279 //____________________________________________________________________
 280 TH2D*
 281 AliCentralMultiplicityTask::Manager::GetSecMapCorrection(UShort_t vtxbin) const
 282 {
 283   if (!fSecmap) {
 284     ::Warning("GetSecMapCorrection","No secondary map defined");
 285     return 0;
 286   }
 287   return fSecmap->GetCorrection(vtxbin);
 288 }
 289 //____________________________________________________________________
 290 TH1D*
 291 AliCentralMultiplicityTask::Manager::GetAcceptanceCorrection(UShort_t vtxbin)
 292   const
 293 {
 294   if (!fAcceptance) {
 295     ::Warning("GetAcceptanceCorrection","No acceptance map defined");
 296     return 0;
 297   }
 298   return fAcceptance->GetCorrection(vtxbin);
 299 }
 300
 301 //____________________________________________________________________
 302 void
 303 AliCentralMultiplicityTask::Manager::Init(UShort_t  sys,
 304                                           UShort_t  sNN,
 305                                           Short_t   field)
 306 {
 307   if(fIsInit) ::Warning("Init","Already initialised - overriding...");
 308
 309   TFile fsec(GetFullFileName(0,sys,sNN,field));
 310   fSecmap =
 311     dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));
 312   if(!fSecmap) {
 313     ::Error("Init", "no central Secondary Map found!") ;
 314     return;
 315   }
 316   TFile facc(GetFullFileName(1,sys,sNN,field));
 317   fAcceptance =
 318     dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
 319   if(!fAcceptance) {
 320     ::Error("Init", "no central Acceptance found!") ;
 321     return;
 322   }
 323
 324   if(fSecmap && fAcceptance) {
 325     fIsInit = kTRUE;
 326     ::Info("Init",
 327            "Central Manager initialised for sys %d, energy %d, field %d",sys,sNN,field);
 328
 329   }
 330
 331 }
 332 //
 333 // EOF
 334 //