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     fUseAcceptance(true),
  41     fFirstEventSeen(false)
  42 {
  43   //
  44   // Constructor
  45   //
  46   DefineOutput(1, TList::Class());
  47   fBranchNames =
  48     "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
  49     "SPDVertex.,PrimaryVertex.";
  50 }
  51 //____________________________________________________________________
  52 AliCentralMultiplicityTask::AliCentralMultiplicityTask()
  53   : AliAnalysisTaskSE(),
  54     fData(0),
  55     fList(0),
  56     fAODCentral(),
  57     fManager(),
  58     fUseSecondary(true),
  59     fUseAcceptance(true),
  60     fFirstEventSeen(false)
  61 {
  62   //
  63   // Constructor
  64   //
  65 }
  66 //____________________________________________________________________
  67 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const AliCentralMultiplicityTask& o)
  68   : AliAnalysisTaskSE(o),
  69     fData(o.fData),
  70     fList(o.fList),
  71     fAODCentral(o.fAODCentral),
  72     fManager(o.fManager),
  73     fUseSecondary(o.fUseSecondary),
  74     fUseAcceptance(o.fUseAcceptance),
  75     fFirstEventSeen(o.fFirstEventSeen)
  76 {
  77   //
  78   // Copy constructor
  79   //
  80 }
  81 //____________________________________________________________________
  82 AliCentralMultiplicityTask&
  83 AliCentralMultiplicityTask::operator=(const AliCentralMultiplicityTask& o)
  84 {
  85   //
  86   // Assignment operator
  87   //
  88   fData           = o.fData;
  89   fList           = o.fList;
  90   fAODCentral     = o.fAODCentral;
  91   fManager        = o.fManager;
  92   fUseSecondary   = o.fUseSecondary;
  93   fUseAcceptance  = o.fUseAcceptance;
  94   fFirstEventSeen = o.fFirstEventSeen;
  95   return *this;
  96 }
  97 //____________________________________________________________________
  98 void AliCentralMultiplicityTask::UserCreateOutputObjects()
  99 {
 100   //
 101   // Create output objects
 102   //
 103   //
 104
 105   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 106   AliAODHandler*      ah =
 107     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 108   if (!ah) AliFatal("No AOD output handler set in analysis manager");
 109
 110
 111   TObject* obj = &fAODCentral;
 112   ah->AddBranch("AliAODCentralMult", &obj);
 113
 114   fList = new TList();
 115   fList->SetOwner();
 116   PostData(1,fList);
 117
 118 }
 119 //____________________________________________________________________
 120 void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/)
 121 {
 122   //
 123   // Process each event
 124   //
 125   // Parameters:
 126   //    option Not used
 127   //
 128
 129   AliESDInputHandler* eventHandler =
 130     dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()
 131                                        ->GetInputEventHandler());
 132   if (!eventHandler) {
 133     AliWarning("No inputhandler found for this event!");
 134     return;}
 135
 136   AliESDEvent* esd = eventHandler->GetEvent();
 137
 138   if(!GetManager().IsInit() && !fFirstEventSeen) {
 139     AliFMDEventInspector inspector;
 140     inspector.ReadRunDetails(esd);
 141     GetManager().Init(inspector.GetCollisionSystem(),
 142                       inspector.GetEnergy(),
 143                       inspector.GetField());
 144
 145     AliInfo("Manager of corrections in AliCentralMultiplicityTask init");
 146     fFirstEventSeen = kTRUE;
 147   }
 148
 149   //Selecting only events with |valid vertex| < 10 cm
 150   const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
 151   if (!vertex) return;
 152   if(!(vertex->GetStatus())) return;
 153   if(vertex->GetNContributors() <= 0) return ;
 154   if(vertex->GetZRes() > 0.1 ) return;
 155   Double_t vertexXYZ[3]={0,0,0};
 156   vertex->GetXYZ(vertexXYZ);
 157   if(TMath::Abs(vertexXYZ[2]) > 10) return;
 158
 159   Double_t delta           = 2 ;
 160   Double_t vertexBinDouble = (vertexXYZ[2] + 10) / delta;
 161   //HHD: The vtxbins are 1-10, not 0-9
 162   Int_t    vtxbin          = Int_t(vertexBinDouble + 1) ;
 163
 164   // Make sure AOD is filled
 165   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 166   AliAODHandler*      ah =
 167     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
 168   if (!ah)
 169     AliFatal("No AOD output handler set in analysis manager");
 170
 171   ah->SetFillAOD(kTRUE);
 172
 173   //Doing analysis
 174   fAODCentral.Clear("");
 175   TH2D *aodHist = &(fAODCentral.GetHistogram());
 176
 177   const AliMultiplicity* spdmult = esd->GetMultiplicity();
 178   //Filling clusters in layer 1 used for tracklets...
 179   for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
 180     aodHist->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
 181
 182   //...and then the unused ones in layer 1
 183   for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++)
 184     aodHist->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),
 185                   spdmult->GetPhiSingle(j));
 186
 187
 188   // Corrections
 189   TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
 190   TH2D* hSecMap     = fManager.GetSecMapCorrection(vtxbin);
 191
 192   if (!hSecMap)     AliFatal("No secondary map!");
 193   if (!hAcceptance) AliFatal("No acceptance!");
 194
 195   if (fUseSecondary && hSecMap) aodHist->Divide(hSecMap);
 196
 197   for(Int_t nx = 1; nx <= aodHist->GetNbinsX(); nx++) {
 198     Float_t accCor = hAcceptance->GetBinContent(nx);
 199     Float_t accErr = hAcceptance->GetBinError(nx);
 200
 201     Bool_t etabinSeen = kFALSE;
 202     for(Int_t ny = 1; ny <= aodHist->GetNbinsY(); ny++) {
 203       // Get currrent value
 204       Float_t aodValue = aodHist->GetBinContent(nx,ny);
 205       Float_t aodErr   = aodHist->GetBinError(nx,ny);
 206
 207       // Set underflow bin
 208       Float_t secCor   = 0;
 209       if(hSecMap) secCor   = hSecMap->GetBinContent(nx,ny);
 210       if (secCor > 0.5) etabinSeen = kTRUE;
 211       if (aodValue < 0.000001) { aodHist->SetBinContent(nx,ny, 0); continue; }
 212
 213       if (!fUseAcceptance) continue;
 214
 215       // Acceptance correction
 216       if (accCor   < 0.000001) accCor = 1;
 217       Float_t aodNew   = aodValue / accCor ;
 218       Float_t error    = aodNew*TMath::Sqrt(TMath::Power(aodErr/aodValue,2) +
 219                                             TMath::Power(accErr/accCor,2) );
 220       aodHist->SetBinContent(nx,ny, aodNew);
 221       //test
 222       aodHist->SetBinError(nx,ny,error);
 223       aodHist->SetBinError(nx,ny,aodErr);
 224
 225     }
 226     //Filling underflow bin if we eta bin is in range
 227     if(etabinSeen) aodHist->SetBinContent(nx,0, 1.);
 228   }
 229
 230   PostData(1,fList);
 231 }
 232 //____________________________________________________________________
 233 void AliCentralMultiplicityTask::Terminate(Option_t* /*option*/)
 234 {
 235   //
 236   // End of job
 237   //
 238   // Parameters:
 239   //    option Not used
 240   //
 241 }
 242 //____________________________________________________________________
 243 void
 244 AliCentralMultiplicityTask::Print(Option_t* /*option*/) const
 245 {
 246   //
 247   // Print information
 248   //
 249   // Parameters:
 250   //    option Not used
 251   //
 252 }
 253 //====================================================================
 254 AliCentralMultiplicityTask::Manager::Manager() :
 255   fAcceptancePath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralAcceptance"),
 256   fSecMapPath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralSecMap"),
 257   fAcceptance(),
 258   fSecmap(),
 259   fAcceptanceName("centralacceptance"),
 260   fSecMapName("centralsecmap"),
 261   fIsInit(kFALSE)
 262 {
 263   //
 264   // Constructor
 265   //
 266 }
 267 //____________________________________________________________________
 268 AliCentralMultiplicityTask::Manager::Manager(const Manager& o)
 269   :fAcceptancePath(o.fAcceptancePath),
 270    fSecMapPath(o.fSecMapPath),
 271    fAcceptance(o.fAcceptance),
 272    fSecmap(o.fSecmap),
 273    fAcceptanceName(o.fAcceptanceName),
 274    fSecMapName(o.fSecMapName),
 275    fIsInit(o.fIsInit)
 276 {
 277   //
 278   // Copy Constructor
 279   //
 280 }
 281 //____________________________________________________________________
 282 AliCentralMultiplicityTask::Manager&
 283 AliCentralMultiplicityTask::Manager::operator=(const Manager& o)
 284 {
 285   //
 286   // Assignment operator
 287   //
 288   fAcceptancePath = o.fAcceptancePath;
 289   fSecMapPath     = o.fSecMapPath;
 290   fAcceptance     = o.fAcceptance;
 291   fSecmap         = o.fSecmap;
 292   fAcceptanceName = o.fAcceptanceName;
 293   fSecMapName     = o.fSecMapName;
 294   fIsInit         = o.fIsInit;
 295   return *this;
 296 }
 297
 298 //____________________________________________________________________
 299 const char*
 300 AliCentralMultiplicityTask::Manager::GetFullFileName(UShort_t what,
 301                                                      UShort_t sys,
 302                                                      UShort_t sNN,
 303                                                      Short_t  field) const
 304 {
 305   //
 306   // Get full path name to object file
 307   //
 308   // Parameters:
 309   //    what   What to get
 310   //    sys    Collision system
 311   //    sNN    Center of mass energy
 312   //    field  Magnetic field
 313   //
 314   // Return:
 315   //
 316   //
 317   return Form("%s/%s",
 318               what == 0 ? GetSecMapPath() : GetAcceptancePath(),
 319               GetFileName(what, sys, sNN, field));
 320 }
 321
 322 //____________________________________________________________________
 323 const char*
 324 AliCentralMultiplicityTask::Manager::GetFileName(UShort_t  what ,
 325                                                  UShort_t  sys,
 326                                                  UShort_t  sNN,
 327                                                  Short_t   field) const
 328 {
 329   //
 330   // Get the full path name
 331   //
 332   // Parameters:
 333   //    what   What to get
 334   //    sys    Collision system
 335   //    sNN    Center of mass energy
 336   //    field  Magnetic field
 337   //
 338   // Return:
 339   //
 340   //
 341   // Must be static - otherwise the data may disappear on return from
 342   // this member function
 343   static TString fname = "";
 344
 345   fname = "";
 346   switch(what) {
 347   case 0:  fname.Append(fSecMapName.Data());     break;
 348   case 1:  fname.Append(fAcceptanceName.Data()); break;
 349   default:
 350     ::Error("GetFileName",
 351             "Invalid indentifier %d for central object, must be 0 or 1!", what);
 352     break;
 353   }
 354   fname.Append(Form("_%s_%04dGeV_%c%1dkG.root",
 355                     AliForwardUtil::CollisionSystemString(sys),
 356                     sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
 357
 358   return fname.Data();
 359 }
 360
 361 //____________________________________________________________________
 362 TH2D*
 363 AliCentralMultiplicityTask::Manager::GetSecMapCorrection(UShort_t vtxbin) const
 364 {
 365   //
 366   // Get the secondary map
 367   //
 368   // Parameters:
 369   //    vtxbin
 370   //
 371   // Return:
 372   //
 373   //
 374   if (!fSecmap) {
 375     ::Warning("GetSecMapCorrection","No secondary map defined");
 376     return 0;
 377   }
 378   return fSecmap->GetCorrection(vtxbin);
 379 }
 380 //____________________________________________________________________
 381 TH1D*
 382 AliCentralMultiplicityTask::Manager::GetAcceptanceCorrection(UShort_t vtxbin)
 383   const
 384 {
 385   //
 386   // Get the acceptance correction
 387   //
 388   // Parameters:
 389   //    vtxbin
 390   //
 391   // Return:
 392   //
 393   //
 394   if (!fAcceptance) {
 395     ::Warning("GetAcceptanceCorrection","No acceptance map defined");
 396     return 0;
 397   }
 398   return fAcceptance->GetCorrection(vtxbin);
 399 }
 400
 401 //____________________________________________________________________
 402 void
 403 AliCentralMultiplicityTask::Manager::Init(UShort_t  sys,
 404                                           UShort_t  sNN,
 405                                           Short_t   field)
 406 {
 407   //
 408   // Initialize
 409   //
 410   // Parameters:
 411   //    sys    Collision system (1: pp, 2: PbPb)
 412   //    sNN    Center of mass energy per nucleon pair [GeV]
 413   //    field  Magnetic field [kG]
 414   //
 415   if(fIsInit) ::Warning("Init","Already initialised - overriding...");
 416
 417   TFile fsec(GetFullFileName(0,sys,sNN,field));
 418   fSecmap =
 419     dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));
 420   if(!fSecmap) {
 421     ::Error("Init", "no central Secondary Map found!") ;
 422     return;
 423   }
 424   TFile facc(GetFullFileName(1,sys,sNN,field));
 425   fAcceptance =
 426     dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
 427   if(!fAcceptance) {
 428     ::Error("Init", "no central Acceptance found!") ;
 429     return;
 430   }
 431
 432   if(fSecmap && fAcceptance) {
 433     fIsInit = kTRUE;
 434     ::Info("Init",
 435            "Central Manager initialised for sys %d, energy %d, field %d",sys,sNN,field);
 436
 437   }
 438
 439 }
 440 //
 441 // EOF
 442 //