The present commit corresponds to an important change in the way the
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibCosmic.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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     Comments to be written here: 
18     1. What do we calibrate.
19     2. How to interpret results
20     3. Simple example
21     4. Analysis using debug streamers.
22
23
24
25     3.Simple example
26     // To make cosmic scan the user interaction neccessary
27     //
28     .x ~/UliStyle.C
29     gSystem->Load("libANALYSIS");
30     gSystem->Load("libTPCcalib");
31     TFile fcalib("CalibObjects.root");
32     TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
33     AliTPCcalibCosmic * cosmic = ( AliTPCcalibCosmic *)array->FindObject("cosmicTPC");
34     
35
36
37 */
38
39
40
41 #include "Riostream.h"
42 #include "TChain.h"
43 #include "TTree.h"
44 #include "TH1F.h"
45 #include "TH2F.h"
46 #include "TList.h"
47 #include "TMath.h"
48 #include "TCanvas.h"
49 #include "TFile.h"
50 #include "TF1.h"
51
52 #include "AliTPCclusterMI.h"
53 #include "AliTPCseed.h"
54 #include "AliESDVertex.h"
55 #include "AliESDEvent.h"
56 #include "AliESDfriend.h"
57 #include "AliESDInputHandler.h"
58 #include "AliAnalysisManager.h"
59
60 #include "AliTracker.h"
61 #include "AliMagF.h"
62 #include "AliTPCCalROC.h"
63
64 #include "AliLog.h"
65
66 #include "AliTPCcalibCosmic.h"
67 #include "AliExternalComparison.h"
68 #include "TTreeStream.h"
69 #include "AliTPCTracklet.h"
70
71 ClassImp(AliTPCcalibCosmic)
72
73
74 AliTPCcalibCosmic::AliTPCcalibCosmic() 
75   :AliTPCcalibBase(),
76    fGainMap(0),
77    fHistNTracks(0),
78    fClusters(0),
79    fModules(0),
80    fHistPt(0),
81    fDeDx(0),
82    fDeDxMIP(0),
83    fMIPvalue(1), 
84    fCutMaxD(5),        // maximal distance in rfi ditection
85    fCutMaxDz(40),      // maximal distance in z ditection
86    fCutTheta(0.03),    // maximal distan theta
87    fCutMinDir(-0.99)   // direction vector products
88 {  
89   AliInfo("Default Constructor");    
90 }
91
92
93 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title) 
94   :AliTPCcalibBase(),
95    fGainMap(0),
96    fHistNTracks(0),
97    fClusters(0),
98    fModules(0),
99    fHistPt(0),
100    fDeDx(0),
101    fDeDxMIP(0),
102    fMIPvalue(1),
103    fCutMaxD(5),        // maximal distance in rfi ditection 
104    fCutMaxDz(40),      // maximal distance in z ditection
105    fCutTheta(0.03),    // maximal distan theta
106    fCutMinDir(-0.99)   // direction vector products
107 {  
108   SetName(name);
109   SetTitle(title);
110
111   fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
112   fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",160,0,160);
113   fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-650,650,600,-700,700);
114   fHistPt = new TH1F("Pt","Pt distribution; p_{T} (GeV); counts",2000,0,50);
115   fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,100.,500,2.,1000);
116   BinLogX(fDeDx);
117   fDeDxMIP =  new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000);
118
119   AliInfo("Non Default Constructor");  
120   //
121 }
122
123 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
124   //
125   //
126   //
127 }
128
129
130
131
132
133 void AliTPCcalibCosmic::Process(AliESDEvent *event) {
134   //
135   //
136   //
137   if (!event) {
138     Printf("ERROR: ESD not available");
139     return;
140   }  
141   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
142   if (!ESDfriend) {
143    Printf("ERROR: ESDfriend not available");
144    return;
145   }
146    
147
148   FindPairs(event); // nearly everything takes place in find pairs...
149
150   if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
151   Int_t ntracks=event->GetNumberOfTracks(); 
152   fHistNTracks->Fill(ntracks);
153   if (ntracks==0) return;
154
155 }
156
157
158
159 void AliTPCcalibCosmic::Analyze() {
160
161   fMIPvalue = CalculateMIPvalue(fDeDxMIP);
162
163   return;
164
165 }
166
167
168
169 void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
170   //
171   // Find cosmic pairs
172   // 
173   // Track0 is choosen in upper TPC part
174   // Track1 is choosen in lower TPC part
175   //
176   if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
177   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
178   Int_t ntracks=event->GetNumberOfTracks(); 
179   TObjArray  tpcSeeds(ntracks);
180   if (ntracks==0) return;
181   Double_t vtxx[3]={0,0,0};
182   Double_t svtxx[3]={0.000001,0.000001,100.};
183   AliESDVertex vtx(vtxx,svtxx);
184   //
185   //track loop
186   //
187   for (Int_t i=0;i<ntracks;++i) {
188    AliESDtrack *track = event->GetTrack(i);
189    fClusters->Fill(track->GetTPCNcls()); 
190   
191    const AliExternalTrackParam * trackIn = track->GetInnerParam();
192    const AliExternalTrackParam * trackOut = track->GetOuterParam();
193    if (!trackIn) continue;
194    if (!trackOut) continue;
195    if (ntracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue;  // filter laser 
196
197
198    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
199    TObject *calibObject;
200    AliTPCseed *seed = 0;
201    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
202      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
203    }
204    if (seed) tpcSeeds.AddAt(seed,i);
205
206    Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
207    if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
208      fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
209      //
210      if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap));
211      //
212      if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159,fGainMap)>300) {
213        TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
214        if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
215        if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
216      }
217
218    }
219
220   }
221
222   if (ntracks<2) return;
223   //
224   // Find pairs
225   //
226   for (Int_t i=0;i<ntracks;++i) {
227     AliESDtrack *track0 = event->GetTrack(i);     
228     // track0 - choosen upper part
229     if (!track0) continue;
230     if (!track0->GetOuterParam()) continue;
231     if (track0->GetOuterParam()->GetAlpha()<0) continue;
232     Double_t dir0[3];
233     track0->GetDirection(dir0);    
234     for (Int_t j=0;j<ntracks;++j) {
235       if (i==j) continue;
236       AliESDtrack *track1 = event->GetTrack(j);   
237       //track 1 lower part
238       if (!track1) continue;
239       if (!track1->GetOuterParam()) continue;
240       if (track1->GetOuterParam()->GetAlpha()>0) continue;
241       //
242       Double_t dir1[3];
243       track1->GetDirection(dir1);
244       
245       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
246       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
247       if (! seed0) continue;
248       if (! seed1) continue;
249       Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
250       Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
251       //
252       Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
253       Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
254       //
255       Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
256       Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
257       //
258       Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
259       Float_t d0  = track0->GetLinearD(0,0);
260       Float_t d1  = track1->GetLinearD(0,0);
261       //
262       // conservative cuts - convergence to be guarantied
263       // applying before track propagation
264       if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
265       if (dir>fCutMinDir) continue;               // direction vector product
266       Float_t bz = AliTracker::GetBz();
267       Float_t dvertex0[2];   //distance to 0,0
268       Float_t dvertex1[2];   //distance to 0,0 
269       track0->GetDZ(0,0,0,bz,dvertex0);
270       track1->GetDZ(0,0,0,bz,dvertex1);
271       if (TMath::Abs(dvertex0[1])>250) continue;
272       if (TMath::Abs(dvertex1[1])>250) continue;
273       //
274       //
275       //
276       Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
277       AliExternalTrackParam param0(*track0);
278       AliExternalTrackParam param1(*track1);
279       //
280       // Propagate using Magnetic field and correct fo material budget
281       //
282       AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
283       AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
284       //
285       // Propagate rest to the 0,0 DCA - z should be ignored
286       //
287       Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
288       Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
289       //      
290       param0.GetDZ(0,0,0,bz,dvertex0);
291       param1.GetDZ(0,0,0,bz,dvertex1);
292       if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
293       //
294       Double_t xyz0[3];//,pxyz0[3];
295       Double_t xyz1[3];//,pxyz1[3];
296       param0.GetXYZ(xyz0);
297       param1.GetXYZ(xyz1);
298       Bool_t isPair = IsPair(&param0,&param1);
299       //
300       if (isPair) FillAcordeHist(track0);
301       //
302       if (fStreamLevel>0){
303         TTreeSRedirector * cstream =  GetDebugStreamer();
304         //printf("My stream=%p\n",(void*)cstream);
305         AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
306         AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
307         AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
308         AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
309         Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0;
310         Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0;
311         Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]);
312         Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]);
313         if (cstream) {
314           (*cstream) << "Track0" <<
315             "run="<<fRun<<              //  run number
316             "event="<<fEvent<<          //  event number
317             "time="<<fTime<<            //  time stamp of event
318             "trigger="<<fTrigger<<      //  trigger
319             "mag="<<fMagF<<             //  magnetic field
320             "dir="<<dir<<               //  direction
321             "OK="<<isPair<<             //  will be accepted
322             "b0="<<b0<<                 //  propagate status
323             "b1="<<b1<<                 //  propagate status
324             "crossI="<<isCrossI<<       //  cross inner
325             "crossO="<<isCrossO<<       //  cross outer
326             //
327             "Orig0.=" << track0 <<      //  original track  0
328             "Orig1.=" << track1 <<      //  original track  1
329             "Tr0.="<<&param0<<          //  track propagated to the DCA 0,0
330             "Tr1.="<<&param1<<          //  track propagated to the DCA 0,0        
331             "Ip0.="<<ip0<<              //  inner param - upper
332             "Ip1.="<<ip1<<              //  inner param - lower
333             "Op0.="<<op0<<              //  outer param - upper
334             "Op1.="<<op1<<              //  outer param - lower
335             //
336             "v00="<<dvertex0[0]<<       //  distance using kalman
337             "v01="<<dvertex0[1]<<       // 
338             "v10="<<dvertex1[0]<<       //
339             "v11="<<dvertex1[1]<<       // 
340             "d0="<<d0<<                 //  linear distance to 0,0
341             "d1="<<d1<<                 //  linear distance to 0,0
342             //
343             //
344             //
345             "x00="<<xyz0[0]<<           // global position close to vertex
346             "x01="<<xyz0[1]<<
347             "x02="<<xyz0[2]<<
348             //
349             "x10="<<xyz1[0]<<           // global position close to vertex
350             "x11="<<xyz1[1]<<
351             "x12="<<xyz1[2]<<
352             //
353             "alpha0="<<alpha0<<
354             "alpha1="<<alpha1<<
355             "dir00="<<dir0[0]<<           // direction upper
356             "dir01="<<dir0[1]<<
357             "dir02="<<dir0[2]<<
358             //
359             "dir10="<<dir1[0]<<           // direction lower
360             "dir11="<<dir1[1]<<
361             "dir12="<<dir1[2]<<
362             //
363             //
364             "Seed0.=" << seed0 <<       //  original seed 0
365             "Seed1.=" << seed1 <<       //  original seed 1
366             //
367             "dedx0="<<dedx0<<           //  dedx0 - all
368             "dedx1="<<dedx1<<           //  dedx1 - all
369             //
370             "dedx0I="<<dedx0I<<         //  dedx0 - inner ROC
371             "dedx1I="<<dedx1I<<         //  dedx1 - inner ROC
372             //
373             "dedx0O="<<dedx0O<<         //  dedx0 - outer ROC
374             "dedx1O="<<dedx1O<<         //  dedx1 - outer ROC
375             "\n";
376         }
377       }      
378     }
379   }  
380 }    
381
382
383
384
385 void  AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
386
387   // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
388   if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
389     
390   const Double_t AcordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
391   const Double_t roof = 210.5;       // distance from x =0 to end of magnet roof
392
393   Double_t r[3];
394   upperTrack->GetXYZ(r);
395   Double_t d[3];
396   upperTrack->GetDirection(d);
397   Double_t x,z;
398   z = r[2] + (d[2]/d[1])*(AcordePlane - r[1]);
399   x = r[0] + (d[0]/d[1])*(AcordePlane - r[1]);
400   
401   if (x > roof) {
402     x = r[0] + (d[0]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
403     z = r[2] + (d[2]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
404   }
405   if (x < -roof) {
406     x = r[0] + (d[0]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);       
407     z = r[2] + (d[2]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
408   } 
409
410   fModules->Fill(z, x);
411  
412 }
413
414
415
416 Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
417
418   TIterator* iter = li->MakeIterator();
419   AliTPCcalibCosmic* cal = 0;
420
421   while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
422     if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
423       //Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
424       return -1;
425     }
426     
427     fHistNTracks->Add(cal->GetHistNTracks());
428     fClusters->Add(cal-> GetHistClusters());
429     fModules->Add(cal->GetHistAcorde());
430     fHistPt->Add(cal->GetHistPt());
431     fDeDx->Add(cal->GetHistDeDx());
432     fDeDxMIP->Add(cal->GetHistMIP());
433   
434   }
435   
436   return 0;
437   
438 }
439
440
441 Bool_t  AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
442   //
443   //
444   /*
445   // 0. Same direction - OPOSITE  - cutDir +cutT    
446   TCut cutDir("cutDir","dir<-0.99")
447   // 1. 
448   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
449   //
450   // 2. The same rphi 
451   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
452   //
453   //
454   //
455   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");  
456   // 1/Pt diff cut
457   */
458   const Double_t *p0 = tr0->GetParameter();
459   const Double_t *p1 = tr1->GetParameter();
460   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
461   if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
462   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
463   
464   Double_t d0[3], d1[3];
465   tr0->GetDirection(d0);    
466   tr1->GetDirection(d1);       
467   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
468   //
469   return kTRUE;  
470 }
471  
472
473
474 Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
475
476   TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
477   funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
478                                 hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
479   hist->Fit(funcDoubleGaus);
480   Double_t MIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
481
482   delete funcDoubleGaus;
483
484   return MIPvalue;
485
486 }
487
488
489
490
491 void AliTPCcalibCosmic::CalculateBetheParams(TH2F */*hist*/, Double_t * /*initialParam*/) {
492   //
493   // Not implemented yet
494   //
495   return;
496
497 }
498
499
500 void AliTPCcalibCosmic::BinLogX(TH1 *h) {
501
502   // Method for the correct logarithmic binning of histograms
503
504   TAxis *axis = h->GetXaxis();
505   int bins = axis->GetNbins();
506
507   Double_t from = axis->GetXmin();
508   Double_t to = axis->GetXmax();
509   Double_t *new_bins = new Double_t[bins + 1];
510    
511   new_bins[0] = from;
512   Double_t factor = pow(to/from, 1./bins);
513   
514   for (int i = 1; i <= bins; i++) {
515    new_bins[i] = factor * new_bins[i-1];
516   }
517   axis->Set(bins, new_bins);
518   delete new_bins;
519   
520 }
521
522 AliExternalTrackParam *AliTPCcalibCosmic::Invert(AliExternalTrackParam *input)
523 {
524   //
525   // Invert paramerameter  - not covariance yet
526   //
527   AliExternalTrackParam *output = new AliExternalTrackParam(*input);
528   Double_t * param = (Double_t*)output->GetParameter();
529   param[0]*=-1;
530   param[3]*=-1;
531   param[4]*=-1;
532   //
533   return output;
534 }
535
536 AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
537   //
538   // 
539   //
540   AliExternalTrackParam *par1R= new AliExternalTrackParam(*track1);
541   par1R->Rotate(track0->GetAlpha());
542   //
543   //
544   Double_t * param = (Double_t*)par1R->GetParameter();
545   Double_t * covar = (Double_t*)par1R->GetCovariance();
546   param[0]*=1;  //OK
547   param[1]*=1;  //OK
548   param[2]*=1;  //?
549   param[3]*=-1; //OK
550   param[4]*=-1; //OK
551   //
552   covar[6] *=-1.; covar[7] *=-1.; covar[8] *=-1.;
553   //covar[10]*=-1.; covar[11]*=-1.; covar[12]*=-1.;
554   covar[13]*=-1.;
555   par1R->PropagateTo(track0->GetX(),0); // bz shold be set -
556   //if (1){
557   //  printf("Print param\n");
558   //  track1->Print();
559   //  par1R->Print();
560   //}
561   return par1R;
562 }
563
564 void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
565   //
566   // Update track 1 with track 2
567   //
568   //
569   //
570   TMatrixD vecXk(5,1);    // X vector
571   TMatrixD covXk(5,5);    // X covariance 
572   TMatrixD matHk(5,5);    // vector to mesurement
573   TMatrixD measR(5,5);    // measurement error 
574   TMatrixD vecZk(5,1);    // measurement
575   //
576   TMatrixD vecYk(5,1);    // Innovation or measurement residual
577   TMatrixD matHkT(5,5);
578   TMatrixD matSk(5,5);    // Innovation (or residual) covariance
579   TMatrixD matKk(5,5);    // Optimal Kalman gain
580   TMatrixD mat1(5,5);     // update covariance matrix
581   TMatrixD covXk2(5,5);   // 
582   TMatrixD covOut(5,5);
583   //
584   Double_t *param1=(Double_t*) track1.GetParameter();
585   Double_t *covar1=(Double_t*) track1.GetCovariance();
586   Double_t *param2=(Double_t*) track2.GetParameter();
587   Double_t *covar2=(Double_t*) track2.GetCovariance();
588   //
589   // copy data to the matrix
590   for (Int_t ipar=0; ipar<5; ipar++){
591     vecXk(ipar,0) = param1[ipar];
592     vecZk(ipar,0) = param2[ipar];
593     for (Int_t jpar=0; jpar<5; jpar++){
594       covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
595       measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
596     }
597   }
598   //
599   //
600   //
601   //
602   matHk(0,0)=1; matHk(1,1)= 1; matHk(2,2)= 1;  
603   matHk(3,3)= 1;    matHk(4,4)= 1;           // vector to measurement
604   //
605   vecYk = vecZk-matHk*vecXk;                 // Innovation or measurement residual
606   matHkT=matHk.T(); matHk.T();
607   matSk = (matHk*(covXk*matHkT))+measR;      // Innovation (or residual) covariance
608   matSk.Invert();
609   matKk = (covXk*matHkT)*matSk;              //  Optimal Kalman gain
610   vecXk += matKk*vecYk;                      //  updated vector 
611   mat1(0,0)=1; mat1(1,1)=1; mat1(2,2)=1; mat1(3,3)=1; mat1(4,4)=1;
612   covXk2 = (mat1-(matKk*matHk));
613   covOut =  covXk2*covXk; 
614   //
615   //
616   //
617   // copy from matrix to parameters
618   if (0) {
619     vecXk.Print();
620     vecZk.Print();
621     //
622     measR.Print();
623     covXk.Print();
624     covOut.Print();
625     //
626     track1.Print();
627     track2.Print();
628   }
629
630   for (Int_t ipar=0; ipar<5; ipar++){
631     param1[ipar]= vecXk(ipar,0) ;
632     for (Int_t jpar=0; jpar<5; jpar++){
633       covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
634     }
635   }
636 }
637
638 void AliTPCcalibCosmic::ProcessTree(TTree * chainTracklet, AliExternalComparison *comp){
639   //
640   // Process the debug streamer tree
641   // Possible to modify selection criteria
642   //
643   TTreeSRedirector * cstream = new TTreeSRedirector("cosmicdump.root");
644   //AliTPCcalibCosmic *cosmic = this;
645   //
646   AliExternalTrackParam * tr0 = 0;
647   AliExternalTrackParam * tr1 = 0;  
648   Int_t npoints =0;
649   {
650     Int_t entries=chainTracklet->GetEntries();
651     for (Int_t i=0; i< entries; i++){
652       chainTracklet->GetBranch("Tr0.")->SetAddress(&tr0);
653       chainTracklet->GetBranch("Tr1.")->SetAddress(&tr1);
654       chainTracklet->GetEntry(i);
655       if (!tr0) continue;
656       if (!tr1) continue;
657       if (tr0->GetY()==0) continue;
658       if (tr1->GetY()==0) continue;
659       // make a local copy
660       AliExternalTrackParam par0(*tr0);
661       AliExternalTrackParam par1(*tr1);
662       AliExternalTrackParam par1R(*tr1);
663       par1R.Rotate(par1.GetAlpha()+TMath::Pi());
664       AliExternalTrackParam *par1T = MakeTrack(tr0,tr1);
665       if (0) {
666         printf("%d\t%d\t\n",i, npoints);
667         par1R.Print();
668         par1T->Print();
669       }
670       AliExternalTrackParam par0U=par0;
671       AliExternalTrackParam par1U=*par1T;
672       //
673       UpdateTrack(par0U,*par1T);
674       UpdateTrack(par1U,par0);
675       //
676       //
677       if (i%100==0) printf("%d\t%d\tt\n",i, npoints);
678       Bool_t accept =   comp->AcceptPair(&par0,par1T);  
679
680       if (1||fStreamLevel>0){
681         (*cstream)<<"Tracklet"<<
682           "accept="<<accept<<
683           "tr0.="<<&par0<<       //original track  up
684           "tr1.="<<&par1<<       //original track  down
685           "tr1R.="<<&par1R<<     //track1 rotated to  0 frame 
686           "tr1T.="<<par1T<<      //track1 transformed to the track 0 frame 
687           //
688           "tr0U.="<<&par0U<<     //track 0 updated with track 1
689           "tr1U.="<<&par1U<<     //track 1 updated with track 0 
690           "\n";
691       }
692       //
693       if (accept) {
694         npoints++;
695         if (comp) comp->Process(&par0,par1T);
696       }
697       delete par1T;
698     }
699   }
700   delete cstream;
701 }
702
703
704
705
706
707
708
709 /*
710
711 void Init(){
712   
713 .x ~/UliStyle.C
714 .x ~/rootlogon.C
715 gSystem->Load("libSTAT.so");
716 gSystem->Load("libANALYSIS");
717 gSystem->Load("libTPCcalib");
718 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
719
720 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
721 AliXRDPROOFtoolkit tool; 
722 TChain * chainCosmic = tool.MakeChain("cosmic.txt","Track0",0,1000000);
723 chainCosmic->Lookup();
724
725 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01");  // OK
726 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2");     // OK
727 TCut cutP1("cutP1","abs(Tr0.fP[1]-Tr1.fP[1])<20");   // OK
728 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
729 TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>50");
730 TCut cutM("cutM","abs(mag)>0.01");
731 TCut cutA=cutT+cutD+cutPt+cutN+cutP1+"trigger!=16";
732
733 TCut cuthpt("abs(Tr0.fP[4])+abs(Tr1.fP[4])<0.2");
734 TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
735
736 //
737 chainCosmic->Draw(">>listELP",cutA,"entryList");
738 //TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
739 //TEntryList *elist = (TEntryList*)gProof->GetOutputList()->At(1);
740 chainCosmic->SetEntryList(elist);
741 //
742 chainCosmic->Draw(">>listV40Z100","abs(d0)<40&&abs(v01)<100","entryList");
743 TEntryList *elistV40Z100 = (TEntryList*)gDirectory->Get("listV40Z100");
744 chainCosmic->SetEntryList(elistV40Z100);
745
746 //
747 // Aliases
748 //
749 chainCosmic->SetAlias("side","(-1+(Tr0.fP[1]>0)*2)");
750 chainCosmic->SetAlias("hpt","abs(Tr0.fP[4])<0.2");
751 chainCosmic->SetAlias("signy","(-1+(Tr0.fP[0]>0)*2)");
752
753 chainCosmic->SetAlias("dy","Tr0.fP[0]+Tr1.fP[0]");
754 chainCosmic->SetAlias("dz","Tr0.fP[1]-Tr1.fP[1]");
755 chainCosmic->SetAlias("d1pt","Tr0.fP[4]+Tr1.fP[4]");
756 chainCosmic->SetAlias("dtheta","(Tr0.fP[3]+Tr1.fP[3])");
757 chainCosmic->SetAlias("dphi","(Tr0.fAlpha-Tr1.fAlpha-pi)");
758
759 chainCosmic->SetAlias("mtheta","(Tr0.fP[3]-Tr1.fP[3])*0.5")
760 chainCosmic->SetAlias("sa","(sin(Tr0.fAlpha+0.))");
761 chainCosmic->SetAlias("ca","(cos(Tr0.fAlpha+0.))");
762
763
764
765 chainCosmic->Draw("dy:sqrt(abs(Tr0.fP[4]))>>hisdyA(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
766 hisdyA->FitSlicesY();
767 hisdyA_2->SetXTitle("#sqrt{1/p_{t}}");
768 hisdyA_2->SetYTitle("#sigma_{y}(cm)");
769 hisdyA_2->SetTitle("Cosmic - Y matching");
770 hisdyA_2->SetMaximum(0.5);
771
772
773 chainCosmic->Draw("dy:sqrt(abs(Tr0.fP[4]))>>hisdyC(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
774 hisdyC->FitSlicesY();
775 hisdyC_2->SetXTitle("#sqrt{1/p_{t}}");
776 hisdyC_2->SetYTitle("#sigma_{y}(cm)");
777 hisdyC_2->SetTitle("Cosmic - Y matching");
778 hisdyC_2->SetMaximum(1);
779 hisdyC_2->SetMinimum(0);
780 hisdyC_2->SetMarkerStyle(22);
781 hisdyA_2->SetMarkerStyle(21);
782 hisdyC_2->SetMarkerSize(1.5);
783 hisdyA_2->SetMarkerSize(1.5);
784 hisdyC_2->Draw();
785 hisdyA_2->Draw("same");
786 gPad->SaveAs("~/Calibration/Cosmic/pic/ymatching.gif")
787
788 chainCosmic->Draw("dz:sqrt(abs(Tr0.fP[4]))>>hisdzA(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
789 hisdzA->FitSlicesY();
790 hisdzA_2->SetXTitle("#sqrt{1/p_{t}}");
791 hisdzA_2->SetYTitle("#sigma_{z}(cm)");
792 hisdzA_2->SetTitle("Cosmic - Z matching - A side ");
793 hisdzA_2->SetMaximum(0.5);
794
795 chainCosmic->Draw("dz:sqrt(abs(Tr0.fP[4]))>>hisdzC(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
796 hisdzC->FitSlicesY();
797 hisdzC_2->SetXTitle("#sqrt{1/p_{t}}");
798 hisdzC_2->SetYTitle("#sigma_{z}(cm)");
799 hisdzC_2->SetTitle("Cosmic - Z matching");
800 hisdzC_2->SetMaximum(0.5);
801 hisdzC_2->SetMarkerStyle(22);
802 hisdzA_2->SetMarkerStyle(21);
803 hisdzC_2->SetMarkerSize(1.5);
804 hisdzA_2->SetMarkerSize(1.5);
805
806 hisdzC_2->Draw();
807 hisdzA_2->Draw("same");
808
809
810 //
811 // PICTURE 1/pt
812 //
813 chainCosmic->Draw("d1pt:sqrt(abs(Tr0.fP[4]))>>hisd1ptA(5,0,1,30,-0.1,0.1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0"+cutM);
814 hisd1ptA->FitSlicesY();
815 hisd1ptA_2->SetXTitle("#sqrt{1/p_{t}}");
816 hisd1ptA_2->SetYTitle("#sigma_{z}(cm)");
817 hisd1ptA_2->SetTitle("Cosmic - Z matching - A side ");
818 hisd1ptA_2->SetMaximum(0.5);
819
820 chainCosmic->Draw("d1pt:sqrt(abs(Tr0.fP[4]))>>hisd1ptC(5,0,1,30,-0.1,0.1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0"+cutM);
821 hisd1ptC->FitSlicesY();
822 hisd1ptC_2->SetXTitle("#sqrt{1/p_{t}}");
823 hisd1ptC_2->SetYTitle("#sigma_{1/pt}(1/GeV)");
824 hisd1ptC_2->SetTitle("Cosmic - 1/pt matching");
825 hisd1ptC_2->SetMaximum(0.05);
826 hisd1ptC_2->SetMarkerStyle(22);
827 hisd1ptA_2->SetMarkerStyle(21);
828 hisd1ptC_2->SetMarkerSize(1.5);
829 hisd1ptA_2->SetMarkerSize(1.5);
830
831 hisd1ptC_2->Draw();
832 hisd1ptA_2->Draw("same");
833 gPad->SaveAs("~/Calibration/Cosmic/pic/1ptmatching.gif")
834
835 //
836 // Theta
837 //
838 chainCosmic->Draw("dtheta:sqrt(abs(Tr0.fP[4]))>>hisdthetaA(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
839 hisdthetaA->FitSlicesY();
840 hisdthetaA_2->SetXTitle("#sqrt{1/p_{t}}");
841 hisdthetaA_2->SetYTitle("#sigma_{#theta}(cm)");
842 hisdthetaA_2->SetTitle("Cosmic - Z matching - A side ");
843 hisdthetaA_2->SetMaximum(0.5);
844
845 chainCosmic->Draw("dtheta:sqrt(abs(Tr0.fP[4]))>>hisdthetaC(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
846 hisdthetaC->FitSlicesY();
847 hisdthetaC_2->SetXTitle("#sqrt{1/p_{t}}");
848 hisdthetaC_2->SetYTitle("#sigma_{#theta}(rad)");
849 hisdthetaC_2->SetTitle("Cosmic - Theta matching");
850 hisdthetaC_2->SetMaximum(0.01);
851 hisdthetaC_2->SetMinimum(0.0);
852 hisdthetaC_2->SetMarkerStyle(22);
853 hisdthetaA_2->SetMarkerStyle(21);
854 hisdthetaC_2->SetMarkerSize(1.5);
855 hisdthetaA_2->SetMarkerSize(1.5);
856
857 hisdthetaC_2->Draw();
858 hisdthetaA_2->Draw("same");
859 gPad->SaveAs("~/Calibration/Cosmic/pic/thetamatching.gif")
860 //
861 // Phi
862 //
863 chainCosmic->Draw("dphi:sqrt(abs(Tr0.fP[4]))>>hisdphiA(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0");
864 hisdphiA->FitSlicesY();
865 hisdphiA_2->SetXTitle("#sqrt{1/p_{t}}");
866 hisdphiA_2->SetYTitle("#sigma_{#phi}(rad)");
867 hisdphiA_2->SetTitle("Cosmic - Z matching - A side ");
868 hisdphiA_2->SetMaximum(0.5);
869
870 chainCosmic->Draw("dphi:sqrt(abs(Tr0.fP[4]))>>hisdphiC(5,0,1,30,-0.01,0.01)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]<0");
871 hisdphiC->FitSlicesY();
872 hisdphiC_2->SetXTitle("#sqrt{1/p_{t}}");
873 hisdphiC_2->SetYTitle("#sigma_{#phi}(rad)");
874 hisdphiC_2->SetTitle("Cosmic - Phi matching");
875 hisdphiC_2->SetMaximum(0.01);
876 hisdphiC_2->SetMinimum(0.0);
877 hisdphiC_2->SetMarkerStyle(22);
878 hisdphiA_2->SetMarkerStyle(21);
879 hisdphiC_2->SetMarkerSize(1.5);
880 hisdphiA_2->SetMarkerSize(1.5);
881
882 hisdphiC_2->Draw();
883 hisdphiA_2->Draw("same");
884 gPad->SaveAs("~/Calibration/Cosmic/pic/phimatching.gif")
885
886
887
888 }
889
890
891 */
892
893
894 /*
895 void MatchTheta(){
896
897 TStatToolkit toolkit;
898 Double_t chi2=0;
899 Int_t    npoints=0;
900 TVectorD fitParamA0;
901 TVectorD fitParamA1;
902 TVectorD fitParamC0;
903 TVectorD fitParamC1;
904 TMatrixD covMatrix;
905
906
907 TString fstring="";
908 // 
909 fstring+="mtheta++";
910 fstring+="ca++";
911 fstring+="sa++";
912 fstring+="ca*mtheta++";
913 fstring+="sa*mtheta++";
914 //
915 fstring+="side++";
916 fstring+="side*mtheta++";
917 fstring+="side*ca++";
918 fstring+="side*sa++";
919 fstring+="side*ca*mtheta++";
920 fstring+="side*sa*mtheta++";
921
922
923 TString *strTheta0 = toolkit.FitPlane(chain,"dtheta",fstring->Data(), "hpt&&!crossI&&!crossO", chi2,npoints,fitParamA0,covMatrix,0.8);
924 chainCosmic->SetAlias("dtheta0",strTheta0.Data());
925 strTheta0->Tokenize("+")->Print();
926
927
928 //fstring+="mtheta++";
929 //fstring+="mtheta^2++";
930 //fstring+="ca*mtheta^2++";
931 //fstring+="sa*mtheta^2++";
932
933
934
935 }
936
937 */
938
939
940
941
942 /*
943  void PosCorrection()
944
945  
946
947  
948  TStatToolkit toolkit;
949  Double_t chi2=0;
950  Int_t    npoints=0;
951  TVectorD fitParam;
952  TMatrixD covMatrix;
953  
954  //Theta
955 chainCosmic->SetAlias("dthe","(Tr0.fP[3]+Tr1.fP[3])");
956 chainCosmic->SetAlias("sign","(-1+(Tr0.fP[1]>0)*2)");
957 chainCosmic->SetAlias("di","(1-abs(Tr0.fP[1])/250)");
958 //
959 //
960 TString strFit="";
961 //
962 strFit+="sign++";                              //1
963 strFit+="Tr0.fP[3]++";                         //2
964 // 
965 strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]>0)++";     //3
966 strFit+="sin(Tr0.fAlpha)*(Tr0.fP[1]<0)++";     //4
967 //
968 strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]>0)++";     //5
969 strFit+="cos(Tr0.fAlpha)*(Tr0.fP[1]<0)++";     //6  
970 //
971 strFit+="sin(Tr0.fAlpha)*Tr0.fP[3]++";         //7
972 strFit+="cos(Tr0.fAlpha)*Tr0.fP[3]++";         //8
973  
974  
975  //                                         
976  TString * thetaParam = toolkit.FitPlane(chain,"dthe", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix,0.8,0,10000)
977  
978  chainCosmic->SetAlias("corTheta",thetaParam->Data());
979  chainCosmic->Draw("dthe:Tr0.fP[1]","","",50000);
980
981
982
983 */
984
985
986
987 /*
988
989 void AliTPCcalibCosmic::dEdxCorrection(){
990   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.01");  // OK
991   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<2");     // OK
992   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
993   TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>100");
994   TCut cutS("cutS","Orig0.fIp.fP[1]*Orig1.fIp.fP[1]>0");
995   TCut cutA=cutT+cutD+cutPt+cutN+cutS;
996
997
998  .x ~/UliStyle.C
999   gSystem->Load("libANALYSIS");
1000   gSystem->Load("libTPCcalib");
1001   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
1002   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
1003  AliXRDPROOFtoolkit tool; 
1004   TChain * chainCosmic = tool.MakeChain("cosmic.txt","Track0",0,1000000);
1005   chainCosmic->Lookup();
1006   
1007   chainCosmic->Draw(">>listEL",cutA,"entryList");
1008   TEntryList *elist = (TEntryList*)gDirectory->Get("listEL");
1009   chainCosmic->SetEntryList(elist);
1010
1011   .x ~/rootlogon.C
1012    gSystem->Load("libSTAT.so");
1013    TStatToolkit toolkit;
1014   Double_t chi2=0;
1015   Int_t    npoints=0;
1016   TVectorD fitParam;
1017   TMatrixD covMatrix;
1018   
1019   chainCosmic->Draw("Tr0.fP[4]+Tr1.fP[4]","OK"+cutA);
1020   
1021   TString strFit;
1022   strFit+="(Tr0.fP[1]/250)++";
1023   strFit+="(Tr0.fP[1]/250)^2++";
1024   strFit+="(Tr0.fP[3])++";
1025   strFit+="(Tr0.fP[3])^2++";
1026
1027   TString * ptParam = TStatToolkit::FitPlane(chain,"Tr0.fP[4]+Tr1.fP[4]", strFit.Data(),"1", chi2,npoints,fitParam,covMatrix) 
1028
1029
1030
1031 */
1032
1033
1034 /*
1035 gSystem->Load("libANALYSIS");
1036 gSystem->Load("libSTAT");
1037 gSystem->Load("libTPCcalib");
1038
1039 TStatToolkit toolkit;
1040 Double_t chi2;
1041 TVectorD fitParam;
1042 TMatrixD covMatrix;
1043 Int_t npoints;
1044 //
1045 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03");  // OK
1046 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5");     // OK
1047 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.2&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
1048 TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>110");
1049 TCut cutA=cutT+cutD+cutPt+cutN;
1050
1051
1052
1053 TTree * chainCosmic = Track0;
1054
1055
1056 chainCosmic->SetAlias("norm","signalTot0.fElements[3]/signalTot1.fElements[3]");
1057 //
1058 chainCosmic->SetAlias("dr1","(signalTot0.fElements[1]/signalTot0.fElements[3])");
1059 chainCosmic->SetAlias("dr2","(signalTot0.fElements[2]/signalTot0.fElements[3])");
1060 chainCosmic->SetAlias("dr4","(signalTot0.fElements[4]/signalTot0.fElements[3])");
1061 chainCosmic->SetAlias("dr5","(signalTot0.fElements[5]/signalTot0.fElements[3])");
1062
1063 TString fstring="";  
1064 fstring+="dr1++";
1065 fstring+="dr2++";
1066 fstring+="dr4++";
1067 fstring+="dr5++";
1068 //
1069 fstring+="dr1*dr2++";
1070 fstring+="dr1*dr4++";
1071 fstring+="dr1*dr5++";
1072 fstring+="dr2*dr4++";
1073 fstring+="dr2*dr5++";
1074 fstring+="dr4*dr5++";
1075
1076
1077
1078 TString *strqdedx = toolkit.FitPlane(chain,"norm",fstring->Data(), cutA, chi2,npoints,fitParam,covMatrix,-1,0,200000);
1079   
1080 chainCosmic->SetAlias("corQT",strqdedx->Data());
1081
1082 */
1083
1084
1085 /*
1086   chainCosmic->SetProof(kTRUE);
1087   chainCosmic->Draw("Seed0.CookdEdxNorm(0,0.6,1,0,159,0,kTRUE,kTRUE):Seed0.CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)",""+cutA,"",100000);
1088
1089
1090 chainCosmic->Draw("Seed0.CookdEdxNorm(0,0.6,1,0,159,0,kTRUE,kTRUE)/Seed1.CookdEdxNorm(0,0.6,1,0,159,0,kTRUE,kTRUE)>>his(100,0.5,1.5)","min(Orig0.fTPCncls,Orig1.fTPCncls)>130"+cutA,"",50000);
1091
1092 */
1093
1094
1095 /*
1096 chainCosmic->Draw("Tr0.fP[1]-Tr1.fP[1]:sqrt(abs(Tr0.fP[4]))>>hisdzA(5,0,1,50,-1,1)","!crossO&&!crossI&&abs(d0)<40&&Tr0.fP[1]>0&&abs(mag)>0.1"+cutA); 
1097
1098 TGraph *grdzA = (TGraph*)gProof->GetOutputList()->At(1)->Clone();
1099
1100
1101
1102  
1103 */
1104
1105
1106
1107