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