Some of the coding violations corrected
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibMaterial.cxx
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
77 #include "AliKFParticle.h"
78 #include "AliKFVertex.h"
79
80 #include "AliESDTZERO.h"
81
82 ClassImp(AliTPCcalibMaterial)
83
84 AliTPCcalibMaterial::AliTPCcalibMaterial():
85   AliTPCcalibBase("calibMaterial","calibMaterial"),
86   fHisMaterial(0),
87   fHisMaterialRPhi(0)
88 {
89   
90 }
91
92 AliTPCcalibMaterial::AliTPCcalibMaterial(const char * name, const char * title):
93   AliTPCcalibBase(name,title),
94   fHisMaterial(0),
95   fHisMaterialRPhi(0)
96 {
97   //
98   //
99   //
100 }
101
102 AliTPCcalibMaterial::~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
114 Long64_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);
127     //
128     if (!fHisMaterial) fHisMaterial=MakeHisto();
129     fHisMaterial->Add(calib->fHisMaterial);
130     fHisMaterialRPhi->Add(calib->fHisMaterialRPhi);
131   }
132   return 0;
133 }
134
135
136
137 void 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;
147   ProcessPairs(event);
148   return;
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();
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
190 THnSparse *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 }
244
245
246
247
248 void  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   //
284   Float_t dca[2]={0};
285   Float_t cov[3]={0};
286   Float_t dca0[2]={0};
287   Float_t dca1[2]={0};
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
546 Bool_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
584 Bool_t AliTPCcalibMaterial::CheckV0(Int_t index, AliESDEvent *event){
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