]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCcalibMaterial.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
TPC DQM update to have 2d occupancy plot and to disable calib plots from
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibMaterial.cxx
CommitLineData
a144a99d 1
2/**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16
17/*
18 // Load libraries
19
20 gSystem->Load("libANALYSIS");
21 gSystem->Load("libTPCcalib");
22
23
24 .x ~/NimStyle.C
25 gSystem->Load("libANALYSIS");
26 gSystem->Load("libTPCcalib");
27
28 // analyze results
29
30 TFile f("CalibObjectsTrain2.root");
31 AliTPCcalibMaterial *calibMaterial = (AliTPCcalibMaterial *)f->Get("alignMaterial");
32
33
34*/
35
36#include "Riostream.h"
37#include "TChain.h"
38#include "TTree.h"
39#include "TH1F.h"
40#include "TH2F.h"
41#include "TH3F.h"
42#include "THnSparse.h"
43#include "TList.h"
44#include "TMath.h"
45#include "TCanvas.h"
46#include "TFile.h"
47#include "TF1.h"
48#include "TVectorD.h"
49#include "TProfile.h"
50#include "TGraphErrors.h"
51#include "TCanvas.h"
52
53#include "AliTPCclusterMI.h"
54#include "AliTPCseed.h"
55#include "AliESDVertex.h"
56#include "AliESDEvent.h"
57#include "AliESDfriend.h"
58#include "AliESDInputHandler.h"
59#include "AliAnalysisManager.h"
60
61#include "AliTracker.h"
62#include "AliMagF.h"
63#include "AliTPCCalROC.h"
64
65#include "AliLog.h"
66
67#include "AliTPCcalibMaterial.h"
68
69#include "TTreeStream.h"
70#include "AliTPCTracklet.h"
71#include "TTimeStamp.h"
72#include "AliTPCcalibDB.h"
73#include "AliTPCcalibLaser.h"
74#include "AliDCSSensorArray.h"
75#include "AliDCSSensor.h"
76
4af75575 77#include "AliKFParticle.h"
78#include "AliKFVertex.h"
79
80#include "AliESDTZERO.h"
81
a144a99d 82ClassImp(AliTPCcalibMaterial)
83
84AliTPCcalibMaterial::AliTPCcalibMaterial():
85 AliTPCcalibBase("calibMaterial","calibMaterial"),
86 fHisMaterial(0),
87 fHisMaterialRPhi(0)
88{
89
90}
91
92AliTPCcalibMaterial::AliTPCcalibMaterial(const char * name, const char * title):
93 AliTPCcalibBase(name,title),
94 fHisMaterial(0),
95 fHisMaterialRPhi(0)
96{
97 //
98 //
99 //
100}
101
102AliTPCcalibMaterial::~AliTPCcalibMaterial(){
103 //
104 // delete histograms
105 // class is owner of all histograms
106 //
107 if (!fHisMaterial) return;
108 delete fHisMaterial;
109 delete fHisMaterialRPhi;
110 fHisMaterial=0;
111}
112
113
114Long64_t AliTPCcalibMaterial::Merge(TCollection *li) {
115 //
116 // Merge histograms
117 //
118 TIterator* iter = li->MakeIterator();
119 AliTPCcalibMaterial* cal = 0;
120
121 while ((cal = (AliTPCcalibMaterial*)iter->Next())) {
122 if (!cal->InheritsFrom(AliTPCcalibMaterial::Class())) {
123 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
124 return -1;
125 }
126 AliTPCcalibMaterial* calib= (AliTPCcalibMaterial*)(cal);
a144a99d 127 //
128 if (!fHisMaterial) fHisMaterial=MakeHisto();
129 fHisMaterial->Add(calib->fHisMaterial);
130 fHisMaterialRPhi->Add(calib->fHisMaterialRPhi);
131 }
132 return 0;
133}
134
135
136
137void AliTPCcalibMaterial::Process(AliESDEvent *event){
138 //
139 //
140 //
141 const Int_t kMinCl=40;
142 const Float_t kMinRatio=0.7;
143 const Float_t kMaxS=0.05;
144 const Float_t kMinDist=5;
145 const Double_t kStep=1.;
146 if (!event) return;
4af75575 147 ProcessPairs(event);
148 return;
a144a99d 149 // TTreeSRedirector * cstream = GetDebugStreamer();
150 //
151 if (!fHisMaterial){
152 MakeHisto();
153 }
154
155 //
156 // fill histogram of track prolongations
157 Float_t dca[2];
158 Int_t ntracks = event->GetNumberOfTracks();
a144a99d 159 for (Int_t itrack=0; itrack<ntracks; itrack++){
160 AliESDtrack *track=event->GetTrack(itrack);
161 if (!track) continue;
162 if (track->GetTPCNcls()<=kMinCl) continue;
163 if ((1.+track->GetTPCNcls())/(1.+track->GetTPCNclsF())<=kMinRatio) continue;
164 if ((1.+track->GetTPCnclsS())/(1.+track->GetTPCNcls())>kMaxS) continue;
165 if (!track->GetInnerParam()) continue;
166 if (track->GetKinkIndex(0)!=0) continue;
167 //
168 track->GetImpactParameters(dca[0],dca[1]);
169 if (TMath::Abs(dca[0])<kMinDist && TMath::Abs(dca[1])<kMinDist) continue;
170 AliExternalTrackParam param(*(track->GetInnerParam()));
171 if (!AliTracker::PropagateTrackTo(&param,90,0.0005,10,kTRUE)) continue;
172 Double_t x[5]={0,0,0,TMath::Sqrt(TMath::Abs(param.GetP()))*param.GetSign(),TMath::Sqrt(TMath::Abs(track->GetTPCsignal()))};
173 //
174 //
175 for (Float_t radius=90; radius>0; radius-=kStep){
176 if (!AliTracker::PropagateTrackTo(&param,radius,0.0005,kStep*0.5,kTRUE)) break;
177 if (TMath::Abs(param.GetSnp())>0.8) break;
178 param.GetXYZ(x);
179 Double_t weight=1./TMath::Sqrt(1.+param.GetSnp()*param.GetSnp()+param.GetTgl()*param.GetTgl());
180 fHisMaterial->Fill(x,weight);
181 Double_t r = TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
182 Double_t phi = TMath::ATan2(x[1],x[0]);
183 x[0]=r;
184 x[1]=phi;
185 fHisMaterialRPhi->Fill(x,weight);
186 }
187 }
188}
189
190THnSparse *AliTPCcalibMaterial::MakeHisto(){
191 //
192 // Make track prolongation histogram
193 //
194 //
195 // gX gY gz p dEdx
196 Int_t bins[5] = {100, 100, 300, 40, 100};
197 Double_t xmin[5] = {-100, -100, -300, -2, 5};
198 Double_t xmax[5] = {100, 100, 300, 2, 33};
199 TString axisName[5]={
200 "gx",
201 "gy",
202 "gz",
203 "p",
204 "dedx"
205 };
206 TString axisTitle[5]={
207 "x (cm)",
208 "y (cm)",
209 "z (cm)",
210 "p (GeV)",
211 "dedx (a.u)"
212 };
213
214 Int_t binsR[5] = {30, 360, 300, 40, 100};
215 Double_t xminR[5] = { 0, -3.14, -300, -2, 5};
216 Double_t xmaxR[5] = {30, 3.14, 300, 2, 33};
217 TString axisNameR[5]={
218 "r",
219 "rphi",
220 "z",
221 "p",
222 "dedx"
223 };
224 TString axisTitleR[5]={
225 "r (cm)",
226 "rphi (cm)",
227 "z (cm)",
228 "p (GeV)",
229 "dedx (a.u)"
230 };
231
232 THnSparse *sparse = new THnSparseF("his_Material", "His Material", 5, bins, xmin, xmax);
233 THnSparse *sparseR = new THnSparseF("his_MaterialRPhi", "His Material Rphi", 5, binsR, xminR, xmaxR);
234 for (Int_t iaxis=0; iaxis<5; iaxis++){
235 sparse->GetAxis(iaxis)->SetName(axisName[iaxis]);
236 sparse->GetAxis(iaxis)->SetTitle(axisTitle[iaxis]);
237 sparseR->GetAxis(iaxis)->SetName(axisNameR[iaxis]);
238 sparseR->GetAxis(iaxis)->SetTitle(axisTitleR[iaxis]);
239 }
240 fHisMaterial=sparse;
241 fHisMaterialRPhi=sparseR;
242 return sparse;
243}
4af75575 244
245
246
247
248void AliTPCcalibMaterial::ProcessPairs(AliESDEvent *event){
249 //
250 // Process pairs of tracks to get a material budget map
251 //
252 //
253 if (!event) return;
254 if (event) AliKFParticle::SetField(event->GetMagneticField()); // set mean magnetic field for KF particles
255 //
256 // 1. Calculate total dEdx for all TPC tracks
257 //
258 const Int_t kMinCl=70;
259 const Double_t kEpsilon=0.000001;
260 const Float_t kMinRatio=0.7;
261 const Float_t kMinDist=1.5;
262 const Float_t kMinDistChi2=8; //
263 const Float_t kMaxDistZ=280; // max distanceZ
264 const Float_t kMaxDistR=250; // max distanceR
265 const Double_t kMaxChi2 =36; // maximal chi2 to define the vertex
266 const Double_t kMaxDistVertexSec=3; // maximal distance to secondary vertex
267
268 if (!event) return;
269 AliESDVertex *vertexSPD = (AliESDVertex *)event->GetPrimaryVertexSPD();
270 AliESDVertex * spdVertex = (AliESDVertex *)event->GetPrimaryVertexSPD();
271 AliESDVertex * trackVertex = (AliESDVertex *)event->GetPrimaryVertexTracks();
272 AliESDVertex * tpcVertex = (AliESDVertex *)event->GetPrimaryVertexTPC();
273 AliESDTZERO * tzero = (AliESDTZERO *)event->GetESDTZERO() ;
274 //
275 Double_t tpcSignalTotPrim=0;
276 Double_t tpcSignalTotSec=0;
277 Int_t ntracksTPC=0;
278 Int_t nTPCPrim=0;
279 Int_t nTPCSec=0;
280 Int_t ntracks=event->GetNumberOfTracks();
281 if ( ntracks<=2 ) return;
282
283 //
ac2d7e53 284 Float_t dca[2]={0};
285 Float_t cov[3]={0};
286 Float_t dca0[2]={0};
287 Float_t dca1[2]={0};
4af75575 288 //
289 //1. Calculate total dEdx for primary and secondary tracks
290 // and count primaries and secondaries
291 Int_t *rejectTrack = new Int_t[ntracks];
292
293 for (Int_t itrack=0; itrack<ntracks; itrack++){
294 AliESDtrack *track=event->GetTrack(itrack);
295 rejectTrack[itrack]=0;
296 if (!track) continue;
297 if (track->GetTPCNcls()<=kMinCl) continue; // skip short tracks
298 ntracksTPC++;
299 if ((1.+track->GetTPCNcls())/(1.+track->GetTPCNclsF())<=kMinRatio) continue;
300 if (!track->GetInnerParam()) continue; // skip not TPC tracks
301 if (track->GetKinkIndex(0)!=0) {rejectTrack[itrack]+=16;continue;} // skip kinks
302 track->GetImpactParameters(dca[0],dca[1]);
303 if (TMath::Abs(dca[0])>kMaxDistR && TMath::Abs(dca[1])>kMaxDistZ) continue;
304 // remove too dip secondaries
305 if (TMath::Abs(dca[0])<kMinDist && TMath::Abs(dca[1])<kMinDist){
306 tpcSignalTotPrim+=track->GetTPCsignal();
307 nTPCPrim++;
308 }else{
309 tpcSignalTotSec+=track->GetTPCsignal();
310 nTPCSec++;
311 };
312 if (cov[0]>kEpsilon &&TMath::Abs(dca[0])>kEpsilon &&TMath::Sqrt(dca[0]*dca[0]/(TMath::Abs(cov[0])))<kMinDistChi2) rejectTrack[itrack]+=1; // primary
313 if (cov[0]>kEpsilon &&TMath::Abs(dca[0])>kEpsilon &&TMath::Abs(dca[0])<kMinDist) rejectTrack[itrack]+=1; // primary
314 if (track->GetTPCsignal()<40) rejectTrack[itrack]+=16;
315 //
316 if (CheckLooper(itrack, event)) rejectTrack[itrack]+=2; // looper
317 if (CheckV0(itrack,event)) rejectTrack[itrack]+=4;
318 }
319
320
321 //
322 // 2. Find secondary vertices - double loop
323 //
324 for (Int_t itrack0=0; itrack0<ntracks; itrack0++){
325 AliESDtrack *track0=event->GetTrack(itrack0);
326 if (!track0) continue;
327 if (track0->GetTPCNcls()<=kMinCl) continue; // skip short tracks
328 if ((1.+track0->GetTPCNcls())/(1.+track0->GetTPCNclsF())<=kMinRatio) continue;
329 if (!track0->GetInnerParam()) continue; // skip not TPC tracks
330 if (track0->GetKinkIndex(0)>0) continue; // skip kinks
331 if (rejectTrack[itrack0]) continue; // skip
332 track0->GetImpactParameters(dca[0],dca[1]);
333 track0->GetImpactParameters(dca0[0],dca0[1]);
334 if (TMath::Abs(dca[0])>kMaxDistR && TMath::Abs(dca[1])>kMaxDistZ) continue;
335 // remove too dip secondaries
336
337 AliKFParticle part0(*track0,211); //assuming pion mass
338 if (track0->Charge()*part0.Q()<0) part0.Q()*=-1; // change sign if opposite
339 //
340 for (Int_t itrack1=itrack0+1; itrack1<ntracks; itrack1++){
341 AliESDtrack *track1=event->GetTrack(itrack1);
342 if (!track1) continue;
343 if (rejectTrack[itrack1]) continue; // skip
344 if (track1->GetTPCNcls()<=kMinCl) continue; // skip short tracks
345 if ((1.+track1->GetTPCNcls())/(1.+track1->GetTPCNclsF())<=kMinRatio) continue;
346 if (!track1->GetInnerParam()) continue; // skip not TPC tracks
347 if (track1->GetKinkIndex(0)!=0) continue; // skip kinks
348 track1->GetImpactParameters(dca1[0],dca1[1]);
349 track1->GetImpactParameters(dca[0],dca[1]);
350 if (TMath::Abs(dca[0])<kMinDist && TMath::Abs(dca[1])<kMinDist) continue;
351 if (TMath::Abs(dca[0])>kMaxDistR && TMath::Abs(dca[1])>kMaxDistZ) continue;
352 AliKFParticle part1(*track1,211); // assuming pion mass
353 if (track1->Charge()*part1.Q()<0) part1.Q()*=-1; // change sign if opposite
354
355 //
356 //
357 AliKFVertex vertex;
358 vertex+=part0;
359 vertex+=part1;
360 if ((vertex.GetChi2()/vertex.GetNDF())> kMaxChi2) continue;
361 if (TMath::Abs(vertex.GetX())>kMaxDistR) continue;
362 if (TMath::Abs(vertex.GetY())>kMaxDistR) continue;
363 if (TMath::Abs(vertex.GetZ())>kMaxDistZ) continue;
364 Double_t errX2=vertex.GetErrX();
365 Double_t errY2=vertex.GetErrY();
366 Double_t errZ2=vertex.GetErrZ();
367 //
368 Double_t err3D=TMath::Sqrt(errX2*errX2+errY2*errY2+errZ2*errZ2/10.);
369 Double_t err2D=TMath::Sqrt(errX2*errX2+errY2*errY2);
370 if (err3D>kMaxDistVertexSec) continue;
371 if (err3D*TMath::Sqrt(vertex.GetChi2()+0.00001)>kMaxDistVertexSec) continue;
372
373 Double_t dvertex=0;
374 dvertex += (vertexSPD->GetX()-vertex.GetX())*(vertexSPD->GetX()-vertex.GetX());
375 dvertex += (vertexSPD->GetY()-vertex.GetY())*(vertexSPD->GetY()-vertex.GetY());
376 dvertex += (vertexSPD->GetZ()-vertex.GetZ())*(vertexSPD->GetZ()-vertex.GetZ());
377 dvertex=TMath::Sqrt(dvertex+0.00000001);
378 if (err3D>0.2*dvertex) continue;
379 if (err3D*TMath::Sqrt(vertex.GetChi2()+0.000001)>0.1*dvertex) continue;
380 Double_t radius = TMath::Sqrt((vertex.GetX()*vertex.GetX()+vertex.GetY()*vertex.GetY()));
381 //
382 AliKFVertex vertex2;
383 vertex2+=part0;
384 vertex2+=part1;
385 //
386
387 for (Int_t itrack2=0; itrack2<ntracks; itrack2++){
388 if (itrack2==itrack0) continue;
389 if (itrack2==itrack1) continue;
390 if (rejectTrack[itrack2]) continue; // skip
391 AliESDtrack *track2=event->GetTrack(itrack2);
392 if (!track2) continue;
393 if (track2->GetTPCNcls()<=kMinCl) continue; // skip short tracks
394 if ((1.+track2->GetTPCNcls())/(1.+track2->GetTPCNclsF())<=kMinRatio) continue;
395 if (!track2->GetInnerParam()) continue; // skip not TPC tracks
396 if (track2->GetKinkIndex(0)>0) continue; // skip kinks
397 track2->GetImpactParameters(dca[0],dca[1]);
398 if (TMath::Abs(dca[0])<kMinDist && TMath::Abs(dca[1])<kMinDist) continue;
399 if (TMath::Abs(dca[0])>kMaxDistR && TMath::Abs(dca[1])>kMaxDistZ) continue;
400 if (TMath::Abs(track2->GetD(vertex.GetX(), vertex.GetY(),event->GetMagneticField()))>kMaxDistVertexSec) continue;
401 Double_t vtxx[3]={vertex2.GetX(),vertex2.GetY(),vertex2.GetZ()};
402 Double_t svtxx[3]={vertex.GetErrX(),vertex.GetErrY(),vertex.GetErrZ()};
403 AliESDVertex vtx(vtxx,svtxx);
404 AliExternalTrackParam param=*track2;
405 Double_t delta[2]={0,0};
406 if (!param.PropagateToDCA(&vtx,event->GetMagneticField(),kMaxDistVertexSec,delta)) continue;
407 if (TMath::Abs(delta[0])>kMaxDistVertexSec) continue;
408 if (TMath::Abs(delta[1])>kMaxDistVertexSec) continue;
409 if (TMath::Abs(delta[0])>6.*TMath::Sqrt(param.GetSigmaY2()+vertex2.GetErrY()*vertex2.GetErrY())+0.1) continue;
410 if (TMath::Abs(delta[1])>6.*TMath::Sqrt(param.GetSigmaZ2()+vertex2.GetErrZ()*vertex2.GetErrZ())+0.5) continue;
411 //
412 AliKFParticle part2(param,211); // assuming pion mass
413 if (track2->Charge()*part2.Q()<0) part2.Q()*=-1; // change sign if opposite
414 vertex2+=part2;
415 rejectTrack[itrack0]+=10; // do noit reuse the track
416 rejectTrack[itrack1]+=10; // do not reuse the track
417 rejectTrack[itrack2]+=10;
418 }
419
420
421 TTreeSRedirector *pcstream = GetDebugStreamer();
422 if (pcstream){
423 //
424 //
425 //
426 Float_t dedx0= track0->GetTPCsignal();
427 Float_t dedx1= track1->GetTPCsignal();
428 AliExternalTrackParam * p0= (AliExternalTrackParam *)track0->GetInnerParam();
429 AliExternalTrackParam * p1= (AliExternalTrackParam *)track1->GetInnerParam();
430 Double_t errX=vertex2.GetErrX();
431 Double_t errY=vertex2.GetErrY();
432 Double_t errZ=vertex2.GetErrZ();
433 Double_t vx = vertex2.GetX();
434 Double_t vy = vertex2.GetY();
435 Double_t vz = vertex2.GetZ();
436 (*pcstream)<<"mapTPC"<<
437 "run="<<fRun<< // run
438 "event="<<fEvent<< // event number
439 "time="<<fTime<< // timeStamp
440 "trigger="<<fTrigger<< // trigger
441 "mag="<<fMagF<< // magnetic field
442 //
443 "spdV.="<<spdVertex<< // spd vertex
444 "trackV.="<<trackVertex<< // track vertex
445 "tpcV.="<<tpcVertex<< // track vertex
446 "tzero.="<<tzero<< // tzero info
447 //
448 "ntracks="<<ntracks<<
449 "ntracksTPC="<<ntracksTPC<<
450 "nPrim="<<nTPCPrim<< // number of primaries
451 "nSec="<<nTPCSec<< // number of secondaries
452 "sigPrim="<<tpcSignalTotPrim<< // total dEdx in primaries
453 "sigSec="<<tpcSignalTotSec<< // total dEdx in secondaries
454 "dedx0="<<dedx0<< // dedx part 0
455 "dedx1="<<dedx1<< // dedx part 1
456 "p0.="<<p0<< // part 0
457 "p1.="<<p1<< //part 1
458 "v.="<<&vertex<< // KF vertex
459 "v2.="<<&vertex2<< // KF vertex all tracks
460 "z0="<<dca0[1]<<
461 "z1="<<dca1[1]<<
462 "rphi0="<<dca0[0]<<
463 "rphi1="<<dca1[0]<<
464 "dvertex="<<dvertex<<
465 "radius="<<radius<<
466 "vx="<<vx<<
467 "vy="<<vy<<
468 "vz="<<vz<<
469 "errX="<<errX<<
470 "errY="<<errY<<
471 "errZ="<<errZ<<
472 "err2D="<<err2D<<
473 "err3D="<<err3D<<
474 "\n";
475 //
476 if (vertex2.GetNDF()>2){
477 (*pcstream)<<"mapVertex"<<
478 "run="<<fRun<< // run
479 "event="<<fEvent<< // event number
480 "time="<<fTime<< // timeStamp
481 "trigger="<<fTrigger<< // trigger
482 "mag="<<fMagF<< // magnetic field
483 //
484 "spdV.="<<spdVertex<< // spd vertex
485 "trackV.="<<trackVertex<< // track vertex
486 "tpcV.="<<tpcVertex<< // track vertex
487 "tzero.="<<tzero<< // tzero info
488 //
489 //
490 "ntracks="<<ntracks<<
491 "ntracksTPC="<<ntracksTPC<<
492 "nPrim="<<nTPCPrim<< // number of primaries
493 "nSec="<<nTPCSec<< // number of secondaries
494 "sigPrim="<<tpcSignalTotPrim<< // total dEdx in primaries
495 "sigSec="<<tpcSignalTotSec<< // total dEdx in secondaries
496 "dedx0="<<dedx0<< // dedx part 0
497 "dedx1="<<dedx1<< // dedx part 1
498 "p0.="<<p0<< // part 0
499 "p1.="<<p1<< //part 1
500 "v.="<<&vertex<< // KF vertex
501 "v2.="<<&vertex2<< // KF vertex all tracks
502 "z0="<<dca0[1]<<
503 "z1="<<dca1[1]<<
504 "rphi0="<<dca0[0]<<
505 "rphi1="<<dca1[0]<<
506 "radius="<<radius<<
507 "vx="<<vx<<
508 "vy="<<vy<<
509 "vz="<<vz<<
510 "errX="<<errX<<
511 "errY="<<errY<<
512 "errZ="<<errZ<<
513 "err2D="<<err2D<<
514 "err3D="<<err3D<<
515 "dvertex="<<dvertex<<
516 "\n";
517 }
518 }
519 }
520 }
521 TTreeSRedirector *pcstream = GetDebugStreamer();
522 if (pcstream){
523 (*pcstream)<<"statTPC"<<
524 "run="<<fRun<< // run
525 "time="<<fTime<< // timeStamp
526 "trigger="<<fTrigger<< // trigger
527 "mag="<<fMagF<< // magnetic field
528 "ntracks="<<ntracks<<
529 "ntracksTPC="<<ntracksTPC<<
530 //
531 "nPrim="<<nTPCPrim<< // number of primaries
532 "nSec="<<nTPCSec<< // number of secondaries
533 "sigPrim="<<tpcSignalTotPrim<< // total dEdx in primaries
534 "sigSec="<<tpcSignalTotSec<< // total dEdx in secondaries
535 //
536 "spdV.="<<spdVertex<< // spd vertex
537 "trackV.="<<trackVertex<< // track vertex
538 "tpcV.="<<tpcVertex<< // track vertex
539 "tzero.="<<tzero<< // tzero info
540 "\n";
541 }
542 delete [] rejectTrack;
543}
544
545
546Bool_t AliTPCcalibMaterial::CheckLooper(Int_t index, AliESDEvent *event){
547 //
548 // check if given track is looper candidate
549 // if looper return kTRUE
550 //
551 Int_t ntracks=event->GetNumberOfTracks();
552 Int_t index1=-1;
553 const Double_t ktglCut=0.03;
554 const Double_t kalphaCut=0.4;
555 static Int_t counter=0;
556 //
557 AliESDtrack * track0 = event->GetTrack(index);
558 AliESDtrack * track1P = 0;
559 for (Int_t itrack1=0; itrack1<ntracks; itrack1++){
560 if (itrack1==index) continue;
561 AliESDtrack *track1=event->GetTrack(itrack1);
562 if (!track1) continue;
563 if (TMath::Abs(TMath::Abs(track1->GetTgl())-TMath::Abs(track0->GetTgl()))>ktglCut) continue;
564 if (TMath::Abs(TMath::Abs(track1->GetAlpha())-TMath::Abs(track0->GetAlpha()))>kalphaCut) continue;
565 index1=index;
566 track1P=track1;
567 }
568 if (index1>=0){
569 TTreeSRedirector *pcstream = GetDebugStreamer();
570 if (pcstream &&counter<100000){
571 counter++;
572 AliExternalTrackParam p0(*track0);
573 AliExternalTrackParam p1(*track1P);
574 (*pcstream)<<"checkLooper"<<
575 "p0.="<<&p0<<
576 "p1.="<<&p1<<
577 "\n";
578 }
579 return kTRUE;
580 }
581 return kFALSE;
582}
583
c57f81bf 584Bool_t AliTPCcalibMaterial::CheckV0(Int_t index, AliESDEvent *event){
4af75575 585 //
586 // check if given track is V0 candidata
587 // if looper return kTRUE
588 //
589 return kFALSE;
590 Int_t ntracks=event->GetNumberOfTracks();
591 Int_t index1=-1;
592 const Double_t kSigmaMass=0.001;
593 const Int_t kChi2Cut=10;
594 static Int_t counter=0;
595 //
596 AliESDtrack * track0 = event->GetTrack(index);
597 AliExternalTrackParam pL(*track0);
598 AliKFParticle part0El(*track0, 11); //assuming mass e
599 AliKFParticle part0Pi(*track0, 211); //assuming mass pi0
600 AliKFParticle part0P(*track0, 2212); //assuming mass proton
601 if (track0->Charge()*part0El.Q()<0) {
602 part0El.Q()*=-1; // change sign if opposite
603 part0Pi.Q()*=-1; // change sign if opposite
604 part0P.Q()*=-1; // change sign if opposite
605 }
606 Bool_t isGamma=0;
607 Bool_t isK0=0;
608
609 for (Int_t itrack1=0; itrack1<ntracks; itrack1++){
610 if (itrack1==index) continue;
611 AliESDtrack *track1=event->GetTrack(itrack1);
612 if (!track1) continue;
613 if (track1->Charge()*track0->Charge()>0) continue;
614 AliKFParticle part1El(*track1, 11); //assuming mass e
615 AliKFParticle part1Pi(*track1, 211); //assuming mass e
616 AliKFParticle part1P(*track1, 2212); //assuming mass e
617 if (track1->Charge()*part1El.Q()<0) {
618 part1El.Q()*=-1; // change sign if opposite
619 part1Pi.Q()*=-1; // change sign if opposite
620 part1P.Q()*=-1; // change sign if opposite
621 }
622 //
623 AliKFVertex vertexG; // gamma conversion candidate
624 vertexG+=part0El;
625 vertexG+=part1El;
626 AliKFVertex vertexGC; // gamma conversion candidate
627 vertexGC+=part0El;
628 vertexGC+=part1El;
629 vertexGC.SetMassConstraint(0,kSigmaMass);
630 AliKFVertex vertexK0; // gamma conversion candidate
631 vertexK0+=part0Pi;
632 vertexK0+=part1Pi;
633 AliKFVertex vertexK0C; // gamma conversion candidate
634 vertexK0C+=part0Pi;
635 vertexK0C+=part1Pi;
636 vertexK0C.SetMassConstraint(0.497614,kSigmaMass);
637 if (vertexGC.GetChi2()<kChi2Cut && vertexG.GetMass()<0.06) isGamma=kTRUE;
638 if (vertexK0C.GetChi2()<kChi2Cut&&TMath::Abs(vertexK0.GetMass()-0.5)<0.06) isK0=kTRUE;
639 if (isGamma||isK0) {
640 index1=index;
641 TTreeSRedirector *pcstream = GetDebugStreamer();
642 if (pcstream&&counter<2000){
643 counter++;
644 AliExternalTrackParam p0(*track0);
645 AliExternalTrackParam p1(*track1);
646 (*pcstream)<<"checkV0"<<
647 "p0.="<<&p0<< //particle 0
648 "p1.="<<&p1<< //particle 1
649 "isGamma="<<isGamma<< // is gamma candidate
650 "isK0s="<<isK0<< // is k0 candidate
651 "vG.="<<&vertexG<< // is gamma candidate
652 "vGC.="<<&vertexGC<< // is gamma candidate
653 "vK.="<<&vertexK0<< // is K0s candidate
654 "vKC.="<<&vertexK0C<< // is K0s candidate
655 "\n";
656 }
657 break;
658 }
659 }
660 if (index1>0) return kTRUE;
661 return kFALSE;
662}
663
664/*
665 //AliXRDPROOFtoolkit::FilterList("mater.list","* mapTPC",1);
666 AliXRDPROOFtoolkit toolkit;
667 TChain *chain = toolkit.MakeChainRandom("mater.list.Good","mapVertex",0,4000);
668 TChain *chainTPC = toolkit.MakeChainRandom("mater.list.Good","mapTPC",0,50000);
669
670
671 TCut cutErr="sqrt(v.fChi2)*err3D<1.0";
672 TCut cutOccu="sqrt(v.fChi2*(errX^2+errY^2))/min(sqrt(v.fP[0]^2+v.fP[1]^2+v.fP[2]^2/20),20)<0.2&&sqrt((errX^2+errY^2))/min(sqrt(v.fP[0]^2+v.fP[1]^2+v.fP[2]^2/20),20)<0.3&&sqrt(v.fChi2)*err3D<1.5";
673 //
674 chainTPC->Draw(">>listTPC",cutOccu,"entryList");
675 TEntryList *elistTPC = (TEntryList*)gDirectory->Get("listTPC");
676 chainTPC->SetEntryList(elistTPC);
677
678 //
679 //
680
681 chainTPC.Draw("v.fP[1]:v.fP[0]","sqrt(v.fP[0]^2+v.fP[1]^2)<100","",10000000);
682
683 TCut cutITS="abs(v.fP[2])-10<sqrt(v.fP[0]^2+v.fP[1]^2)";
684 chainTPC.Draw("v.fP[1]:v.fP[0]",cutITS+"(sqrt(v.fP[0]^2+v.fP[1]^2)<50)&&err2D<1&&err3D/dvertex<0.05","",500000);
685
686
687 chainTPC.Draw("v.fP[1]:v.fP[0]>>his(400,-100,100,400,-100,100)","err2D<1&&err2D/dvertex<0.2","colz",50000);
688
689
690 chainTPC.Draw("radius:vz>>his(400,-100,100,400,0,100)","err2D<1&&err2D/dvertex<0.05","colz",5000);
691
692 TTree * tree = chain->CloneTree();
693 TFile f("material.root","recreate");
694 tree->Write();
695 f.Close();
696
697
698
699 TCut cutChi2="v.fChi2<2&&sqrt(errX^2+errY^2)<2.";
700 TFile f("material.root")
701 TTree * tree = (TTree*)f.Get("mapVertex");
702
703 tree->Draw("v.fChi2",cutChi2,"",1000)
704
705
706
707
708*/
709
710