*** V interface for TPCCalibTasks ***
[u/mrichter/AliRoot.git] / TPC / Calib / AliTPCcalibCosmic.cxx
1
2
3 /**************************************************************************
4  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5  *                                                                        *
6  * Author: The ALICE Off-line Project.                                    *
7  * Contributors are mentioned in the code where appropriate.              *
8  *                                                                        *
9  * Permission to use, copy, modify and distribute this software and its   *
10  * documentation strictly for non-commercial purposes is hereby granted   *
11  * without fee, provided that the above copyright notice appears in all   *
12  * copies and that both the copyright notice and this permission notice   *
13  * appear in the supporting documentation. The authors make no claims     *
14  * about the suitability of this software for any purpose. It is          *
15  * provided "as is" without express or implied warranty.                  *
16  **************************************************************************/
17
18 /*
19     Comments to be written here: 
20     1. What do we calibrate.
21     2. How to interpret results
22     3. Simple example
23     4. Analysis using debug streamers.
24
25
26
27     3.Simple example
28     // To make cosmic scan the user interaction neccessary
29     //
30      
31   */
32
33
34
35 #include "Riostream.h"
36 #include "TChain.h"
37 #include "TTree.h"
38 #include "TH1F.h"
39 #include "TH2F.h"
40 #include "TList.h"
41 #include "TMath.h"
42 #include "TCanvas.h"
43 #include "TFile.h"
44 #include "TF1.h"
45 #include "THnSparse.h"
46 #include "TDatabasePDG.h"
47
48 #include "AliTPCclusterMI.h"
49 #include "AliTPCseed.h"
50 #include "AliESDVertex.h"
51 #include "AliESDEvent.h"
52 #include "AliESDtrack.h"
53 #include "AliESDfriend.h"
54 #include "AliESDInputHandler.h"
55 #include "AliAnalysisManager.h"
56
57 #include "AliVEvent.h"
58 #include "AliVTrack.h"
59 #include "AliVfriendEvent.h"
60 #include "AliVfriendTrack.h"
61
62 #include "AliTracker.h"
63 #include "AliMagF.h"
64 #include "AliTPCCalROC.h"
65 #include "AliTPCParam.h"
66 #include "AliLog.h"
67
68 #include "AliTPCcalibCosmic.h"
69 #include "TTreeStream.h"
70 #include "AliTPCTracklet.h"
71 //#include "AliESDcosmic.h"
72 #include "AliRecoParam.h"
73 #include "AliMultiplicity.h"
74 #include "AliTPCTransform.h"
75 #include "AliTPCcalibDB.h"
76 #include "AliTPCseed.h"
77 #include "AliGRPObject.h"
78 #include "AliTPCCorrection.h"
79 ClassImp(AliTPCcalibCosmic)
80
81
82 AliTPCcalibCosmic::AliTPCcalibCosmic() 
83   :AliTPCcalibBase(),
84    fHistNTracks(0),
85    fClusters(0),
86    fModules(0),
87    fHistPt(0),
88    fDeDx(0),
89    fDeDxMIP(0),
90    fMIPvalue(1), 
91    fCutMaxD(5),        // maximal distance in rfi ditection
92    fCutMaxDz(40),      // maximal distance in z ditection
93    fCutTheta(0.03),    // maximal distan theta
94    fCutMinDir(-0.99),   // direction vector products
95    fCosmicTree(0)      // tree with cosmic data
96 {  
97   //
98   // CONSTRUCTOR - SEE COMMENTS ABOVE
99   //
100   AliInfo("Default Constructor");    
101   for (Int_t ihis=0; ihis<6;ihis++){
102     fHistoDelta[ihis]=0;
103     fHistoPull[ihis]=0;
104   }
105   for (Int_t ihis=0; ihis<4;ihis++){
106     fHistodEdxMax[ihis]    =0;
107     fHistodEdxTot[ihis]    =0;
108   }
109 }
110
111
112 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title) 
113   :AliTPCcalibBase(),
114    fHistNTracks(0),
115    fClusters(0),
116    fModules(0),
117    fHistPt(0),
118    fDeDx(0),
119    fDeDxMIP(0),
120    fMIPvalue(1),
121    fCutMaxD(5),        // maximal distance in rfi ditection 
122    fCutMaxDz(40),      // maximal distance in z ditection
123    fCutTheta(0.03),    // maximal distan theta
124    fCutMinDir(-0.99),  // direction vector products
125    fCosmicTree(0)      // tree with cosmic data
126 {  
127   //
128   // cONSTRUCTOR - SEE COMENTS ABOVE
129   //
130   SetName(name);
131   SetTitle(title);
132
133   fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
134   fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",160,0,160);
135   fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-650,650,600,-700,700);
136   fHistPt = new TH1F("Pt","Pt distribution; p_{T} (GeV); counts",2000,0,50);
137   fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,100.,500,2.,1000);
138   BinLogX(fDeDx);
139   fDeDxMIP =  new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000);
140   Init();
141   AliInfo("Non Default Constructor");  
142   //
143 }
144
145 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
146   //
147   // destructor
148   //
149   for (Int_t ihis=0; ihis<6;ihis++){
150     delete fHistoDelta[ihis];
151     delete fHistoPull[ihis];
152   }
153   for (Int_t ihis=0; ihis<4;ihis++){
154     delete fHistodEdxTot[ihis];
155     delete fHistodEdxMax[ihis];
156   }
157
158   delete fHistNTracks;            //  histogram showing number of ESD tracks per event
159   delete fClusters;               //  histogram showing the number of clusters per track
160   delete fModules;                //  2d histogram of tracks which are propagated to the ACORDE scintillator array
161   delete fHistPt;                 //  Pt histogram of reconstructed tracks
162   delete fDeDx;                   //  dEdx spectrum showing the different particle types
163   delete fDeDxMIP;                //  TPC signal close to the MIP region of muons 0.4 < p < 0.45 GeV
164 }
165
166
167 void AliTPCcalibCosmic::Init(){
168   //
169   // init component
170   // Make performance histograms
171   //
172
173   // tracking performance bins
174   // 0 - delta of interest
175   // 1 - min (track0, track1) number of clusters
176   // 2 - R  - vertex radius
177   // 3 - P1 - mean z
178   // 4 - P2 - snp(phi)    at inner wall of TPC
179   // 5 - P3 - tan(theta)  at inner wall of TPC
180   // 6 - P4 - 1/pt mean
181   // 7 - pt - pt mean
182   // 8 - alpha
183   // 9 - is corss indicator
184   Int_t ndim=10;
185   Double_t xminTrack[10], xmaxTrack[10];
186   Int_t binsTrack[10];
187   TString axisName[10];
188   //
189   binsTrack[0] =100;
190   axisName[0]  ="#Delta";
191   //
192   binsTrack[1] =8;
193   xminTrack[1] =80; xmaxTrack[1]=160;
194   axisName[1]  ="N_{cl}";
195   //
196   binsTrack[2] =10;
197   xminTrack[2] =0; xmaxTrack[2]=90;  // 
198   axisName[2]  ="dca_{r} (cm)";
199   //
200   binsTrack[3] =25;
201   xminTrack[3] =-250; xmaxTrack[3]=250;  // 
202   axisName[3]  ="z (cm)";
203   //
204   binsTrack[4] =10;
205   xminTrack[4] =-0.8; xmaxTrack[4]=0.8;  // 
206   axisName[4]  ="sin(#phi)";
207   //
208   binsTrack[5] =10;
209   xminTrack[5] =-1; xmaxTrack[5]=1;  // 
210   axisName[5]  ="tan(#theta)";
211   //
212   binsTrack[6] =40;
213   xminTrack[6] =-2; xmaxTrack[6]=2;  // 
214   axisName[6]  ="1/pt (1/GeV)";
215   //
216   binsTrack[7] =50;
217   xminTrack[7] =1; xmaxTrack[7]=1000;  // 
218   axisName[7]  ="pt (GeV)";
219   //
220   binsTrack[8] =18;
221   xminTrack[8] =0; xmaxTrack[8]=TMath::Pi();  // 
222   axisName[8]  ="alpha";
223   //
224   binsTrack[9] =3;
225   xminTrack[9] =-0.1; xmaxTrack[9]=2.1;  // 
226   axisName[9]  ="cross";
227   //
228   // delta y
229   xminTrack[0] =-1; xmaxTrack[0]=1;  // 
230   fHistoDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", ndim, binsTrack,xminTrack, xmaxTrack);
231   xminTrack[0] =-5; xmaxTrack[0]=5;  // 
232   fHistoPull[0] = new THnSparseS("#Delta_{Y} (unit)","#Delta_{Y} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
233   //
234   // delta z
235   xminTrack[0] =-1; xmaxTrack[0]=1;  // 
236   fHistoDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", ndim, binsTrack,xminTrack, xmaxTrack);
237   xminTrack[0] =-5; xmaxTrack[0]=5;  // 
238   fHistoPull[1] = new THnSparseS("#Delta_{Z} (unit)","#Delta_{Z} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
239   //
240   // delta P2
241   xminTrack[0] =-10; xmaxTrack[0]=10;  // 
242   fHistoDelta[2] = new THnSparseS("#Delta_{#phi} (mrad)","#Delta_{#phi} (mrad)", ndim, binsTrack,xminTrack, xmaxTrack);
243   xminTrack[0] =-5; xmaxTrack[0]=5;  // 
244   fHistoPull[2] = new THnSparseS("#Delta_{#phi} (unit)","#Delta_{#phi} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
245   //
246   // delta P3
247   xminTrack[0] =-10; xmaxTrack[0]=10;  // 
248   fHistoDelta[3] = new THnSparseS("#Delta_{#theta} (mrad)","#Delta_{#theta} (mrad)", ndim, binsTrack,xminTrack, xmaxTrack);
249   xminTrack[0] =-5; xmaxTrack[0]=5;  // 
250   fHistoPull[3] = new THnSparseS("#Delta_{#theta} (unit)","#Delta_{#theta} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
251   //
252   // delta P4
253   xminTrack[0] =-0.2; xmaxTrack[0]=0.2;  // 
254   fHistoDelta[4] = new THnSparseS("#Delta_{1/pt} (1/GeV)","#Delta_{1/pt} (1/GeV)", ndim, binsTrack,xminTrack, xmaxTrack);
255   xminTrack[0] =-5; xmaxTrack[0]=5;  // 
256   fHistoPull[4] = new THnSparseS("#Delta_{1/pt} (unit)","#Delta_{1/pt} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
257   
258   //
259   // delta Pt
260   xminTrack[0] =-0.5; xmaxTrack[0]=0.5;  // 
261   fHistoDelta[5] = new THnSparseS("#Delta_{pt}/p_{t}","#Delta_{pt}/p_{t}", ndim, binsTrack,xminTrack, xmaxTrack);
262   xminTrack[0] =-5; xmaxTrack[0]=5;  // 
263   fHistoPull[5] = new THnSparseS("#Delta_{pt}/p_{t} (unit)","#Delta_{pt}/p_{t} (unit)", ndim, binsTrack,xminTrack, xmaxTrack);
264   //
265
266   for (Int_t idedx=0;idedx<4;idedx++){
267     xminTrack[0] =0.5; xmaxTrack[0]=1.5;  // 
268     binsTrack[1] =40;
269     xminTrack[1] =10; xmaxTrack[1]=160;
270
271     fHistodEdxMax[idedx] = new THnSparseS(Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
272                                           Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx), 
273                                           ndim, binsTrack,xminTrack, xmaxTrack);
274     fHistodEdxTot[idedx] = new THnSparseS(Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
275                                           Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx), 
276                                           ndim, binsTrack,xminTrack, xmaxTrack);
277   }
278   
279
280
281   for (Int_t ivar=0;ivar<6;ivar++){
282     for (Int_t ivar2=0;ivar2<ndim;ivar2++){      
283       fHistoDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
284       fHistoDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
285       fHistoPull[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
286       fHistoPull[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
287       BinLogX(fHistoDelta[ivar],7);
288       BinLogX(fHistoPull[ivar],7);
289       if (ivar<4){
290         fHistodEdxMax[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
291         fHistodEdxMax[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
292         fHistodEdxTot[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
293         fHistodEdxTot[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
294         BinLogX(fHistodEdxMax[ivar],7);
295         BinLogX(fHistodEdxTot[ivar],7);
296       }
297     }
298   }
299 }
300
301
302 void AliTPCcalibCosmic::Add(const AliTPCcalibCosmic* cosmic){
303   //
304   // merge the content of the cosmic componentnts
305   //
306   for (Int_t ivar=0; ivar<6;ivar++){
307     if (fHistoDelta[ivar] && cosmic->fHistoDelta[ivar]){
308       fHistoDelta[ivar]->Add(cosmic->fHistoDelta[ivar]);
309     }
310     if (fHistoPull[ivar] && cosmic->fHistoPull[ivar]){
311       fHistoPull[ivar]->Add(cosmic->fHistoPull[ivar]);
312     }
313   }
314   for (Int_t ivar=0; ivar<4;ivar++){
315     if (fHistodEdxMax[ivar] && cosmic->fHistodEdxMax[ivar]){
316       fHistodEdxMax[ivar]->Add(cosmic->fHistodEdxMax[ivar]);
317     }
318     if (fHistodEdxTot[ivar] && cosmic->fHistodEdxTot[ivar]){
319       fHistodEdxTot[ivar]->Add(cosmic->fHistodEdxTot[ivar]);
320     }
321   }
322   if (cosmic->fCosmicTree){
323     if (!fCosmicTree) {
324       fCosmicTree = new TTree("pairs","pairs");
325       fCosmicTree->SetDirectory(0);
326     }
327     AliTPCcalibCosmic::AddTree(fCosmicTree,cosmic->fCosmicTree);
328   }
329 }
330
331
332
333
334 void AliTPCcalibCosmic::Process(AliVEvent *event) {
335   //
336   // Process of the ESD event  - fill calibration components
337   //
338   if (!event) {
339     Printf("ERROR: event not available");
340     return;
341   }  
342    
343   //
344   //Int_t isOK=kTRUE;
345   // COSMIC not signed properly
346   //  UInt_t specie = event->GetEventSpecie();  // select only cosmic events
347   //if (specie==AliRecoParam::kCosmic || specie==AliRecoParam::kCalib) {
348   //  isOK = kTRUE;
349   //}
350   //if (!isOK) return;
351   // Work around
352   FindCosmicPairs(event);
353   //const AliMultiplicity *multiplicity = event->GetMultiplicity();
354   //  Int_t ntracklets = multiplicity->GetNumberOfTracklets();
355   //if (ntracklets>6) return; // filter out "normal" event with high multiplicity
356   //const TString &trigger = event->GetFiredTriggerClasses();
357   //if (trigger.Contains("C0OB0")==0) return;
358
359
360   FindPairs(event); // nearly everything takes place in find pairs...
361
362   if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
363   Int_t ntracks=event->GetNumberOfTracks(); 
364   fHistNTracks->Fill(ntracks);
365   
366 }
367
368
369 void AliTPCcalibCosmic::FillHistoPerformance(const AliExternalTrackParam *par0, const AliExternalTrackParam *par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam */*inner1*/, AliTPCseed *seed0,  AliTPCseed *seed1, const AliExternalTrackParam *param0Combined , Int_t cross){
370   //
371   // par0,par1       - parameter of tracks at DCA 0
372   // inner0,inner1   - parameter of tracks at the TPC entrance
373   // seed0, seed1    - detailed track information
374   // param0Combined  - Use combined track parameters for binning
375   //
376   Int_t kMinCldEdx =20;
377   Int_t ncl0 = seed0->GetNumberOfClusters();
378   Int_t ncl1 = seed1->GetNumberOfClusters();
379   const Double_t kpullCut    = 10;
380   Double_t x[10];
381   Double_t xyz0[3];
382   Double_t xyz1[3];
383   par0->GetXYZ(xyz0);
384   par1->GetXYZ(xyz1);
385   Double_t radius0 = TMath::Sqrt(xyz0[0]*xyz0[0]+xyz0[1]*xyz0[1]);
386   Double_t radius1 = TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
387   inner0->GetXYZ(xyz0);
388   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
389   // bin parameters
390   x[1] = TMath::Min(ncl0,ncl1);
391   x[2] = (radius0+radius1)*0.5;
392   x[3] = param0Combined->GetZ();
393   x[4] = inner0->GetSnp();
394   x[5] = param0Combined->GetTgl();
395   x[6] = param0Combined->GetSigned1Pt();
396   x[7] = param0Combined->Pt();
397   x[8] = alpha;
398   x[9] = cross;
399   // deltas
400   Double_t delta[6];
401   Double_t sigma[6];
402   delta[0] = (par0->GetY()+par1->GetY());
403   delta[1] = (par0->GetZ()-par1->GetZ());
404   delta[2] = (par0->GetAlpha()-par1->GetAlpha()-TMath::Pi());
405   delta[3] = (par0->GetTgl()+par1->GetTgl());
406   delta[4] = (par0->GetParameter()[4]+par1->GetParameter()[4]);
407   delta[5] = (par0->Pt()-par1->Pt())/((par0->Pt()+par1->Pt())*0.5);
408   //
409   sigma[0] = TMath::Sqrt(par0->GetSigmaY2()+par1->GetSigmaY2());
410   sigma[1] = TMath::Sqrt(par0->GetSigmaZ2()+par1->GetSigmaZ2());
411   sigma[2] = TMath::Sqrt(par0->GetSigmaSnp2()+par1->GetSigmaSnp2());
412   sigma[3] = TMath::Sqrt(par0->GetSigmaTgl2()+par1->GetSigmaTgl2());
413   sigma[4] = TMath::Sqrt(par0->GetSigma1Pt2()+par1->GetSigma1Pt2());
414   sigma[5] = sigma[4]*((par0->Pt()+par1->Pt())*0.5);
415   //
416   Bool_t isOK = kTRUE;
417   for (Int_t ivar=0;ivar<6;ivar++){
418     if (sigma[ivar]==0) isOK=kFALSE;
419     x[0]= delta[ivar]/sigma[ivar];
420     if (TMath::Abs(x[0])>kpullCut) isOK = kFALSE;
421   }
422   //
423
424   if (isOK) for (Int_t ivar=0;ivar<6;ivar++){
425     x[0]= delta[ivar];    // Modifiation 10.10 use not normalized deltas
426     if (ivar==2 || ivar ==3) x[0]*=1000;  // angles in radian
427     fHistoDelta[ivar]->Fill(x);
428     if (sigma[ivar]>0){
429       x[0]= delta[ivar]/sigma[ivar];
430       fHistoPull[ivar]->Fill(x);
431     }
432   }
433
434   //                                            
435   // Fill dedx performance
436   //
437   for (Int_t ipad=0; ipad<4;ipad++){
438     //
439     //
440     //
441     Int_t row0=0;
442     Int_t row1=160;
443     if (ipad==0) row1=63;
444     if (ipad==1) {row0=63; row1=63+64;}
445     if (ipad==2) {row0=128;}
446     Int_t   nclUp       = TMath::Nint(seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
447     Int_t   nclDown     = TMath::Nint(seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
448     Int_t   minCl       = TMath::Min(nclUp,nclDown);
449     if (minCl<kMinCldEdx) continue;
450     x[1] = minCl;
451     //
452     Float_t dEdxTotUp   = seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
453     Float_t dEdxTotDown = seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
454     Float_t dEdxMaxUp   = seed0->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
455     Float_t dEdxMaxDown = seed1->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
456     //
457     if (dEdxTotDown<=0) continue;
458     if (dEdxMaxDown<=0) continue;
459     x[0]=dEdxTotUp/dEdxTotDown;
460     fHistodEdxTot[ipad]->Fill(x);
461     x[0]=dEdxMaxUp/dEdxMaxDown;
462     fHistodEdxMax[ipad]->Fill(x);
463   }
464
465
466   
467 }
468
469 void AliTPCcalibCosmic::FindPairs(const AliVEvent *event){
470   //
471   // Find cosmic pairs
472   // 
473   // Track0 is choosen in upper TPC part
474   // Track1 is choosen in lower TPC part
475   //
476   if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
477   AliVfriendEvent *friendEvent=event->FindFriend();
478   Int_t ntracks=event->GetNumberOfTracks(); 
479   TObjArray  tpcSeeds(ntracks);
480   if (ntracks==0) return;
481   Double_t vtxx[3]={0,0,0};
482   Double_t svtxx[3]={0.000001,0.000001,100.};
483   AliESDVertex vtx(vtxx,svtxx);
484   //
485   //track loop
486   //
487   for (Int_t i=0;i<ntracks;++i) {
488    AliVTrack *track = event->GetVTrack(i);
489    fClusters->Fill(track->GetTPCNcls()); 
490   
491    const AliExternalTrackParam * trackIn = track->GetInnerParam();
492    const AliExternalTrackParam * trackOut = track->GetOuterParam();
493    if (!trackIn) continue;
494    if (!trackOut) continue;
495    if (ntracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue;  // filter laser 
496
497
498    const AliVfriendTrack *friendTrack = friendEvent->GetTrack(i);
499    if (!friendTrack) continue;
500    TObject *calibObject;
501    AliTPCseed *seed = 0;   
502    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
503      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
504    }
505    if (seed) tpcSeeds.AddAt(seed,i);
506
507    Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
508    if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
509      fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159));
510      //
511      if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159));
512      //
513      // if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159)>300) {
514 //        //TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
515 //        //if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
516 //        // if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
517 //      }
518
519    }
520
521   }
522
523   if (ntracks<2) return;
524   //
525   // Find pairs
526   //
527   for (Int_t i=0;i<ntracks;++i) {
528       AliVTrack *track0 = event->GetVTrack(i);
529     // track0 - choosen upper part
530     if (!track0) continue;
531     if (!track0->GetOuterParam()) continue;
532     if (track0->GetOuterParam()->GetAlpha()<0) continue;
533     Double_t dir0[3];
534     track0->GetDirection(dir0);    
535     for (Int_t j=0;j<ntracks;++j) {
536       if (i==j) continue;
537       AliVTrack *track1 = event->GetVTrack(j);
538       //track 1 lower part
539       if (!track1) continue;
540       if (!track1->GetOuterParam()) continue;
541       if (track1->GetOuterParam()->GetAlpha()>0) continue;
542       //
543       Double_t dir1[3];
544       track1->GetDirection(dir1);
545       
546       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
547       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
548       if (! seed0) continue;
549       if (! seed1) continue;
550       Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159);
551       Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159);
552       //
553       Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63);
554       Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63);
555       //
556       Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159);
557       Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159);
558       //
559       Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
560       Float_t d0  = track0->GetLinearD(0,0);
561       Float_t d1  = track1->GetLinearD(0,0);
562       //
563       // conservative cuts - convergence to be guarantied
564       // applying before track propagation
565       if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
566       if (dir>fCutMinDir) continue;               // direction vector product
567       Float_t bz = AliTracker::GetBz();
568       Float_t dvertex0[2];   //distance to 0,0
569       Float_t dvertex1[2];   //distance to 0,0 
570       track0->GetDZ(0,0,0,bz,dvertex0);
571       track1->GetDZ(0,0,0,bz,dvertex1);
572       if (TMath::Abs(dvertex0[1])>250) continue;
573       if (TMath::Abs(dvertex1[1])>250) continue;
574       //
575       //
576       //
577       Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
578       AliExternalTrackParam param0;
579       param0.CopyFromVTrack(track0);
580
581       AliExternalTrackParam param1;
582       param1.CopyFromVTrack(track1);
583       //
584       // Propagate using Magnetic field and correct fo material budget
585       //
586       Double_t sign0=-1;
587       Double_t sign1=1;
588       Double_t maxsnp=0.90;
589       AliTracker::PropagateTrackToBxByBz(&param0,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign0);
590       AliTracker::PropagateTrackToBxByBz(&param1,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign1);
591       //
592       // Propagate rest to the 0,0 DCA - z should be ignored
593       //
594       Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
595       Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
596       //      
597       param0.GetDZ(0,0,0,bz,dvertex0);
598       param1.GetDZ(0,0,0,bz,dvertex1);
599       if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
600       //
601       Double_t xyz0[3];//,pxyz0[3];
602       Double_t xyz1[3];//,pxyz1[3];
603       param0.GetXYZ(xyz0);
604       param1.GetXYZ(xyz1);
605       Bool_t isPair = IsPair(&param0,&param1);
606       //
607       if (isPair) FillAcordeHist(track0);
608       if (isPair &&param0.Pt()>1) {
609         const TString &trigger = event->GetFiredTriggerClasses(); 
610         UInt_t specie = event->GetEventSpecie();
611         printf("COSMIC ?\t%s\t%d\t%f\t%f\n", trigger.Data(),specie, param0.GetZ(), param1.GetZ()); 
612       }
613       //
614       // combined track params 
615       //
616       AliExternalTrackParam *par0U=MakeCombinedTrack(&param0,&param1);
617       AliExternalTrackParam *par1U=MakeCombinedTrack(&param1,&param0);
618       
619       //
620       if (fStreamLevel>0){
621         TTreeSRedirector * cstream =  GetDebugStreamer();
622         //printf("My stream=%p\n",(void*)cstream);
623         AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
624         AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
625         AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
626         AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
627         Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0;
628         Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0;
629         Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]);
630         Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]);
631         //
632         //
633         //
634         Int_t cross =0;  // 0 no cross, 2 cross on both sides
635         if (isCrossI) cross+=1; 
636         if (isCrossO) cross+=1; 
637         FillHistoPerformance(&param0, &param1, ip0, ip1, seed0, seed1,par0U, cross);
638         if (cstream) {
639           (*cstream) << "Track0" <<
640             "run="<<fRun<<              //  run number
641             "event="<<fEvent<<          //  event number
642             "time="<<fTime<<            //  time stamp of event
643             "trigger="<<fTrigger<<      //  trigger
644             "triggerClass="<<&fTriggerClass<<      //  trigger
645             "mag="<<fMagF<<             //  magnetic field
646             "dir="<<dir<<               //  direction
647             "OK="<<isPair<<             //  will be accepted
648             "b0="<<b0<<                 //  propagate status
649             "b1="<<b1<<                 //  propagate status
650             "crossI="<<isCrossI<<       //  cross inner
651             "crossO="<<isCrossO<<       //  cross outer
652             //
653             "Orig0.=" << track0 <<      //  original track  0
654             "Orig1.=" << track1 <<      //  original track  1
655             "Tr0.="<<&param0<<          //  track propagated to the DCA 0,0
656             "Tr1.="<<&param1<<          //  track propagated to the DCA 0,0        
657             "Ip0.="<<ip0<<              //  inner param - upper
658             "Ip1.="<<ip1<<              //  inner param - lower
659             "Op0.="<<op0<<              //  outer param - upper
660             "Op1.="<<op1<<              //  outer param - lower
661             "Up0.="<<par0U<<           //  combined track 0
662             "Up1.="<<par1U<<           //  combined track 1
663             //
664             "v00="<<dvertex0[0]<<       //  distance using kalman
665             "v01="<<dvertex0[1]<<       // 
666             "v10="<<dvertex1[0]<<       //
667             "v11="<<dvertex1[1]<<       // 
668             "d0="<<d0<<                 //  linear distance to 0,0
669             "d1="<<d1<<                 //  linear distance to 0,0
670             //
671             //
672             //
673             "x00="<<xyz0[0]<<           // global position close to vertex
674             "x01="<<xyz0[1]<<
675             "x02="<<xyz0[2]<<
676             //
677             "x10="<<xyz1[0]<<           // global position close to vertex
678             "x11="<<xyz1[1]<<
679             "x12="<<xyz1[2]<<
680             //
681             "alpha0="<<alpha0<<
682             "alpha1="<<alpha1<<
683             "dir00="<<dir0[0]<<           // direction upper
684             "dir01="<<dir0[1]<<
685             "dir02="<<dir0[2]<<
686             //
687             "dir10="<<dir1[0]<<           // direction lower
688             "dir11="<<dir1[1]<<
689             "dir12="<<dir1[2]<<
690             //
691             //
692             "Seed0.=" << seed0 <<       //  original seed 0
693             "Seed1.=" << seed1 <<       //  original seed 1
694             //
695             "dedx0="<<dedx0<<           //  dedx0 - all
696             "dedx1="<<dedx1<<           //  dedx1 - all
697             //
698             "dedx0I="<<dedx0I<<         //  dedx0 - inner ROC
699             "dedx1I="<<dedx1I<<         //  dedx1 - inner ROC
700             //
701             "dedx0O="<<dedx0O<<         //  dedx0 - outer ROC
702             "dedx1O="<<dedx1O<<         //  dedx1 - outer ROC
703             "\n";
704         }
705       }
706       delete par0U;
707       delete par1U;
708     }
709   }  
710 }    
711
712
713
714
715 void  AliTPCcalibCosmic::FillAcordeHist(AliVTrack *upperTrack) {
716
717   // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
718   if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
719     
720   const Double_t acordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
721   const Double_t roof = 210.5;       // distance from x =0 to end of magnet roof
722
723   Double_t r[3];
724   upperTrack->GetXYZ(r);
725   Double_t d[3];
726   upperTrack->GetDirection(d);
727   Double_t x,z;
728   z = r[2] + (d[2]/d[1])*(acordePlane - r[1]);
729   x = r[0] + (d[0]/d[1])*(acordePlane - r[1]);
730   
731   if (x > roof) {
732     x = r[0] + (d[0]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
733     z = r[2] + (d[2]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
734   }
735   if (x < -roof) {
736     x = r[0] + (d[0]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);       
737     z = r[2] + (d[2]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);
738   } 
739
740   fModules->Fill(z, x);
741  
742 }
743
744
745
746 Long64_t AliTPCcalibCosmic::Merge(TCollection *const li) {
747   //
748   // component merging
749   //
750
751   TIterator* iter = li->MakeIterator();
752   AliTPCcalibCosmic* cal = 0;
753
754   while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
755     if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
756       //Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
757       return -1;
758     }
759     
760     fHistNTracks->Add(cal->GetHistNTracks());
761     fClusters->Add(cal-> GetHistClusters());
762     fModules->Add(cal->GetHistAcorde());
763     fHistPt->Add(cal->GetHistPt());
764     fDeDx->Add(cal->GetHistDeDx());
765     fDeDxMIP->Add(cal->GetHistMIP());
766     Add(cal);
767   }
768   return 0;
769   
770 }
771
772
773 Bool_t  AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1) const{
774   //
775   //
776   /*
777   // 0. Same direction - OPOSITE  - cutDir +cutT    
778   TCut cutDir("cutDir","dir<-0.99")
779   // 1. 
780   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
781   //
782   // 2. The same rphi 
783   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
784   //
785   //
786   //
787   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");  
788   // 1/Pt diff cut
789   */
790   const Double_t *p0 = tr0->GetParameter();
791   const Double_t *p1 = tr1->GetParameter();
792   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
793   if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
794   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
795   
796   Double_t d0[3], d1[3];
797   tr0->GetDirection(d0);    
798   tr1->GetDirection(d1);       
799   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
800   //
801   return kTRUE;  
802 }
803  
804
805
806 Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
807   //
808   // Calculate the MIP value - gaussian fit used
809   //
810
811   TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
812   funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
813                                 hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
814   hist->Fit(funcDoubleGaus);
815   Double_t aMIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
816
817   delete funcDoubleGaus;
818
819   return aMIPvalue;
820
821 }
822
823
824
825
826 void AliTPCcalibCosmic::CalculateBetheParams(TH2F */*hist*/, Double_t * /*initialParam*/) {
827   //
828   // Not implemented yet
829   //
830   return;
831
832 }
833
834
835 void AliTPCcalibCosmic::BinLogX(THnSparse *const h, Int_t axisDim) {
836
837   // Method for the correct logarithmic binning of histograms
838
839   TAxis *axis = h->GetAxis(axisDim);
840   int bins = axis->GetNbins();
841
842   Double_t from = axis->GetXmin();
843   Double_t to = axis->GetXmax();
844   Double_t *newBins = new Double_t[bins + 1];
845
846   newBins[0] = from;
847   Double_t factor = pow(to/from, 1./bins);
848
849   for (int i = 1; i <= bins; i++) {
850    newBins[i] = factor * newBins[i-1];
851   }
852   axis->Set(bins, newBins);
853   delete [] newBins;
854
855 }
856
857
858 void AliTPCcalibCosmic::BinLogX(TH1 *const h) {
859
860   // Method for the correct logarithmic binning of histograms
861
862   TAxis *axis = h->GetXaxis();
863   int bins = axis->GetNbins();
864
865   Double_t from = axis->GetXmin();
866   Double_t to = axis->GetXmax();
867   Double_t *newBins = new Double_t[bins + 1];
868    
869   newBins[0] = from;
870   Double_t factor = pow(to/from, 1./bins);
871   
872   for (int i = 1; i <= bins; i++) {
873    newBins[i] = factor * newBins[i-1];
874   }
875   axis->Set(bins, newBins);
876   delete [] newBins;
877   
878 }
879
880
881 AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
882   //
883   // Make a atrack using the kalman update of track0 and track1
884   //
885   AliExternalTrackParam *par1R= new AliExternalTrackParam(*track1);
886   par1R->Rotate(track0->GetAlpha());
887   par1R->PropagateTo(track0->GetX(),AliTracker::GetBz()); 
888   //
889   //
890   Double_t * param = (Double_t*)par1R->GetParameter();
891   Double_t * covar = (Double_t*)par1R->GetCovariance();
892
893   param[0]*=1;  //OK
894   param[1]*=1;  //OK
895   param[2]*=1;  //?
896   param[3]*=-1; //OK
897   param[4]*=-1; //OK
898   //
899   covar[6] *=-1.; covar[7] *=-1.; covar[8] *=-1.;
900   //covar[10]*=-1.; covar[11]*=-1.; covar[12]*=-1.;
901   covar[13]*=-1.;
902   return par1R;
903 }
904
905 AliExternalTrackParam *AliTPCcalibCosmic::MakeCombinedTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
906   //
907   // Make combined track
908   //
909   //
910   AliExternalTrackParam * par1T = MakeTrack(track0,track1);
911   AliExternalTrackParam * par0U = new AliExternalTrackParam(*track0);
912   //
913   UpdateTrack(*par0U,*par1T);
914   delete par1T;
915   return par0U;
916 }
917
918
919 void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
920   //
921   // Update track 1 with track 2
922   //
923   //
924   //
925   TMatrixD vecXk(5,1);    // X vector
926   TMatrixD covXk(5,5);    // X covariance 
927   TMatrixD matHk(5,5);    // vector to mesurement
928   TMatrixD measR(5,5);    // measurement error 
929   TMatrixD vecZk(5,1);    // measurement
930   //
931   TMatrixD vecYk(5,1);    // Innovation or measurement residual
932   TMatrixD matHkT(5,5);
933   TMatrixD matSk(5,5);    // Innovation (or residual) covariance
934   TMatrixD matKk(5,5);    // Optimal Kalman gain
935   TMatrixD mat1(5,5);     // update covariance matrix
936   TMatrixD covXk2(5,5);   // 
937   TMatrixD covOut(5,5);
938   //
939   Double_t *param1=(Double_t*) track1.GetParameter();
940   Double_t *covar1=(Double_t*) track1.GetCovariance();
941   Double_t *param2=(Double_t*) track2.GetParameter();
942   Double_t *covar2=(Double_t*) track2.GetCovariance();
943   //
944   // copy data to the matrix
945   for (Int_t ipar=0; ipar<5; ipar++){
946     for (Int_t jpar=0; jpar<5; jpar++){
947       covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
948       measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
949       matHk(ipar,jpar)=0;
950       mat1(ipar,jpar)=0;
951     }
952     vecXk(ipar,0) = param1[ipar];
953     vecZk(ipar,0) = param2[ipar];
954     matHk(ipar,ipar)=1;
955     mat1(ipar,ipar)=0;
956   }
957   //
958   //
959   //
960   //
961   //
962   vecYk = vecZk-matHk*vecXk;                 // Innovation or measurement residual
963   matHkT=matHk.T(); matHk.T();
964   matSk = (matHk*(covXk*matHkT))+measR;      // Innovation (or residual) covariance
965   matSk.Invert();
966   matKk = (covXk*matHkT)*matSk;              //  Optimal Kalman gain
967   vecXk += matKk*vecYk;                      //  updated vector 
968   covXk2 = (mat1-(matKk*matHk));
969   covOut =  covXk2*covXk; 
970   //
971   //
972   //
973   // copy from matrix to parameters
974   if (0) {
975     vecXk.Print();
976     vecZk.Print();
977     //
978     measR.Print();
979     covXk.Print();
980     covOut.Print();
981     //
982     track1.Print();
983     track2.Print();
984   }
985
986   for (Int_t ipar=0; ipar<5; ipar++){
987     param1[ipar]= vecXk(ipar,0) ;
988     for (Int_t jpar=0; jpar<5; jpar++){
989       covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
990     }
991   }
992 }
993
994
995
996 void AliTPCcalibCosmic::FindCosmicPairs(const AliVEvent *event) {
997   //
998   // find cosmic pairs trigger by random trigger
999   //
1000   //
1001   AliESDVertex vtxSPD;
1002   event->GetPrimaryVertexSPD(vtxSPD);
1003   AliESDVertex *vertexSPD=&vtxSPD;
1004
1005   AliESDVertex vtxTPC;
1006   event->GetPrimaryVertexTPC(vtxTPC);
1007   AliESDVertex *vertexTPC=&vtxTPC;
1008
1009   AliVfriendEvent *friendEvent=event->FindFriend();
1010   const Double_t kMinPt=1;
1011   const Double_t kMinPtMax=0.8;
1012   const Double_t kMinNcl=50;
1013   const Double_t kMaxDelta[5]={2,600,0.02,0.02,0.1};
1014   Int_t ntracks=event->GetNumberOfTracks(); 
1015   //  Float_t dcaTPC[2]={0,0};
1016   // Float_t covTPC[3]={0,0,0};
1017
1018   UInt_t specie = event->GetEventSpecie();  // skip laser events
1019   if (specie==AliRecoParam::kCalib) return;
1020   
1021
1022
1023   for (Int_t itrack0=0;itrack0<ntracks;itrack0++) {
1024     AliVTrack *track0 = event->GetVTrack(itrack0);
1025     if (!track0) continue;
1026     if (!track0->IsOn(AliVTrack::kTPCrefit)) continue;
1027
1028     if (TMath::Abs(AliTracker::GetBz())>1&&track0->Pt()<kMinPt) continue;
1029     if (track0->GetTPCncls()<kMinNcl) continue;
1030     if (TMath::Abs(track0->GetY())<kMaxDelta[0]) continue; 
1031     if (track0->GetKinkIndex(0)>0) continue;
1032     const Double_t * par0=track0->GetParameter(); //track param at rhe DCA
1033     //rm primaries
1034     //
1035     //track0->GetImpactParametersTPC(dcaTPC,covTPC);
1036     //if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
1037     //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
1038     //    const AliExternalTrackParam * trackIn0 = track0->GetInnerParam();
1039     for (Int_t itrack1=itrack0+1;itrack1<ntracks;itrack1++) {
1040       AliVTrack *track1 = event->GetVTrack(itrack1);
1041       if (!track1) continue;  
1042       if (!track1->IsOn(AliVTrack::kTPCrefit)) continue;
1043       if (track1->GetKinkIndex(0)>0) continue;
1044       if (TMath::Abs(AliTracker::GetBz())>1&&track1->Pt()<kMinPt) continue;
1045       if (track1->GetTPCncls()<kMinNcl) continue;
1046       if (TMath::Abs(AliTracker::GetBz())>1&&TMath::Max(track1->Pt(), track0->Pt())<kMinPtMax) continue;
1047       if (TMath::Abs(track1->GetY())<kMaxDelta[0]) continue;
1048       //track1->GetImpactParametersTPC(dcaTPC,covTPC);
1049       //      if (TMath::Abs(dcaTPC[0])<kMaxDelta[0]) continue;
1050       //if (TMath::Abs(dcaTPC[1])<kMaxDelta[0]*2) continue;
1051       //
1052       const Double_t* par1=track1->GetParameter(); //track param at rhe DCA
1053       //
1054       Bool_t isPair=kTRUE;
1055       for (Int_t ipar=0; ipar<5; ipar++){
1056         if (ipar==4&&TMath::Abs(AliTracker::GetBz())<1) continue; // 1/pt not defined for B field off
1057         if (TMath::Abs(TMath::Abs(par0[ipar])-TMath::Abs(par1[ipar]))>kMaxDelta[ipar]) isPair=kFALSE;
1058       }
1059       if (!isPair) continue;
1060       if (TMath::Abs(TMath::Abs(track0->GetAlpha()-track1->GetAlpha())-TMath::Pi())>kMaxDelta[2]) isPair=kFALSE;
1061       //delta with correct sign
1062       /*
1063         TCut cut0="abs(t1.fP[0]+t0.fP[0])<2"
1064         TCut cut3="abs(t1.fP[3]+t0.fP[3])<0.02"
1065         TCut cut4="abs(t1.fP[4]+t0.fP[4])<0.2"
1066       */
1067       if  (TMath::Abs(par0[0]+par1[0])>kMaxDelta[0]) isPair=kFALSE; //delta y   opposite sign
1068       if  (TMath::Abs(par0[3]+par1[3])>kMaxDelta[3]) isPair=kFALSE; //delta tgl opposite sign
1069       if  (TMath::Abs(AliTracker::GetBz())>1 && TMath::Abs(par0[4]+par1[4])>kMaxDelta[4]) isPair=kFALSE; //delta 1/pt opposite sign
1070       if (!isPair) continue;
1071       TString filename(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1072       Int_t eventNumber = event->GetEventNumberInFile(); 
1073       Bool_t hasFriend=(friendEvent) ? (friendEvent->GetTrack(itrack0)!=0):0;
1074       Bool_t hasITS=(track0->GetNcls(0)+track1->GetNcls(0)>4);
1075       printf("DUMPHPTCosmic:%s|%f|%d|%d|%d\n",filename.Data(),(TMath::Min(track0->Pt(),track1->Pt())), eventNumber,hasFriend,hasITS);
1076
1077
1078       //      const AliExternalTrackParam * trackIn1 = track1->GetInnerParam();      
1079       //
1080       //       
1081       TTreeSRedirector * pcstream =  GetDebugStreamer();
1082       Int_t ntracksSPD = vertexSPD->GetNContributors();
1083       Int_t ntracksTPC = vertexTPC->GetNContributors();
1084       //
1085       const AliVfriendTrack *friendTrack0 = friendEvent->GetTrack(itrack0);
1086       if (!friendTrack0) continue;
1087       const AliVfriendTrack *friendTrack1 = friendEvent->GetTrack(itrack1);
1088       if (!friendTrack1) continue;
1089       TObject *calibObject;
1090       AliTPCseed *seed0 = 0;   
1091       AliTPCseed *seed1 = 0;
1092       //
1093       for (Int_t l=0;(calibObject=friendTrack0->GetCalibObject(l));++l) {
1094         if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
1095       }
1096       for (Int_t l=0;(calibObject=friendTrack1->GetCalibObject(l));++l) {
1097         if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
1098       }
1099
1100       //
1101       if (pcstream){
1102         (*pcstream)<<"pairs"<<
1103           "run="<<fRun<<              //  run number
1104           "event="<<fEvent<<          //  event number
1105           "time="<<fTime<<            //  time stamp of event
1106           "trigger="<<fTrigger<<      //  trigger
1107           "triggerClass="<<&fTriggerClass<<      //  trigger
1108           "mag="<<fMagF<<             //  magnetic field
1109           //
1110           "nSPD="<<ntracksSPD<<
1111           "nTPC="<<ntracksTPC<<
1112           "vSPD.="<<vertexSPD<<         //primary vertex -SPD
1113           "vTPC.="<<vertexTPC<<         //primary vertex -TPC
1114           "t0.="<<track0<<              //track0
1115           "t1.="<<track1<<              //track1
1116       //"ft0.="<<dummyfriendTrack0<<       //track0
1117       //"ft1.="<<dummyfriendTrack1<<       //track1
1118           "s0.="<<seed0<<               //track0
1119           "s1.="<<seed1<<               //track1
1120           "\n";      
1121       }
1122
1123       //**********************TEMPORARY!!*******************************************
1124       // more investigation is needed with Tree ///!!!
1125       //all dummy stuff here is just for code to compile and work with ESD
1126
1127       AliESDfriendTrack *dummyfriendTrack0 = (AliESDfriendTrack*)friendTrack0;
1128       AliESDfriendTrack *dummyfriendTrack1 = (AliESDfriendTrack*)friendTrack1;
1129
1130       AliESDtrack *dummytrack0 = (AliESDtrack*)track0;
1131       AliESDtrack *dummytrack1 = (AliESDtrack*)track1;
1132
1133       if ((pcstream)&&(dummyfriendTrack0)){
1134           (*pcstream)<<"ft0.="<<dummyfriendTrack0<<"\n";
1135       }
1136       if ((pcstream)&&(dummyfriendTrack1)){
1137           (*pcstream)<<"ft1.="<<dummyfriendTrack1<<"\n";
1138       }
1139
1140       if (!fCosmicTree) {
1141         fCosmicTree = new TTree("pairs","pairs");
1142         fCosmicTree->SetDirectory(0);
1143       }
1144       if (fCosmicTree->GetEntries()==0){
1145         //
1146         fCosmicTree->SetDirectory(0);
1147     fCosmicTree->Branch("t0.",&dummytrack0);
1148     fCosmicTree->Branch("t1.",&dummytrack1);
1149     fCosmicTree->Branch("ft0.",&dummyfriendTrack0);
1150     fCosmicTree->Branch("ft1.",&dummyfriendTrack1);
1151       }else{
1152     fCosmicTree->SetBranchAddress("t0.",&dummytrack0);
1153     fCosmicTree->SetBranchAddress("t1.",&dummytrack1);
1154     fCosmicTree->SetBranchAddress("ft0.",&dummyfriendTrack0);
1155     fCosmicTree->SetBranchAddress("ft1.",&dummyfriendTrack1);
1156       }
1157       fCosmicTree->Fill();
1158     }
1159   }
1160 }
1161
1162
1163 void  AliTPCcalibCosmic::Terminate(){
1164   //
1165   // copy the cosmic tree to memory resident tree
1166   //
1167   static Int_t counter=0;
1168   printf("AliTPCcalibCosmic::Terminate\t%d\n",counter);
1169   counter++;
1170   AliTPCcalibBase::Terminate();
1171 }
1172
1173
1174 void AliTPCcalibCosmic::AddTree(TTree* treeOutput, TTree * treeInput){
1175   //
1176   // Add the content of tree: 
1177   // Notice automatic copy of tree in ROOT does not work for such complicated tree
1178   //  
1179   return;
1180   //if (TMath::Abs(fMagF)<0.1) return; // work around - otherwise crashes 
1181   AliESDtrack *track0=new AliESDtrack;     ///!!!
1182   AliESDtrack *track1=new AliESDtrack;
1183   AliESDfriendTrack *ftrack0=new AliESDfriendTrack;
1184   AliESDfriendTrack *ftrack1=new AliESDfriendTrack;
1185   treeInput->SetBranchAddress("t0.",&track0);   
1186   treeInput->SetBranchAddress("t1.",&track1);
1187   treeInput->SetBranchAddress("ft0.",&ftrack0); 
1188   treeInput->SetBranchAddress("ft1.",&ftrack1);
1189   treeOutput->SetDirectory(0);
1190   //
1191   Int_t entries= treeInput->GetEntries();
1192   Int_t step=1+Int_t(TMath::Log(1+treeOutput->GetEntries()/10.));
1193   for (Int_t i=0; i<entries; i+=step){
1194     treeInput->SetBranchAddress("t0.",&track0); 
1195     treeInput->SetBranchAddress("t1.",&track1);
1196     treeInput->SetBranchAddress("ft0.",&ftrack0);       
1197     treeInput->SetBranchAddress("ft1.",&ftrack1);
1198     treeInput->GetEntry(i);
1199     if (!track0) continue;
1200     if (!track1) continue;
1201     if (!ftrack0) continue;
1202     if (!ftrack1) continue;
1203     if (track0->GetTPCncls()<=0) continue;
1204     if (track1->GetTPCncls()<=0) continue;
1205     if (!track0->GetInnerParam()) continue;
1206     if (!track1->GetInnerParam()) continue;
1207     if (!track0->GetTPCInnerParam()) continue;
1208     if (!track1->GetTPCInnerParam()) continue;
1209     //track0
1210     treeOutput->SetBranchAddress("t0.",&track0);        
1211     treeOutput->SetBranchAddress("t1.",&track1);
1212     treeOutput->SetBranchAddress("ft0.",&ftrack0);      
1213     treeOutput->SetBranchAddress("ft1.",&ftrack1);    
1214     treeOutput->Fill();
1215     delete track0;
1216     delete track1;
1217     delete ftrack0;
1218     delete ftrack1;
1219     track0=0;
1220     track1=0;
1221     ftrack0=0;
1222     ftrack1=0;
1223   }
1224 }
1225
1226
1227
1228 void AliTPCcalibCosmic::MakeFitTree(TTree * treeInput, TTreeSRedirector *pcstream, const TObjArray * corrArray, Int_t step, Int_t run){
1229   //
1230   // Make fit tree
1231   // refit the tracks with original points + corrected points for each correction
1232   // Input:
1233   //   treeInput - tree with cosmic tracks
1234   //   pcstream  - debug output
1235
1236   // Algorithm:
1237   // Loop over pair of cosmic tracks:
1238   //   1. Find trigger offset between cosmic event and "physic" trigger
1239   //   2. Refit tracks with current transformation
1240   //   3. Refit tracks using additional "primitive" distortion on top
1241   // Best correction estimated as linear combination of corrections 
1242   // minimizing the observed distortions
1243   // Observed distortions - matching in the y,z, snp, theta and 1/pt
1244   //
1245   const Double_t kResetCov=20.;
1246   const Double_t kMaxDelta[5]={1,40,0.03,0.01,0.2};
1247   const Double_t kSigma=2.;    
1248   const Double_t kMaxTime=1050;
1249   const Double_t kMaxSnp=0.8;
1250   Int_t ncorr=corrArray->GetEntries();
1251   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1252   AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
1253   AliGRPObject*  grp = AliTPCcalibDB::Instance()->GetGRP(run);
1254   Double_t time=0.5*(grp->GetTimeStart() +grp->GetTimeEnd()); 
1255   transform->SetCurrentRun(run);
1256   transform->SetCurrentTimeStamp(TMath::Nint(time));
1257   Double_t covar[15];
1258   for (Int_t i=0;i<15;i++) covar[i]=0;
1259   covar[0]=kSigma*kSigma;
1260   covar[2]=kSigma*kSigma;
1261   covar[5]=kSigma*kSigma/Float_t(150*150);
1262   covar[9]=kSigma*kSigma/Float_t(150*150);
1263   covar[14]=0.2*0.2;
1264   Double_t *distortions = new Double_t[ncorr+1];
1265
1266   AliESDtrack *track0=new AliESDtrack;
1267   AliESDtrack *track1=new AliESDtrack;
1268   AliESDfriendTrack *ftrack0=new AliESDfriendTrack;
1269   AliESDfriendTrack *ftrack1=new AliESDfriendTrack;
1270   treeInput->SetBranchAddress("t0.",&track0);   
1271   treeInput->SetBranchAddress("t1.",&track1);
1272   treeInput->SetBranchAddress("ft0.",&ftrack0); 
1273   treeInput->SetBranchAddress("ft1.",&ftrack1);
1274   Int_t entries= treeInput->GetEntries();
1275   for (Int_t i=0; i<entries; i+=step){    
1276     treeInput->GetEntry(i);
1277     if (i%20==0) printf("%d\n",i);
1278     if (!ftrack0->GetTPCOut()) continue;
1279     if (!ftrack1->GetTPCOut()) continue;
1280     AliTPCseed *seed0=0;
1281     AliTPCseed *seed1=0;
1282     TObject *calibObject;
1283     for (Int_t l=0;(calibObject=ftrack0->GetCalibObject(l));++l) {
1284       if ((seed0=dynamic_cast<AliTPCseed*>(calibObject))) break;
1285     }
1286     for (Int_t l=0;(calibObject=ftrack1->GetCalibObject(l));++l) {
1287       if ((seed1=dynamic_cast<AliTPCseed*>(calibObject))) break;
1288     }
1289     if (!seed0) continue;
1290     if (!seed1) continue;
1291     if (TMath::Abs(seed0->GetSnp())>kMaxSnp) continue;
1292     if (TMath::Abs(seed1->GetSnp())>kMaxSnp) continue;
1293     //
1294     //
1295     Int_t nclA0=0, nclC0=0;     // number of clusters
1296     Int_t nclA1=0, nclC1=0;     // number of clusters
1297     Int_t ncl0=0,ncl1=0;
1298     Double_t rmin0=300, rmax0=-300;  // variables to estimate the time0 - trigger offset
1299     Double_t rmin1=300, rmax1=-300;
1300     Double_t tmin0=2000, tmax0=-2000;
1301     Double_t tmin1=2000, tmax1=-2000;
1302     //
1303     //
1304     // calculate trigger offset usig "missing clusters"
1305     for (Int_t irow=0; irow<159; irow++){
1306       AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
1307       if (cluster0 &&cluster0->GetX()>10){
1308         if (cluster0->GetX()<rmin0) { rmin0=cluster0->GetX(); tmin0=cluster0->GetTimeBin();}
1309         if (cluster0->GetX()>rmax0) { rmax0=cluster0->GetX(); tmax0=cluster0->GetTimeBin();}
1310         ncl0++;
1311         if (cluster0->GetDetector()%36<18) nclA0++;
1312         if (cluster0->GetDetector()%36>=18) nclC0++;
1313       }  
1314       AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
1315       if (cluster1&&cluster1->GetX()>10){
1316         if (cluster1->GetX()<rmin1) { rmin1=cluster1->GetX();  tmin1=cluster1->GetTimeBin();}
1317         if (cluster1->GetX()>rmax1) { rmax1=cluster1->GetX(); tmax1=cluster1->GetTimeBin();}
1318         ncl1++;
1319         if (cluster1->GetDetector()%36<18) nclA1++;
1320         if (cluster1->GetDetector()%36>=18) nclC1++;
1321       }
1322     }
1323     Int_t cosmicType=0;  // types of cosmic topology
1324     if ((nclA0>nclC0) && (nclA1>nclC1)) cosmicType=0; // AA side
1325     if ((nclA0<nclC0) && (nclA1<nclC1)) cosmicType=1; // CC side
1326     if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=2; // AC side
1327     if ((nclA0<nclC0) && (nclA1>nclC1)) cosmicType=3; // CA side
1328     //if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=6; // AC side out of time
1329     //if ((nclA0>nclC0) && (nclA1<nclC1)) cosmicType=7; // CA side out of time
1330     //
1331     Double_t deltaTime = 0;   // correction for trigger not in time - equalizing the time arival
1332     if ((cosmicType>1)&&TMath::Abs(track1->GetZ()-track0->GetZ())>4){
1333       cosmicType+=4;
1334       deltaTime=0.5*(track1->GetZ()-track0->GetZ())/param->GetZWidth();
1335       if (nclA0>nclC0) deltaTime*=-1; // if A side track
1336     }
1337     //
1338     TVectorD vectorDT(8);
1339     Int_t crossCounter=0;
1340     Double_t deltaTimeCross = AliTPCcalibCosmic::GetDeltaTime(rmin0, rmax0, rmin1, rmax1, tmin0, tmax0, tmin1, tmax1, TMath::Abs(track0->GetY()),vectorDT);
1341     Bool_t isOKTrigger=kTRUE;
1342     for (Int_t ic=0; ic<6;ic++) {
1343       if (TMath::Abs(vectorDT[ic])>0) {
1344         if (vectorDT[ic]+vectorDT[6]<0) isOKTrigger=kFALSE;
1345         if (vectorDT[ic]+vectorDT[7]>kMaxTime) isOKTrigger=kFALSE;
1346         if (isOKTrigger){
1347           crossCounter++; 
1348         }
1349       }
1350     }
1351     Double_t deltaTimeCluster=deltaTime;
1352     if ((cosmicType==0 || cosmicType==1) && crossCounter>0){
1353       deltaTimeCluster=deltaTimeCross;
1354       cosmicType+=8;
1355     }
1356     if (nclA0*nclC0>0 || nclA1*nclC1>0) cosmicType+=16;  // mixed A side C side - bad for visualization
1357     //
1358     // Apply current transformation
1359     //
1360     //
1361     for (Int_t irow=0; irow<159; irow++){
1362       AliTPCclusterMI *cluster0=seed0->GetClusterPointer(irow);
1363       if (cluster0 &&cluster0->GetX()>10){
1364         Double_t x0[3]={ static_cast<Double_t>(cluster0->GetRow()),cluster0->GetPad(),cluster0->GetTimeBin()+deltaTimeCluster};
1365         Int_t index0[1]={cluster0->GetDetector()};
1366         transform->Transform(x0,index0,0,1);  
1367         cluster0->SetX(x0[0]);
1368         cluster0->SetY(x0[1]);
1369         cluster0->SetZ(x0[2]);
1370         //
1371       }
1372       AliTPCclusterMI *cluster1=seed1->GetClusterPointer(irow);
1373       if (cluster1&&cluster1->GetX()>10){
1374         Double_t x1[3]={ static_cast<Double_t>(cluster1->GetRow()),cluster1->GetPad(),cluster1->GetTimeBin()+deltaTimeCluster};
1375         Int_t index1[1]={cluster1->GetDetector()};
1376         transform->Transform(x1,index1,0,1);  
1377         cluster1->SetX(x1[0]);
1378         cluster1->SetY(x1[1]);
1379         cluster1->SetZ(x1[2]);
1380       }
1381     }
1382     //
1383     //
1384     Double_t alpha=track0->GetAlpha();   // rotation frame
1385     Double_t cos = TMath::Cos(alpha);
1386     Double_t sin = TMath::Sin(alpha);
1387     Double_t mass =  TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
1388     AliExternalTrackParam  btrack0=*(ftrack0->GetTPCOut());
1389     AliExternalTrackParam  btrack1=*(ftrack1->GetTPCOut());
1390     btrack0.Rotate(alpha);
1391     btrack1.Rotate(alpha);
1392     // change the sign for track 1
1393     Double_t* par1=(Double_t*)btrack0.GetParameter();
1394     par1[3]*=-1;
1395     par1[4]*=-1;
1396     btrack0.AddCovariance(covar);
1397     btrack1.AddCovariance(covar);
1398     btrack0.ResetCovariance(kResetCov);
1399     btrack1.ResetCovariance(kResetCov);
1400     Bool_t isOK=kTRUE;
1401     Bool_t isOKT=kTRUE;
1402     TObjArray tracks0(ncorr+1);
1403     TObjArray tracks1(ncorr+1);
1404     //    
1405     Double_t dEdx0Tot=seed0->CookdEdxAnalytical(0.02,0.6,kTRUE);
1406     Double_t dEdx1Tot=seed1->CookdEdxAnalytical(0.02,0.6,kTRUE);
1407     Double_t dEdx0Max=seed0->CookdEdxAnalytical(0.02,0.6,kFALSE);
1408     Double_t dEdx1Max=seed1->CookdEdxAnalytical(0.02,0.6,kFALSE);
1409     //if (TMath::Abs((dEdx0Max+1)/(dEdx0Tot+1)-1.)>0.1) isOK=kFALSE;
1410     //if (TMath::Abs((dEdx1Max+1)/(dEdx1Tot+1)-1.)>0.1) isOK=kFALSE;
1411     ncl0=0; ncl1=0;
1412     for (Int_t icorr=-1; icorr<ncorr; icorr++){
1413       AliExternalTrackParam  rtrack0=btrack0;
1414       AliExternalTrackParam  rtrack1=btrack1;
1415       AliTPCCorrection *corr = 0;
1416       if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
1417       //
1418       for (Int_t irow=159; irow>0; irow--){ 
1419         AliTPCclusterMI *cluster=seed0->GetClusterPointer(irow);
1420         if (!cluster) continue;
1421         if (!isOKT) break;
1422         Double_t rD[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1423         transform->RotatedGlobal2Global(cluster->GetDetector()%36,rD);  // transform to global
1424         Float_t  r[3]={static_cast<Float_t>(rD[0]),static_cast<Float_t>(rD[1]),static_cast<Float_t>(rD[2])};
1425         if (corr){
1426           corr->DistortPoint(r, cluster->GetDetector());
1427         }
1428         //
1429         Double_t cov[3]={0.01,0.,0.01}; 
1430         Double_t lx =cos*r[0]+sin*r[1];      
1431         Double_t ly =-sin*r[0]+cos*r[1];
1432         rD[1]=ly; rD[0]=lx; rD[2]=r[2];  //transform to track local
1433         if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, lx,mass,1.,kFALSE)) isOKT=kFALSE;;
1434         if (!rtrack0.Update(&rD[1],cov)) isOKT =kFALSE;
1435         if (icorr<0) ncl0++;
1436       }
1437       //
1438       for (Int_t irow=159; irow>0; irow--){ 
1439         AliTPCclusterMI *cluster=seed1->GetClusterPointer(irow);
1440         if (!cluster) continue;
1441         if (!isOKT) break;
1442         Double_t rD[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1443         transform->RotatedGlobal2Global(cluster->GetDetector()%36,rD);
1444         Float_t  r[3]={static_cast<Float_t>(rD[0]),static_cast<Float_t>(rD[1]),static_cast<Float_t>(rD[2])};
1445         if (corr){
1446           corr->DistortPoint(r, cluster->GetDetector());
1447         }
1448         //
1449         Double_t cov[3]={0.01,0.,0.01}; 
1450         Double_t lx =cos*r[0]+sin*r[1];      
1451         Double_t ly =-sin*r[0]+cos*r[1];
1452         rD[1]=ly; rD[0]=lx; rD[2]=r[2];
1453         if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, lx,mass,1.,kFALSE)) isOKT=kFALSE;
1454         if (!rtrack1.Update(&rD[1],cov)) isOKT=kFALSE;
1455         if (icorr<0) ncl1++;
1456       }
1457       if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, 0,mass,10.,kFALSE)) isOKT=kFALSE;
1458       if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, 0,mass,10.,kFALSE)) isOKT=kFALSE;
1459       if (!AliTracker::PropagateTrackToBxByBz(&rtrack0, 0,mass,1.,kFALSE))  isOKT=kFALSE;
1460       if (!AliTracker::PropagateTrackToBxByBz(&rtrack1, 0,mass,1.,kFALSE))  isOKT=kFALSE;
1461       const Double_t *param0=rtrack0.GetParameter();
1462       const Double_t *param1=rtrack1.GetParameter();
1463       for (Int_t ipar=0; ipar<4; ipar++){
1464         if (TMath::Abs(param1[ipar]-param0[ipar])>kMaxDelta[ipar]) isOK=kFALSE;
1465       }
1466       tracks0.AddAt(rtrack0.Clone(), icorr+1);
1467       tracks1.AddAt(rtrack1.Clone(), icorr+1);
1468       AliExternalTrackParam out0=*(ftrack0->GetTPCOut());
1469       AliExternalTrackParam out1=*(ftrack1->GetTPCOut());
1470       Int_t nentries=TMath::Min(ncl0,ncl1);
1471
1472       if (icorr<0) {
1473         (*pcstream)<<"cosmic"<<
1474           "isOK="<<isOK<<              // correct all propagation update and also residuals
1475           "isOKT="<<isOKT<<            // correct all propagation update
1476           "isOKTrigger="<<isOKTrigger<< // correct? estimate of trigger offset
1477           "id="<<cosmicType<<
1478           //
1479           //
1480           "cross="<<crossCounter<<
1481           "vDT.="<<&vectorDT<<
1482           //
1483           "dTime="<<deltaTime<<        // delta time using the A-c side cross
1484           "dTimeCross="<<deltaTimeCross<< // delta time using missing clusters
1485           //
1486           "dEdx0Max="<<dEdx0Max<<
1487           "dEdx0Tot="<<dEdx0Tot<<
1488           "dEdx1Max="<<dEdx1Max<<
1489           "dEdx1Tot="<<dEdx1Tot<<
1490           //
1491           "track0.="<<track0<<         // original track 0
1492           "track1.="<<track1<<         // original track 1
1493           "out0.="<<&out0<<             // outer track  0
1494           "out1.="<<&out1<<             // outer track  1
1495           "rtrack0.="<<&rtrack0<<      // refitted track with current transform
1496           "rtrack1.="<<&rtrack1<<     //          
1497           "ncl0="<<ncl0<<
1498           "ncl1="<<ncl1<<
1499           "entries="<<nentries<<       // number of clusters
1500           "\n";
1501       }
1502     }
1503     //
1504
1505     if (isOK){        
1506       Int_t nentries=TMath::Min(ncl0,ncl1);    
1507       for (Int_t ipar=0; ipar<5; ipar++){
1508         for (Int_t icorr=-1; icorr<ncorr; icorr++){
1509           AliTPCCorrection *corr = 0;
1510           if (icorr>=0) corr = (AliTPCCorrection*)corrArray->At(icorr);
1511           //
1512           AliExternalTrackParam *param0=(AliExternalTrackParam *) tracks0.At(icorr+1);
1513           AliExternalTrackParam *param1=(AliExternalTrackParam *) tracks1.At(icorr+1);
1514           distortions[icorr+1]=param1->GetParameter()[ipar]-param0->GetParameter()[ipar];
1515           if (icorr>=0){
1516             distortions[icorr+1]-=distortions[0];
1517           }
1518           //
1519           if (icorr<0){
1520             Double_t bz=AliTrackerBase::GetBz();
1521             Double_t gxyz[3];
1522             param0->GetXYZ(gxyz);
1523             Int_t dtype=20;
1524             Double_t theta=param0->GetParameter()[3];
1525             Double_t phi = alpha;
1526             Double_t snp = track0->GetInnerParam()->GetSnp();
1527             Double_t mean= distortions[0];
1528             Int_t index = param0->GetIndex(ipar,ipar);
1529             Double_t rms=TMath::Sqrt(param1->GetCovariance()[index]+param1->GetCovariance()[index]);
1530             if (crossCounter<1) rms*=1;
1531             Double_t sector=9*phi/TMath::Pi();
1532             Double_t dsec   = sector-TMath::Nint(sector+0.5);
1533             Double_t gx=gxyz[0],gy=gxyz[1],gz=gxyz[2];
1534             Double_t refX=TMath::Sqrt(gx*gx+gy*gy);
1535             Double_t dRrec=0;
1536             //      Double_t pt=(param0->GetSigned1Pt()+param1->GetSigned1Pt())*0.5;
1537             Double_t pt=(param0->GetSigned1Pt()+param1->GetSigned1Pt())*0.5;
1538
1539             (*pcstream)<<"fit"<<  // dump valus for fit
1540               "run="<<run<<       //run number
1541               "bz="<<bz<<         // magnetic filed used
1542               "dtype="<<dtype<<   // detector match type 20
1543               "ptype="<<ipar<<   // parameter type
1544               "theta="<<theta<<   // theta
1545               "phi="<<phi<<       // phi 
1546               "snp="<<snp<<       // snp
1547               "mean="<<mean<<     // mean dist value
1548               "rms="<<rms<<       // rms
1549               "sector="<<sector<<
1550               "dsec="<<dsec<<
1551               //
1552               "refX="<<refX<<      // reference radius
1553               "gx="<<gx<<         // global position
1554               "gy="<<gy<<         // global position
1555               "gz="<<gz<<         // global position
1556               "dRrec="<<dRrec<<      // delta Radius in reconstruction
1557               "pt="<<pt<<            //1/pt
1558               "id="<<cosmicType<<     //type of the cosmic used      
1559               "entries="<<nentries;// number of entries in bin
1560             (*pcstream)<<"cosmicDebug"<<
1561               "p0.="<<param0<<    // dump distorted track 0
1562               "p1.="<<param1;    // dump distorted track 1
1563           }
1564           if (icorr>=0){
1565             (*pcstream)<<"fit"<<
1566               Form("%s=",corr->GetName())<<distortions[icorr+1];    // dump correction value
1567             (*pcstream)<<"cosmicDebug"<<
1568               Form("%s=",corr->GetName())<<distortions[icorr+1]<<    // dump correction value
1569               Form("%sp0.=",corr->GetName())<<param0<<    // dump distorted track 0
1570               Form("%sp1.=",corr->GetName())<<param1;    // dump distorted track 1
1571           }
1572         } //loop corrections      
1573         (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n";
1574         (*pcstream)<<"cosmicDebug"<<"isOK="<<isOK<<"\n";
1575       } //loop over parameters
1576     } // dump results
1577   }//loop tracks
1578   delete [] distortions;
1579 }
1580
1581
1582
1583 Double_t AliTPCcalibCosmic::GetDeltaTime(Double_t rmin0, Double_t rmax0, Double_t rmin1, Double_t rmax1, Double_t tmin0, Double_t tmax0, Double_t tmin1, Double_t tmax1, Double_t dcaR, TVectorD &vectorDT)
1584 {
1585   //
1586   // Estimate trigger offset between random cosmic event and "physics" trigger
1587   // Efficiency about 50 % of cases:
1588   // Cases:
1589   // 0. Tracks crossing A side and C side - no match in z - 30 % of cases
1590   // 1. Track only on one side and  dissapear at small or at high radiuses - 50 % of cases
1591   //    1.a) rmax<Rc    && tmax<Tcmax && tmax>tmin    => deltaT=Tcmax-tmax 
1592   //    1.b) rmin>Rcmin && tmin<Tcmax && tmin>tmax    => deltaT=Tcmax-tmin  
1593   //                      // case the z matching gives proper time
1594   //    1.c) rmax<Rc    && tmax>Tcmin && tmax<tmin    => deltaT=-tmax
1595   //
1596   // check algorithm:
1597   //    TCut cutStop = "(min(rmax0,rmax1)<235||abs(rmin0-rmin1)>10)"; // tracks not registered for full time
1598
1599   // Combinations:
1600   // 0-1 - forbidden
1601   // 0-2 - forbidden
1602   // 0-3 - occur - wrong correlation
1603   // 1-2 - occur - wrong correlation
1604   // 1-3 - forbidden
1605   // 2-3 - occur - small number of outlyers -20%
1606   // Frequency:
1607   //   0 - 106
1608   //   1 - 265
1609   //   2 - 206
1610   //   3 - 367
1611   //
1612   const Double_t kMaxRCut=235;  // max radius
1613   const Double_t kMinRCut=TMath::Max(dcaR,90.);  // min radius
1614   const Double_t kMaxDCut=30;   // max distance for minimal radius
1615   const Double_t kMinTime=110;
1616   const Double_t kMaxTime=950;  
1617   Double_t deltaT=0;
1618   Int_t counter=0;
1619   vectorDT[6]=TMath::Min(TMath::Min(tmin0,tmin1),TMath::Min(tmax0,tmax1));
1620   vectorDT[7]=TMath::Max(TMath::Max(tmin0,tmin1),TMath::Max(tmax0,tmax1));
1621   if (TMath::Min(rmax0,rmax1)<kMaxRCut){
1622     // max cross - deltaT>0
1623     if (rmax0<kMaxRCut && tmax0 <kMaxTime && tmax0>tmin0) vectorDT[0]=kMaxTime-tmax0; // disapear at CE
1624     if (rmax1<kMaxRCut && tmax1 <kMaxTime && tmax1>tmin1) vectorDT[1]=kMaxTime-tmax1; // disapear at CE
1625     // min cross - deltaT<0 - OK they are correlated
1626     if (rmax0<kMaxRCut && tmax0 >kMinTime && tmax0<tmin0) vectorDT[2]=-tmax0;         // disapear at ROC
1627     if (rmax1<kMaxRCut && tmax1 >kMinTime && tmax1<tmin1) vectorDT[3]=-tmax1;         // disapear at ROC
1628   }  
1629   if (rmin0> kMinRCut+kMaxDCut && tmin0 <kMaxTime && tmin0>tmax0) vectorDT[4]=kMaxTime-tmin0;
1630   if (rmin1> kMinRCut+kMaxDCut && tmin1 <kMaxTime && tmin1>tmax1) vectorDT[5]=kMaxTime-tmin1;
1631   Bool_t isOK=kTRUE;
1632   for (Int_t i=0; i<6;i++) {
1633     if (TMath::Abs(vectorDT[i])>0) {
1634       counter++; 
1635       if (vectorDT[i]+vectorDT[6]<0) isOK=kFALSE;
1636       if (vectorDT[i]+vectorDT[7]>kMaxTime) isOK=kFALSE;
1637       if (isOK) deltaT=vectorDT[i];
1638     }
1639   }
1640   return deltaT;  
1641 }
1642