PWGLF/SPECTRA/ChargedHadrons/dNdPt/AliPtResolAnalysisPbPb.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 //------------------------------------------------------------------------------
  16 // AliPtResolAnalysisPbPb class.
  17 //
  18 // a. functionality:
  19 // - fills analysis control histograms
  20 //
  21 // b. data members:
  22 // - control histograms
  23 //
  24 // Author: J.Otwinowski 04/11/2008
  25 //------------------------------------------------------------------------------
  26
  27 #include "TH1.h"
  28 #include "TH2.h"
  29 #include "TCanvas.h"
  30 #include "THnSparse.h"
  31
  32 #include "AliHeader.h"
  33 #include "AliInputEventHandler.h"
  34 #include "AliAnalysisManager.h"
  35 #include "AliGenEventHeader.h"
  36 #include "AliStack.h"
  37 #include "AliESDEvent.h"
  38 #include "AliMCEvent.h"
  39 #include "AliESDtrackCuts.h"
  40 #include "AliLog.h"
  41 #include "AliMultiplicity.h"
  42 #include "AliTracker.h"
  43 #include "AliCentrality.h"
  44
  45 #include "AlidNdPtEventCuts.h"
  46 #include "AlidNdPtAcceptanceCuts.h"
  47 #include "AliPhysicsSelection.h"
  48 #include "AliTriggerAnalysis.h"
  49
  50 #include "AliPWG0Helper.h"
  51 #include "AlidNdPtHelper.h"
  52 #include "AliPtResolAnalysisPbPb.h"
  53
  54
  55 using namespace std;
  56
  57 ClassImp(AliPtResolAnalysisPbPb)
  58
  59 //_____________________________________________________________________________
  60   AliPtResolAnalysisPbPb::AliPtResolAnalysisPbPb(): AlidNdPt(),
  61   fAnalysisFolder(0),
  62   fTrackParamHist(0),
  63   fTrackParamHist2(0),
  64   fCentralityEstimator(0)
  65 {
  66   // default constructor
  67   Init();
  68 }
  69
  70 //_____________________________________________________________________________
  71 AliPtResolAnalysisPbPb::AliPtResolAnalysisPbPb(Char_t* name, Char_t* title): AlidNdPt(name,title),
  72   fAnalysisFolder(0),
  73   fTrackParamHist(0),
  74   fTrackParamHist2(0),
  75   fCentralityEstimator(0)
  76 {
  77   Init();
  78 }
  79
  80 //_____________________________________________________________________________
  81 AliPtResolAnalysisPbPb::~AliPtResolAnalysisPbPb() {
  82   //
  83   // destructor
  84   //
  85   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
  86   if(fTrackParamHist) delete fTrackParamHist; fTrackParamHist=0;
  87   if(fTrackParamHist2) delete fTrackParamHist2; fTrackParamHist2=0;
  88 }
  89
  90 //_____________________________________________________________________________
  91 void AliPtResolAnalysisPbPb::Init(){
  92   //
  93   // Generic histograms to be corrected
  94   //
  95   //1/pT:#sigma(1/pT):centrality
  96   Int_t binsTrackParamHist[3]={400,300,11};
  97   Double_t minTrackParamHist[3]={0,0,0};
  98   Double_t maxTrackParamHist[3]={1,0.015,100};
  99
 100   Double_t centrBins[12] = {0.0,5.,10.,20.,30.,40.,50.,60.,70.,80.,90.,100};
 101
 102
 103   fTrackParamHist = new THnSparseF("fTrackParamHist","1/pT:#sigma(1/pT):centrality",3,binsTrackParamHist,minTrackParamHist,maxTrackParamHist);
 104   fTrackParamHist->SetBinEdges(2,centrBins);
 105   fTrackParamHist->GetAxis(0)->SetTitle("1/pT (GeV/c)^{-1}");
 106   fTrackParamHist->GetAxis(1)->SetTitle("#sigma(1/pT)");
 107   fTrackParamHist->GetAxis(2)->SetTitle("centrality");
 108   fTrackParamHist->Sumw2();
 109
 110   //pt:sigma(1/pT)*pT:centrality
 111   const Int_t ptNbins = 73;
 112   Double_t bins[74] = {0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0 };
 113
 114   Int_t binsTrackParamHist2[3]={ptNbins,200,11};
 115   Double_t minTrackParamHist2[3]={0,0,0};
 116   Double_t maxTrackParamHist2[3]={100,0.2,100};
 117
 118   fTrackParamHist2 = new THnSparseF("fTrackParamHist2","pT:#sigma(1/pT)*pT:centrality",3,binsTrackParamHist2,minTrackParamHist2,maxTrackParamHist2);
 119   fTrackParamHist2->SetBinEdges(0,bins);
 120   fTrackParamHist2->GetAxis(0)->SetTitle("pT (GeV/c)");
 121   fTrackParamHist2->GetAxis(1)->SetTitle("#sigma(1/pT)*pT");
 122   fTrackParamHist2->GetAxis(2)->SetTitle("centrality");
 123   fTrackParamHist2->Sumw2();
 124
 125   // init folder
 126   fAnalysisFolder = CreateFolder("folderdNdPt","Analysis dNdPt Folder");
 127
 128 }
 129
 130 //_____________________________________________________________________________
 131 void AliPtResolAnalysisPbPb::Process(AliESDEvent *const esdEvent, AliMCEvent *const mcEvent)
 132 {
 133   //
 134   // Process real and/or simulated events
 135   //
 136   if(!esdEvent) {
 137     AliDebug(AliLog::kError, "esdEvent not available");
 138     return;
 139   }
 140
 141   // get selection cuts
 142   AlidNdPtEventCuts *evtCuts = GetEventCuts();
 143   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts();
 144   AliESDtrackCuts *esdTrackCuts = GetTrackCuts();
 145
 146   if(!evtCuts || !accCuts  || !esdTrackCuts) {
 147     AliDebug(AliLog::kError, "cuts not available");
 148     return;
 149   }
 150
 151   // trigger selection
 152   Bool_t isEventTriggered = kTRUE;
 153   AliPhysicsSelection *physicsSelection = NULL;
 154   AliTriggerAnalysis* triggerAnalysis = NULL;
 155
 156   //
 157   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 158   if (!inputHandler)
 159   {
 160     Printf("ERROR: Could not receive input handler");
 161     return;
 162   }
 163
 164   if(evtCuts->IsTriggerRequired())
 165   {
 166     // always MB
 167     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
 168
 169     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
 170     if(!physicsSelection) return;
 171     //SetPhysicsTriggerSelection(physicsSelection);
 172
 173     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
 174       // set trigger (V0AND)
 175       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
 176       if(!triggerAnalysis) return;
 177       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
 178     }
 179
 180
 181    // centrality determination
 182    Float_t centralityF = -1.;
 183    AliCentrality *esdCentrality = esdEvent->GetCentrality();
 184    centralityF = esdCentrality->GetCentralityPercentile(fCentralityEstimator.Data());
 185
 186   // get reconstructed vertex
 187   const AliESDVertex* vtxESD = 0;
 188   Bool_t isRecVertex = kFALSE;
 189   if(evtCuts->IsRecVertexRequired())
 190   {
 191     Bool_t bRedoTPCVertex = evtCuts->IsRedoTPCVertex();
 192     Bool_t bUseConstraints = evtCuts->IsUseBeamSpotConstraint();
 193     vtxESD = AlidNdPtHelper::GetVertex(esdEvent,evtCuts,accCuts,esdTrackCuts,GetAnalysisMode(),kFALSE,bRedoTPCVertex,bUseConstraints);
 194     isRecVertex = AlidNdPtHelper::TestRecVertex(vtxESD, esdEvent->GetPrimaryVertexSPD(), GetAnalysisMode(), kFALSE);
 195   }
 196
 197   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD) && isRecVertex;
 198   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
 199   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
 200
 201   //
 202   TObjArray *allChargedTracks=0;
 203   // check event cuts
 204   if(isEventOK && isEventTriggered)
 205   {
 206     // get all charged tracks
 207     allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,GetAnalysisMode());
 208     if(!allChargedTracks) return;
 209
 210     Int_t entries = allChargedTracks->GetEntries();
 211     //printf("entries %d \n",entries);
 212
 213     // fill histograms
 214     for(Int_t i=0; i<entries;++i)
 215     {
 216       AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);
 217       if(!track) continue;
 218       if(track->Charge()==0) continue;
 219
 220       // only postive charged
 221       if(GetParticleMode() == AlidNdPtHelper::kPlus && track->Charge() < 0)
 222         continue;
 223
 224       // only negative charged
 225       if(GetParticleMode() == AlidNdPtHelper::kMinus && track->Charge() > 0)
 226         continue;
 227
 228       if(esdTrackCuts->AcceptTrack(track))
 229       {
 230         if(accCuts->AcceptTrack(track))
 231         {
 232           //Double_t x, par[5], cov[15];
 233           //track->GetExternalParameters(x, p);
 234           //track->GetExternalCovariance(cov);
 235
 236           Double_t v[3] = {track->OneOverPt(),TMath::Sqrt(track->GetSigma1Pt2()),centralityF};
 237           fTrackParamHist->Fill(v);
 238
 239           Double_t v2[3] = {track->Pt(),track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),centralityF};
 240           fTrackParamHist2->Fill(v2);
 241         }
 242       }
 243     }
 244   }
 245  }
 246 }
 247
 248 //_____________________________________________________________________________
 249 Long64_t AliPtResolAnalysisPbPb::Merge(TCollection* const list)
 250 {
 251   // Merge list of objects (needed by PROOF)
 252
 253   if (!list)
 254   return 0;
 255
 256   if (list->IsEmpty())
 257   return 1;
 258
 259   TIterator* iter = list->MakeIterator();
 260   TObject* obj = 0;
 261
 262   //
 263   //TList *collPhysSelection = new TList;
 264
 265   // collection of generated histograms
 266
 267   Int_t count=0;
 268   while((obj = iter->Next()) != 0) {
 269     AliPtResolAnalysisPbPb* entry = dynamic_cast<AliPtResolAnalysisPbPb*>(obj);
 270     if (entry == 0) continue;
 271
 272     //
 273     fTrackParamHist->Add(entry->fTrackParamHist);
 274     fTrackParamHist2->Add(entry->fTrackParamHist2);
 275   }
 276
 277 return count;
 278 }
 279
 280 //_____________________________________________________________________________
 281 void AliPtResolAnalysisPbPb::Analyse()
 282 {
 283   // Analyse histograms
 284   //
 285   TH1::AddDirectory(kFALSE);
 286   TObjArray *aFolderObj = new TObjArray;
 287   if(!aFolderObj) return;
 288
 289   //
 290   // Reconstructed event vertex
 291   //
 292
 293   // export objects to analysis folder
 294   fAnalysisFolder = ExportToFolder(aFolderObj);
 295   if(!fAnalysisFolder) {
 296     if(aFolderObj) delete aFolderObj;
 297     return;
 298   }
 299
 300   // delete only TObjArray
 301   if(aFolderObj) delete aFolderObj;
 302 }
 303
 304 //_____________________________________________________________________________
 305 TFolder* AliPtResolAnalysisPbPb::ExportToFolder(TObjArray * const array)
 306 {
 307   // recreate folder avery time and export objects to new one
 308   //
 309   if(!array) return NULL;
 310
 311   AliPtResolAnalysisPbPb * comp=this;
 312   TFolder *folder = comp->GetAnalysisFolder();
 313
 314   TString name, title;
 315   TFolder *newFolder = 0;
 316   Int_t i = 0;
 317   Int_t size = array->GetSize();
 318
 319   if(folder) {
 320      // get name and title from old folder
 321      name = folder->GetName();
 322      title = folder->GetTitle();
 323
 324          // delete old one
 325      delete folder;
 326
 327          // create new one
 328      newFolder = CreateFolder(name.Data(),title.Data());
 329      newFolder->SetOwner();
 330
 331          // add objects to folder
 332      while(i < size) {
 333            newFolder->Add(array->At(i));
 334            i++;
 335          }
 336   }
 337
 338 return newFolder;
 339 }
 340
 341 //_____________________________________________________________________________
 342 TFolder* AliPtResolAnalysisPbPb::CreateFolder(TString name,TString title) {
 343 // create folder for analysed histograms
 344 //
 345 TFolder *folder = 0;
 346   folder = new TFolder(name.Data(),title.Data());
 347
 348   return folder;
 349 }