remove obsolate classes
[u/mrichter/AliRoot.git] / PWG1 / ITS / AliAnalysisTaskITSsaTracks.cxx
1 #include "AliAnalysisTaskSE.h"
2 #include "AliTrackPointArray.h"
3 #include "AliAnalysisManager.h"
4 #include "AliAnalysisDataContainer.h"
5 #include "AliESDEvent.h"
6 #include "AliESDfriend.h"
7 #include "AliStack.h"
8 #include "AliPID.h"
9 #include "AliITSPIDResponse.h"
10 #include "AliMCEventHandler.h"
11 #include "AliMCEvent.h"
12 #include <TParticle.h>
13 #include <TSystem.h>
14 #include <TTree.h>
15 #include <TNtuple.h>
16 #include <TH1F.h>
17 #include <TH2F.h>
18 #include <TChain.h>
19 #include "AliESDInputHandlerRP.h"
20 #include "AliAnalysisTaskITSsaTracks.h"
21
22 /**************************************************************************
23  * Copyright(c) 1998-2012, ALICE Experiment at CERN, All rights reserved. *
24  *                                                                        *
25  * Author: The ALICE Off-line Project.                                    *
26  * Contributors are mentioned in the code where appropriate.              *
27  *                                                                        *
28  * Permission to use, copy, modify and distribute this software and its   *
29  * documentation strictly for non-commercial purposes is hereby granted   *
30  * without fee, provided that the above copyright notice appears in all   *
31  * copies and that both the copyright notice and this permission notice   *
32  * appear in the supporting documentation. The authors make no claims     *
33  * about the suitability of this software for any purpose. It is          *
34  * provided "as is" without express or implied warranty.                  *
35  **************************************************************************/
36
37 //*************************************************************************
38 // Implementation of class AliAnalysisTaskITSsaTracks
39 // AliAnalysisTaskSE to extract QA and performance histos for ITS standalone tracks
40 // 
41 //
42 // Authors: L. Milano, milano@to.infn.it
43 //          F. Prino, prino@to.infn.it
44 //          
45 //*************************************************************************
46
47 ClassImp(AliAnalysisTaskITSsaTracks)
48 //______________________________________________________________________________
49 AliAnalysisTaskITSsaTracks::AliAnalysisTaskITSsaTracks() : AliAnalysisTaskSE("ITSsa resolution"), 
50   fOutput(0),
51   fHistNEvents(0),
52   fHistMCPhiResid(0),
53   fHistPhiResid(0),
54   fNtupleTracks(0),
55   fNPtBins(kMaxPtBins),
56   fMinITSpts(4),
57   fMinSPDpts(1),
58   fMinPtsforPid(3),
59   fMinTPCpts(50),
60   fMaxITSChi2Clu(100.),
61   fFillNtuple(kFALSE),
62   fReadMC(kFALSE),
63   fUseMCId(kFALSE)
64 {
65   //
66   for(Int_t ilay=0; ilay<6;ilay++) fRequirePoint[ilay]=kFALSE;
67   DefineInput(0, TChain::Class());
68   DefineOutput(1, TList::Class());
69   const Int_t nbins = 29;
70   Double_t xbins[nbins+1]={0.06,0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,
71                            0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,
72                            0.85,0.90,0.95,1.00,1.20,1.40,1.60,1.80,1.90,2.00};
73   SetPtBins(nbins,xbins);
74 }
75
76
77 //___________________________________________________________________________
78 AliAnalysisTaskITSsaTracks::~AliAnalysisTaskITSsaTracks(){
79   //
80   delete fHistNEvents;
81   for(Int_t iType=0; iType<kNtrackTypes; iType++){
82     delete fHistPt[iType];
83     delete fHistPtGood[iType];
84     delete fHistPtFake[iType];
85     delete fHistEtaPhi[iType];
86     delete fHistEtaPhiGood[iType];
87     delete fHistEtaPhiFake[iType];
88     delete fHistChi2[iType];
89     delete fHistChi2Good[iType];
90     delete fHistChi2Fake[iType];
91     delete fHistNclu[iType];
92     delete fHistNcluGood[iType];
93     delete fHistNcluFake[iType];
94     delete fHistdedxvsPt3cls[iType];
95     delete fHistdedxvsPt4cls[iType];
96   }
97   for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
98     delete fHistPtTPCITS[iSpec];
99     delete fHistPtITSsa[iSpec];
100     delete fHistPtITSpureSA[iSpec];
101     delete fHistEtaPhiTPCITS[iSpec];
102     delete fHistEtaPhiITSsa[iSpec];
103     delete fHistEtaPhiITSpureSA[iSpec];
104     delete fHistNcluTPCITS[iSpec];
105     delete fHistNcluITSsa[iSpec];
106     delete fHistNcluITSpureSA[iSpec];
107     delete fHistd0rphiITSpureSA[iSpec];
108     delete fHistd0zITSpureSA[iSpec];
109     delete fHistCluInLayTPCITS[iSpec];
110     delete fHistCluInLayITSsa[iSpec];
111     delete fHistCluInLayITSpureSA[iSpec];
112     delete fHistOuterLayITSpureSA[iSpec];
113     delete fHistPtResid[iSpec];
114     delete fHistPtRelResid[iSpec];
115     delete fHistInvPtResid[iSpec];
116     delete fHistInvPtRelResid[iSpec];
117     delete fHistMCPtResid[iSpec];
118     delete fHistMCPtRelResid[iSpec];
119     delete fHistMCInvPtResid[iSpec];
120     delete fHistMCInvPtRelResid[iSpec];
121   } 
122
123   delete fHistMCPhiResid;
124   delete fHistPhiResid;
125   delete fNtupleTracks;
126   if (fOutput) {
127     delete fOutput;
128     fOutput = 0;
129   }
130 }
131    
132 //________________________________________________________________________
133 void AliAnalysisTaskITSsaTracks::SetPtBins(Int_t n, Double_t* lim){
134   // define pt bins for analysis
135   if(n>kMaxPtBins){
136     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
137     return;
138   }else{
139     fNPtBins=n;
140     for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
141     for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
142   }
143 }
144 //___________________________________________________________________________
145 void AliAnalysisTaskITSsaTracks::UserCreateOutputObjects() {
146   // create output histos
147
148   fOutput = new TList();
149   fOutput->SetOwner();
150   fOutput->SetName("OutputHistos");
151
152   fNtupleTracks= new TNtuple("ntupleTracks","ntupleTracks","pt:px:py:pz:eta:phi:d0xy:d0z:dedx3:dedx4:dedx5:dedx6:truncmean:chi2:status:ITSrefit:TPCin:TPCrefit:ITSpureSA:nITSclu:nTPCclu:clumap:spdin:spdout:sddin:sddout:ssdin:ssdout:label:ptgen:pxgen:pygen:pzgen:etagen:phigen:ntracks");
153   // kinematics: pt, p eta,phi        ->  6 variables
154   // impact parameters: d0xy, d0z     ->  2 variables
155   // dE/dx: 4 Layers + trunc mean     ->  5 variables
156   // chi2 and track status:           ->  6 variables
157   // cluster infos:                   ->  9 variables
158   // MC info:                         ->  7 variables
159   // Multiplicity                     ->  1 variable
160   // Total:                              36
161   fOutput->Add(fNtupleTracks);
162
163   fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-0.5,2.5);
164   fHistNEvents->Sumw2();
165   fHistNEvents->SetMinimum(0);
166   fOutput->Add(fHistNEvents);
167
168   //binning for the dedx histogram
169   const Int_t hnbinsdedx=400;
170   Double_t hxmindedx = 0.01;
171   Double_t hxmaxdedx = 10;
172   Double_t hlogxmindedx = TMath::Log10(hxmindedx);
173   Double_t hlogxmaxdedx = TMath::Log10(hxmaxdedx);
174   Double_t hbinwidthdedx = (hlogxmaxdedx-hlogxmindedx)/hnbinsdedx;
175   Double_t hxbinsdedx[hnbinsdedx+1];
176   hxbinsdedx[0] = 0.01; 
177   for (Int_t i=1;i<=hnbinsdedx;i++) {
178     hxbinsdedx[i] = hxmindedx + TMath::Power(10,hlogxmindedx+i*hbinwidthdedx);
179   }
180   
181   TString tit[kNtrackTypes]={"TPCITS","ITSsa","ITSpureSA"};
182   for(Int_t iType=0; iType<kNtrackTypes; iType++){
183     fHistPt[iType]=new TH1F(Form("hPt%s",tit[iType].Data()),"",100,0.,2.);
184     fHistPt[iType]->Sumw2();
185     fOutput->Add(fHistPt[iType]);
186     fHistPtGood[iType]=new TH1F(Form("hPtGood%s",tit[iType].Data()),"",100,0.,2.);
187     fHistPtGood[iType]->Sumw2();
188     fOutput->Add(fHistPtGood[iType]);
189     fHistPtFake[iType]=new TH1F(Form("hPtFake%s",tit[iType].Data()),"",100,0.,2.);
190     fHistPtFake[iType]->Sumw2();
191     fOutput->Add(fHistPtFake[iType]);
192
193     fHistEtaPhi[iType] = new TH2F(Form("hEtaPhi%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
194     fHistEtaPhi[iType]->Sumw2();
195     fOutput->Add(fHistEtaPhi[iType]);
196     fHistEtaPhiGood[iType] = new TH2F(Form("hEtaPhiGood%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
197     fHistEtaPhiGood[iType]->Sumw2();
198     fOutput->Add(fHistEtaPhiGood[iType]);
199     fHistEtaPhiFake[iType] = new TH2F(Form("hEtaPhiFake%s",tit[iType].Data()),"",50,-1,1.,100,0.,2.*TMath::Pi());
200     fHistEtaPhiFake[iType]->Sumw2();
201     fOutput->Add(fHistEtaPhiFake[iType]);
202
203     fHistChi2[iType]=new TH1F(Form("hChi2%s",tit[iType].Data()),"",100,0.,10.);
204     fHistChi2[iType]->Sumw2();
205     fOutput->Add(fHistChi2[iType]);
206     fHistChi2Good[iType]=new TH1F(Form("hChi2Good%s",tit[iType].Data()),"",100,0.,10.);
207     fHistChi2Good[iType]->Sumw2();
208     fOutput->Add(fHistChi2Good[iType]);
209     fHistChi2Fake[iType]=new TH1F(Form("hChi2Fake%s",tit[iType].Data()),"",100,0.,10.);
210     fHistChi2Fake[iType]->Sumw2();
211     fOutput->Add(fHistChi2Fake[iType]);
212
213     fHistNclu[iType]=new TH1F(Form("hNclu%s",tit[iType].Data()),"",7,-0.5,6.5);
214     fHistNclu[iType]->Sumw2();
215     fOutput->Add(fHistNclu[iType]);
216     fHistNcluGood[iType]=new TH1F(Form("hNcluGood%s",tit[iType].Data()),"",7,-0.5,6.5);
217     fHistNcluGood[iType]->Sumw2();
218     fOutput->Add(fHistNcluGood[iType]);
219     fHistNcluFake[iType]=new TH1F(Form("hNcluFake%s",tit[iType].Data()),"",7,-0.5,6.5);
220     fHistNcluFake[iType]->Sumw2();
221     fOutput->Add(fHistNcluFake[iType]);
222
223     fHistdedxvsPt3cls[iType] = new TH2F(Form("hdedxvsPt3cls%s",tit[iType].Data()),"",hnbinsdedx,hxbinsdedx,900,0,1000);
224     fHistdedxvsPt3cls[iType]->Sumw2();
225     fOutput->Add(fHistdedxvsPt3cls[iType]);
226
227     fHistdedxvsPt4cls[iType] = new TH2F(Form("hdedxvsPt4cls%s",tit[iType].Data()),"",hnbinsdedx,hxbinsdedx,900,0,1000);
228     fHistdedxvsPt4cls[iType]->Sumw2();
229     fOutput->Add(fHistdedxvsPt4cls[iType]);
230   }
231
232   //-----------------------------------------------------------
233
234   TString spname[3]={"Pion","Kaon","Proton"};
235   TString hisname;
236   const Int_t nbins = fNPtBins;
237   Double_t xbins[nbins+1];
238   for(Int_t ibin=0; ibin<=nbins; ibin++) xbins[ibin]=fPtLimits[ibin];
239
240   for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
241
242     hisname.Form("hPtTPCITS%s",spname[iSpec].Data());
243     fHistPtTPCITS[iSpec] = new TH1F(hisname.Data(),"",100,0.,2.);
244     fHistPtTPCITS[iSpec]->Sumw2();
245     fOutput->Add(fHistPtTPCITS[iSpec]);
246
247     hisname.Form("hPtITSsa%s",spname[iSpec].Data());
248     fHistPtITSsa[iSpec] = new TH1F(hisname.Data(),"",100,0.,2.);
249     fHistPtITSsa[iSpec]->Sumw2();
250     fOutput->Add(fHistPtITSsa[iSpec]);
251
252     hisname.Form("hPtITSpureSA%s",spname[iSpec].Data());
253     fHistPtITSpureSA[iSpec] = new TH1F(hisname.Data(),"",100,0.,2.);
254     fHistPtITSpureSA[iSpec]->Sumw2();
255     fOutput->Add(fHistPtITSpureSA[iSpec]);
256
257     //---
258
259     hisname.Form("hEtaPhiTPCITS%s",spname[iSpec].Data());
260     fHistEtaPhiTPCITS[iSpec] = new TH2F(hisname.Data(),"",50,-1,1.,50,0.,2.*TMath::Pi());
261     fHistEtaPhiTPCITS[iSpec]->Sumw2();
262     fOutput->Add(fHistEtaPhiTPCITS[iSpec]);
263
264     hisname.Form("hEtaPhiITSsa%s",spname[iSpec].Data());
265     fHistEtaPhiITSsa[iSpec] = new TH2F(hisname.Data(),"",50,-1,1.,50,0.,2.*TMath::Pi());
266     fHistEtaPhiITSsa[iSpec]->Sumw2();
267     fOutput->Add(fHistEtaPhiITSsa[iSpec]);
268
269     hisname.Form("hEtaPhiITSpureSA%s",spname[iSpec].Data());
270     fHistEtaPhiITSpureSA[iSpec] = new TH2F(hisname.Data(),"",50,-1,1.,50,0.,2.*TMath::Pi());
271     fHistEtaPhiITSpureSA[iSpec]->Sumw2();
272     fOutput->Add(fHistEtaPhiITSpureSA[iSpec]);
273
274     //---
275
276     hisname.Form("hNcluTPCITS%s",spname[iSpec].Data());
277     fHistNcluTPCITS[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-0.5,6.5);
278     fHistNcluTPCITS[iSpec]->Sumw2();
279     fOutput->Add(fHistNcluTPCITS[iSpec]);
280
281     hisname.Form("hNcluITSsa%s",spname[iSpec].Data());
282     fHistNcluITSsa[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-0.5,6.5);
283     fHistNcluITSsa[iSpec]->Sumw2();
284     fOutput->Add(fHistNcluITSsa[iSpec]);
285
286     hisname.Form("hNcluITSpureSA%s",spname[iSpec].Data());
287     fHistNcluITSpureSA[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-0.5,6.5);
288     fHistNcluITSpureSA[iSpec]->Sumw2();
289     fOutput->Add(fHistNcluITSpureSA[iSpec]);
290
291     //---
292
293     hisname.Form("hd0rphiITSpureSA%s",spname[iSpec].Data());
294     fHistd0rphiITSpureSA[iSpec] = new TH2F(hisname.Data(),"",nbins,xbins,2000,-1,1);
295     fHistd0rphiITSpureSA[iSpec]->Sumw2();
296     fOutput->Add(fHistd0rphiITSpureSA[iSpec]);
297
298     hisname.Form("hd0zITSpureSA%s",spname[iSpec].Data());
299     fHistd0zITSpureSA[iSpec] = new TH2F(hisname.Data(),"",nbins,xbins,2000,-1,1);
300     fHistd0zITSpureSA[iSpec]->Sumw2();
301     fOutput->Add(fHistd0zITSpureSA[iSpec]);
302
303     //---
304
305     hisname.Form("hCluInLayTPCITS%s",spname[iSpec].Data());
306     fHistCluInLayTPCITS[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
307     fHistCluInLayTPCITS[iSpec]->Sumw2();
308     fOutput->Add(fHistCluInLayTPCITS[iSpec]);
309     
310     hisname.Form("hCluInLayITSsa%s",spname[iSpec].Data());
311     fHistCluInLayITSsa[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
312     fHistCluInLayITSsa[iSpec]->Sumw2();
313     fOutput->Add(fHistCluInLayITSsa[iSpec]);
314
315     hisname.Form("hCluInLayITSpureSA%s",spname[iSpec].Data());
316     fHistCluInLayITSpureSA[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
317     fHistCluInLayITSpureSA[iSpec]->Sumw2();
318     fOutput->Add(fHistCluInLayITSpureSA[iSpec]);
319     
320     hisname.Form("hOuterLayITSpureSA%s",spname[iSpec].Data());
321     fHistOuterLayITSpureSA[iSpec] = new TH2F(hisname.Data(),"",100,0.,2.,7,-1.5,5.5);
322     fHistOuterLayITSpureSA[iSpec]->Sumw2();
323     fOutput->Add(fHistOuterLayITSpureSA[iSpec]);
324     
325     //---
326
327     hisname.Form("hPtResid%s",spname[iSpec].Data());
328     fHistPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
329     fHistPtResid[iSpec]->Sumw2();
330     fOutput->Add(fHistPtResid[iSpec]);
331     hisname.Form("hPtRelResid%s",spname[iSpec].Data());
332     fHistPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
333     fHistPtRelResid[iSpec]->Sumw2();
334     fOutput->Add(fHistPtRelResid[iSpec]);
335     
336     hisname.Form("hInvPtResid%s",spname[iSpec].Data());
337     fHistInvPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
338     fHistInvPtResid[iSpec]->Sumw2();
339     fOutput->Add(fHistInvPtResid[iSpec]);
340     hisname.Form("hInvPtRelResid%s",spname[iSpec].Data());
341     fHistInvPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
342     fHistInvPtRelResid[iSpec]->Sumw2();
343     fOutput->Add(fHistInvPtRelResid[iSpec]);
344     
345     hisname.Form("hMCPtResid%s",spname[iSpec].Data());
346     fHistMCPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
347     fHistMCPtResid[iSpec]->Sumw2();
348     fOutput->Add(fHistMCPtResid[iSpec]);
349     hisname.Form("hMCPtRelResid%s",spname[iSpec].Data());
350     fHistMCPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
351     fHistMCPtRelResid[iSpec]->Sumw2();
352     fOutput->Add(fHistMCPtRelResid[iSpec]);
353     
354     hisname.Form("hMCInvPtResid%s",spname[iSpec].Data());
355     fHistMCInvPtResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-1.,1.);
356     fHistMCInvPtResid[iSpec]->Sumw2();
357     fOutput->Add(fHistMCInvPtResid[iSpec]);
358     hisname.Form("hMCInvPtRelResid%s",spname[iSpec].Data());
359     fHistMCInvPtRelResid[iSpec]=new TH2F(hisname.Data(),hisname.Data(),nbins,xbins,100,-0.5,0.5);
360     fHistMCInvPtRelResid[iSpec]->Sumw2();
361     fOutput->Add(fHistMCInvPtRelResid[iSpec]);
362     
363   }
364
365   fHistMCPhiResid=new TH2F("hMCPhiResid","",nbins,xbins,100,-0.5,0.5);
366   fHistMCPhiResid->Sumw2();
367   fOutput->Add(fHistMCPhiResid);
368   fHistPhiResid=new TH2F("hPhiResid","",nbins,xbins,100,-0.5,0.5);
369   fHistPhiResid->Sumw2();
370   fOutput->Add(fHistPhiResid);
371
372   PostData(1,fOutput);
373
374 }
375 //______________________________________________________________________________
376 void AliAnalysisTaskITSsaTracks::UserExec(Option_t *)
377 {
378   //
379
380   AliESDEvent *esd = (AliESDEvent*) (InputEvent());
381
382
383   if(!esd) {
384     printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
385     return;
386   } 
387
388
389   if(!ESDfriend()) {
390     printf("AliAnalysisTaskSDDRP::Exec(): bad ESDfriend\n");
391     return;
392   }
393   
394   AliStack* stack=0;
395
396   if(fReadMC){
397     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
398     if (!eventHandler) {
399       Printf("ERROR: Could not retrieve MC event handler");
400       return;
401     }
402     AliMCEvent* mcEvent = eventHandler->MCEvent();
403     if (!mcEvent) {
404       Printf("ERROR: Could not retrieve MC event");
405       return;
406     }
407     stack = mcEvent->Stack();
408     if (!stack) {
409       Printf("ERROR: stack not available");
410       return;
411     }
412   }
413
414
415   fHistNEvents->Fill(0);
416   const AliESDVertex *spdv=esd->GetPrimaryVertexSPD();
417   if(spdv->GetNContributors()<=0) return;
418   fHistNEvents->Fill(1);
419   
420   const Int_t ntSize=36;
421   Float_t xnt[ntSize];
422   
423   Int_t ntracks = esd->GetNumberOfTracks();
424   for (Int_t iTrack=0; iTrack < ntracks; iTrack++) {
425     AliESDtrack * track = esd->GetTrack(iTrack);
426     if (!track) continue;
427     Float_t pttrack=track->Pt();
428     Float_t ptrack=track->P();
429     Float_t pxtrack=track->Px();
430     Float_t pytrack=track->Py();
431     Float_t pztrack=track->Pz();
432     Float_t impactXY=-999, impactZ=-999;
433     track->GetImpactParameters(impactXY, impactZ);
434     Double_t dedxlay[4];
435     track->GetITSdEdxSamples(dedxlay);
436     Float_t dedx=track->GetITSsignal();
437     Float_t chi2=track->GetITSchi2();
438     Int_t status=track->GetStatus();
439     Bool_t isITSrefit=kFALSE;
440     Bool_t isTPCin=kFALSE;
441     Bool_t isTPCrefit=kFALSE;
442     Bool_t isPureSA=kFALSE;
443     if(status&AliESDtrack::kITSrefit) isITSrefit=kTRUE; 
444     if(status&AliESDtrack::kTPCin) isTPCin=kTRUE; 
445     if(status&AliESDtrack::kTPCrefit) isTPCrefit=kTRUE; 
446     if(status&AliESDtrack::kITSpureSA) isPureSA=kTRUE;
447     Bool_t isSA=kTRUE;
448     if(status&AliESDtrack::kTPCin) isSA=kFALSE;
449     Int_t nTPCclus=track->GetNcls(1);
450     Int_t nITSclus=track->GetNcls(0);
451     UChar_t clumap=track->GetITSClusterMap();
452     Int_t idet;
453     Float_t xloc,zloc;
454     Int_t statusLay[6];
455     for(Int_t iLay=0; iLay<6;iLay++) track->GetITSModuleIndexInfo(iLay,idet,statusLay[iLay],xloc,zloc);
456     Int_t trlabel=track->GetLabel();
457     Float_t ptgen=-999.;
458     Float_t pgen=-999.;
459     Float_t pxgen=-999.;
460     Float_t pygen=-999.;
461     Float_t pzgen=-999.;
462     Float_t etagen=-999.;
463     Float_t phigen=-999.;
464     if(fReadMC){
465       TParticle* part = stack->Particle(TMath::Abs(trlabel));
466       ptgen=part->Pt();
467       pgen=part->P();
468       pxgen=part->Px();
469       pygen=part->Py();
470       pzgen=part->Pz();
471       etagen=part->Eta();
472       phigen=part->Phi();
473     }
474
475     Int_t indexn=0;
476     xnt[indexn++]=pttrack;
477     xnt[indexn++]=pxtrack;
478     xnt[indexn++]=pytrack;
479     xnt[indexn++]=pztrack;
480     xnt[indexn++]=track->Eta();
481     xnt[indexn++]=track->Phi();
482     xnt[indexn++]=impactXY;
483     xnt[indexn++]=impactZ;
484     for(Int_t iLay=0; iLay<4; iLay++) xnt[indexn++]=dedxlay[iLay];
485     xnt[indexn++]=dedx;
486     xnt[indexn++]=chi2;
487     xnt[indexn++]=status;
488     xnt[indexn++]=isITSrefit;
489     xnt[indexn++]=isTPCin;
490     xnt[indexn++]=isTPCrefit;
491     xnt[indexn++]=isPureSA;
492     xnt[indexn++]=nITSclus;
493     xnt[indexn++]=nTPCclus;
494     xnt[indexn++]=clumap;
495     for(Int_t iLay=0; iLay<6; iLay++) xnt[indexn++]=(Float_t)statusLay[iLay];
496     xnt[indexn++]=trlabel;
497     xnt[indexn++]=ptgen;
498     xnt[indexn++]=pxgen;
499     xnt[indexn++]=pygen;
500     xnt[indexn++]=pzgen;
501     xnt[indexn++]=etagen;
502     xnt[indexn++]=phigen;
503     xnt[indexn++]=ntracks;
504
505     if(indexn>ntSize) printf("AliAnalysisTaskITSsaTracks: ERROR ntuple insexout of range\n");
506     if(fFillNtuple) fNtupleTracks->Fill(xnt);
507
508     if(!(status&AliESDtrack::kITSrefit)) continue;
509     if(nITSclus<fMinITSpts)continue;
510     if((chi2/nITSclus) > fMaxITSChi2Clu) continue;
511
512     Int_t nSPDPoints=0;
513     for(Int_t i=0; i<2; i++){
514       if(clumap&(1<<i)) ++nSPDPoints;
515     }
516     if(nSPDPoints<fMinSPDpts) continue;
517
518     Int_t nPointsForPid=0;
519     for(Int_t i=2; i<6; i++){
520       if(clumap&(1<<i)) ++nPointsForPid;
521     }
522
523     Int_t iTrackType=-1;
524     if(status&AliESDtrack::kTPCin){
525       iTrackType=kTypeTPCITS;
526     }else{
527       if(status&AliESDtrack::kITSpureSA) iTrackType=kTypeITSpureSA;
528       else iTrackType=kTypeITSsa;
529     }
530     if(iTrackType<0 || iTrackType>=kNtrackTypes) continue;
531     fHistPt[iTrackType]->Fill(pttrack);
532     fHistEtaPhi[iTrackType]->Fill(track->Eta(),track->Phi());
533     fHistChi2[iTrackType]->Fill(chi2/nITSclus);
534     fHistNclu[iTrackType]->Fill(nITSclus);
535     if(nPointsForPid==3){
536       fHistdedxvsPt3cls[iTrackType]->Fill(pttrack,dedx);
537     }
538     if(nPointsForPid==4){
539       fHistdedxvsPt4cls[iTrackType]->Fill(pttrack,dedx);
540     }
541     if(fReadMC){
542       if(trlabel<0){
543         fHistPtFake[iTrackType]->Fill(pttrack);
544         fHistEtaPhiFake[iTrackType]->Fill(track->Eta(),track->Phi());
545         fHistChi2Fake[iTrackType]->Fill(chi2/nITSclus);
546         fHistNcluFake[iTrackType]->Fill(nITSclus);
547       }else{
548         fHistPtGood[iTrackType]->Fill(pttrack);
549         fHistEtaPhiGood[iTrackType]->Fill(track->Eta(),track->Phi());
550         fHistChi2Good[iTrackType]->Fill(chi2/nITSclus);
551         fHistNcluGood[iTrackType]->Fill(nITSclus);
552       }
553     }
554
555     Int_t hadronSpecie=-1;
556     if(fReadMC && fUseMCId){
557       if(trlabel>=0){
558         TParticle* part = stack->Particle(trlabel);
559         Int_t pdg=TMath::Abs(part->GetPdgCode());
560         if(pdg==211) hadronSpecie=kPion;
561         else if(pdg==321) hadronSpecie=kKaon;
562         else if(pdg==2212) hadronSpecie=kProton;
563       }
564     }else{
565       if(nPointsForPid<fMinPtsforPid)continue;
566       AliITSPIDResponse pidits(fReadMC);
567       Float_t nSigmaPion=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kPion,nPointsForPid,isSA);
568       if(nSigmaPion>-2. && nSigmaPion<2.){
569         hadronSpecie=kPion;
570       }else{
571         Float_t nSigmaKaon=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kKaon,nPointsForPid,isSA);
572         if(nSigmaKaon>-2. && nSigmaKaon<2.){
573           hadronSpecie=kKaon;
574         }else{
575           Float_t nSigmaProton=pidits.GetNumberOfSigmas(ptrack,dedx,AliPID::kProton,nPointsForPid,isSA);
576           if(nSigmaProton>-2. && nSigmaProton<2.){
577             hadronSpecie=kProton;
578           }
579         }
580       } 
581     }
582     if(hadronSpecie<0) continue;
583     if(iTrackType==kTypeTPCITS){ // TPC+ITS tracks
584       fHistPtTPCITS[hadronSpecie]->Fill(pttrack);
585       fHistEtaPhiTPCITS[hadronSpecie]->Fill(track->Eta(),track->Phi());
586       fHistNcluTPCITS[hadronSpecie]->Fill(pttrack,nITSclus);
587       fHistCluInLayTPCITS[hadronSpecie]->Fill(-1.);
588       for(Int_t iBit=0; iBit<6; iBit++){
589         if(clumap&(1<<iBit)) fHistCluInLayTPCITS[hadronSpecie]->Fill(pttrack,iBit);
590       }
591     }else if(iTrackType==kTypeITSpureSA){ // TPC+ITS tracks
592       fHistPtITSpureSA[hadronSpecie]->Fill(pttrack);
593       fHistEtaPhiITSpureSA[hadronSpecie]->Fill(track->Eta(),track->Phi());
594       fHistNcluITSpureSA[hadronSpecie]->Fill(pttrack,nITSclus);
595       fHistd0rphiITSpureSA[hadronSpecie]->Fill(pttrack,impactXY);
596       fHistd0zITSpureSA[hadronSpecie]->Fill(pttrack,impactZ);
597       fHistCluInLayITSpureSA[hadronSpecie]->Fill(-1.);
598       Int_t outerLay=-1;
599       for(Int_t iBit=0; iBit<6; iBit++){
600         if(clumap&(1<<iBit)){
601           fHistCluInLayITSpureSA[hadronSpecie]->Fill(pttrack,iBit);
602           if(iBit>outerLay) outerLay=iBit;
603         }
604         }
605       fHistOuterLayITSpureSA[hadronSpecie]->Fill(pttrack,outerLay);     
606       
607         
608       if(fReadMC){  
609         if(trlabel>=0){
610           TParticle* part = stack->Particle(trlabel);
611           ptgen=part->Pt();
612           Float_t invpttrack=track->OneOverPt();
613           Float_t invptgen=0.;
614           if(ptgen>0.) invptgen=1./ptgen;
615           fHistMCPtResid[hadronSpecie]->Fill(pttrack,pttrack-ptgen);
616           fHistMCPtRelResid[hadronSpecie]->Fill(pttrack,(pttrack-ptgen)/ptgen);
617           fHistMCInvPtResid[hadronSpecie]->Fill(pttrack,invpttrack-invptgen);
618           fHistMCInvPtRelResid[hadronSpecie]->Fill(pttrack,(invpttrack-invptgen)/invptgen);       
619           Float_t deltaphi=track->Phi()-part->Phi();
620           if(deltaphi<-TMath::Pi()) deltaphi+=2*TMath::Pi();
621           if(deltaphi>TMath::Pi()) deltaphi-=2*TMath::Pi();
622           fHistMCPhiResid->Fill(pttrack,deltaphi);
623         }
624       }
625     }else if(iTrackType==kTypeITSsa){ // TPC+ITS tracks
626       fHistPtITSsa[hadronSpecie]->Fill(pttrack);
627       fHistEtaPhiITSsa[hadronSpecie]->Fill(track->Eta(),track->Phi());
628       fHistNcluITSsa[hadronSpecie]->Fill(pttrack,nITSclus);
629       fHistCluInLayITSsa[hadronSpecie]->Fill(-1.);
630       for(Int_t iBit=0; iBit<6; iBit++){
631         if(clumap&(1<<iBit)) fHistCluInLayITSsa[hadronSpecie]->Fill(pttrack,iBit);
632       }
633     }
634   
635     // match TPCITS with ITSpureSA
636     if(nITSclus<fMinITSpts || nTPCclus<fMinTPCpts) continue;      
637     Int_t idxMI[12],idxSA[12];
638     for(Int_t icl=0; icl<12; icl++){ 
639       idxMI[icl]=-1; 
640       idxSA[icl]=-1;
641     }
642     Int_t ncls=track->GetClusters(0,idxMI);
643     if(fMinITSpts<6){
644       Bool_t accept=kTRUE;
645       for(Int_t ilay=0; ilay<6; ilay++){
646         if(fRequirePoint[ilay]){
647           Int_t mask = 1<<ilay;
648           if(!(clumap & mask)){ 
649             accept=kFALSE;
650             break;
651           }
652         }
653       }
654       if(!accept) continue;
655     }
656     // Sort
657     for(Int_t i=0;i<12;i++){
658       for(Int_t j=i+1;j<12;j++){
659         if(idxMI[j]>idxMI[i]){
660           Int_t tmp=idxMI[j];
661           idxMI[j]=idxMI[i];
662           idxMI[i]=tmp;
663         }
664       }
665     }
666     //    for(Int_t i=0; i<12; i++) printf("%d ",idxMI[i]);
667     //    printf("\n");
668     if(idxMI[0]<0 && idxMI[0]==idxMI[1]) continue;
669     Bool_t matched=kFALSE;
670     for (Int_t iTrack2 = 0; iTrack2 < ntracks; iTrack2++) {
671       if(matched) break;
672       if(iTrack2==iTrack) continue;
673       AliESDtrack* track2 = esd->GetTrack(iTrack2);
674       Int_t status2=track2->GetStatus();
675       if(!(status2&AliESDtrack::kITSrefit)) continue;
676       if(!(status2&AliESDtrack::kITSpureSA)) continue;      
677       Int_t clumap2=track2->GetITSClusterMap();
678       Int_t nITSclus2=track2->GetNcls(0);
679       Int_t nTPCclus2=track2->GetNcls(1);
680       if(nITSclus2<fMinITSpts || nTPCclus2>0) continue; 
681       Int_t ncls2=track2->GetClusters(0,idxSA);
682       if(ncls2!=ncls) continue;
683       if(fMinITSpts<6){
684         Bool_t accept=kTRUE;
685         for(Int_t ilay=0; ilay<6; ilay++){
686           if(fRequirePoint[ilay]){
687             Int_t mask = 1<<ilay;
688             if(!(clumap2 & mask)){ 
689               accept=kFALSE;
690               break;
691             }
692           }
693         }
694         if(!accept) continue;
695       }
696       // Sort
697       for(Int_t i=0;i<12;i++){
698         for(Int_t j=i+1;j<12;j++){
699           if(idxSA[j]>idxSA[i]){
700             Int_t tmp=idxSA[j];
701             idxSA[j]=idxSA[i];
702             idxSA[i]=tmp;
703           }
704         }
705       }
706       Int_t match=0;
707       for(Int_t icl=0; icl<ncls; icl++){
708         if(idxSA[icl]!=idxMI[icl]){
709           match=0; 
710           break;
711         }
712         else match++;
713       }
714       if(match==ncls && match>0){
715         matched=kTRUE;
716         Float_t pt1=track->Pt();
717         Float_t pt2=track2->Pt();
718         Float_t ptm1=track->OneOverPt();
719         Float_t ptm2=track2->OneOverPt();
720         fHistPtResid[hadronSpecie]->Fill(pt1,pt2-pt1);
721         fHistPtRelResid[hadronSpecie]->Fill(pt1,(pt2-pt1)/pt1);
722         fHistInvPtResid[hadronSpecie]->Fill(pt1,ptm2-ptm1);
723         fHistInvPtRelResid[hadronSpecie]->Fill(pt1,(ptm2-ptm1)/ptm1);
724         Float_t deltaphi=track->Phi()-track2->Phi();
725         if(deltaphi<-TMath::Pi()) deltaphi+=2*TMath::Pi();
726         if(deltaphi>TMath::Pi()) deltaphi-=2*TMath::Pi();
727         fHistPhiResid->Fill(pt1,deltaphi);
728       }
729     }
730   }
731
732   PostData(1,fOutput);
733   
734 }
735 //______________________________________________________________________________
736 void AliAnalysisTaskITSsaTracks::Terminate(Option_t */*option*/)
737 {
738   // Terminate analysis
739   fOutput = dynamic_cast<TList*> (GetOutputData(1));
740   if (!fOutput) {     
741     printf("ERROR: fOutput not available\n");
742     return;
743   }
744   fNtupleTracks = dynamic_cast<TNtuple*>(fOutput->FindObject("ntupleTracks"));
745   fHistNEvents= dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
746
747   TString spname[3]={"Pion","Kaon","Proton"};
748
749   for(Int_t iSpec=0; iSpec<kNspecies; iSpec++){
750     fHistPtTPCITS[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtTPCITS%s",spname[iSpec].Data())));
751     fHistPtITSsa[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtITSsa%s",spname[iSpec].Data())));
752     fHistPtITSpureSA[iSpec]= dynamic_cast<TH1F*>(fOutput->FindObject(Form("hPtITSpureSA%s",spname[iSpec].Data())));
753
754     fHistEtaPhiTPCITS[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hEtaPhiTPCITS%s",spname[iSpec].Data())));
755     fHistEtaPhiITSsa[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hEtaPhiITSsa%s",spname[iSpec].Data())));
756     fHistEtaPhiITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hEtaPhiITSpureSA%s",spname[iSpec].Data())));
757     
758     fHistNcluTPCITS[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hNcluTPCITS%s",spname[iSpec].Data())));
759     fHistNcluITSsa[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hNcluITSsa%s",spname[iSpec].Data())));
760     fHistNcluITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hNcluITSpureSA%s",spname[iSpec].Data())));
761     
762     fHistd0rphiITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hd0rphiITSpureSA%s",spname[iSpec].Data())));
763     fHistd0zITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hd0zITSpureSA%s",spname[iSpec].Data())));
764     
765     fHistCluInLayTPCITS[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hCluInLayTPCITS%s",spname[iSpec].Data())));
766     fHistCluInLayITSsa[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hCluInLayITSsa%s",spname[iSpec].Data())));
767     fHistCluInLayITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hCluInLayITSpureSA%s",spname[iSpec].Data())));
768     
769     fHistOuterLayITSpureSA[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hOuterLayITSpureSA%s",spname[iSpec].Data())));
770   
771     fHistPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hPtResid%s",spname[iSpec].Data())));
772     fHistPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hPtRelResid%s",spname[iSpec].Data())));
773     fHistInvPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hInvPtResid%s",spname[iSpec].Data())));
774     fHistInvPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hInvPtRelResid%s",spname[iSpec].Data())));
775     fHistMCPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCPtResid%s",spname[iSpec].Data())));
776     fHistMCPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCPtRelResid%s",spname[iSpec].Data())));
777     fHistMCInvPtResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCInvPtResid%s",spname[iSpec].Data())));
778     fHistMCInvPtRelResid[iSpec]= dynamic_cast<TH2F*>(fOutput->FindObject(Form("hMCInvPtRelResid%s",spname[iSpec].Data())));
779   
780   }
781   fHistMCPhiResid= dynamic_cast<TH2F*>(fOutput->FindObject("hMCPhiResid"));
782   fHistPhiResid= dynamic_cast<TH2F*>(fOutput->FindObject("hPhiResid"));
783   return;
784 }
785
786
787
788
789