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