856dfb93f1e7e6fbb1d7265d05581e502a806ce6
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / AliAnalysisTaskDiHadronPID.cxx
1 // ----------------------------------------------------------------------------
2 // This class makes di-hadron correlations, with TOF and TPC signals for
3 // the associated particles. It runs on AOD049 and AOD73 productions.
4 // ----------------------------------------------------------------------------
5 // Author: Misha Veldhoen (misha.veldhoen@cern.ch)
6 // Start: July 21st, 2011.
7 // Last edit: Mar 2nd 2012. (v 8.00)
8 // ----------------------------------------------------------------------------
9 //
10
11 #include <iostream>
12 #include "TChain.h"
13 #include "TTree.h"
14 #include "TH1F.h"
15 #include "TH2F.h"
16 #include "TH3F.h"
17 #include "TCanvas.h"
18 #include "TFile.h"
19
20 #include "AliAODTrack.h"
21 #include "AliAODEvent.h"
22 #include "AliAODInputHandler.h"
23 #include "AliAODVertex.h"
24 //#include "AliAODPid.h"
25
26 #include "AliAnalysisManager.h"
27 #include "AliInputEventHandler.h"
28 #include "AliPIDResponse.h"
29 #include "AliTPCPIDResponse.h"
30 //#include "AliTOFPIDResponse.h"
31
32 #include "AliAnalysisTaskDiHadronPID.h"
33
34 using namespace std;
35
36 ClassImp(AliAnalysisTaskDiHadronPID);
37
38 //_____________________________________________________________________________
39 AliAnalysisTaskDiHadronPID::AliAnalysisTaskDiHadronPID():
40         AliAnalysisTaskSE(),
41         fPIDResponse(0x0),
42         fAODEvent(0x0),
43         fAODHeader(0x0),
44         fAODVertex(0x0),
45         fAODTrack(0x0),
46         fPIDPartners(0x0),
47         fCentrality(0x0),
48         fVertexZ(0x0),
49     fTrackCuts(0x0),
50         fPtSpectrum(0x0),
51         fAssociatedDistribution(0x0),
52         fTPCnSigmaProton(0x0),
53         fTOFnSigmaProton(0x0),
54         fTPCnSigmaPion(0x0),
55         fTOFnSigmaPion(0x0),
56         fTPCnSigmaKaon(0x0),
57         fTOFnSigmaKaon(0x0),
58         fTPCSignal(0x0),
59         fTOFSignal(0x0),
60         fDiHadron(0x0),
61         fMixedEvents(0x0),
62         fHistoList(0x0),
63         fVerbose(kFALSE),
64         fMask(0),       
65         fCalculateMixedEvents(kFALSE),
66         fTrigBufferIndex(0),    
67         fTrigBufferSize(0)
68
69 {
70         //
71         // Default Constructor.
72         //
73         
74         // Trigger buffer.
75         for(Int_t i=0; i<25000; i++) {
76                 for(Int_t j=0; j<4; j++) {
77                         fTrigBuffer[i][j]=0;
78                 }                               
79         }       
80
81         
82         fMask = 1<<7;
83         
84     // The identified di-hadron correlations.
85         for (Int_t i = 0; i < 3; i++) {
86                 for (Int_t j = 0; j < 10; j++) {
87             fDiHadronTPC[i][j]=0x0;
88             fDiHadronTOF[i][j]=0x0;
89                         fDiHadronTPCTOF[i][j]=0x0;
90                 }
91         }
92     
93 }
94
95 //_____________________________________________________________________________
96 AliAnalysisTaskDiHadronPID::AliAnalysisTaskDiHadronPID(const char *name):
97         AliAnalysisTaskSE(name),
98         fPIDResponse(0x0),
99         fAODEvent(0x0),
100         fAODHeader(0x0),
101         fAODVertex(0x0),
102         fAODTrack(0x0),
103         fPIDPartners(0x0),
104         fCentrality(0x0),
105         fVertexZ(0x0),
106     fTrackCuts(0x0),
107         fPtSpectrum(0x0),
108     fAssociatedDistribution(0x0),
109         fTPCnSigmaProton(0x0),
110         fTOFnSigmaProton(0x0),
111         fTPCnSigmaPion(0x0),
112         fTOFnSigmaPion(0x0),
113         fTPCnSigmaKaon(0x0),
114         fTOFnSigmaKaon(0x0),
115         fTPCSignal(0x0),
116         fTOFSignal(0x0),
117         fDiHadron(0x0),
118         fMixedEvents(0x0),
119         fHistoList(0x0),
120         fVerbose(kFALSE),
121         fMask(0),
122         fCalculateMixedEvents(kFALSE),
123         fTrigBufferIndex(0),    
124         fTrigBufferSize(0)
125
126 {
127         //
128         // Named Constructor.
129         //
130         
131         DefineInput(0, TChain::Class());
132         DefineOutput(1, TList::Class());
133         
134         // Trigger buffer.
135         for(Int_t i=0; i<25000; i++) {
136                 for(Int_t j=0; j<4; j++) {
137                         fTrigBuffer[i][j]=0;
138                 }                               
139         }       
140         
141         fMask = 1<<7;
142     
143     // The identified di-hadron correlations.
144         for (Int_t i = 0; i < 3; i++) {
145                 for (Int_t j = 0; j < 10; j++) {
146             fDiHadronTPC[i][j]=0x0;
147             fDiHadronTOF[i][j]=0x0;
148                         fDiHadronTPCTOF[i][j]=0x0;
149                 }
150         }
151
152 }
153
154 //_____________________________________________________________________________
155 AliAnalysisTaskDiHadronPID::~AliAnalysisTaskDiHadronPID() {
156
157         //
158         // Destructor.
159         //
160         
161     if(fPIDPartners) {
162         delete fPIDPartners;
163         fPIDPartners=0;
164     }
165
166 }
167
168 //_____________________________________________________________________________
169 void AliAnalysisTaskDiHadronPID::UserCreateOutputObjects() 
170
171 {
172         //
173         // Output objects.
174         //
175     
176         // The Analysis Manager.
177         AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
178         AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*> (manager->GetInputEventHandler());
179
180         // Pointers to PID Response objects.    
181         fPIDResponse = inputHandler->GetPIDResponse(); 
182         cout << "PID Response object: " << fPIDResponse << endl;
183
184         // Create the output of the task.       
185         fHistoList = new TList();
186         fHistoList->SetOwner(kTRUE); 
187         
188         // Ranges in dPhi, dEta, and the number of bins for di-hadron correlations.
189         Int_t binsHisto[4] = {36,25};
190         Double_t minHisto[4] = {-TMath::Pi()/2.,-1.8};
191         Double_t maxHisto[4] = {3.*TMath::Pi()/2,1.8};
192         
193     // --- EVENT SAMPLE PLOTS (PER EVENT) ---
194     
195         // Centrality Histogram.
196         fCentrality = new TH1F("fCentrality","Centrality;Centrality;N_{evt}",10,0,100);
197         fHistoList->Add(fCentrality);
198         
199         // Vertex Z Histogram.
200         fVertexZ = new TH1F("fVertexZ","Vertex Z position;z (cm);N_{evt}",30,-15,15);
201         fHistoList->Add(fVertexZ);
202         
203         // --- UNIDENTIFIED SPECTRA ---
204         
205         fPtSpectrum = new TH1F("fPtSpectrum","p_{T} Spectrum Unidentified;p_{T}",100,0,50);
206         fHistoList->Add(fPtSpectrum);
207         
208         fAssociatedDistribution = new TH3F("fAssociatedDistribution","Associated Distribution;#phi;#eta;p_{T_assoc}",binsHisto[0],0.,2.*TMath::Pi(),binsHisto[1],-0.9,0.9,10,0.,5.);
209     fAssociatedDistribution->Sumw2();
210     fHistoList->Add(fAssociatedDistribution);
211     
212     // --- TRACK CUTS ---
213     
214     fTrackCuts = new TH1F("fTrackCuts","Track Cuts;Cut;Count",4,0,4);
215     (fTrackCuts->GetXaxis())->SetBinLabel(1,"Standard Cuts");
216     (fTrackCuts->GetXaxis())->SetBinLabel(2,"+ TPCpid");
217     (fTrackCuts->GetXaxis())->SetBinLabel(3,"+ TOFpid");
218     (fTrackCuts->GetXaxis())->SetBinLabel(4,"+ no TOFmismatch");
219     fHistoList->Add(fTrackCuts);
220         
221     // --- QA PLOTS PID ---
222     
223         fTPCnSigmaProton = new TH2F("fTPCnSigmaProton","n#sigma plot for the TPC (Protons);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
224         fHistoList->Add(fTPCnSigmaProton);
225         
226         fTOFnSigmaProton = new TH2F("fTOFnSigmaProton","n#sigma plot for the TOF (Protons);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
227         fHistoList->Add(fTOFnSigmaProton);
228
229         fTPCnSigmaPion = new TH2F("fTPCnSigmaPion","n#sigma plot for the TPC (Pions);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
230         fHistoList->Add(fTPCnSigmaPion);
231         
232         fTOFnSigmaPion = new TH2F("fTOFnSigmaPion","n#sigma plot for the TOF (Pions);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
233         fHistoList->Add(fTOFnSigmaPion);
234         
235         fTPCnSigmaKaon = new TH2F("fTPCnSigmaKaon","n#sigma plot for the TPC (Kaons);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
236         fHistoList->Add(fTPCnSigmaKaon);
237         
238         fTOFnSigmaKaon = new TH2F("fTOFnSigmaKaon","n#sigma plot for the TOF (Kaons);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
239         fHistoList->Add(fTOFnSigmaKaon);        
240         
241         // --- PID SIGNALS ---
242         
243         fTPCSignal = new TH3F("fTPCSignal","TPC Signal;#eta;p_{T} GeV/c;dE/dx",25,-.9,.9,100.,0.,5.,150,0.,300.);
244         fHistoList->Add(fTPCSignal);
245         
246         fTOFSignal = new TH3F("fTOFSignal","TOF Signal;#eta;p_{T} GeV/c;t",25,-.9,.9,100.,0.,5.,150,10000.,25000.);
247         fHistoList->Add(fTOFSignal);
248         
249         // --- UNIDENTIFIED DI-HADRON CORRELATIONS & MIXED EVENTS ---
250                 
251         // Di Hadron Correlations, unidentified. (Dphi,Deta,p_T_assoc)
252         fDiHadron = new TH3F("fDiHadron","Di-Hadron Correlations;#Delta#phi;#Delta#eta;p_{T,assoc}",binsHisto[0],minHisto[0],maxHisto[0],binsHisto[1],minHisto[1],maxHisto[1],10,0.,5.);
253         fHistoList->Add(fDiHadron);
254                 
255     if (fCalculateMixedEvents) {
256         fMixedEvents = new TH3F("fMixedEvents","Mixed Events;#Delta#phi;#Delta#eta;p_{T_assoc}",binsHisto[0],minHisto[0],maxHisto[0],binsHisto[1],minHisto[1],maxHisto[1],10,0.,5.);
257         fMixedEvents->Sumw2();
258         fHistoList->Add(fMixedEvents);
259     }
260                 
261     // --- DI-HADRON CORRELATIONS WITH TPC AND TOF SIGNALS ---
262     
263         // Di Hadron Correlations with two PID signals, for each p_T bin and particle species separately. (i.e. 30 histo's) 
264         // Axes: {Dphi, Deta, TPC signal, TOF signal}
265     TString basenameTPC("fDiHadronTPC");
266     TString basenameTOF("fDiHadronTOF");
267         TString basenameTPCTOF("fDiHadronTPCTOF");
268         TString basetitle("Di-Hadron correlation");
269         TString finalname, finaltitle;
270         TString species[3] = {"Pion","Kaon","Proton"};
271         TString ptbins[11] = {"0.0","0.5","1.0","1.5","2.0","2.5","3.0","3.5","4.0","4.5","5.0"};
272     
273     // Unzoomed pictures.
274         /*
275         Int_t binsTPC[3][10] = {{100,100,100,100,100,100,100,100,100,100},
276                             {100,100,100,100,100,100,100,100,100,100},
277                             {100,100,100,100,100,100,100,100,100,100}};
278                 
279         Double_t minTPC[3][10] =    {{   -50., -50.,-30.,-30.,-30.,-30.,-30.,-30.,-30.,-30.},
280                                  {      -100., -50.,-20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.},
281                                  {  -200.,-150.,-50.,-30.,-20.,-20.,-20.,-20.,-20.,-20.}};
282         Double_t maxTPC[3][10] =    {{ 150., 100., 50., 30., 25., 25., 25., 25., 25., 25.},
283                                  {  50., 100., 40., 30., 30., 30., 30., 30., 30., 30.},
284                                  { 100.,  50., 50., 30., 30., 30., 30., 30., 30., 30.}};
285         
286         Int_t binsTOF[3][10] = {{100,100,100,100,100,100,100,100,100,100},
287                             {100,100,100,100,100,100,100,100,100,100},
288                             {100,100,100,100,100,100,100,100,100,100}};
289         
290         Double_t minTOF[3][10] =    {{-2000.,-2000.,-1000., -500., -500., -500., -500., -500., -500., -500.},
291                                  {-2000.,-2000.,-2000.,-1000.,-1000.,-1000.,-1000.,-1000.,-1000.,-1000.},
292                                  {-2000.,-2000.,-6000.,-3000.,-2000.,-1500.,-1500.,-1500.,-1500.,-1500.}};
293         Double_t maxTOF[3][10] =    {{ 2000., 2000., 5000., 3000., 2000., 1500., 1500., 1500., 1500., 1500.},
294                                  { 2000., 2000., 5000., 2000., 1500., 1500., 1500., 1500., 1500., 1500.},
295                                  { 2000., 2000., 2000., 1000., 1000., 1000., 1000., 1000., 1000., 1000.}};
296         */
297     
298     // Zoomed pictures.
299     Int_t binsTPC[3][10] = {{100,100,100,100,100,100,100,100,100,100},
300                             {100,100,100,100,100,100,100,100,100,100},
301                             {100,100,100,100,100,100,100,100,100,100}};
302     
303         Double_t minTPC[3][10] =    {{   -20., -20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.},
304                                  {       -20., -20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.},
305                                  {   -40., -20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.}};
306         Double_t maxTPC[3][10] =    {{  20.,  20., 20., 20., 20., 20., 20., 20., 20., 20.},
307                                  {  20.,  20., 20., 20., 20., 20., 20., 20., 20., 20.},
308                                  {  40.,  20., 20., 20., 20., 30., 30., 30., 30., 30.}};
309         
310         Int_t binsTOF[3][10] = {{100,100,100,100,100,100,100,100,100,100},
311                             {100,100,100,100,100,100,100,100,100,100},
312                             {100,100,100,100,100,100,100,100,100,100}};
313         
314         Double_t minTOF[3][10] =    {{-1000.,-1000.,-500.,  -500., -500., -400., -400., -400., -400., -400.},
315                                  { -800., -800., -800., -800., -800., -600., -500., -500., -400., -400.},
316                                  {-1000.,-1000.,-1000.,-1000., -800.,-1000.,-1000., -800., -700., -700.}};
317         Double_t maxTOF[3][10] =    {{ 1000., 1000., 1000., 1000., 1000., 1000., 1000.,  900.,  800.,  700.},
318                                  { 1000., 1000.,  500.,  500.,  500.,  900.,  700.,  700.,  600.,  500.},
319                                  { 1000., 1000., 1000.,  500.,  500.,  500.,  500.,  500.,  500.,  500.}};
320     
321         // Recall that AliPID::kPion = 2, AliPID::kKaon = 3, AliPID::kProton = 4.
322         for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
323                                 
324                 for (Int_t iPtBin = 0; iPtBin < 10; iPtBin++) {
325                         
326                         // setting the variables for each histogram.
327                         finaltitle = basetitle;
328                         (((((finaltitle += " (") += species[iSpecies]) += ") ") += ptbins[iPtBin]) += " < P_t < ") += ptbins[iPtBin+1];
329             finaltitle+=";#Delta#phi;#Delta#eta;TPC signal;TOF signal;";
330             
331                         binsHisto[2] = binsTPC[iSpecies][iPtBin];
332                         binsHisto[3] = binsTOF[iSpecies][iPtBin];
333                         minHisto[2] = minTPC[iSpecies][iPtBin];
334                         minHisto[3] = minTOF[iSpecies][iPtBin];
335                         maxHisto[2] = maxTPC[iSpecies][iPtBin];
336                         maxHisto[3] = maxTOF[iSpecies][iPtBin];
337                 
338             // Make the di-hadron correlations with different pid signals.
339             finalname = basenameTPC;
340                         (((finalname += "_") += species[iSpecies]) += "_") += iPtBin;
341             fDiHadronTPC[iSpecies][iPtBin] = new TH3F(finalname,finaltitle,binsHisto[0],minHisto[0],maxHisto[0],binsHisto[1],minHisto[1],maxHisto[1],binsHisto[2],minHisto[2],maxHisto[2]);
342             fHistoList->Add(fDiHadronTPC[iSpecies][iPtBin]);
343
344             finalname = basenameTOF;
345                         (((finalname += "_") += species[iSpecies]) += "_") += iPtBin;
346             fDiHadronTOF[iSpecies][iPtBin] = new TH3F(finalname,finaltitle,binsHisto[0],minHisto[0],maxHisto[0],binsHisto[1],minHisto[1],maxHisto[1],binsHisto[3],minHisto[3],maxHisto[3]);
347             fHistoList->Add(fDiHadronTOF[iSpecies][iPtBin]);
348
349             finalname = basenameTPCTOF;
350                         (((finalname += "_") += species[iSpecies]) += "_") += iPtBin;
351             fDiHadronTPCTOF[iSpecies][iPtBin] = new THnSparseF(finalname,finaltitle,4,binsHisto,minHisto,maxHisto);
352             fHistoList->Add(fDiHadronTPCTOF[iSpecies][iPtBin]);
353             
354                 }
355         }
356     
357         PostData(1, fHistoList);
358         
359 }
360
361 //_____________________________________________________________________________
362 Bool_t AliAnalysisTaskDiHadronPID::SelectTrack(AliAODTrack *track, Int_t cuts)
363
364 {
365         
366     // This selects tracks with three different kinds of track cuts:
367     //
368     //  Case 0: 
369     //   - Tracks must pass the filterbit,
370     //   - Tracks must have eta < 0.9.
371     //  Case 1:
372     //   - TPCpid hit is demanded.
373     //  Case 2:
374     //   - TOFpid hit is demanded.
375     //  Case 3:
376     //   - no TOFmismatch is demanded.
377     //
378     
379     ULong_t status;
380     
381     switch (cuts) {
382         case 0:
383             //cout<<"test0"<<endl;
384             if (!(track->TestFilterMask(fMask))) {return kFALSE;}
385             if (!(TMath::Abs(track->Eta())<0.9)) {return kFALSE;}
386             break;
387         case 1:
388             //cout<<"test1"<<endl;
389             status=GetTrackPartner(track)->GetStatus();  
390             if (!((status&AliAODTrack::kTPCpid)==AliAODTrack::kTPCpid)) {return kFALSE;}
391             break;
392         case 2:
393             //cout<<"test2"<<endl;
394             status=GetTrackPartner(track)->GetStatus();  
395             if (!((status&AliAODTrack::kTOFpid)==AliAODTrack::kTOFpid)) {return kFALSE;}
396             break;
397         case 3:
398             //cout<<"test3"<<endl;
399             status=GetTrackPartner(track)->GetStatus();  
400             if ((status&AliAODTrack::kTOFmismatch)==AliAODTrack::kTOFmismatch) {return kFALSE;}
401             break;
402         default:
403             return kFALSE;
404             break;
405     }
406
407     return kTRUE;
408     
409 }
410
411 //____________________________________________________________________________
412 Bool_t AliAnalysisTaskDiHadronPID::SelectEvent(AliAODVertex* vertex)
413
414 {
415         
416         //
417         // Event Selection.
418         //
419         
420         Double_t primVtx[3];
421         vertex->GetXYZ(primVtx);
422         if (TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2])>10.) {
423                 if (fVerbose) cout << "AliAnalysisTaskDiHadronPID::SelectEvent: Vertex Out of Range." << endl;
424         if (fVerbose) cout << "AliAnalysisTaskDiHadronPID::SelectEvent: Event not selected." << endl;
425                 return kFALSE;
426         }
427         if (fVerbose) cout << "AliAnalysisTaskDiHadronPID::SelectEvent: Vertex is OK." << endl;
428     
429     // We also wish to make a 0-10% centrality cut.
430     if (fAODHeader->GetCentrality()>10.) {
431         if (fVerbose) cout << "AliAnalysisTaskDiHadronPID::SelectEvent: Non-central event." << endl;
432         if (fVerbose) cout << "AliAnalysisTaskDiHadronPID::SelectEvent: Event not selected." << endl;
433                 return kFALSE;
434     }
435         if (fVerbose) cout << "AliAnalysisTaskDiHadronPID::SelectEvent: Central Event." << endl;
436     
437     if (fVerbose) cout << "AliAnalysisTaskDiHadronPID::SelectEvent: Event selected." << endl;
438         return kTRUE;
439         
440 }
441
442 //_____________________________________________________________________________
443 void AliAnalysisTaskDiHadronPID::FillPIDPartnersArray() {
444
445         // Initialize the mapping for corresponding PID tracks. (see
446         // GetTrackPartner(AliAODTrack* track)). 
447         //
448         
449         if (!fAODEvent) {
450         cout << "ERROR in CreatePIDPartersArray(): fAODEvent not set." << endl;
451                 return;
452         }
453         
454         if (!fMask) {
455                 cout << "ERROR in CreatePIDPartersArray(): fMask not set." << endl;
456                 return;
457         }
458         
459         fPIDPartners = new TObjArray();
460         AliAODTrack* track = 0x0;
461                 
462         for (Int_t iTrack = 0; iTrack < fAODEvent->GetNumberOfTracks(); iTrack++) {
463                 
464                 track = fAODEvent->GetTrack(iTrack);
465                 
466                 // cout << "Track: " << iTrack << " FilterMaskPass: "<< track->TestFilterMask(fMask) << " Track ID: " << track->GetID() << endl;
467                 
468                 // I.e., if it does NOT pass the filtermask.
469                 if (!(track->TestFilterMask(fMask))) {
470             fPIDPartners->AddAtAndExpand(track,track->GetID());
471             if (track->GetID()<1) cout<<"Track ID: "<<track->GetID()<<" Partner ID: "<<(-track->GetID()-1)<<endl;
472                 }
473         
474         }
475         
476 }
477
478 //_____________________________________________________________________________
479 AliAODTrack* AliAnalysisTaskDiHadronPID::GetTrackPartner(AliAODTrack* track) {
480         
481         //
482         // Returns the "parner track" of track, which contains the pid information.
483         //
484         
485         AliAODTrack* partner = 0x0;
486             
487     partner = (AliAODTrack*)(fPIDPartners->At(-track->GetID()-1));
488             
489         if (!partner&&fVerbose) cout<<"GetTrackPartner: No Partner found!"<<endl;
490         
491         return partner;
492         
493 }
494
495 //_____________________________________________________________________________
496 Double_t AliAnalysisTaskDiHadronPID::PhiRange(Double_t DPhi)
497
498 {
499         //
500         // Puts the argument in the range [-pi/2,3 pi/2].
501         //
502         
503         if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();
504         if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi();      
505
506         return DPhi;
507         
508 }
509
510 //_____________________________________________________________________________
511 void AliAnalysisTaskDiHadronPID::UserExec(Option_t *)
512
513 {
514         //
515         // UserExec.
516         //
517         
518         // Input the event.
519         fAODEvent = dynamic_cast<AliAODEvent*>(InputEvent());
520         
521         if (!fAODEvent) {
522                 cout << "ERROR: No AliAODEvent pointer could be created." << endl;
523                 return;
524         }
525         
526         // Get the event header.
527         fAODHeader = fAODEvent->GetHeader();
528         
529         if (!fAODHeader) {
530                 cout << "ERROR: No AliAODHeader pointer could be created."<<endl;
531                 return;
532         }
533         
534         // Get the event vertex.
535         fAODVertex = fAODEvent->GetPrimaryVertex();
536         
537         if (!fAODVertex) {
538                 cout << "ERROR: No AliAODVertex pointer could be created." << endl;
539                 return;
540         }
541                 
542         // Display basic event information.
543         if (fVerbose) cout << endl;
544         if (fVerbose) cout << "Event centrality: " << fAODHeader->GetCentrality() << endl;
545         if (fVerbose) cout << "Event Vertex Z: " << fAODVertex->GetZ() << endl;
546         if (fVerbose) cout << "Event tracks in AOD: " << fAODEvent->GetNumberOfTracks() << endl;
547         if (fVerbose) cout << endl;
548
549     if (!SelectEvent(fAODVertex)) return;
550         // Fill the TObjArray which holds PID partners.
551         FillPIDPartnersArray();
552         
553         // Filling Centrality/ VertexZ Hisogram.
554         fCentrality->Fill(fAODHeader->GetCentrality());
555         fVertexZ->Fill(fAODVertex->GetZ());
556         
557         TObjArray *triggers             = new TObjArray();
558         TObjArray *associateds  = new TObjArray();
559         
560         // 1. Identifying sets of triggers and associateds. OK!
561         // 2. Filling Pt distribution (UnID). OK!
562         // 3. Making the nSigma plots for TPC and TOF.
563             
564         for (Int_t iTrack = 0; iTrack < fAODEvent->GetNumberOfTracks(); iTrack++) {
565         
566                 fAODTrack = fAODEvent->GetTrack(iTrack);
567
568         // Skip track if there is no pointer, or if it doesn't pass the strandard track cuts.
569                 if ((!fAODTrack)||(!SelectTrack(fAODTrack,0))) continue;
570         
571         if (fAODTrack->GetID()>-1) cout<<"Track ID: "<<fAODTrack->GetID()<<endl;
572         
573         // Make the pT spectrum with only standard track cuts.
574         fPtSpectrum->Fill(fAODTrack->Pt());  
575         
576                 if (fAODTrack->Pt() > 5.0) {
577                         if (fVerbose) cout << "AliAnalysisTaskDiHadronPID::UserExec: Trigger found!" << endl;
578                         triggers->AddLast(fAODTrack);
579                 }
580                 
581                 if (fAODTrack->Pt() < 5.0) {
582         
583             fTrackCuts->AddBinContent(1);
584         
585             // Now we demand a TPC hit.            
586             if (!SelectTrack(fAODTrack,1)) continue;
587             //cout<<"Passed Track cuts 1"<<endl;
588             fTrackCuts->AddBinContent(2);
589                         
590                         // the associateds array contains tracks pT < 5.0 GeV/c && TPC hit.
591                         associateds->AddLast(fAODTrack);
592             if (fVerbose&&(fAODTrack->GetID()>-1)) cout<<"Assoc. Track ID: "<<fAODTrack->GetID()<<endl;
593
594                         // Make the nSigma plots.
595                         Double_t mom, nSigma;
596
597                         mom = GetTrackPartner(fAODTrack)->GetTPCmomentum();
598                         nSigma = fPIDResponse->NumberOfSigmasTPC(GetTrackPartner(fAODTrack),AliPID::kProton);
599                         fTPCnSigmaProton->Fill(mom,nSigma);
600                         nSigma = fPIDResponse->NumberOfSigmasTPC(GetTrackPartner(fAODTrack),AliPID::kPion);
601                         fTPCnSigmaPion->Fill(mom,nSigma);
602             nSigma = fPIDResponse->NumberOfSigmasTPC(GetTrackPartner(fAODTrack),AliPID::kKaon);
603                         fTPCnSigmaKaon->Fill(mom,nSigma);
604
605                         fTPCSignal->Fill(fAODTrack->Eta(),fAODTrack->Pt(),GetTrackPartner(fAODTrack)->GetTPCsignal());
606                         
607             if (!SelectTrack(fAODTrack,2)) continue;
608             fTrackCuts->AddBinContent(3);
609             if (!SelectTrack(fAODTrack,3)) continue;
610             fTrackCuts->AddBinContent(4);
611             
612                         mom = GetTrackPartner(fAODTrack)->P();
613                         nSigma = fPIDResponse->NumberOfSigmasTOF(GetTrackPartner(fAODTrack),AliPID::kProton);
614                         fTOFnSigmaProton->Fill(mom,nSigma);                     
615             nSigma = fPIDResponse->NumberOfSigmasTOF(GetTrackPartner(fAODTrack),AliPID::kPion);
616                         fTOFnSigmaPion->Fill(mom,nSigma);       
617             nSigma = fPIDResponse->NumberOfSigmasTOF(GetTrackPartner(fAODTrack),AliPID::kKaon);
618                         fTOFnSigmaKaon->Fill(mom,nSigma);
619             
620                         fTOFSignal->Fill(fAODTrack->Eta(),fAODTrack->Pt(),GetTrackPartner(fAODTrack)->GetTOFsignal());
621
622                 }
623                                 
624         }
625         
626     //cout<<"Done with the first loop."<<endl;
627     
628         // 1. Making the di-hadron correlations. (to do: set this up a bit nicer!)
629     // {Dphi, Deta, TPC signal, TOF signal}
630         Double_t histoFill[4];
631         AliAODTrack* currentTrigger = 0x0;
632         AliAODTrack* currentAssociated = 0x0;
633         AliAODTrack* currentPIDPartner = 0x0;
634         //AliAODPid* currentPIDObject = 0x0;
635         
636         AliTPCPIDResponse& TPCPIDResponse = fPIDResponse->GetTPCResponse();
637         //AliTOFPIDResponse& TOFPIDResponse = fPIDResponse->GetTOFResponse();
638         
639         for (Int_t iTrig = 0; iTrig < triggers->GetEntriesFast(); iTrig++){
640         
641                 currentTrigger = (AliAODTrack*)(triggers->At(iTrig));
642                 
643                 for (Int_t iAssoc = 0; iAssoc < associateds->GetEntriesFast(); iAssoc++) {
644                 
645                         currentAssociated = (AliAODTrack*)(associateds->At(iAssoc));
646                                                 
647                         histoFill[0] = PhiRange(currentTrigger->Phi() - currentAssociated->Phi());
648                         histoFill[1] = currentTrigger->Eta() - currentAssociated->Eta();
649                         Double_t pt = currentAssociated->Pt();
650
651                         // Be aware that there may be a caveat here when Pt = 5.00000000
652                         const Int_t ptbin = (Int_t)(2*currentAssociated->Pt());
653                         
654             fDiHadron->Fill(histoFill[0],histoFill[1],pt);
655                         
656                         currentPIDPartner = GetTrackPartner(currentAssociated);
657             //currentPIDObject = currentPIDPartner->GetDetPid();
658                         
659                         if (currentPIDPartner/*&&currentPIDObject*/) {
660                                 
661                                 Double_t TPCmom = currentPIDPartner->GetTPCmomentum();
662                 Double_t TPCsignal = currentPIDPartner->GetTPCsignal();
663                 Double_t TOFsignal = -999.;
664                 Double_t expectedTOFsignalKaon=0,expectedTOFsignalPion=0,expectedTOFsignalProton=0;
665                 Double_t times[AliPID::kSPECIES]; // For the expected time per particle species. 
666             
667                                 if (SelectTrack(currentAssociated,2)) {
668                     TOFsignal = currentPIDPartner->GetTOFsignal();
669                     
670                     //currentPIDObject->GetIntegratedTimes(times);
671                     currentPIDPartner->GetIntegratedTimes(times);
672                     
673                     expectedTOFsignalPion = times[AliPID::kPion];
674                     expectedTOFsignalKaon = times[AliPID::kKaon];
675                     expectedTOFsignalProton = times[AliPID::kProton]; 
676                 }
677                 
678                                 Double_t expectedTPCsignalPion = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kPion);
679                 Double_t expectedTPCsignalKaon = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kKaon);
680                                 Double_t expectedTPCsignalProton = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kProton);
681                                 
682                 histoFill[2] = TPCsignal - expectedTPCsignalPion;
683                 fDiHadronTPC[0][ptbin]->Fill(histoFill[0],histoFill[1],histoFill[2]);                
684                 if (SelectTrack(currentAssociated,2)&&SelectTrack(currentAssociated,3)) {
685                     histoFill[3] = TOFsignal - expectedTOFsignalPion;
686                     fDiHadronTOF[0][ptbin]->Fill(histoFill[0],histoFill[1],histoFill[3]);
687                     fDiHadronTPCTOF[0][ptbin]->Fill(histoFill);
688                 }
689
690                 histoFill[2] = TPCsignal - expectedTPCsignalKaon;
691                 fDiHadronTPC[1][ptbin]->Fill(histoFill[0],histoFill[1],histoFill[2]);                
692                 if (SelectTrack(currentAssociated,2)&&SelectTrack(currentAssociated,3)) {
693                     histoFill[3] = TOFsignal - expectedTOFsignalKaon;
694                     fDiHadronTOF[1][ptbin]->Fill(histoFill[0],histoFill[1],histoFill[3]);
695                     fDiHadronTPCTOF[1][ptbin]->Fill(histoFill);
696                 }
697                 
698                 histoFill[2] = TPCsignal - expectedTPCsignalProton;
699                 fDiHadronTPC[2][ptbin]->Fill(histoFill[0],histoFill[1],histoFill[2]);                
700                 if (SelectTrack(currentAssociated,2)&&SelectTrack(currentAssociated,3)) {
701                     histoFill[3] = TOFsignal - expectedTOFsignalProton;
702                     fDiHadronTOF[2][ptbin]->Fill(histoFill[0],histoFill[1],histoFill[3]);
703                     fDiHadronTPCTOF[2][ptbin]->Fill(histoFill);
704                 }
705                                 
706                 fAssociatedDistribution->Fill(currentAssociated->Phi(),currentAssociated->Eta(),currentAssociated->Pt());
707                         }
708                 }
709         }
710         
711     if (fCalculateMixedEvents) {
712         
713         
714         // Loop over the trigger buffer.
715         if (fVerbose) cout << "AliAnalysisTaskDiHadronPID::UserExec: Mixing the events with "<<fTrigBufferSize<<" triggers from the buffer." <<endl;
716         if (fVerbose) cout << "AliAnalysisTaskDiHadronPID::UserExec: Buffer size: "<<fTrigBufferIndex<<endl;
717         
718         for (Int_t iTrig=0;iTrig<fTrigBufferSize;iTrig++) {
719             
720             // Check if the trigger and the associated have a reconstructed
721             // vertext no further than 2cm apart.
722             
723             // fTrigBuffer[i][0] = z
724             // fTrigBuffer[i][1] = phi
725             // fTrigBuffer[i][2] = eta
726             // fTrigBuffer[i][3] = p_t
727             
728             if (TMath::Abs(fTrigBuffer[iTrig][0]-fAODVertex->GetZ())<2.) {
729                 
730                 if (fVerbose) cout<<"AliAnalysisTaskDiHadronPID::UserExec: Mixing with trigger Z: "<<fTrigBuffer[iTrig][0]<<", Pt: "<<fTrigBuffer[iTrig][3]<<endl;
731                 
732                 for (Int_t iAssoc = 0; iAssoc < associateds->GetEntriesFast(); iAssoc++) {
733                     
734                     currentAssociated = (AliAODTrack*)(associateds->At(iAssoc));
735                     currentPIDPartner = GetTrackPartner(currentAssociated);
736                     //currentPIDObject = currentPIDPartner->GetDetPid();
737                     
738                     if (currentPIDPartner/*&&currentPIDObject*/) {
739                         
740                         Double_t DPhi = PhiRange(fTrigBuffer[iTrig][1] - currentAssociated->Phi());
741                         Double_t DEta = fTrigBuffer[iTrig][2] - currentAssociated->Eta();
742                         Double_t Ptassoc = currentAssociated->Pt();
743                         
744                         fMixedEvents->Fill(DPhi,DEta,Ptassoc);
745                     }
746                 }
747             }
748         }
749     
750         // Copy the triggers from the current event into the buffer.
751         if (fAODVertex->GetZ()<10.) {
752         
753             if (fVerbose) cout<<"AliAnalysisTaskDiHadronPID::UserExec: Copying "<<triggers->GetEntriesFast()<<" triggers to the buffer with vertex z = "<<fAODVertex->GetZ()<<endl;
754         
755             for (Int_t iTrig = 0; iTrig<triggers->GetEntriesFast(); iTrig++) {
756             
757                 currentTrigger = (AliAODTrack*)(triggers->At(iTrig));
758                 if (fVerbose) cout<<"AliAnalysisTaskDiHadronPID::UserExec: Trigger pt = "<<currentTrigger->Pt()<<endl;
759             
760                 fTrigBuffer[fTrigBufferIndex][0] = fAODVertex->GetZ();
761                 fTrigBuffer[fTrigBufferIndex][1] = currentTrigger->Phi();
762                 fTrigBuffer[fTrigBufferIndex][2] = currentTrigger->Eta();
763                 fTrigBuffer[fTrigBufferIndex][3] = currentTrigger->Pt();
764                 fTrigBufferIndex++;
765                 if (fTrigBufferSize<25000) {fTrigBufferSize++;}
766                 if (fTrigBufferIndex==25000) {fTrigBufferIndex=0;}
767             }
768         }
769     }        
770         
771     if (fVerbose) cout<<"AliAnalysisTaskDiHadronPID::UserExec: Trigger buffer index: "<<fTrigBufferIndex<<", and size: "<<fTrigBufferSize<<endl;
772     
773         delete triggers;
774         delete associateds;
775          
776         PostData(1,fHistoList);
777         
778 }
779
780 //_____________________________________________________________________________
781 void AliAnalysisTaskDiHadronPID::Terminate(Option_t *)
782
783 {
784         //
785         // Terminate.
786     //
787     
788 }
789