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