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