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