end-of-line normalization
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / dNdPt / AliPtResolAnalysis.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 // AliPtResolAnalysis 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
44 #include "AlidNdPtEventCuts.h"
45 #include "AlidNdPtAcceptanceCuts.h"
46 #include "AliPhysicsSelection.h"
47 #include "AliTriggerAnalysis.h"
48
49 #include "AliPWG0Helper.h"
50 #include "AlidNdPtHelper.h"
51 #include "AliPtResolAnalysis.h"
52
53 using namespace std;
54
55 ClassImp(AliPtResolAnalysis)
56
57 //_____________________________________________________________________________
58   AliPtResolAnalysis::AliPtResolAnalysis(): AlidNdPt(),
59   fAnalysisFolder(0),
60   fTrackParamHist(0),
61   fTrackParamHist2(0)
62 {
63   // default constructor
64   Init();
65 }
66
67 //_____________________________________________________________________________
68 AliPtResolAnalysis::AliPtResolAnalysis(Char_t* name, Char_t* title): AlidNdPt(name,title),
69   fAnalysisFolder(0),
70   fTrackParamHist(0),
71   fTrackParamHist2(0)
72 {
73   Init();
74 }
75
76 //_____________________________________________________________________________
77 AliPtResolAnalysis::~AliPtResolAnalysis() {
78   //
79   // destructor
80   //
81   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
82   if(fTrackParamHist) delete fTrackParamHist; fTrackParamHist=0;
83   if(fTrackParamHist2) delete fTrackParamHist2; fTrackParamHist2=0;
84 }
85
86 //_____________________________________________________________________________
87 void AliPtResolAnalysis::Init(){
88   //
89   // Generic histograms to be corrected
90   //
91   //1/pT:#sigma(1/pT)
92   Int_t binsTrackParamHist[2]={400,300};
93   Double_t minTrackParamHist[2]={0,0}; 
94   Double_t maxTrackParamHist[2]={1,0.015};
95
96   fTrackParamHist = new THnSparseF("fTrackParamHist","1/pT:#sigma(1/pT)",2,binsTrackParamHist,minTrackParamHist,maxTrackParamHist);
97   fTrackParamHist->GetAxis(0)->SetTitle("1/pT (GeV/c)^{-1}");
98   fTrackParamHist->GetAxis(1)->SetTitle("#sigma(1/pT)");
99   fTrackParamHist->Sumw2();
100   
101   //pt:sigma(1/pT)*pT
102   const Int_t ptNbins = 73;
103   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 };
104
105   Int_t binsTrackParamHist2[2]={ptNbins,200};
106   Double_t minTrackParamHist2[2]={0,0}; 
107   Double_t maxTrackParamHist2[2]={100,0.2};
108
109   fTrackParamHist2 = new THnSparseF("fTrackParamHist2","pT:#sigma(1/pT)*pT",2,binsTrackParamHist2,minTrackParamHist2,maxTrackParamHist2);
110   fTrackParamHist2->SetBinEdges(0,bins);
111   fTrackParamHist2->GetAxis(0)->SetTitle("pT (GeV/c)");
112   fTrackParamHist2->GetAxis(1)->SetTitle("#sigma(1/pT)*pT");
113   fTrackParamHist2->Sumw2();
114
115   // init folder
116   fAnalysisFolder = CreateFolder("folderdNdPt","Analysis dNdPt Folder");
117 }
118
119 //_____________________________________________________________________________
120 void AliPtResolAnalysis::Process(AliESDEvent *const esdEvent, AliMCEvent *const mcEvent)
121 {
122   //
123   // Process real and/or simulated events
124   //
125   if(!esdEvent) {
126     AliDebug(AliLog::kError, "esdEvent not available");
127     return;
128   }
129
130   // get selection cuts
131   AlidNdPtEventCuts *evtCuts = GetEventCuts(); 
132   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); 
133   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); 
134
135   if(!evtCuts || !accCuts  || !esdTrackCuts) {
136     AliDebug(AliLog::kError, "cuts not available");
137     return;
138   }
139
140   // trigger selection
141   Bool_t isEventTriggered = kTRUE;
142   AliPhysicsSelection *physicsSelection = NULL;
143   AliTriggerAnalysis* triggerAnalysis = NULL;
144
145   // 
146   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
147   if (!inputHandler)
148   {
149     Printf("ERROR: Could not receive input handler");
150     return;
151   }
152
153   if(evtCuts->IsTriggerRequired())  
154   {
155     // always MB
156     isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
157
158     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
159     if(!physicsSelection) return;
160     //SetPhysicsTriggerSelection(physicsSelection);
161
162     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
163       // set trigger (V0AND)
164       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
165       if(!triggerAnalysis) return;
166       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
167     }
168
169   // get reconstructed vertex  
170   const AliESDVertex* vtxESD = 0; 
171   Bool_t isRecVertex = kFALSE;
172   if(evtCuts->IsRecVertexRequired()) 
173   {
174     Bool_t bRedoTPCVertex = evtCuts->IsRedoTPCVertex();
175     Bool_t bUseConstraints = evtCuts->IsUseBeamSpotConstraint();
176     vtxESD = AlidNdPtHelper::GetVertex(esdEvent,evtCuts,accCuts,esdTrackCuts,GetAnalysisMode(),kFALSE,bRedoTPCVertex,bUseConstraints); 
177     isRecVertex = AlidNdPtHelper::TestRecVertex(vtxESD, esdEvent->GetPrimaryVertexSPD(), GetAnalysisMode(), kFALSE);
178   }
179
180   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD) && isRecVertex; 
181   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
182   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
183
184   //
185   TObjArray *allChargedTracks=0;
186   // check event cuts
187   if(isEventOK && isEventTriggered)
188   {
189     // get all charged tracks
190     allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,GetAnalysisMode());
191     if(!allChargedTracks) return;
192
193     Int_t entries = allChargedTracks->GetEntries();
194     //printf("entries %d \n",entries);
195
196     // fill histograms
197     for(Int_t i=0; i<entries;++i) 
198     {
199       AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);
200       if(!track) continue;
201       if(track->Charge()==0) continue;
202
203       // only postive charged 
204       if(GetParticleMode() == AlidNdPtHelper::kPlus && track->Charge() < 0) 
205         continue;
206       
207       // only negative charged 
208       if(GetParticleMode() == AlidNdPtHelper::kMinus && track->Charge() > 0) 
209         continue;
210
211       if(esdTrackCuts->AcceptTrack(track)) 
212       {
213         if(accCuts->AcceptTrack(track)) 
214         {
215           //Double_t x, par[5], cov[15];
216           //track->GetExternalParameters(x, p);
217           //track->GetExternalCovariance(cov);
218
219           Double_t v1[2] = {track->OneOverPt(),TMath::Sqrt(track->GetSigma1Pt2())};
220           fTrackParamHist->Fill(v1);
221
222           Double_t v2[2] = {track->Pt(),track->Pt()*TMath::Sqrt(track->GetSigma1Pt2())};
223           fTrackParamHist2->Fill(v2);
224         }
225       }  
226     }
227   }
228  }
229 }
230
231 //_____________________________________________________________________________
232 Long64_t AliPtResolAnalysis::Merge(TCollection* const list) 
233 {
234   // Merge list of objects (needed by PROOF)
235
236   if (!list)
237   return 0;
238
239   if (list->IsEmpty())
240   return 1;
241
242   TIterator* iter = list->MakeIterator();
243   TObject* obj = 0;
244
245   //
246   //TList *collPhysSelection = new TList;
247
248   // collection of generated histograms
249
250   Int_t count=0;
251   while((obj = iter->Next()) != 0) {
252     AliPtResolAnalysis* entry = dynamic_cast<AliPtResolAnalysis*>(obj);
253     if (entry == 0) continue; 
254     
255     //
256     fTrackParamHist->Add(entry->fTrackParamHist);
257     fTrackParamHist2->Add(entry->fTrackParamHist2);
258   }
259
260 return count;
261 }
262
263 //_____________________________________________________________________________
264 void AliPtResolAnalysis::Analyse() 
265 {
266   // Analyse histograms
267   //
268   TH1::AddDirectory(kFALSE);
269   TObjArray *aFolderObj = new TObjArray;
270   if(!aFolderObj) return;
271   
272   //
273   // Reconstructed event vertex
274   //
275   
276   // export objects to analysis folder
277   fAnalysisFolder = ExportToFolder(aFolderObj);
278   if(!fAnalysisFolder) { 
279     if(aFolderObj) delete aFolderObj;
280     return;
281   }
282
283   // delete only TObjArray
284   if(aFolderObj) delete aFolderObj;
285 }
286
287 //_____________________________________________________________________________
288 TFolder* AliPtResolAnalysis::ExportToFolder(TObjArray * const array) 
289 {
290   // recreate folder every time and export objects to new one
291   //
292   if(!array) return NULL;
293
294   AliPtResolAnalysis * comp=this;
295   TFolder *folder = comp->GetAnalysisFolder();
296
297   TString name, title;
298   TFolder *newFolder = 0;
299   Int_t i = 0;
300   Int_t size = array->GetSize();
301
302   if(folder) { 
303      // get name and title from old folder
304      name = folder->GetName();  
305      title = folder->GetTitle();  
306
307          // delete old one
308      delete folder;
309
310          // create new one
311      newFolder = CreateFolder(name.Data(),title.Data());
312      newFolder->SetOwner();
313
314          // add objects to folder
315      while(i < size) {
316            newFolder->Add(array->At(i));
317            i++;
318          }
319   }
320
321 return newFolder;
322 }
323
324 //_____________________________________________________________________________
325 TFolder* AliPtResolAnalysis::CreateFolder(TString name,TString title) { 
326 // create folder for analysed histograms
327 //
328 TFolder *folder = 0;
329   folder = new TFolder(name.Data(),title.Data());
330
331   return folder;
332 }