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