- coding rule violations corrected
[u/mrichter/AliRoot.git] / PWG1 / cosmic / AliAnalysisTaskCosmic.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2007, 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 // Analysis Task for the Quality Assurance of Cosmic Data
19 // Two track segments in the are matched in angle and charged. T
20 // The quality of the matching in is checked by comparing 
21 // the tarnsverse momenta and starting points of the track segments
22 //
23 // Author
24 // Andreas Morsch
25 // andreas.morsch@cern.ch
26
27 #include "TChain.h"
28 #include "TTree.h"
29
30 #include "TH1F.h"
31 #include "TH2F.h"
32 #include "TProfile.h"
33 #include "TList.h"
34
35 #include "AliESDEvent.h"
36 #include "AliESDtrack.h"
37 #include "AliAnalysisTaskCosmic.h"
38
39 ClassImp(AliAnalysisTaskCosmic)
40
41 //________________________________________________________________________
42 AliAnalysisTaskCosmic::AliAnalysisTaskCosmic(const char *name) 
43     : AliAnalysisTaskSE(name),
44       fHists(0),
45       fhDZvsZ(0),
46       fhDZvsPhi(0),
47       fhCh1Ch2(0),
48       fhPh1Ph2(0),
49       fhCl1Cl2G(0),
50       fhCl1Cl2B(0)
51 {
52   //
53   // Constructor
54   //
55     
56     for (Int_t i = 0; i < 6; i++) {
57         fhPt[i]      = 0;
58         fhTheta[i]   = 0;
59         fhPhi[i]     = 0;
60         fhDPhi[i]    = 0;
61         fhDTheta[i]  = 0;
62         fhDZ[i]      = 0;
63         fhDX[i]      = 0;
64         fhDY[i]      = 0;
65         fhDPt[i]     = 0;
66         fhD1ovPt[i]  = 0;
67         fpDPt[i]     = 0;
68         fhDPtovPt[i] = 0;
69         fpDPtS[i]    = 0;
70     }
71     
72   DefineOutput(1,  TList::Class());
73 }
74
75 AliAnalysisTaskCosmic::AliAnalysisTaskCosmic(const AliAnalysisTaskCosmic& task)
76     : AliAnalysisTaskSE(task),
77       fHists(0),
78       fhDZvsZ(0),
79       fhDZvsPhi(0),
80       fhCh1Ch2(0),
81       fhPh1Ph2(0),
82       fhCl1Cl2G(0),
83       fhCl1Cl2B(0)
84 {
85   // Copy constructor
86     fHists       = new TList();
87     fhDZvsZ      = (TH2F*)    ((task.fhDZvsZ)   ->Clone()); 
88     fhDZvsPhi    = (TH2F*)    ((task.fhDZvsPhi) ->Clone());
89     fhCh1Ch2     = (TH2F*)    ((task.fhCh1Ch2)  ->Clone());
90     fhPh1Ph2     = (TH2F*)    ((task.fhPh1Ph2)  ->Clone());
91     fhCl1Cl2G    = (TH2F*)    ((task.fhCl1Cl2G) ->Clone());
92     fhCl1Cl2B    = (TH2F*)    ((task.fhCl1Cl2B) ->Clone());
93     for (Int_t i = 0; i < 6; i++) {
94         fhPt[i]       =    (TH1F*) ((task.fhPt[i])        ->Clone());
95         fhTheta[i]    =    (TH1F*) ((task.fhTheta[i])     ->Clone());
96         fhPhi[i]      =    (TH1F*) ((task.fhPhi[i])       ->Clone());
97         fhDPhi[i]     =    (TH1F*) ((task.fhDPhi[i])      ->Clone());
98         fhDTheta[i]   =    (TH1F*) ((task.fhDTheta[i])    ->Clone());
99         fhDZ[i]       =    (TH1F*) ((task.fhDZ[i])        ->Clone());
100         fhDX[i]       =    (TH1F*) ((task.fhDX[i])        ->Clone());
101         fhDY[i]       =    (TH1F*) ((task.fhDY[i])        ->Clone());
102         fhDPt[i]      =    (TH1F*) ((task.fhDPt[i])       ->Clone());
103         fhD1ovPt[i]   =    (TH1F*) ((task.fhD1ovPt[i])    ->Clone());
104         fhDPtovPt[i]  =    (TH1F*) ((task.fhDPtovPt[i])   ->Clone());
105
106         fpDPt[i]      =    (TProfile*) ((task.fpDPt[i])       ->Clone());
107         fpDPtS[i]     =    (TProfile*) ((task.fpDPtS[i])      ->Clone());
108     }
109 }
110
111 //______________________________________________________________________________
112 AliAnalysisTaskCosmic& AliAnalysisTaskCosmic::operator=(const AliAnalysisTaskCosmic& task)
113 {
114   // Assignment operator
115   if(this!=&task) {
116
117     AliAnalysisTaskSE::operator=(task);
118
119     fHists       = new TList();
120     fhDZvsZ      = (TH2F*)    ((task.fhDZvsZ)   ->Clone()); 
121     fhDZvsPhi    = (TH2F*)    ((task.fhDZvsPhi) ->Clone());
122     fhCh1Ch2     = (TH2F*)    ((task.fhCh1Ch2)  ->Clone());
123     fhPh1Ph2     = (TH2F*)    ((task.fhPh1Ph2)  ->Clone());
124     fhCl1Cl2G    = (TH2F*)    ((task.fhCl1Cl2G) ->Clone());
125     fhCl1Cl2B    = (TH2F*)    ((task.fhCl1Cl2B) ->Clone());
126     for (Int_t i = 0; i < 6; i++) {
127         fhPt[i]       =    (TH1F*) ((task.fhPt[i])        ->Clone());
128         fhTheta[i]    =    (TH1F*) ((task.fhTheta[i])     ->Clone());
129         fhPhi[i]      =    (TH1F*) ((task.fhPhi[i])       ->Clone());
130         fhDPhi[i]     =    (TH1F*) ((task.fhDPhi[i])      ->Clone());
131         fhDTheta[i]   =    (TH1F*) ((task.fhDTheta[i])    ->Clone());
132         fhDZ[i]       =    (TH1F*) ((task.fhDZ[i])        ->Clone());
133         fhDX[i]       =    (TH1F*) ((task.fhDX[i])        ->Clone());
134         fhDY[i]       =    (TH1F*) ((task.fhDY[i])        ->Clone());
135         fhDPt[i]      =    (TH1F*) ((task.fhDPt[i])       ->Clone());
136         fhD1ovPt[i]   =    (TH1F*) ((task.fhD1ovPt[i])    ->Clone());
137         fhDPtovPt[i]  =    (TH1F*) ((task.fhDPtovPt[i])   ->Clone());
138
139         fpDPt[i]      =    (TProfile*) ((task.fpDPt[i])       ->Clone());
140         fpDPtS[i]     =    (TProfile*) ((task.fpDPtS[i])      ->Clone());
141     }
142   }
143   return *this;
144 }
145
146
147 //________________________________________________________________________
148 void AliAnalysisTaskCosmic::UserCreateOutputObjects()
149 {
150   // Create histograms
151   // Called once
152     TString ext[6] = {"PC", "NC", "PZ", "NZ", "_Good", "_Bad"};
153     
154     char name[12];
155     
156     fHists = new TList();
157
158     for (Int_t i = 0; i < 6; i++) {
159         // Pt
160         sprintf(name, "fhPt%2s", ext[i].Data());
161         fhPt[i]   = new TH1F(name,   " pT distribution",        800, 0., 200.);
162         fhPt[i]->SetXTitle("p_{T} [Gev]");
163         // Phi
164         sprintf(name, "fhPhi%2s", ext[i].Data());
165         fhPhi[i]   = new TH1F(name,   "Phi distribution",        62,  0., 2. * TMath::Pi());
166         fhPhi[i]->SetXTitle("#phi [rad]");
167         // Theta
168         sprintf(name, "fhTheta%2s", ext[i].Data());
169         fhTheta[i]   = new TH1F(name,   "Theta distribution",    62,  0., TMath::Pi());    
170         fhTheta[i]->SetXTitle("#theta [rad]");
171         // Delta Phi
172         sprintf(name, "fhDPhi%2s", ext[i].Data());
173         fhDPhi[i]   = new TH1F(name,   "DeltaPhi distribution",        320, -0.4, 0.4);
174         fhDPhi[i]->SetXTitle("#Delta#phi [rad]");
175         // Delta Theta
176         sprintf(name, "fhDTheta%2s", ext[i].Data());
177         fhDTheta[i]   = new TH1F(name,   "DeltaTheta distribution",    320, -0.4, 0.4);    
178         fhDTheta[i]->SetXTitle("#Delta#theta [rad]");
179         // Delta Z
180         sprintf(name, "fhDZ%2s", ext[i].Data());
181         fhDZ[i]   = new TH1F(name,   "DeltaZ distribution",        200, -10., 10.);
182         fhDZ[i]->SetXTitle("#DeltaZ [cm]");
183         // Delta X
184         sprintf(name, "fhDX%2s", ext[i].Data());
185         fhDX[i]   = new TH1F(name,   "DeltaX distribution",        200, -10., 10.);
186         fhDX[i]->SetXTitle("#DeltaX [cm]");
187         // Delta Y
188         sprintf(name, "fhDY%2s", ext[i].Data());
189         fhDY[i]   = new TH1F(name,   "DeltaY distribution",        200, -10, 10.);
190         fhDY[i]->SetXTitle("#DeltaY [cm]");
191         // Delta Pt
192         sprintf(name, "fhDPt%2s", ext[i].Data());
193         fhDPt[i]   = new TH1F(name,   "DeltaPt distribution",        200, -20., 20.);
194         fhDPt[i]->SetXTitle("#Delta p_{T}  [GeV]");
195
196         // Delta 1/Pt
197         sprintf(name, "fhD1ovPt%2s", ext[i].Data());
198         fhD1ovPt[i]   = new TH1F(name,   "Delta 1/Pt distribution",        200, -1., 1.);
199         fhD1ovPt[i]->SetXTitle("#Delta 1/Pt");
200
201         // Delta Pt over Pt
202         sprintf(name, "fhDPtovPt%2s", ext[i].Data());
203         fhDPtovPt[i]   = new TH1F(name,   "DeltaPt/Pt distribution",        200, -2., 2.);
204         fhDPtovPt[i]->SetXTitle("#DeltaPt/Pt");
205
206
207
208         // Delta Pt/ Pt vs Pt
209         sprintf(name, "fpDPt%2s", ext[i].Data());
210         fpDPt[i] = new TProfile(name, "#Delta Pt / Pt", 20, 0., 20., -1, 1., "S");
211         fpDPt[i]->SetXTitle("p_{T} [GeV]");
212         fpDPt[i]->SetYTitle("#Delta 1/p_{T} [GeV^{-1}]");
213         // Delta Pt error
214         sprintf(name, "fpDPtS%2s", ext[i].Data());
215         fpDPtS[i] = new TProfile(name, "#Delta Pt / <sigma>", 20, 0., 20., 0., 10.);
216         fpDPtS[i]->SetXTitle("p_{T}");
217         fpDPtS[i]->SetYTitle("#Delta p_{T} / <#sigma_{p_{T}}>");
218
219
220         fHists->Add(fhPt[i]);
221         fHists->Add(fhPhi[i]);
222         fHists->Add(fhTheta[i]);
223         fHists->Add(fhDPhi[i]);
224         fHists->Add(fhDPt[i]);
225         fHists->Add(fhD1ovPt[i]);
226         fHists->Add(fhDPtovPt[i]);
227         fHists->Add(fhDTheta[i]);
228         fHists->Add(fhDZ[i]);
229         fHists->Add(fhDX[i]);
230         fHists->Add(fhDY[i]);
231         fHists->Add(fpDPt[i]);
232         fHists->Add(fpDPtS[i]);
233     }
234
235     fhDZvsZ      = new TH2F("fhDZvsZ",    "dz vs z", 500, -250., 250., 100, -10., 10.);
236     fhDZvsZ->SetXTitle("z_{in} * sign(z_{in}) * sign(z_{out}) [cm]");
237     fhDZvsZ->SetYTitle("#Delta z [cm]");
238     
239     
240     fhDZvsPhi    = new TH2F("fhDZvsPhi", "dz vs phi", 36, 0., TMath::Pi(),  50, -2., 2.);
241
242     fhCh1Ch2   = new TH2F("fCh1Ch2",   "ch1 vs ch2", 8, -2., 2., 8, -2., 2.);
243     fhPh1Ph2   = new TH2F("fPh1Ph2",   "ph1 vs ph2", 128, 0., 2. * TMath::Pi(), 128, 0., 2. * TMath::Pi());
244     fhCl1Cl2G  = new TH2F("fCl1Cl2G",   "#cl vs #cl", 200, 0., 200., 200, 0., 200);
245     fhCl1Cl2B  = new TH2F("fCl1Cl2B",   "#cl vs #cl", 200, 0., 200., 200, 0., 200);    
246
247     fHists->Add(fhDZvsZ);
248     fHists->Add(fhDZvsPhi);
249     fHists->Add(fhCh1Ch2);
250     fHists->Add(fhPh1Ph2);
251     fHists->Add(fhCl1Cl2G);
252     fHists->Add(fhCl1Cl2B);
253     fHists->SetOwner();
254
255 }
256
257 //________________________________________________________________________
258 void AliAnalysisTaskCosmic::UserExec(Option_t *) 
259 {
260   // Main loop
261   // Called for each event
262
263   if (!fInputEvent) {
264     Printf("ERROR: fESD not available");
265     return;
266   }
267   
268   AliESDEvent* esdE = (AliESDEvent*)fInputEvent;
269   
270   if (esdE->GetNumberOfTracks() != 2) {
271       PostData(1, fHists);
272       return;
273   }
274   
275   
276   for (Int_t iTracks = 0; iTracks < esdE->GetNumberOfTracks(); iTracks++) {
277       AliESDtrack* track = esdE->GetTrack(iTracks);
278
279
280       
281 //      const AliExternalTrackParam * track = trackG->GetTPCInnerParam();
282 //      if (!track) continue;
283       // Pt cut
284       Float_t pt     = track->Pt();
285       Short_t charge = track->Charge();
286       Float_t phi    = track->Phi();
287       if (phi > 0. && phi  < TMath::Pi()) charge*=(-1);
288
289
290
291       // TPC track
292       UInt_t status = track->GetStatus();
293       if ((status&AliESDtrack::kTPCrefit) ==0) continue;
294
295       Int_t nClustersTPC = track->GetTPCclusters(0);
296       if (nClustersTPC < 50) continue;
297       
298       // 
299
300       Float_t z      = track->GetZ();
301       Float_t x      = track->Xv();
302       Float_t y      = track->Yv();
303       Float_t theta  = track->Theta();
304
305
306       const AliExternalTrackParam * trackOut = track->GetOuterParam();
307       Float_t zOut = 0.;
308       if (trackOut)zOut = trackOut->Zv();
309       
310
311       if (charge > 0) {
312           fhPt[kPosC]   ->Fill(pt);
313       } else {
314           fhPt[kNegC]   ->Fill(pt);
315       }
316
317       //      if ((TMath::Abs(esdE->GetCurrentL3()) > 1.e-3) && (pt < 1. || pt > 50.)) continue;
318
319       
320       //
321       if (charge > 0) {
322           fhPhi[kPosC]  ->Fill(phi);
323           fhTheta[kPosC]->Fill(theta);
324       } else {
325           fhPhi[kNegC]  ->Fill(phi);
326           fhTheta[kNegC]->Fill(theta);
327       }
328
329       
330       if (z > 0) {
331           fhPt[kPosZ]   ->Fill(pt);
332           fhPhi[kPosZ]  ->Fill(phi);
333           fhTheta[kPosZ]->Fill(theta);
334       } else {
335           fhPt[kNegZ]   ->Fill(pt);
336           fhPhi[kNegZ]  ->Fill(phi);
337           fhTheta[kNegZ]->Fill(theta);
338       }
339
340
341
342       // Tracks coming from above
343       if (phi > 0. &&  phi < TMath::Pi()) {
344           
345 //        printf("Track#1 %5d %5d %13.3f %13.3f \n", Int_t(Entry()), iTracks, phi, zOut);
346           
347           Float_t dphiMin = 999.;
348           Float_t rMin    = 999.;
349           Int_t   jMin    = -1;
350           // Search for a matching track (in dphi-dtheta)
351           for (Int_t jTracks = 0; jTracks < esdE->GetNumberOfTracks(); jTracks++) {
352               if (jTracks == iTracks) continue;
353
354               AliESDtrack* track2 = esdE->GetTrack(jTracks);
355
356               UInt_t status = track2->GetStatus();
357               if ((status&AliESDtrack::kTPCrefit) ==0) continue;
358
359               Int_t nClustersTPC2 = track2->GetTPCclusters(0);
360               if (nClustersTPC2 < 50) continue;
361               //              if ((TMath::Abs(esdE->GetCurrentL3()) > 1.e-3) && (track2->Pt() < 1. || track2->Pt() > 50.)) continue;
362 //            const AliExternalTrackParam * track2 = trackG2->GetTPCInnerParam();
363 //            if (!track2) continue;
364               
365               Float_t phi2   = track2->Phi() - TMath::Pi();
366               Float_t theta2 = TMath::Pi()   - track2->Theta();
367               Float_t dphi   = phi2   - phi;
368               Float_t dtheta = theta2 - theta;
369               
370               if (dphi >  TMath::Pi()) dphi -= 2. * TMath::Pi();
371               if (dphi < -TMath::Pi()) dphi += 2. * TMath::Pi();
372
373               Float_t dR = TMath::Sqrt(dphi * dphi + dtheta * dtheta);
374
375               if (dR < rMin) {
376                   rMin    = dR;
377                   dphiMin = dphi;
378                   jMin = jTracks;
379               }
380
381           } // tracks 2
382           if (jMin != -1) {
383               // we found a matching track candidate ...
384               AliESDtrack* track2 = esdE->GetTrack(jMin);
385               const AliExternalTrackParam * trackOut2 = track2->GetOuterParam();
386               Float_t zOut2 = 0.;
387               if (trackOut2) zOut2 = trackOut2->Zv();
388
389
390               Float_t theta2 = - track2->Theta() + TMath::Pi();
391               Float_t z2 = track2->GetZ();
392               Float_t x2 = track2->Xv();
393               Float_t y2 = track2->Yv();
394               Short_t charge2 = track2->Charge();
395               Float_t dz  = z2 - z;
396               Float_t dx  = x2 - x;
397               Float_t dy  = y2 - y;
398               Float_t pt1 = track->Pt();
399               Float_t pt2 = track2->Pt();
400               Float_t dpt = pt2 - pt1;
401               Float_t d1pt = 1./pt2 - 1./pt1;
402               
403               Float_t ptm = 0.5 * (pt1 + pt2);
404
405               Int_t nClustersTPC2 = track2->GetTPCclusters(0);
406
407               if (charge > 0.) {
408                   fhDPhi[kPosC]  ->Fill(dphiMin);
409                   fhDTheta[kPosC]->Fill(theta2 - theta);
410               } else {
411                   fhDPhi[kNegC]  ->Fill(dphiMin);
412                   fhDTheta[kNegC]->Fill(theta2 - theta);
413               }
414
415               if (z > 0.) {
416                   fhDPhi[kPosZ]  ->Fill(dphiMin);
417                   fhDTheta[kPosZ]->Fill(theta2 - theta);
418               } else {
419                   fhDPhi[kNegZ]  ->Fill(dphiMin);
420                   fhDTheta[kNegZ]->Fill(theta2 - theta);
421               }
422
423               if (TMath::Abs(dpt)/ptm > 0.5) {
424                   fhDPhi[kBad]  ->Fill(dphiMin);
425                   fhDTheta[kBad]->Fill(theta2 - theta);
426               } else {
427                   fhDPhi[kGood]  ->Fill(dphiMin);
428                   fhDTheta[kGood]->Fill(theta2 - theta);
429               }
430               // Good matches ...             
431 //            if (TMath::Abs(rMin < 0.04) && (charge == charge2) && TMath::Abs(dz) < 0.5 &&  TMath::Abs(dy) && TMath::Abs(dx) < 0.5 ) 
432 //            if (TMath::Abs(rMin < 0.04) && (charge == charge2) && (zOut * zOut2) < 0.) 
433               if (TMath::Abs(rMin < 0.04) && (charge == charge2)) 
434               {
435
436                   if (TMath::Abs(dpt)/ptm  < .1) {
437                       fhCl1Cl2G->Fill(nClustersTPC, nClustersTPC2);
438                   }
439                   
440                   if (TMath::Abs(dpt)/ptm  > .5) fhCl1Cl2B->Fill(nClustersTPC, nClustersTPC2);
441                   fhPh1Ph2->Fill(track->Phi(), track2->Phi());
442
443
444                   Double_t sigmaPt1 = TMath::Sqrt(track->GetSigma1Pt2());
445                   Double_t sigmaPt2 = TMath::Sqrt(track->GetSigma1Pt2());
446                   Double_t sigmaPt = 0.5 * TMath::Sqrt(sigmaPt1 * sigmaPt1 + sigmaPt2 * sigmaPt2);
447                   if (TMath::Abs(dpt)/ptm > 0.2 && pt > 10. &&  charge > 0.)
448 //                    printf("Track#2 %5d %5d %13.3f %13.3f %13.3f %13.3f\n", Entry(), jMin, track2->Phi(), dz, dpt, zOut2);                      
449                   if (TMath::Abs(zOut) > 1.e-10 && TMath::Abs(zOut2) > 1.e-10) fhDZvsZ->Fill(TMath::Abs(zOut) * zOut/zOut2, dz);
450                   if (zOut * zOut2 > 0. && zOut < 0) fhDZvsPhi->Fill(phi, dz);
451                   
452                   fhCh1Ch2->Fill(charge, track2->Charge());
453                   
454                   
455                   if (charge > 0.) {
456                       fhDPt[kPosC]->Fill(dpt);
457                       fhD1ovPt[kPosC]->Fill(d1pt);
458                       fpDPt[kPosC]->Fill(ptm, d1pt);
459                       if (ptm > 5. && ptm < 15.) fhDPtovPt[kPosC]->Fill(dpt/ptm);
460                       fhDZ[kPosC]->Fill(dz);
461                       fhDX[kPosC]->Fill(dx);
462                       fhDY[kPosC]->Fill(dy);
463                       fpDPtS[kPosC]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
464                   } else {
465                       fhDPt[kNegC]->Fill(dpt);
466                       fhD1ovPt[kNegC]->Fill(d1pt);
467                       fpDPt[kNegC]->Fill(ptm, d1pt);
468                       if (ptm > 5. && ptm < 15.) fhDPtovPt[kNegC]->Fill(dpt/ptm);
469                       fhDZ[kNegC]->Fill(dz);
470                       fhDX[kNegC]->Fill(dx);
471                       fhDY[kNegC]->Fill(dy);
472                       fpDPtS[kNegC]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
473                   }
474
475                   if (z > 0.) {
476                       fhDPt[kPosZ]->Fill(dpt);
477                       fhD1ovPt[kPosZ]->Fill(d1pt);
478                       fpDPt[kPosZ]->Fill(ptm, d1pt);
479                       fhDZ[kPosZ]->Fill(dz);
480                       fhDX[kPosZ]->Fill(dx);
481                       fhDY[kPosZ]->Fill(dy);
482                       fpDPtS[kPosZ]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
483                   } else {
484                       fhDPt[kNegZ]->Fill(dpt);
485                       fhD1ovPt[kNegZ]->Fill(d1pt);
486                       fpDPt[kNegZ]->Fill(ptm, d1pt);
487                       fhDZ[kNegZ]->Fill(dz);
488                       fhDX[kNegZ]->Fill(dx);
489                       fhDY[kNegZ]->Fill(dy);
490                       fpDPtS[kNegZ]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
491                   } 
492
493
494                   if (TMath::Abs(dpt)/ptm > 0.5) {
495                       fhDPt[kBad]->Fill(dpt);
496                       fpDPt[kBad]->Fill(ptm, TMath::Abs(dpt)/ptm);
497                       fhDZ[kBad]->Fill(dz);
498                       fhDX[kBad]->Fill(dx);
499                       fhDY[kBad]->Fill(dy);
500                       fpDPtS[kBad]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
501                   } else {
502                       fhDPt[kGood]->Fill(dpt);
503                       fpDPt[kGood]->Fill(ptm, TMath::Abs(dpt)/ptm);
504                       fhDZ[kGood]->Fill(dz);
505                       fhDX[kGood]->Fill(dx);
506                       fhDY[kGood]->Fill(dy);
507                       fpDPtS[kGood]->Fill(ptm, 2. * TMath::Abs(1./pt1 - 1./pt2)/sigmaPt);
508                   } 
509                   
510               } // good matches
511           } // found possible match
512       } // upper
513   } // tracks 1
514   PostData(1, fHists);
515 }      
516
517 //________________________________________________________________________
518 void AliAnalysisTaskCosmic::Terminate(Option_t *) 
519 {
520 }