]>
Commit | Line | Data |
---|---|---|
8d63376d | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-2009, 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 | ///////////////////////////////////////////////////////////// | |
17 | // | |
18 | // AliAnalysisTask to extract from ESD tracks the information | |
19 | // on ITS tracking efficiency and resolutions. | |
20 | // | |
21 | // Author: A.Dainese, andrea.dainese@pd.infn.it | |
22 | ///////////////////////////////////////////////////////////// | |
23 | ||
24 | #include <TChain.h> | |
25 | #include <TTree.h> | |
26 | #include <TNtuple.h> | |
27 | #include <TBranch.h> | |
28 | #include <TClonesArray.h> | |
29 | #include <TObjArray.h> | |
30 | #include <TH1F.h> | |
31 | #include <TH2F.h> | |
32 | #include <TCanvas.h> | |
33 | #include <TParticle.h> | |
34 | ||
35 | #include "AliAnalysisTask.h" | |
36 | #include "AliAnalysisManager.h" | |
37 | ||
38 | #include "AliMultiplicity.h" | |
39 | #include "AliVertexerTracks.h" | |
40 | #include "AliESDtrack.h" | |
41 | #include "AliExternalTrackParam.h" | |
42 | #include "AliESDVertex.h" | |
43 | #include "AliESDEvent.h" | |
44 | #include "AliESDfriend.h" | |
45 | #include "AliESDInputHandler.h" | |
46 | #include "AliESDInputHandlerRP.h" | |
47 | #include "AliTrackPointArray.h" | |
48 | #include "AliITSRecPoint.h" | |
49 | ||
50 | #include "AliMCEventHandler.h" | |
51 | #include "AliMCEvent.h" | |
52 | #include "AliStack.h" | |
53 | #include "AliLog.h" | |
54 | ||
55 | #include "AliGenEventHeader.h" | |
56 | #include "AliAnalysisTaskITSTrackingCheck.h" | |
57 | ||
58 | ||
59 | ClassImp(AliAnalysisTaskITSTrackingCheck) | |
60 | ||
61 | //________________________________________________________________________ | |
62 | AliAnalysisTaskITSTrackingCheck::AliAnalysisTaskITSTrackingCheck(const char *name) : | |
63 | AliAnalysisTask(name, "ITSTrackingCheckTask"), | |
64 | fReadMC(kTRUE), | |
65 | fReadRPLabels(kFALSE), | |
66 | fESD(0), | |
67 | fESDfriend(0), | |
68 | fOutput(0), | |
69 | fHistNtracks(0), | |
70 | fHistNclsITSMI(0), | |
71 | fHistNclsITSSA(0), | |
72 | fHistClusterMapITSMI(0), | |
73 | fHistClusterMapITSMIok(0), | |
74 | fHistClusterMapITSMIbad(0), | |
75 | fHistClusterMapITSMIskipped(0), | |
76 | fHistClusterMapITSMIoutinz(0), | |
77 | fHistClusterMapITSMInorefit(0), | |
78 | fHistClusterMapITSMInocls(0), | |
79 | fHistClusterMapITSSA(0), | |
80 | fHistClusterMapITSSAok(0), | |
81 | fHistClusterMapITSSAbad(0), | |
82 | fHistClusterMapITSSAskipped(0), | |
83 | fHistClusterMapITSSAoutinz(0), | |
84 | fHistClusterMapITSSAnorefit(0), | |
85 | fHistClusterMapITSSAnocls(0), | |
86 | fHistPhiTPC(0), | |
87 | fHistPtTPC(0), | |
88 | fHistPtITSMI2(0), | |
89 | fHistPtITSMI3(0), | |
90 | fHistPtITSMI4(0), | |
91 | fHistPtITSMI5(0), | |
92 | fHistPtITSMI6(0), | |
93 | fHistPtITSMISPD(0), | |
94 | fNtupleESDTracks(0), | |
95 | fNtupleITSAlignExtra(0), | |
96 | fNtupleITSAlignSPDTracklets(0) | |
97 | { | |
98 | // Constructor | |
99 | ||
100 | for(Int_t i=0; i<10; i++) fCountsPerPtBin[i]=0; | |
101 | ||
102 | // Define input and output slots here | |
103 | // Input slot #0 works with a TChain | |
104 | DefineInput(0, TChain::Class()); | |
105 | // Output slot #0 writes into a TList container | |
106 | DefineOutput(0, TList::Class()); //My private output | |
107 | } | |
108 | //________________________________________________________________________ | |
109 | AliAnalysisTaskITSTrackingCheck::~AliAnalysisTaskITSTrackingCheck() | |
110 | { | |
111 | // Destructor | |
112 | ||
113 | // histograms are in the output list and deleted when the output | |
114 | // list is deleted by the TSelector dtor | |
115 | ||
116 | if (fOutput) { | |
117 | delete fOutput; | |
118 | fOutput = 0; | |
119 | } | |
120 | } | |
121 | ||
122 | //________________________________________________________________________ | |
123 | void AliAnalysisTaskITSTrackingCheck::ConnectInputData(Option_t *) | |
124 | { | |
125 | // Connect ESD or AOD here | |
126 | // Called once | |
127 | ||
128 | TTree* tree = dynamic_cast<TTree*> (GetInputData(0)); | |
129 | if (!tree) { | |
130 | Printf("ERROR: Could not read chain from input slot 0"); | |
131 | } else { | |
132 | // Disable all branches and enable only the needed ones | |
133 | // The next two lines are different when data produced as AliESDEvent is read | |
134 | ||
135 | tree->SetBranchStatus("ESDfriend*", 1); | |
136 | tree->SetBranchAddress("ESDfriend.",&fESDfriend); | |
137 | ||
138 | AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); | |
139 | ||
140 | ||
141 | if (!esdH) { | |
142 | Printf("ERROR: Could not get ESDInputHandler"); | |
143 | } else { | |
144 | fESD = esdH->GetEvent(); | |
145 | ||
146 | } | |
147 | } | |
148 | ||
149 | return; | |
150 | } | |
151 | ||
152 | //________________________________________________________________________ | |
153 | void AliAnalysisTaskITSTrackingCheck::CreateOutputObjects() | |
154 | { | |
155 | // Create histograms | |
156 | // Called once | |
157 | ||
158 | for(Int_t i=0; i<10; i++) fCountsPerPtBin[i]=0; | |
159 | ||
160 | // Several histograms are more conveniently managed in a TList | |
161 | fOutput = new TList; | |
162 | fOutput->SetOwner(); | |
163 | ||
164 | fHistNtracks = new TH1F("fHistNtracks", "N ESD tracks; N tracks; Events",5000, -0.5, 4999.5); | |
165 | fHistNtracks->Sumw2(); | |
166 | fHistNtracks->SetMinimum(0); | |
167 | fOutput->Add(fHistNtracks); | |
168 | ||
169 | fHistNclsITSMI = new TH1F("fHistNclsITSMI", "N ITS clusters per track (MI); N clusters; Counts",7, -0.5, 6.5); | |
170 | fHistNclsITSMI->Sumw2(); | |
171 | fHistNclsITSMI->SetMinimum(0); | |
172 | fOutput->Add(fHistNclsITSMI); | |
173 | ||
174 | fHistNclsITSSA = new TH1F("fHistNclsITSSA", "N ITS clusters per track (SA); N clusters; Counts",7, -0.5, 6.5); | |
175 | fHistNclsITSSA->Sumw2(); | |
176 | fHistNclsITSSA->SetMinimum(0); | |
177 | fOutput->Add(fHistNclsITSSA); | |
178 | ||
179 | fHistClusterMapITSMI = new TH1F("fHistClusterMapITSMI", "N tracks with point on Layer (MI); Layer; N tracks",6, -0.5, 5.5); | |
180 | fHistClusterMapITSMI->Sumw2(); | |
181 | fHistClusterMapITSMI->SetMinimum(0); | |
182 | fOutput->Add(fHistClusterMapITSMI); | |
183 | ||
184 | fHistClusterMapITSSA = new TH1F("fHistClusterMapITSSA", "N tracks with point on Layer (SA); Layer; N tracks",6, -0.5, 5.5); | |
185 | fHistClusterMapITSSA->Sumw2(); | |
186 | fHistClusterMapITSSA->SetMinimum(0); | |
187 | fOutput->Add(fHistClusterMapITSSA); | |
188 | ||
189 | fHistClusterMapITSMIok = new TH1F("fHistClusterMapITSMIok", "N tracks with ok on Layer (MI); Layer; N tracks",6, -0.5, 5.5); | |
190 | fHistClusterMapITSMIok->Sumw2(); | |
191 | fHistClusterMapITSMIok->SetMinimum(0); | |
192 | fOutput->Add(fHistClusterMapITSMIok); | |
193 | ||
194 | fHistClusterMapITSSAok = new TH1F("fHistClusterMapITSSAok", "N tracks with ok on Layer (SA); Layer; N tracks",6, -0.5, 5.5); | |
195 | fHistClusterMapITSSAok->Sumw2(); | |
196 | fHistClusterMapITSSAok->SetMinimum(0); | |
197 | fOutput->Add(fHistClusterMapITSSAok); | |
198 | ||
199 | fHistClusterMapITSMIbad = new TH1F("fHistClusterMapITSMIbad", "N tracks with bad on Layer (MI); Layer; N tracks",6, -0.5, 5.5); | |
200 | fHistClusterMapITSMIbad->Sumw2(); | |
201 | fHistClusterMapITSMIbad->SetMinimum(0); | |
202 | fOutput->Add(fHistClusterMapITSMIbad); | |
203 | ||
204 | fHistClusterMapITSSAbad = new TH1F("fHistClusterMapITSSAbad", "N tracks with bad on Layer (SA); Layer; N tracks",6, -0.5, 5.5); | |
205 | fHistClusterMapITSSAbad->Sumw2(); | |
206 | fHistClusterMapITSSAbad->SetMinimum(0); | |
207 | fOutput->Add(fHistClusterMapITSSAbad); | |
208 | ||
209 | fHistClusterMapITSMIskipped = new TH1F("fHistClusterMapITSMIskipped", "N tracks with skip on Layer (MI); Layer; N tracks",6, -0.5, 5.5); | |
210 | fHistClusterMapITSMIskipped->Sumw2(); | |
211 | fHistClusterMapITSMIskipped->SetMinimum(0); | |
212 | fOutput->Add(fHistClusterMapITSMIskipped); | |
213 | ||
214 | fHistClusterMapITSSAskipped = new TH1F("fHistClusterMapITSSAskipped", "N tracks with skip on Layer (SA); Layer; N tracks",6, -0.5, 5.5); | |
215 | fHistClusterMapITSSAskipped->Sumw2(); | |
216 | fHistClusterMapITSSAskipped->SetMinimum(0); | |
217 | fOutput->Add(fHistClusterMapITSSAskipped); | |
218 | ||
219 | fHistClusterMapITSMIoutinz = new TH1F("fHistClusterMapITSMIoutinz", "N tracks out in z on Layer (MI); Layer; N tracks",6, -0.5, 5.5); | |
220 | fHistClusterMapITSMIoutinz->Sumw2(); | |
221 | fHistClusterMapITSMIoutinz->SetMinimum(0); | |
222 | fOutput->Add(fHistClusterMapITSMIoutinz); | |
223 | ||
224 | fHistClusterMapITSSAoutinz = new TH1F("fHistClusterMapITSSAoutinz", "N tracks with out in z on Layer (SA); Layer; N tracks",6, -0.5, 5.5); | |
225 | fHistClusterMapITSSAoutinz->Sumw2(); | |
226 | fHistClusterMapITSSAoutinz->SetMinimum(0); | |
227 | fOutput->Add(fHistClusterMapITSSAoutinz); | |
228 | ||
229 | fHistClusterMapITSMInorefit = new TH1F("fHistClusterMapITSMInorefit", "N tracks with norefit on Layer (MI); Layer; N tracks",6, -0.5, 5.5); | |
230 | fHistClusterMapITSMInorefit->Sumw2(); | |
231 | fHistClusterMapITSMInorefit->SetMinimum(0); | |
232 | fOutput->Add(fHistClusterMapITSMInorefit); | |
233 | ||
234 | fHistClusterMapITSSAnorefit = new TH1F("fHistClusterMapITSSAnorefit", "N tracks with norefit on Layer (SA); Layer; N tracks",6, -0.5, 5.5); | |
235 | fHistClusterMapITSSAnorefit->Sumw2(); | |
236 | fHistClusterMapITSSAnorefit->SetMinimum(0); | |
237 | fOutput->Add(fHistClusterMapITSSAnorefit); | |
238 | ||
239 | fHistClusterMapITSMInocls = new TH1F("fHistClusterMapITSMInocls", "N tracks with nocls on Layer (MI); Layer; N tracks",6, -0.5, 5.5); | |
240 | fHistClusterMapITSMInocls->Sumw2(); | |
241 | fHistClusterMapITSMInocls->SetMinimum(0); | |
242 | fOutput->Add(fHistClusterMapITSMInocls); | |
243 | ||
244 | fHistClusterMapITSSAnocls = new TH1F("fHistClusterMapITSSAnocls", "N tracks with nocls on Layer (SA); Layer; N tracks",6, -0.5, 5.5); | |
245 | fHistClusterMapITSSAnocls->Sumw2(); | |
246 | fHistClusterMapITSSAnocls->SetMinimum(0); | |
247 | fOutput->Add(fHistClusterMapITSSAnocls); | |
248 | ||
249 | fHistPhiTPC = new TH1F("fHistPhiTPC","Azimuthal distribution of TPC tracks; #phi; N tracks",100, -3.1415, 3.1415); | |
250 | fHistPhiTPC->Sumw2(); | |
251 | fHistPhiTPC->SetMinimum(0); | |
252 | fOutput->Add(fHistPhiTPC); | |
253 | ||
254 | fHistPtTPC = new TH1F("fHistPtTPC","pt distribution of TPC tracks; p_{t} [GeV/c]; N tracks",80, 0, 40); | |
255 | fHistPtTPC->Sumw2(); | |
256 | fHistPtTPC->SetMinimum(0); | |
257 | fOutput->Add(fHistPtTPC); | |
258 | ||
259 | fHistPtITSMI6 = new TH1F("fHistPtITSMI6","pt distribution of ITSMI6 tracks; p_{t} [GeV/c]; N tracks",80, 0, 40); | |
260 | fHistPtITSMI6->Sumw2(); | |
261 | fHistPtITSMI6->SetMinimum(0); | |
262 | fOutput->Add(fHistPtITSMI6); | |
263 | ||
264 | fHistPtITSMI5 = new TH1F("fHistPtITSMI5","pt distribution of ITSMI5 tracks; p_{t} [GeV/c]; N tracks",80, 0, 40); | |
265 | fHistPtITSMI5->Sumw2(); | |
266 | fHistPtITSMI5->SetMinimum(0); | |
267 | fOutput->Add(fHistPtITSMI5); | |
268 | ||
269 | fHistPtITSMI4 = new TH1F("fHistPtITSMI4","pt distribution of ITSMI4 tracks; p_{t} [GeV/c]; N tracks",80, 0, 40); | |
270 | fHistPtITSMI4->Sumw2(); | |
271 | fHistPtITSMI4->SetMinimum(0); | |
272 | fOutput->Add(fHistPtITSMI4); | |
273 | ||
274 | fHistPtITSMI3 = new TH1F("fHistPtITSMI3","pt distribution of ITSMI3 tracks; p_{t} [GeV/c]; N tracks",80, 0, 40); | |
275 | fHistPtITSMI3->Sumw2(); | |
276 | fHistPtITSMI3->SetMinimum(0); | |
277 | fOutput->Add(fHistPtITSMI3); | |
278 | ||
279 | fHistPtITSMI2 = new TH1F("fHistPtITSMI2","pt distribution of ITSMI2 tracks; p_{t} [GeV/c]; N tracks",80, 0, 40); | |
280 | fHistPtITSMI2->Sumw2(); | |
281 | fHistPtITSMI2->SetMinimum(0); | |
282 | fOutput->Add(fHistPtITSMI2); | |
283 | ||
284 | fHistPtITSMISPD = new TH1F("fHistPtITSMISPD","pt distribution of ITSMISPD tracks; p_{t} [GeV/c]; N tracks",80, 0, 40); | |
285 | fHistPtITSMISPD->Sumw2(); | |
286 | fHistPtITSMISPD->SetMinimum(0); | |
287 | fOutput->Add(fHistPtITSMISPD); | |
288 | ||
289 | ||
290 | ||
291 | // ntuples | |
292 | // | |
293 | fNtupleESDTracks = new TNtuple("fNtupleESDTracks","tracks","pt:eta:phi:d0:z0:sigmad0:sigmaz0:ptMC:pdgMC:d0MC:d0MCv:z0MCv:sigmad0MCv:sigmaz0MCv:ITSflag"); | |
294 | fOutput->Add(fNtupleESDTracks); | |
295 | ||
296 | fNtupleITSAlignExtra = new TNtuple("fNtupleITSAlignExtra","ITS alignment checks: extra clusters","layer:x:y:z:dxy:dz:xloc:zloc:npoints"); | |
297 | fOutput->Add(fNtupleITSAlignExtra); | |
298 | ||
299 | fNtupleITSAlignSPDTracklets = new TNtuple("fNtupleITSAlignSPDTracklets","ITS alignment checks: SPD tracklets wrt SPD vertex","phi:theta:z:dxy:dz"); | |
300 | fOutput->Add(fNtupleITSAlignSPDTracklets); | |
301 | ||
302 | return; | |
303 | } | |
304 | ||
305 | //________________________________________________________________________ | |
306 | void AliAnalysisTaskITSTrackingCheck::Exec(Option_t *) | |
307 | { | |
308 | // Main loop | |
309 | // Called for each event | |
310 | ||
311 | if (!fESD) { | |
312 | Printf("ERROR: fESD not available"); | |
313 | return; | |
314 | } | |
315 | ||
316 | ||
317 | const AliESDVertex *vertexESD = fESD->GetPrimaryVertexTracks(); | |
318 | if(!vertexESD) return; | |
319 | if(!(vertexESD->GetStatus())) return; | |
320 | ||
321 | // *********** MC info *************** | |
322 | TArrayF mcVertex(3); | |
323 | mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.; | |
324 | Float_t dNchdy=-999.; | |
325 | ||
326 | TParticle *part=0; | |
327 | AliESDVertex *vertexMC=0; | |
328 | AliStack *stack=0; | |
329 | if (fReadMC) { | |
330 | AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); | |
331 | if (!eventHandler) { | |
332 | Printf("ERROR: Could not retrieve MC event handler"); | |
333 | return; | |
334 | } | |
335 | ||
336 | AliMCEvent* mcEvent = eventHandler->MCEvent(); | |
337 | if (!mcEvent) { | |
338 | Printf("ERROR: Could not retrieve MC event"); | |
339 | return; | |
340 | } | |
341 | ||
342 | stack = mcEvent->Stack(); | |
343 | if (!stack) { | |
344 | AliDebug(AliLog::kError, "Stack not available"); | |
345 | return; | |
346 | } | |
347 | ||
348 | AliHeader* header = mcEvent->Header(); | |
349 | if (!header) { | |
350 | AliDebug(AliLog::kError, "Header not available"); | |
351 | return; | |
352 | } | |
353 | AliGenEventHeader* genHeader = header->GenEventHeader(); | |
354 | genHeader->PrimaryVertex(mcVertex); | |
355 | ||
356 | ||
357 | Int_t ngenpart = (Int_t)stack->GetNtrack(); | |
358 | //printf("# generated particles = %d\n",ngenpart); | |
359 | dNchdy=0; | |
360 | for(Int_t ip=0; ip<ngenpart; ip++) { | |
361 | part = (TParticle*)stack->Particle(ip); | |
362 | // keep only electrons, muons, pions, kaons and protons | |
363 | Int_t apdg = TMath::Abs(part->GetPdgCode()); | |
364 | if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue; | |
365 | // reject secondaries | |
366 | if(TMath::Sqrt((part->Vx()-mcVertex[0])*(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])*(part->Vy()-mcVertex[1]))>0.0010) continue; | |
367 | // reject incoming protons | |
368 | Double_t energy = part->Energy(); | |
369 | if(energy>900.) continue; | |
370 | Double_t pz = part->Pz(); | |
371 | Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13)); | |
372 | if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1 | |
373 | } | |
374 | //printf("# primary particles = %7.1f\n",dNchdy); | |
375 | } | |
376 | // *********** MC info *************** | |
377 | Double_t mcVtxPos[3]={mcVertex[0],mcVertex[1],mcVertex[2]},mcVtxSigma[3]={0,0,0}; | |
378 | vertexMC = new AliESDVertex(mcVtxPos,mcVtxSigma); | |
379 | ||
380 | // *********** ESD friends *********** | |
381 | fESD->SetESDfriend(fESDfriend); //Attach the friend to the ESD | |
382 | // *********** ESD friends *********** | |
8d63376d | 383 | if(!fESDfriend) printf("no ESD friend\n"); |
384 | ||
385 | // | |
386 | ||
387 | // ********** Trigger ***************** | |
388 | ULong64_t triggerMask; | |
389 | ULong64_t spdFO = (1 << 14); | |
390 | ULong64_t v0left = (1 << 11); | |
391 | ULong64_t v0right = (1 << 12); | |
392 | ||
393 | triggerMask=fESD->GetTriggerMask(); | |
394 | // MB1: SPDFO || V0L || V0R | |
395 | Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right))); | |
396 | //MB2: GFO && V0R | |
397 | //triggerMask & spdFO && ((triggerMask&v0left) || (triggerMask&v0right)) | |
398 | // ************ Trigger ****************** | |
399 | if(!eventTriggered) return; | |
400 | ||
401 | // SPD vertex | |
402 | const AliESDVertex *spdv=fESD->GetPrimaryVertexSPD(); | |
403 | ||
404 | ||
405 | Int_t ntracks = fESD->GetNumberOfTracks(); | |
406 | //printf("Tracks # = %d\n",fESD->GetNumberOfTracks()); | |
407 | ||
408 | fHistNtracks->Fill(ntracks); | |
409 | // Post the data already here | |
410 | PostData(0, fOutput); | |
411 | ||
412 | Int_t idet,status; Float_t xloc,zloc; | |
413 | ||
414 | // loop on tracks | |
415 | for(Int_t itr=0; itr<ntracks; itr++) { | |
416 | AliESDtrack *track = fESD->GetTrack(itr); | |
417 | ||
418 | if(track->GetKinkIndex(0)>0) continue; | |
419 | ||
420 | Int_t trkLabel = TMath::Abs(track->GetLabel()); | |
421 | Int_t nclsITS = track->GetNcls(0); | |
422 | ||
423 | ||
424 | Bool_t itsrefit=kFALSE,tpcin=kFALSE,itsfindableAcc=kFALSE; | |
425 | if ((track->GetStatus() & AliESDtrack::kITSrefit)) itsrefit=kTRUE; | |
426 | if ((track->GetStatus() & AliESDtrack::kTPCin)) tpcin=kTRUE; | |
427 | ||
428 | if(tpcin) { | |
429 | fHistNclsITSMI->Fill(nclsITS); | |
430 | } else { | |
431 | fHistNclsITSSA->Fill(nclsITS); | |
432 | } | |
433 | ||
434 | // TPC track in ITS acceptance | |
435 | if(tpcin && TMath::Abs(track->Eta())<0.9 && track->GetNcls(1)>50) { | |
436 | itsfindableAcc=kTRUE; | |
437 | fHistPtTPC->Fill(track->Pt()); | |
438 | fHistPhiTPC->Fill(track->Phi()); | |
439 | } | |
440 | ||
441 | ||
442 | for(Int_t layer=0; layer<6; layer++) { | |
443 | if(TESTBIT(track->GetITSClusterMap(),layer)) { | |
444 | if(tpcin) { | |
445 | fHistClusterMapITSMI->Fill(layer); | |
446 | } else { | |
447 | fHistClusterMapITSSA->Fill(layer); | |
448 | } | |
449 | } | |
450 | track->GetITSModuleIndexInfo(layer,idet,status,xloc,zloc); | |
451 | if(tpcin) { | |
452 | if(status==1) fHistClusterMapITSMIok->Fill(layer); | |
453 | if(status==2) fHistClusterMapITSMIbad->Fill(layer); | |
454 | if(status==3) fHistClusterMapITSMIskipped->Fill(layer); | |
455 | if(status==4) fHistClusterMapITSMIoutinz->Fill(layer); | |
456 | if(status==5) fHistClusterMapITSMInocls->Fill(layer); | |
457 | if(status==6) fHistClusterMapITSMInorefit->Fill(layer); | |
458 | } else { | |
459 | if(status==1) fHistClusterMapITSSAok->Fill(layer); | |
460 | if(status==2) fHistClusterMapITSSAbad->Fill(layer); | |
461 | if(status==3) fHistClusterMapITSSAskipped->Fill(layer); | |
462 | if(status==4) fHistClusterMapITSSAoutinz->Fill(layer); | |
463 | if(status==5) fHistClusterMapITSSAnocls->Fill(layer); | |
464 | if(status==6) fHistClusterMapITSSAnorefit->Fill(layer); | |
465 | } | |
466 | } | |
467 | ||
468 | if(itsfindableAcc) { | |
469 | if(nclsITS==6) fHistPtITSMI6->Fill(track->Pt()); | |
470 | if(nclsITS==5) fHistPtITSMI5->Fill(track->Pt()); | |
471 | if(nclsITS==4) fHistPtITSMI4->Fill(track->Pt()); | |
472 | if(nclsITS==3) fHistPtITSMI3->Fill(track->Pt()); | |
473 | if(nclsITS==2) fHistPtITSMI2->Fill(track->Pt()); | |
474 | if(track->HasPointOnITSLayer(0) && track->HasPointOnITSLayer(1)) | |
475 | fHistPtITSMISPD->Fill(track->Pt()); | |
476 | } | |
477 | ||
478 | ||
479 | Int_t iITSflag=0; //ITSflag takes the value 0 if the track has no cluster assigned in the SPDs, 1 (2) if one cluster is assigned in SPD1(2), 3 if two clusters are present. Then the same adding 10,20 or 30 for SDD and 100,200 or 300 for SSD | |
480 | ||
481 | if(track->HasPointOnITSLayer(0)) iITSflag+=1; | |
482 | if(track->HasPointOnITSLayer(1)) iITSflag+=2; | |
483 | if(track->HasPointOnITSLayer(2)) iITSflag+=10; | |
484 | if(track->HasPointOnITSLayer(3)) iITSflag+=20; | |
485 | if(track->HasPointOnITSLayer(4)) iITSflag+=100; | |
486 | if(track->HasPointOnITSLayer(5)) iITSflag+=200; | |
487 | ||
488 | if(iITSflag==333 && track->GetNcls(0)<6) | |
489 | printf(" ERROR %d %d\n",track->GetNcls(0),track->GetLabel()); | |
490 | ||
491 | // number of associated ITS clusters | |
492 | iITSflag += 1000*track->GetNcls(0); | |
493 | ||
ab846928 | 494 | // number of associated TPC clusters |
495 | iITSflag += 100000*track->GetNcls(1); | |
496 | ||
8d63376d | 497 | // if MC info and is available |
498 | // write the number of ITS clusters produced by this track | |
499 | Int_t nITSclsMC=0; | |
500 | if(fReadMC && fReadRPLabels) { | |
501 | nITSclsMC = NumberOfITSClustersMC(trkLabel); | |
502 | if(nITSclsMC>=0) iITSflag += 10000*nITSclsMC; | |
503 | } | |
504 | ||
505 | // impact parameter to VertexTracks | |
506 | Double_t d0z0[2],covd0z0[3]; | |
507 | track->PropagateToDCA(vertexESD,fESD->GetMagneticField(),100.,d0z0,covd0z0); | |
508 | if(covd0z0[0]<0. || covd0z0[2]<0.) continue; | |
509 | ||
510 | // if MC info is available: get particle properties | |
511 | Float_t ptMC=-999.,pdgMC=-999.,d0MC=-999.; | |
512 | Double_t d0z0MCv[2]={-999.,-999.},covd0z0MCv[3]={1.,1.,1.}; | |
513 | if(fReadMC) { | |
514 | part = (TParticle*)stack->Particle(trkLabel); | |
515 | ptMC=part->Pt(); | |
516 | pdgMC=part->GetPdgCode(); | |
517 | d0MC=ParticleImpParMC(part,vertexMC,0.1*fESD->GetMagneticField()); | |
518 | track->PropagateToDCA(vertexMC,fESD->GetMagneticField(),100.,d0z0MCv,covd0z0MCv); | |
519 | if(covd0z0MCv[0]<0. || covd0z0MCv[2]<0.) continue; | |
ab846928 | 520 | // flag fake tracks |
521 | if(track->GetLabel()<0) iITSflag *= -1; | |
8d63376d | 522 | } |
523 | ||
ab846928 | 524 | Double_t sigmad0MCv=TMath::Sqrt(covd0z0MCv[0]); |
525 | if(!itsrefit) sigmad0MCv *= -1.; | |
526 | ||
8d63376d | 527 | // fill ntuple with track properties |
528 | if(SelectPt(track->Pt())) { | |
ab846928 | 529 | Float_t fillArray[19]={track->Pt(),track->Eta(),track->Phi(),d0z0[0],d0z0[1],TMath::Sqrt(covd0z0[0]),TMath::Sqrt(covd0z0[2]),ptMC,pdgMC,d0MC,d0z0MCv[0],d0z0MCv[1],sigmad0MCv,TMath::Sqrt(covd0z0MCv[2]),(Float_t)iITSflag}; |
8d63376d | 530 | fNtupleESDTracks->Fill(fillArray); |
531 | } | |
532 | ||
533 | //--------------------------------------------- | |
534 | // AliTrackPoints: alignment checks | |
535 | // | |
ab846928 | 536 | if(!fESDfriend) continue; |
537 | ||
538 | ||
8d63376d | 539 | const AliTrackPointArray *array = track->GetTrackPointArray(); |
540 | AliTrackPoint point; | |
541 | Int_t pointOnLayer[6]={0,0,0,0,0,0}; | |
542 | Int_t indexAssociated[6]={-1,-1,-1,-1,-1,-1},indexExtra=-1; | |
543 | Bool_t extra=kFALSE; | |
544 | Int_t layerId,layerExtra=-1; | |
545 | for(Int_t ipt=0; ipt<array->GetNPoints(); ipt++) { | |
546 | array->GetPoint(point,ipt); | |
547 | Float_t r = TMath::Sqrt(point.GetX()*point.GetX()+point.GetY()*point.GetY()); | |
548 | ||
549 | if(r>3 && r<6) { | |
550 | layerId = 0; | |
551 | } else if(r>6 && r<8) { | |
552 | layerId = 1; | |
553 | } else if(r>8 && r<18) { | |
554 | layerId = 2; | |
555 | } else if(r>18 && r<25) { | |
556 | layerId = 3; | |
557 | } else if(r>25 && r<38) { | |
558 | layerId = 4; | |
559 | } else if(r>38 && r<50) { | |
560 | layerId = 5; | |
561 | } else { | |
562 | layerId=100; | |
563 | } | |
564 | ||
565 | // only ITS points | |
566 | if(layerId>5) continue; | |
567 | ||
568 | if(!point.IsExtra()) { | |
569 | pointOnLayer[layerId]++; | |
570 | indexAssociated[layerId]=ipt; | |
571 | } else { | |
572 | // this is an extra cluster | |
573 | extra=kTRUE; | |
574 | layerExtra=layerId; | |
575 | indexExtra=ipt; | |
576 | } | |
577 | } // end loop on AliTrackPoints | |
578 | ||
579 | TString vtitle = spdv->GetTitle(); | |
580 | if(!vtitle.Contains("3D")) continue; | |
581 | ||
582 | // SPD tracklet | |
583 | if(indexAssociated[0]>=0 && indexAssociated[1]>=0) { | |
584 | AliTrackPoint pointSPD1,pointSPD2; | |
585 | array->GetPoint(pointSPD1,indexAssociated[0]); | |
586 | array->GetPoint(pointSPD2,indexAssociated[1]); | |
587 | Float_t phi=TMath::ATan2(pointSPD2.GetY()-pointSPD1.GetY(),pointSPD2.GetX()-pointSPD1.GetX()); | |
588 | Float_t lambda=TMath::ATan((pointSPD2.GetZ()-pointSPD1.GetZ())/TMath::Sqrt((pointSPD2.GetX()-pointSPD1.GetX())*(pointSPD2.GetX()-pointSPD1.GetX())+(pointSPD2.GetY()-pointSPD1.GetY())*(pointSPD2.GetY()-pointSPD1.GetY()))); | |
589 | Float_t theta=0.5*TMath::Pi()-lambda; | |
590 | TParticle particle(211,0,0,0,0,0,TMath::Cos(phi),TMath::Sin(phi),TMath::Tan(lambda),10.,pointSPD1.GetX(),pointSPD1.GetY(),pointSPD1.GetZ(),0); | |
591 | AliESDtrack tracklet(&particle); | |
592 | Float_t dz[2]; | |
593 | // distance to primary SPD (only if 3D and high multiplicity) | |
594 | if(spdv->GetNContributors()>10) { | |
595 | tracklet.GetDZ(spdv->GetXv(),spdv->GetYv(),spdv->GetZv(),0,dz); | |
596 | fNtupleITSAlignSPDTracklets->Fill(phi,theta,0.5*(pointSPD1.GetZ()+pointSPD2.GetZ()),dz[0],dz[1]); | |
597 | } | |
598 | } | |
599 | ||
600 | // distance to extra | |
601 | if(extra && spdv->GetNContributors()>4) { | |
602 | AliTrackPoint pointExtra,pointAssociated; | |
603 | array->GetPoint(pointAssociated,indexAssociated[layerExtra]); | |
604 | array->GetPoint(pointExtra,indexExtra); | |
605 | Float_t phiExtra = TMath::ATan2(pointExtra.GetY()-spdv->GetYv(),pointExtra.GetX()-spdv->GetXv()); | |
606 | Float_t phiAssociated = TMath::ATan2(pointAssociated.GetY()-spdv->GetYv(),pointAssociated.GetX()-spdv->GetXv()); | |
607 | Float_t rExtra = TMath::Sqrt((pointExtra.GetX()-spdv->GetXv())*(pointExtra.GetX()-spdv->GetXv())+(pointExtra.GetY()-spdv->GetYv())*(pointExtra.GetY()-spdv->GetYv())); | |
608 | Float_t rAssociated = TMath::Sqrt((pointAssociated.GetX()-spdv->GetXv())*(pointAssociated.GetX()-spdv->GetXv())+(pointAssociated.GetY()-spdv->GetYv())*(pointAssociated.GetY()-spdv->GetYv())); | |
609 | Float_t dzExtra[2]; | |
610 | dzExtra[0] = (phiExtra-phiAssociated)*0.5*(rExtra+rAssociated); | |
611 | dzExtra[1] = pointExtra.GetZ()-pointAssociated.GetZ()-(rExtra-rAssociated)*(pointAssociated.GetZ()-spdv->GetZv())/rAssociated; | |
612 | Float_t xlocExtra=-100.,zlocExtra=-100.; | |
613 | fNtupleITSAlignExtra->Fill(layerExtra,pointExtra.GetX(),pointExtra.GetY(),pointExtra.GetZ(),xlocExtra,zlocExtra,dzExtra[0],dzExtra[1],nclsITS); | |
614 | } | |
615 | ||
616 | ||
617 | } // end loop on tracks | |
618 | ||
619 | if(vertexMC) { delete vertexMC; vertexMC=0; } | |
620 | ||
621 | return; | |
622 | } | |
623 | ||
624 | //________________________________________________________________________ | |
625 | void AliAnalysisTaskITSTrackingCheck::Terminate(Option_t *) | |
626 | { | |
627 | // Draw result to the screen | |
628 | // Called once at the end of the query | |
629 | fOutput = dynamic_cast<TList*> (GetOutputData(0)); | |
630 | if (!fOutput) { | |
631 | Printf("ERROR: fOutput not available"); | |
632 | return; | |
633 | } | |
634 | ||
635 | fHistNtracks = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNtracks")); | |
636 | fHistNclsITSMI = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNclsITSMI")); | |
637 | fHistNclsITSSA = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNclsITSSA")); | |
638 | fHistClusterMapITSMI = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSMI")); | |
639 | fHistClusterMapITSMIbad = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSMIbad")); | |
640 | fHistClusterMapITSMIskipped = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSMIskipped")); | |
641 | fHistClusterMapITSMIoutinz = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSMIoutinz")); | |
642 | fHistClusterMapITSMInorefit = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSMInorefit")); | |
643 | fHistClusterMapITSMInocls = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSMInocls")); | |
644 | fHistClusterMapITSSA = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSSA")); | |
645 | fHistClusterMapITSSAbad = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSSAbad")); | |
646 | fHistClusterMapITSSAskipped = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSSAskipped")); | |
647 | fHistClusterMapITSSAoutinz = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSSAoutinz")); | |
648 | fHistClusterMapITSSAnorefit = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSSAnorefit")); | |
649 | fHistClusterMapITSSAnocls = dynamic_cast<TH1F*>(fOutput->FindObject("fHistClusterMapITSSAnocls")); | |
650 | fNtupleESDTracks = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleESDTracks")); | |
651 | fNtupleITSAlignExtra = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleITSAlignExtra")); | |
652 | fNtupleITSAlignSPDTracklets = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleITSAlignSPDTracklets")); | |
653 | ||
654 | return; | |
655 | } | |
656 | //--------------------------------------------------------------------------- | |
657 | Int_t AliAnalysisTaskITSTrackingCheck::NumberOfITSClustersMC(Int_t label) const | |
658 | { | |
659 | // | |
660 | // Return number of ITS clusters produced by MC particle with given label | |
661 | // | |
662 | ||
663 | AliESDInputHandlerRP *esdHRP = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); | |
664 | if(!esdHRP) return -1; | |
665 | TTree *cTree = (TTree*)esdHRP->GetTreeR("ITS"); | |
666 | if(!cTree) return -1; | |
667 | TClonesArray *clusters=0; // new TClonesArray("AliITSRecPoint",10000); | |
668 | cTree->SetBranchAddress("ITSRecPoints",&clusters); | |
669 | if(!clusters) return -1; | |
670 | ||
671 | AliITSRecPoint *c=0; | |
672 | Int_t i,n,icl,lay,ilab; | |
673 | Int_t ncls[6]={0,0,0,0,0,0}; | |
674 | Int_t nclstot=0; | |
675 | ||
676 | for(i=0; i<2198; i++) { | |
677 | cTree->GetEvent(i); | |
678 | n=clusters->GetEntriesFast(); | |
679 | for (icl=0; icl<n; icl++) { | |
680 | c=(AliITSRecPoint*)clusters->UncheckedAt(icl); | |
681 | lay=c->GetLayer(); | |
682 | for(ilab=0;ilab<3;ilab++) { | |
683 | if(c->GetLabel(ilab)==label) ncls[lay]++; | |
684 | } | |
685 | } | |
686 | } | |
687 | for(i=0;i<6;i++) { if(ncls[i]) nclstot++; } | |
688 | ||
689 | return nclstot; | |
690 | } | |
691 | //--------------------------------------------------------------------------- | |
692 | Double_t AliAnalysisTaskITSTrackingCheck::ParticleImpParMC(TParticle *part, | |
693 | AliESDVertex *vert, | |
694 | Double_t bzT) const | |
695 | { | |
696 | // | |
697 | // Return the MC value of the impact parameter | |
698 | // | |
699 | ||
700 | Double_t vx=part->Vx()-vert->GetX(); | |
701 | Double_t vy=part->Vy()-vert->GetY(); | |
702 | ||
703 | Double_t pt=part->Pt(); | |
704 | Double_t px=part->Px(); | |
705 | Double_t py=part->Py(); | |
706 | Double_t charge = (part->GetPdgCode()>0. ? 1. : -1.); | |
707 | if(TMath::Abs(part->GetPdgCode())<100) charge*=-1.; | |
708 | ||
709 | if(px<0.000001) px=0.000001; | |
710 | Double_t rAnd=((10./2.99792458)*pt/bzT)*100.; | |
711 | Double_t center[3],d0; | |
712 | center[0]=vx-(1./charge)*rAnd*(py/pt); | |
713 | center[1]=vy+(1./charge)*rAnd*(px/pt); | |
714 | center[2]=TMath::Sqrt(center[0]*center[0]+center[1]*center[1]); | |
715 | d0 = -center[2]+rAnd; | |
716 | ||
717 | return d0; | |
718 | } | |
719 | //--------------------------------------------------------------------------- | |
720 | Bool_t AliAnalysisTaskITSTrackingCheck::SelectPt(Double_t pt) | |
721 | { | |
722 | // | |
723 | // Keep only tracks in given pt bins | |
724 | // | |
ab846928 | 725 | Double_t ptlower[10]={0.29,0.49,0.75,0.9,1.9,3.5,6.5, 9.,19.,27.}; |
8d63376d | 726 | Double_t ptupper[10]={0.31,0.51,0.85,1.1,2.1,4.5,7.5,11.,21.,33.}; |
727 | ||
728 | for(Int_t i=0; i<10; i++) { | |
ab846928 | 729 | if(pt>ptlower[i] && pt<ptupper[i]) { |
8d63376d | 730 | fCountsPerPtBin[i]++; |
731 | return kTRUE; | |
732 | } | |
733 | } | |
734 | return kFALSE; | |
735 | } | |
736 | ||
737 | ||
738 | ||
739 |