07422ae9b1adb32f4662f4df1ca5a9bdee726aa2
[u/mrichter/AliRoot.git] / prod / acrcaf / qa_pp / AliAnalysisTaskQASym.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 /* $Id$ */
17
18 //------------------------------
19 // Analysis task for quality-assurance of central tarcking
20 // mainly based on fundamental symmetries 
21 //
22 // contact eva.sicking@cern.ch
23 // authors 
24 // Authors: Jan Fiete Grosse-Oetringhaus, Christian Klein-Boesing,
25 //          Andreas Morsch, Eva Sicking
26
27
28 #include "TChain.h"
29 #include "TTree.h"
30 #include "TH1F.h"
31 #include "TH2F.h"
32 #include "TCanvas.h"
33 #include "TList.h"
34 #include "TParticle.h"
35 #include "TParticlePDG.h"
36 #include "TProfile.h"
37 #include "TNtuple.h"
38 #include "TFile.h"
39
40 #include "AliAnalysisTask.h"
41 #include "AliAnalysisManager.h"
42
43 #include "AliESDEvent.h"
44 #include "AliLog.h"
45 #include "AliESDVertex.h"
46 #include "AliESDInputHandler.h"
47 #include "AliESDtrackCuts.h"
48 #include "AliMultiplicity.h"
49
50
51 #include "AliAnalysisTaskQASym.h"
52 #include "AliExternalTrackParam.h"
53 #include "AliTrackReference.h"
54
55 #include "AliHeader.h"
56 #include "AliGenEventHeader.h"
57 #include "AliGenDPMjetEventHeader.h"
58
59 ClassImp(AliAnalysisTaskQASym)
60
61 //________________________________________________________________________
62 AliAnalysisTaskQASym::AliAnalysisTaskQASym(const char *name) 
63   : AliAnalysisTaskSE(name) 
64     ,fFieldOn(kTRUE)
65     ,fHists(0)
66     ,fHistRECpt(0)
67     ,fEta(0)
68     ,fEtaPhi(0)
69     ,fEtaPt(0)
70     ,fQPt(0)
71     ,fDca(0)
72     ,fqRec(0)
73     ,fsigmaPt(0)
74   
75     ,fRecPtPos(0)
76     ,fRecPtNeg(0)
77     ,fRecPhiPos(0)
78     ,fRecPhiNeg(0)
79     ,fRecEtaPos(0)
80     ,fRecEtaNeg(0)
81     ,fRecEtaPtPos(0)
82     ,fRecEtaPtNeg(0)
83     ,fRecDcaPos(0)
84     ,fRecDcaNeg(0)
85     ,fRecDcaNegInv(0)
86     ,fRecDPos(0)
87     ,fRecDNeg(0)
88
89
90     ,fRecQPtPosEta(0)
91     ,fRecQPtNegEta(0)
92     ,fRecPtPosEta(0)
93     ,fRecPtNegEta(0)
94     ,fRecPhiPosEta(0)
95     ,fRecPhiNegEta(0)
96     ,fRecDcaPosEta(0)
97     ,fRecDcaNegEta(0)
98     ,fRecDPosEta(0)
99     ,fRecDNegEta(0)
100
101     ,fRecPtPosVz(0)
102     ,fRecPtNegVz(0)
103     ,fRecEtaPosVz(0)
104     ,fRecEtaNegVz(0)
105     ,fRecPhiPosVz(0)
106     ,fRecPhiNegVz(0)
107     ,fSignedDcaPosVz(0)
108     ,fSignedDcaNegVz(0)
109     ,fRecQPtPosEtaVz(0)
110     ,fRecQPtNegEtaVz(0)
111     ,fRecEtaPtPosVz(0)
112     ,fRecEtaPtNegVz(0)
113
114
115     ,fDeltaPhiAll(0)
116     ,fDeltaPhiLeading(0) 
117     ,fDiffDcaD(0)
118     ,fPhiRec(0)
119     ,fThetaRec(0)
120     ,fNumber(0)
121     ,fVx(0)
122     ,fVy(0)
123     ,fVz(0)
124     ,fCuts(0)
125   
126 {
127     // Constructor
128     //
129   for(Int_t i = 0;i<18;++i){
130     fRecPtTpcSector[i] = 0;
131     fRecEtaTpcSector[i] = 0;
132     fSignedDcaTpcSector[i] = 0;
133     fRecQPtTpcSector[i] = 0;
134     fRecEtaPtTpcSector[i] = 0;
135   }
136
137   for(Int_t i = 0;i< 7;++i){
138     fRecPtPosLadder[i] = 0;
139     fRecPtNegLadder[i] = 0;
140     fRecPhiPosLadder[i] = 0;
141     fRecPhiNegLadder[i] = 0;
142     fRecEtaPosLadder[i] = 0;
143     fRecEtaNegLadder[i] = 0;
144     fSignDcaPos[i] = 0;
145     fSignDcaNeg[i] = 0;
146     fSignDcaNegInv[i] = 0;
147     fPtSigmaPos[i] =0;
148     fPtSigmaNeg[i] =0;
149     fqPtRec[i] =0;
150     fDcaSigmaPos[i] =0;
151     fDcaSigmaNeg[i] =0;
152   }
153
154   DefineOutput(1,  TList::Class()); 
155
156   
157   
158 }
159
160
161 //________________________________________________________________________
162 void AliAnalysisTaskQASym::UserCreateOutputObjects()
163 {
164   // Create histograms
165   // Called once
166
167   Double_t range = 1.;
168   Double_t pt = 20.;
169
170   fHists = new TList();
171   fHistRECpt   = new TH1F("fHistRECpt", 
172                           " p_{T}",
173                           100, 0., pt);
174   fEta   = new TH1F("fEta", 
175                           " #eta",
176                           200, -2., 2.);
177   fEtaPhi   = new TH2F("fEtaPhi", 
178                           " #eta - #phi",
179                           200, -2., 2., 128, 0., 2. * TMath::Pi());
180   
181   fThetaRec   = new TH1F("fThetaRec", 
182                           " #theta",
183                           180, 0., TMath::Pi());
184   fPhiRec   = new TH1F("fPhiRec", 
185                           " #phi",
186                           180, 0., 2*TMath::Pi());
187   fNumber   = new TH1F("fNumber", 
188                        "number of tracks per event",
189                        50, 0.5, 49.5);
190   fVx   = new TH1F("fVx", 
191                    "X of vertex",
192                    100, -5., 5.);
193   fVy   = new TH1F("fVy", 
194                    "Y of vertex",
195                    100, -5., 5.);
196   fVz   = new TH1F("fVz", 
197                    "Z of vertex",
198                    500, -50., 50.);
199
200   fEtaPt   = new TH1F("fEtaPt", 
201                           " #eta/p_{T} ",
202                           100, -1., 1.);
203
204   fQPt   = new TH1F("fQPt", 
205                           " charge/p_{T} ",
206                           100, -1., 1.);
207
208   fDca   = new TH1F("fDca", 
209                           " dca ",
210                           200, -range, range);
211
212
213   fqRec    = new TH1F("fqRec",   
214                           " charge all reconstructed particle",
215                           21, -9.5, 10.5);
216   
217   fsigmaPt    = new TH1F("fsigmaPt",   
218                           "Log_{10}(#sigma_{p_{T}})",
219                           200, -2., 8.);
220
221
222
223
224   //------------
225   for(Int_t ITSlayer_case=0;ITSlayer_case<7;ITSlayer_case++){
226
227     fSignDcaPos[ITSlayer_case]   = new TH1F(Form("fSignDcaPos%d", ITSlayer_case),  
228                                    " Signed dca", 
229                                    200, -range, range);
230     fSignDcaPos[ITSlayer_case]->GetXaxis()->SetTitle("dca");
231     fSignDcaPos[ITSlayer_case]->GetYaxis()->SetTitle("");
232    
233  
234     fSignDcaNeg[ITSlayer_case]   = new TH1F(Form("fSignDcaNeg%d", ITSlayer_case),  
235                                    " Signed dcas",
236                                    200, -range, range);
237     fSignDcaNeg[ITSlayer_case]->GetXaxis()->SetTitle("dca");
238     fSignDcaNeg[ITSlayer_case]->GetYaxis()->SetTitle("");
239
240     fSignDcaNegInv[ITSlayer_case]   = new TH1F(Form("fSignDcaNegInv%d", ITSlayer_case),  
241                                    " inverse Signed dca ",
242                                    200, -range, range);
243     fSignDcaNegInv[ITSlayer_case]->GetXaxis()->SetTitle("-dca");
244     fSignDcaNegInv[ITSlayer_case]->GetYaxis()->SetTitle("");
245
246
247
248
249     fPtSigmaPos[ITSlayer_case]   = new TH1F(Form("fPtSigmaPos%d", ITSlayer_case),  
250                                    " #sigma_{pT} ",
251                                    208, -2., 8.);
252     fPtSigmaPos[ITSlayer_case]->GetXaxis()->SetTitle("Log_{10}(#sigma_{pT})");
253     fPtSigmaPos[ITSlayer_case]->GetYaxis()->SetTitle("");
254     
255     
256     fPtSigmaNeg[ITSlayer_case]   = new TH1F(Form("fPtSigmaNeg%d",ITSlayer_case),  
257                                   " #sigma_{pT}",
258                                    208, -2., 8.);
259     fPtSigmaNeg[ITSlayer_case]->GetXaxis()->SetTitle("Log_{10}(#sigma_{pT})");
260     fPtSigmaNeg[ITSlayer_case]->GetYaxis()->SetTitle("");
261
262
263
264
265
266     fqPtRec[ITSlayer_case]   = new TH1F(Form("fqPtRec%d",ITSlayer_case),  
267                                   "q/ p_{T}",
268                                    200, -100., 100.);
269     fqPtRec[ITSlayer_case]->GetXaxis()->SetTitle("q_{tr}/p_{T, tr} (GeV/c)");
270     fqPtRec[ITSlayer_case]->GetYaxis()->SetTitle("");
271
272   
273    
274
275
276     fDcaSigmaPos[ITSlayer_case]   = new TH2F(Form("fDcaSigmaPos%d", ITSlayer_case),  
277                                    " p_{T} shift vs #sigma_{pT} ",
278                                    200, -range, range,200, -4., 4. );
279     fDcaSigmaPos[ITSlayer_case]->GetXaxis()->SetTitle("signed DCA)");
280     fDcaSigmaPos[ITSlayer_case]->GetYaxis()->SetTitle("log_{10}(#sigma_{pT})");
281     
282     
283     fDcaSigmaNeg[ITSlayer_case]   = new TH2F(Form("fDcaSigmaNeg%d", ITSlayer_case),  
284                                    " p_{T} shift vs #sigma_{pT} ",
285                                    200, -range, range,200, -4., 4. );
286     fDcaSigmaNeg[ITSlayer_case]->GetXaxis()->SetTitle("signed DCA");
287     fDcaSigmaNeg[ITSlayer_case]->GetYaxis()->SetTitle("log_{10}(#sigma_{pT})");
288
289
290     
291     
292  
293     
294     // YIELDs---------- positive and negative particles
295     
296     fRecPtPos   = new TH1F("fRecPtPos", 
297                            " p_{T}",
298                            100, 0.,pt);
299     fRecPtPos->GetXaxis()->SetTitle("p_{T} (GeV/c)");
300     fRecPtNeg   = new TH1F("fRecPtNeg", 
301                            " p_{T} ",
302                            100, 0., pt);
303     fRecPtNeg->GetXaxis()->SetTitle("p_{T} (GeV/c)");
304
305     
306     fRecPhiPos   = new TH1F("fRecPhiPos", 
307                             "#phi",
308                             361, 0., 360.);
309     fRecPhiPos->GetXaxis()->SetTitle("#phi (deg)");
310   
311     fRecPhiNeg   = new TH1F("fRecPhiNeg", 
312                             "#phi ",
313                             361, 0., 360.);
314     fRecPhiNeg->GetXaxis()->SetTitle("#phi (deg)");
315     
316     fRecEtaPos   = new TH1F("fRecEtaPos", 
317                             "#eta",
318                             200, -2., 2.);
319     fRecEtaPos->GetXaxis()->SetTitle("#eta");
320
321     fRecEtaNeg   = new TH1F("fRecEtaNeg", 
322                             "#eta",
323                             200, -2., 2.);
324     fRecEtaNeg->GetXaxis()->SetTitle("#eta");
325     
326     fRecEtaPtPos   = new TH1F("fRecEtaPtPos", 
327                               "#eta/p_{T}",
328                               200, -0.1, .1);
329     fRecEtaPtPos->GetXaxis()->SetTitle("#eta/p_{T}");
330
331     fRecEtaPtNeg   = new TH1F("fRecEtaPtNeg", 
332                               "#eta/p_{T}",
333                               200, -.1, .1);
334     fRecEtaPtNeg->GetXaxis()->SetTitle("#eta/p_{T}");
335
336     fRecDcaPos   = new TH1F("fRecDcaPos", 
337                          " dca",
338                            100, -range, range);
339     fRecDcaPos->GetXaxis()->SetTitle("dca (cm)");
340     fRecDcaNeg   = new TH1F("fRecDcaNeg", 
341                            " dca",
342                            100, -range, range);
343     fRecDcaNeg->GetXaxis()->SetTitle("dca (cm)");
344
345     fRecDcaNegInv   = new TH1F("fRecDcaNegInv", 
346                            " dca",
347                            100, -range, range);
348     fRecDcaNegInv->GetXaxis()->SetTitle("dca (cm)");
349
350
351     fRecDPos   = new TH1F("fRecDPos", 
352                          " d",
353                            100, -range, range);
354     fRecDPos->GetXaxis()->SetTitle("d (cm)");
355     fRecDNeg   = new TH1F("fRecDNeg", 
356                            "d",
357                            100, -range, range);
358     fRecDNeg->GetXaxis()->SetTitle("d (cm)");
359
360
361     //  YIELDs ---------------- positive and negative eta
362     
363     
364     fRecQPtPosEta   = new TH1F("fRecQPtPosEta", 
365                                "q/p_{T}",
366                                200, -0.5, 0.5);
367     fRecQPtPosEta->GetXaxis()->SetTitle("q/p_{T} ");
368
369     fRecQPtNegEta   = new TH1F("fRecQPtNegEta", 
370                                "q/p_{T}",
371                                200, -0.5, 0.5);
372     fRecQPtNegEta->GetXaxis()->SetTitle("q/p_{T}");
373     
374     fRecPtPosEta   = new TH1F("fRecPtPosEta", 
375                               " p_{T} ",
376                               100, 0., pt);
377     fRecPtPosEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
378
379     fRecPtNegEta   = new TH1F("fRecPtNegEta", 
380                               " p_{T}",
381                               100, 0., pt);
382     fRecPtNegEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
383     
384     fRecPhiPosEta   = new TH1F("fRecPhiPosEta", 
385                             "#phi",
386                             361, 0., 360);
387     fRecPhiPosEta->GetXaxis()->SetTitle("#phi (deg)");
388
389     fRecPhiNegEta   = new TH1F("fRecPhiNegEta", 
390                             "#phi ",
391                             361, 0, 360);
392     fRecPhiNegEta->GetXaxis()->SetTitle("#phi (deg)");
393
394     fRecDcaPosEta   = new TH1F("fRecDcaPosEta", 
395                          " dca ",
396                            100, -range, range);
397     fRecDcaPosEta->GetXaxis()->SetTitle("dca (cm)");
398     fRecDcaNegEta   = new TH1F("fRecDcaNegEta", 
399                            " dca",
400                            100, -range, range);
401     fRecDcaNegEta->GetXaxis()->SetTitle("dca (cm)");
402
403     fRecDPosEta   = new TH1F("fRecDPosEta", 
404                          " d",
405                            100, -range, range);
406     fRecDPosEta->GetXaxis()->SetTitle("d (cm)");
407     fRecDNegEta   = new TH1F("fRecDNegEta", 
408                            "d",
409                            100, -5., 5.);
410     fRecDNegEta->GetXaxis()->SetTitle("d (cm)");
411
412
413     
414     //  YIELDs ---------------- for TPC sectors
415     for(Int_t sector=0; sector<18;sector++){
416       
417
418       fRecPtTpcSector[sector]   = new TH1F(Form("fRecPtTpcSector%02d",sector), 
419                                            Form("p_{T} distribution: TPC sector %d",
420                                                 sector),100, 0., pt);
421       fRecPtTpcSector[sector]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
422
423       fRecEtaTpcSector[sector]   = new TH1F(Form("fRecEtaTpcSector%02d",sector), 
424                                            Form("#eta distribution: TPC sector %d",
425                                                 sector),200, -2., 2.);
426       fRecEtaTpcSector[sector]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
427      
428
429       fSignedDcaTpcSector[sector]   = new TH1F(Form("fSignedDcaTpcSector%02d",sector), 
430                                            Form("dca distribution: TPC sector %d",
431                                                 sector),200, -range, range );
432       fSignedDcaTpcSector[sector]->GetXaxis()->SetTitle("dca");
433
434       fRecQPtTpcSector[sector]   = new TH1F(Form("fRecQPtTpcSector%02d",sector), 
435                                            Form("Q/ p_{T} distribution: TPC sector %d",
436                                                 sector),100, -1., 1.);
437       fRecQPtTpcSector[sector]->GetXaxis()->SetTitle("Q/p_{T} (GeV/c)");
438
439       fRecEtaPtTpcSector[sector]   = new TH1F(Form("fRecEtaPtTpcSector%02d",sector), 
440                                            Form("#eta/ p_{T} distribution: TPC sector %d",
441                                                 sector),100, -1., 1.);
442       fRecEtaPtTpcSector[sector]->GetXaxis()->SetTitle("#eta/p_{T} (GeV/c)");
443  
444     }
445     // YIELDS ITS ladder
446     for(Int_t i=0;i<7;i++){
447       fRecPtPosLadder[i]   = new TH1F(Form("fRecPtPosLadder%d", i), 
448                              " p_{T} distribution",
449                              100, 0., pt);
450       fRecPtPosLadder[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
451       fRecPtNegLadder[i]   = new TH1F(Form("fRecPtNegLadder%d",i), 
452                              " p_{T} distribution ",
453                              100, 0., pt);
454       fRecPtNegLadder[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
455
456
457       fRecPhiPosLadder[i]   = new TH1F(Form("fRecPhiPosLadder%d",i), 
458                                  "#phi distribution: all pos eta",
459                                  361, 0., 360);
460       fRecPhiPosLadder[i]->GetXaxis()->SetTitle("#phi (deg)");
461       
462       fRecPhiNegLadder[i]   = new TH1F(Form("fRecPhiNegLadder%d", i),
463                                  "#phi distribution: all neg eta",
464                                  361, 0, 360);
465       fRecPhiNegLadder[i]->GetXaxis()->SetTitle("#phi (deg)");
466
467
468
469       fRecEtaPosLadder[i]   = new TH1F(Form("fRecEtaPosLadder%d",i), 
470                                        "#eta distribution",
471                                  200, -2., 2.);
472       fRecEtaPosLadder[i]->GetXaxis()->SetTitle("#eta)");
473       
474       fRecEtaNegLadder[i]   = new TH1F(Form("fRecEtaNegLadder%d", i),
475                                  "#eta distribution",
476                                  200, -2., 2.);
477       fRecEtaNegLadder[i]->GetXaxis()->SetTitle("#eta");
478     }
479
480     Double_t vzmax = 30;
481
482     fRecPtPosVz = new TH2F("fRecPtPosVz", 
483                            "p_{T} distribution vs Vz()",
484                            100, -1., 2., 200,-vzmax,vzmax);
485     fRecPtPosVz->GetXaxis()->SetTitle("log_{10}(p_{T})");
486     
487     fRecPtNegVz = new TH2F("fRecPtNegVz",
488                            "p_{T} distribution vs Vz()",
489                            100, -1., 2.,200,-vzmax,vzmax);
490     fRecPtNegVz->GetXaxis()->SetTitle("Log_{10}(p_{T})");
491     
492    
493     fRecEtaPosVz= new TH2F("fRecEtaPosVz", 
494                           "#eta distribution vs Vz()",
495                           100, -2., 2., 200,-vzmax,vzmax);
496     fRecEtaPosVz->GetXaxis()->SetTitle("#eta");
497     fRecEtaNegVz = new TH2F("fRecEtaNegVz",
498                            "#eta distribution vs Vz()",
499                            100, -2., 2.,200,-vzmax,vzmax);
500     fRecEtaNegVz->GetXaxis()->SetTitle("#eta");
501
502     fRecPhiPosVz= new TH2F("fRecPhiPosVz", 
503                           "#eta distribution vs Vz()",
504                           361, 0., 360., 200,-vzmax,vzmax);
505     fRecPhiPosVz->GetXaxis()->SetTitle("#phi (deg)");
506     fRecPhiNegVz = new TH2F("fRecPhiNegVz",
507                            "dca vs Vz()",
508                            361, 0., 360.,200,-vzmax,vzmax);
509     fRecPhiNegVz->GetXaxis()->SetTitle("#phi (deg)");
510
511     fSignedDcaPosVz= new TH2F("fSignedDcaPosVz", 
512                           "#eta distribution vs Vz()",
513                           200, -range, range, 200,-vzmax,vzmax);
514     fSignedDcaPosVz->GetXaxis()->SetTitle("dca (cm)");
515     fSignedDcaNegVz = new TH2F("fSignedDcaNegVz",
516                            "dca vs Vz()",
517                            200, -range, range,200,-vzmax,vzmax);
518     fSignedDcaNegVz->GetXaxis()->SetTitle("dca (cm)");
519
520     fRecQPtPosEtaVz= new TH2F("fRecQPtPosEtaVz",
521                           " Q/p_{T} distribution vs Vz()",
522                           100, -1., 1., 200,-vzmax,vzmax);
523     fRecQPtPosEtaVz->GetXaxis()->SetTitle("Q/p_{T}");
524     fRecQPtNegEtaVz = new TH2F("fRecQPtNegEtaVz",
525                            " Q/p_{T} distribution vs Vz()",
526                            100, -1., 1.,200,-vzmax,vzmax);
527     fRecQPtNegEtaVz->GetXaxis()->SetTitle("Q/p_{T}");
528
529  
530     fRecEtaPtPosVz= new TH2F("fRecEtaPtPosVz",
531                           " #eta/p_{T} distribution vs Vz()",
532                           100, -1., 1., 200,-vzmax,vzmax);
533     fRecEtaPtPosVz->GetXaxis()->SetTitle("#eta/p_{T");
534     fRecEtaPtNegVz = new TH2F("fRecEtaPtNegVz",
535                            " #eta/p_{T} distribution vs Vz()",
536                            100, -1., 1.,200,-vzmax,vzmax);
537     fRecEtaPtNegVz->GetXaxis()->SetTitle("#eta/p_{T}");
538
539     //new
540     fDeltaPhiAll = new TH1F("fDeltaPhiAll",
541                            " #Delta #phi",200,-360,360);
542     fDeltaPhiAll->GetXaxis()->SetTitle("#Delta #phi");
543
544
545     fDeltaPhiLeading = new TH2F("fDeltaPhiLeading",
546                            " #Delta #phi",361,-360,360, 361,0, 360);
547     fDeltaPhiLeading->GetXaxis()->SetTitle("#Delta #phi (deg.)");
548     fDeltaPhiLeading->GetYaxis()->SetTitle("#phi_{leading particle} (deg.)");
549
550     fDiffDcaD    = new TH1F("fDiffDcaD",   
551                             "dca-d",
552                             200, -5., 5.);
553     
554   }
555
556   fHists->SetOwner();
557
558   fHists->Add(fHistRECpt);
559   fHists->Add(fEta);
560   fHists->Add(fEtaPhi);
561   fHists->Add(fThetaRec);
562   fHists->Add(fPhiRec);
563   fHists->Add(fNumber);
564   fHists->Add(fVx);
565   fHists->Add(fVy);
566   fHists->Add(fVz);
567
568   fHists->Add(fEtaPt);
569   fHists->Add(fQPt);
570   fHists->Add(fDca);
571
572   fHists->Add(fDeltaPhiAll);
573   fHists->Add(fDeltaPhiLeading);
574   fHists->Add(fDiffDcaD);
575
576   fHists->Add(fqRec);
577   fHists->Add(fsigmaPt);
578
579   fHists->Add(fRecPtPos);
580   fHists->Add(fRecPtNeg);
581   fHists->Add(fRecPhiPos);
582   fHists->Add(fRecPhiNeg);
583   fHists->Add(fRecEtaPos);
584   fHists->Add(fRecEtaNeg);
585   fHists->Add(fRecEtaPtPos);
586   fHists->Add(fRecEtaPtNeg);
587   fHists->Add(fRecDcaPos);
588   fHists->Add(fRecDcaNeg);
589   fHists->Add(fRecDcaNegInv);
590   fHists->Add(fRecDPos);
591   fHists->Add(fRecDNeg);
592
593
594   fHists->Add(fRecQPtPosEta);
595   fHists->Add(fRecQPtNegEta);
596   fHists->Add(fRecPtPosEta);
597   fHists->Add(fRecPtNegEta);
598   fHists->Add(fRecPhiPosEta);
599   fHists->Add(fRecPhiNegEta);
600   fHists->Add(fRecDcaPosEta);
601   fHists->Add(fRecDcaNegEta);
602   fHists->Add(fRecDPosEta);
603   fHists->Add(fRecDNegEta);
604
605
606   for(Int_t i=0;i<18;i++){
607     fHists->Add(fRecPtTpcSector[i]);
608     fHists->Add(fRecEtaTpcSector[i]);
609     fHists->Add(fSignedDcaTpcSector[i]);
610     fHists->Add(fRecQPtTpcSector[i]);
611     fHists->Add(fRecEtaPtTpcSector[i]);
612   }
613
614   for(Int_t i=0;i<7;i++){
615     fHists->Add(fRecPtPosLadder[i]);
616     fHists->Add(fRecPtNegLadder[i]);
617     fHists->Add(fRecPhiPosLadder[i]);
618     fHists->Add(fRecPhiNegLadder[i]);
619     fHists->Add(fRecEtaPosLadder[i]);
620     fHists->Add(fRecEtaNegLadder[i]);
621   }
622
623   fHists->Add(fRecPtPosVz);
624   fHists->Add(fRecPtNegVz);
625   fHists->Add(fRecEtaPosVz);
626   fHists->Add(fRecEtaNegVz);
627   fHists->Add(fRecPhiPosVz);
628   fHists->Add(fRecPhiNegVz);
629   fHists->Add(fSignedDcaPosVz);
630   fHists->Add(fSignedDcaNegVz);
631   fHists->Add(fRecQPtPosEtaVz);
632   fHists->Add(fRecQPtNegEtaVz);
633   fHists->Add(fRecEtaPtPosVz);
634   fHists->Add(fRecEtaPtNegVz);
635
636
637   for(Int_t i=0;i<7;i++){
638     fHists->Add(fSignDcaPos[i]);
639     fHists->Add(fSignDcaNeg[i]);
640     fHists->Add(fSignDcaNegInv[i]);
641  
642     fHists->Add(fPtSigmaPos[i]);
643     fHists->Add(fPtSigmaNeg[i]);
644     fHists->Add(fqPtRec[i]);
645   
646     fHists->Add(fDcaSigmaPos[i]);
647     fHists->Add(fDcaSigmaNeg[i]);
648  
649
650   } 
651   
652     
653     
654   for (Int_t i=0; i<fHists->GetEntries(); ++i) {
655     TH1 *h1 = dynamic_cast<TH1*>(fHists->At(i));
656     if (h1){
657       // Printf("%s ",h1->GetName());
658       h1->Sumw2();
659     }
660   }
661     // BKC
662
663
664 }
665
666 //__________________________________________________________
667
668 void AliAnalysisTaskQASym::UserExec(Option_t *) 
669 {
670   printf("I'm here \n");
671   AliVEvent *event = InputEvent();
672   if (!event) {
673      Printf("ERROR: Could not retrieve event");
674      return;
675   }
676
677
678   if(Entry()==0){
679     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
680     if(esd){
681       Printf("We are reading from ESD");
682     }
683    
684   }
685
686   Printf("There are %d tracks in this event", event->GetNumberOfTracks());
687
688   
689   Int_t   leadingTrack  =   0;
690   Float_t leadingEnergy = -20.;
691   Float_t leadingPhi    =   0;//TMath::Pi();
692
693
694   if(event->GetNumberOfTracks()!=0) fNumber->Fill(event->GetNumberOfTracks());
695
696   const AliVVertex* vertex = event->GetPrimaryVertex();
697   Float_t vz = vertex->GetZ();
698   if (TMath::Abs(vz) > 10.) return;
699
700   for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) {
701     
702
703     AliVParticle *track = event->GetTrack(iTrack);
704     AliESDtrack *esdtrack =  dynamic_cast<AliESDtrack*>(track);
705     if (!track) {
706       Printf("ERROR: Could not receive track %d", iTrack);
707       continue;
708     }
709     
710     //if (!fCuts->AcceptTrack(esdtrack)) continue;
711     const AliExternalTrackParam * tpcPSO = esdtrack->GetTPCInnerParam();
712     const AliExternalTrackParam *  tpcP = esdtrack;
713     if (!tpcP) continue;
714    
715
716     if (tpcP->Pt() > 50. && tpcPSO) {
717       printf("High Pt %5d %5d %13.3f %13.3f \n", event->GetPeriodNumber(), event->GetOrbitNumber(),
718              tpcPSO->Pt(), esdtrack->Pt());
719       // AliFatal("Jet");
720     } 
721     
722 //    if (tpcPSO) fRecPt12->Fill(tpcPSO->Pt(), esdtrack->Pt());
723
724     if(tpcP->E()>leadingEnergy){
725       leadingTrack=iTrack;
726       leadingEnergy=tpcP->E();
727       leadingPhi=tpcP->Phi();
728     }
729    
730
731
732     
733     //propagate to dca
734     esdtrack->PropagateToDCA(event->GetPrimaryVertex(),
735                              event->GetMagneticField(), 10000.);
736     
737     
738     // if(tpcP->Pt()<2.)continue;
739
740
741     fqRec->Fill(tpcP->Charge());
742   
743
744     Double_t sigmapt = tpcP->GetSigma1Pt2();
745     sigmapt= sqrt(sigmapt);
746     sigmapt= sigmapt *(tpcP->Pt()*tpcP->Pt()); 
747
748     if(sigmapt < 1.e-12) continue;
749     fsigmaPt->Fill(TMath::Log10(sigmapt));
750    
751
752     // hits in ITS layer
753     Int_t cas=-1;
754     if(esdtrack->HasPointOnITSLayer(0)) 
755       cas=0;
756     else if(!esdtrack->HasPointOnITSLayer(0)
757             &&  esdtrack->HasPointOnITSLayer(1)) 
758       cas=1;
759     else if(!esdtrack->HasPointOnITSLayer(0)
760             && !esdtrack->HasPointOnITSLayer(1) 
761             &&  esdtrack->HasPointOnITSLayer(2)) 
762       cas=2;
763     else if(!esdtrack->HasPointOnITSLayer(0)
764             && !esdtrack->HasPointOnITSLayer(1) 
765             && !esdtrack->HasPointOnITSLayer(2)
766             &&  esdtrack->HasPointOnITSLayer(3)) 
767       cas=3;
768     else if(!esdtrack->HasPointOnITSLayer(0)
769             && !esdtrack->HasPointOnITSLayer(1) 
770             && !esdtrack->HasPointOnITSLayer(2)
771             && !esdtrack->HasPointOnITSLayer(3)
772             &&  esdtrack->HasPointOnITSLayer(4)) 
773       cas=4;
774     else if(   !esdtrack->HasPointOnITSLayer(0)
775             && !esdtrack->HasPointOnITSLayer(1)
776             && !esdtrack->HasPointOnITSLayer(2)
777             && !esdtrack->HasPointOnITSLayer(3)
778             && !esdtrack->HasPointOnITSLayer(4) 
779             &&  esdtrack->HasPointOnITSLayer(5)) 
780       cas=5;
781     else 
782       cas=6;
783   
784    
785    
786     //------------------- 
787
788     Double_t sdcatr = (tpcP->Py()*tpcP->Xv()
789                - tpcP->Px()*tpcP->Yv())/tpcP->Pt();
790   
791
792
793     fqPtRec[cas]->Fill(tpcP->Charge()/tpcP->Pt());
794     
795     fHistRECpt->Fill(tpcP->Pt());
796     fEta->Fill(tpcP->Eta());
797     fThetaRec->Fill(tpcP->Theta());
798     fPhiRec->Fill(tpcP->Phi());
799     fVx->Fill(tpcP->Xv());
800     fVy->Fill(tpcP->Yv());
801     fVz->Fill(tpcP->Zv());
802   
803
804     fEtaPt->Fill(tpcP->Eta()/tpcP->Pt());
805     fQPt->Fill(tpcP->Charge()/tpcP->Pt());
806     fDca->Fill(sdcatr);
807
808
809
810     Float_t xy, z;
811     esdtrack->GetImpactParameters(xy,z);
812     fDiffDcaD->Fill(sdcatr+xy);
813    
814     
815     //for positive particles
816
817     if(tpcP->Charge()>0){
818       fRecPtPos->Fill(tpcP->Pt());
819       fRecPtPosLadder[cas]->Fill(tpcP->Pt());
820       fRecPtPosVz->Fill(TMath::Log10(tpcP->Pt()),tpcP->Zv());
821       fRecPhiPos->Fill(TMath::RadToDeg()*tpcP->Phi());
822     
823      
824       fRecPhiPosLadder[cas]->Fill(TMath::RadToDeg()*tpcP->Phi());
825       fRecPhiPosVz->Fill(TMath::RadToDeg()*tpcP->Phi(),tpcP->Zv());
826       fSignedDcaPosVz->Fill(sdcatr,tpcP->Zv());
827
828       fRecEtaPos->Fill(tpcP->Eta());
829       fRecEtaPosLadder[cas]->Fill(tpcP->Eta());
830       fRecEtaPtPos->Fill(tpcP->Eta()/tpcP->Pt());
831       fRecEtaPosVz->Fill(tpcP->Eta(),tpcP->Zv());
832       fRecEtaPtPosVz->Fill(tpcP->Eta()/tpcP->Pt(),tpcP->Zv());
833      
834       fRecDcaPos->Fill(sdcatr);
835       fRecDPos->Fill(xy);
836       fSignDcaPos[cas]->Fill(sdcatr);
837     
838      
839       fDcaSigmaPos[cas]->Fill(sdcatr, TMath::Log10(sigmapt));
840     
841       fPtSigmaPos[cas]->Fill(TMath::Log10(sigmapt));
842     }
843     //and negative particles
844     else {
845       fRecPtNeg->Fill(tpcP->Pt());
846       fRecPtNegLadder[cas]->Fill(tpcP->Pt());
847       fRecPtNegVz->Fill(TMath::Log10(tpcP->Pt()),tpcP->Zv());
848            
849       fRecPhiNeg->Fill(TMath::RadToDeg()*tpcP->Phi());
850       fRecPhiNegLadder[cas]->Fill(TMath::RadToDeg()*tpcP->Phi());
851       fRecPhiNegVz->Fill(TMath::RadToDeg()*tpcP->Phi(),tpcP->Zv());
852       fSignedDcaNegVz->Fill(sdcatr,tpcP->Zv());
853       fRecEtaPtNegVz->Fill(tpcP->Eta()/tpcP->Pt(),tpcP->Zv());
854
855       fRecEtaNeg->Fill(tpcP->Eta());
856       fRecEtaNegLadder[cas]->Fill(tpcP->Eta());
857       fRecEtaPtNeg->Fill(tpcP->Eta()/tpcP->Pt());
858       fRecEtaNegVz->Fill(tpcP->Eta(),tpcP->Zv());
859      
860       fRecDcaNeg->Fill(sdcatr);
861       fRecDcaNegInv->Fill(-sdcatr);
862       fRecDNeg->Fill(xy);
863       fSignDcaNeg[cas]->Fill(sdcatr);
864       fSignDcaNegInv[cas]->Fill(-sdcatr);
865      
866      
867       fDcaSigmaNeg[cas]->Fill(sdcatr,TMath::Log10(sigmapt));
868    
869       fPtSigmaNeg[cas]->Fill(TMath::Log10(sigmapt));
870     }
871     
872
873
874     //all particles with positive eta
875     if(tpcP->Eta()>0){
876       fRecQPtPosEta->Fill(tpcP->Charge()/tpcP->Pt());
877       fRecPtPosEta->Fill(tpcP->Pt());
878       fRecPhiPosEta->Fill(TMath::RadToDeg()*tpcP->Phi());
879       fRecQPtPosEtaVz->Fill(tpcP->Charge()/tpcP->Pt(),tpcP->Zv());
880       fRecDcaPosEta->Fill(sdcatr);
881       fRecDPosEta->Fill(xy);
882     }
883     //all particles with negative eta (and eta==0)
884     else{
885       fRecQPtNegEta->Fill(tpcP->Charge()/tpcP->Pt());
886       fRecPtNegEta->Fill(tpcP->Pt());
887       fRecPhiNegEta->Fill(TMath::RadToDeg()*tpcP->Phi());
888       fRecQPtNegEtaVz->Fill(tpcP->Charge()/tpcP->Pt(),tpcP->Zv());
889       fRecDcaNegEta->Fill(sdcatr);
890       fRecDNegEta->Fill(xy);
891
892     }
893      
894
895     //spectren detected by TPC sectors
896     //pt cut on 1 GeV/c ?!
897     // if(tpcP->Pt()<1.) continue;
898     fRecPtTpcSector[Int_t(tpcP->Phi()*
899                           TMath::RadToDeg()/20)]->Fill(tpcP->Pt());
900     fRecEtaTpcSector[Int_t(tpcP->Phi()*
901                           TMath::RadToDeg()/20)]->Fill(tpcP->Eta());
902     fSignedDcaTpcSector[Int_t(tpcP->Phi()*
903                           TMath::RadToDeg()/20)]->Fill(sdcatr); 
904     fRecQPtTpcSector[Int_t(tpcP->Phi()*
905                           TMath::RadToDeg()/20)]->Fill(tpcP->Charge()/tpcP->Pt());
906     fRecEtaPtTpcSector[Int_t(tpcP->Phi()*
907                           TMath::RadToDeg()/20)]->Fill(tpcP->Eta()/tpcP->Pt());
908      
909
910
911
912
913
914   // another track loop
915     for (Int_t iTrack2 = 0; iTrack2 < event->GetNumberOfTracks(); iTrack2++) {
916       
917       if(leadingTrack==iTrack2) continue;
918
919       AliVParticle *track2 = event->GetTrack(iTrack2);
920       AliESDtrack* esdtrack2 =  dynamic_cast<AliESDtrack*>(track2);
921       if (!track2) {
922         Printf("ERROR: Could not receive track %d", iTrack);
923         continue;
924       }
925       if (!fCuts->AcceptTrack(esdtrack2)) continue;
926       //propagate to dca
927       esdtrack2->PropagateToDCA(event->GetPrimaryVertex(),
928                                event->GetMagneticField(), 10000.);
929  
930       fDeltaPhiLeading->Fill((leadingPhi-esdtrack2->Phi())*TMath::RadToDeg(),
931                              leadingPhi*TMath::RadToDeg() );
932
933      
934
935     }//second track loop
936   }//first track loop
937
938   
939  
940
941   // Post output data.
942   // PostData(1, fHistPt);
943   PostData(1, fHists);
944 }      
945
946
947
948
949
950 //________________________________________________________________________
951 void AliAnalysisTaskQASym::Terminate(Option_t *) 
952 {
953   // Terminate
954
955 }  
956
957
958
959
960