Adding visualization of the tree aliases
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibTime.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
19 1. What do we calibrate.
20
21   Time dependence of gain and drift velocity in order to account for changes in: temperature, pressure, gas composition.
22
23   AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime",0, 1213.9e+06, 1213.96e+06, 0.04e+04, 0.04e+04);
24
25 2. How to interpret results
26
27 3. Simple example
28
29   a) determine the required time range:
30
31   AliXRDPROOFtoolkit tool;
32   TChain * chain = tool.MakeChain("pass2.txt","esdTree",0,6000);
33   chain->Draw("GetTimeStamp()")
34
35   b) analyse calibration object on Proof in calibration train 
36
37   AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime", StartTimeStamp, EndTimeStamp, IntegrationTimeVdrift);
38
39   c) plot results
40   .x ~/NimStyle.C
41   gSystem->Load("libANALYSIS");
42   gSystem->Load("libTPCcalib");
43
44   TFile f("CalibObjectsTrain1.root");
45   AliTPCcalibTime *calib = (AliTPCcalibTime *)f->Get("calibTime");
46   calib->GetHistoDrift("all")->Projection(2,0)->Draw()
47   calib->GetFitDrift("all")->Draw("lp")
48
49 4. Analysis using debug streamers.    
50
51   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
52   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
53   AliXRDPROOFtoolkit tool;
54   AliXRDPROOFtoolkit::FilterList("timeitstpc.txt","* itstpc",1) 
55   AliXRDPROOFtoolkit::FilterList("timetoftpc.txt","* pointMatch",1) 
56   AliXRDPROOFtoolkit::FilterList("time.txt","* trackInfo",1) 
57   AliXRDPROOFtoolkit::FilterList("timelaser.txt","* laserInfo",1) 
58
59   TChain * chainTPCITS = tool.MakeChainRandom("timeitstpc.txt.Good","itstpc",0,10000); 
60   TChain * chainTPCTOF = tool.MakeChainRandom("timetoftpc.txt","pointMatch",0,500); 
61   TChain * chainTime = tool.MakeChainRandom("time.txt.Good","trackInfo",0,10000);
62   TChain * chainLaser = tool.MakeChainRandom("timelaser.txt.Good","laserInfo",0,10000);
63   chainTime->Lookup();
64   chainLaser->Lookup();
65 */
66
67 #include "Riostream.h"
68 #include "TChain.h"
69 #include "TTree.h"
70 #include "TH1F.h"
71 #include "TH2F.h"
72 #include "TH3F.h"
73 #include "THnSparse.h"
74 #include "TList.h"
75 #include "TMath.h"
76 #include "TCanvas.h"
77 #include "TFile.h"
78 #include "TF1.h"
79 #include "TVectorD.h"
80 #include "TProfile.h"
81 #include "TGraphErrors.h"
82 #include "TCanvas.h"
83 #include "AliTPCclusterMI.h"
84 #include "AliTPCseed.h"
85 #include "AliESDVertex.h"
86 #include "AliESDEvent.h"
87 #include "AliESDfriend.h"
88 #include "AliESDInputHandler.h"
89 #include "AliAnalysisManager.h"
90
91 #include "AliTracker.h"
92 #include "AliMagF.h"
93 #include "AliTPCCalROC.h"
94
95 #include "AliLog.h"
96
97 #include "AliTPCcalibTime.h"
98 #include "AliRelAlignerKalman.h"
99
100 #include "TTreeStream.h"
101 #include "AliTPCTracklet.h"
102 #include "TTimeStamp.h"
103 #include "AliTPCcalibDB.h"
104 #include "AliTPCcalibLaser.h"
105 #include "AliDCSSensorArray.h"
106 #include "AliDCSSensor.h"
107
108 #include "TDatabasePDG.h"
109 #include "AliTrackPointArray.h"
110
111 ClassImp(AliTPCcalibTime)
112
113
114 AliTPCcalibTime::AliTPCcalibTime() 
115   :AliTPCcalibBase(), 
116    fLaser(0),       // pointer to laser calibration
117    fDz(0),          // current delta z
118    fCutMaxD(3),        // maximal distance in rfi ditection
119    fCutMaxDz(25),      // maximal distance in rfi ditection
120    fCutTheta(0.03),    // maximal distan theta
121    fCutMinDir(-0.99),  // direction vector products
122    fCutTracks(10),
123    fArrayDz(0),          //NEW! Tmap of V drifts for different triggers
124    fAlignITSTPC(0),      //alignemnt array ITS TPC match
125    fAlignTRDTPC(0),      //alignemnt array TRD TPC match 
126    fAlignTOFTPC(0),      //alignemnt array TOF TPC match
127    fTimeBins(0),
128    fTimeStart(0),
129    fTimeEnd(0),
130    fPtBins(0),
131    fPtStart(0),
132    fPtEnd(0),
133    fVdriftBins(0),
134    fVdriftStart(0),
135    fVdriftEnd(0),
136    fRunBins(0),
137    fRunStart(0),
138    fRunEnd(0)
139 //   fBinsVdrift(fTimeBins,fPtBins,fVdriftBins),
140 //   fXminVdrift(fTimeStart,fPtStart,fVdriftStart),
141 //   fXmaxVdrift(fTimeEnd,fPtEnd,fVdriftEnd)
142 {  
143   AliInfo("Default Constructor");  
144   for (Int_t i=0;i<3;i++) {
145     fHistVdriftLaserA[i]=0;
146     fHistVdriftLaserC[i]=0;
147   }
148   for (Int_t i=0;i<10;i++) {
149     fCosmiMatchingHisto[i]=0;
150   }
151 }
152
153 AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeVdrift)
154   :AliTPCcalibBase(),
155    fLaser(0),            // pointer to laser calibration
156    fDz(0),               // current delta z
157    fCutMaxD(5*0.5356),   // maximal distance in rfi ditection
158    fCutMaxDz(40),   // maximal distance in rfi ditection
159    fCutTheta(5*0.004644),// maximal distan theta
160    fCutMinDir(-0.99),    // direction vector products
161    fCutTracks(10),
162    fArrayDz(0),            //Tmap of V drifts for different triggers
163    fAlignITSTPC(0),      //alignemnt array ITS TPC match
164    fAlignTRDTPC(0),      //alignemnt array TRD TPC match 
165    fAlignTOFTPC(0),      //alignemnt array TOF TPC match
166    fTimeBins(0),
167    fTimeStart(0),
168    fTimeEnd(0),
169    fPtBins(0),
170    fPtStart(0),
171    fPtEnd(0),
172    fVdriftBins(0),
173    fVdriftStart(0),
174    fVdriftEnd(0),
175    fRunBins(0),
176    fRunStart(0),
177    fRunEnd(0)
178 {
179   SetName(name);
180   SetTitle(title);
181   for (Int_t i=0;i<3;i++) {
182     fHistVdriftLaserA[i]=0;
183     fHistVdriftLaserC[i]=0;
184   }
185
186   AliInfo("Non Default Constructor");
187   fTimeBins   =(EndTime-StartTime)/deltaIntegrationTimeVdrift;
188   fTimeStart  =StartTime; //(((TObjString*)(mapGRP->GetValue("fAliceStartTime")))->GetString()).Atoi();
189   fTimeEnd    =EndTime;   //(((TObjString*)(mapGRP->GetValue("fAliceStopTime")))->GetString()).Atoi();
190   fPtBins     = 400;
191   fPtStart    = -0.04;
192   fPtEnd      =  0.04;
193   fVdriftBins = 500;
194   fVdriftStart= -0.1;
195   fVdriftEnd  =  0.1;
196   fRunBins    = 100001;
197   fRunStart   = -1.5;
198   fRunEnd     = 99999.5;
199
200   Int_t    binsVdriftLaser[4] = {fTimeBins , fPtBins , fVdriftBins*20, fRunBins };
201   Double_t xminVdriftLaser[4] = {fTimeStart, fPtStart, fVdriftStart  , fRunStart};
202   Double_t xmaxVdriftLaser[4] = {fTimeEnd  , fPtEnd  , fVdriftEnd    , fRunEnd  };
203   TString axisTitle[4]={
204     "T",
205     "#delta_{P/T}",
206     "value",
207     "run"
208   };
209   TString histoName[3]={
210     "Loffset",
211     "Lcorr",
212     "Lgy"
213   };
214
215   
216   for (Int_t i=0;i<3;i++) {
217     fHistVdriftLaserA[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
218     fHistVdriftLaserC[i] = new THnSparseF("HistVdriftLaser","HistVdriftLaser;time;p/T ratio;Vdrift;run",4,binsVdriftLaser,xminVdriftLaser,xmaxVdriftLaser);
219     fHistVdriftLaserA[i]->SetName(histoName[i]);
220     fHistVdriftLaserC[i]->SetName(histoName[i]);
221     for (Int_t iaxis=0; iaxis<4;iaxis++){
222       fHistVdriftLaserA[i]->GetAxis(iaxis)->SetName(axisTitle[iaxis]);
223       fHistVdriftLaserC[i]->GetAxis(iaxis)->SetName(axisTitle[iaxis]);
224     }
225   }
226   fBinsVdrift[0] = fTimeBins;
227   fBinsVdrift[1] = fPtBins;
228   fBinsVdrift[2] = fVdriftBins;
229   fBinsVdrift[3] = fRunBins;
230   fXminVdrift[0] = fTimeStart;
231   fXminVdrift[1] = fPtStart;
232   fXminVdrift[2] = fVdriftStart;
233   fXminVdrift[3] = fRunStart;
234   fXmaxVdrift[0] = fTimeEnd;
235   fXmaxVdrift[1] = fPtEnd;
236   fXmaxVdrift[2] = fVdriftEnd;
237   fXmaxVdrift[3] = fRunEnd;
238
239   fArrayDz=new TObjArray();
240   fAlignITSTPC = new TObjArray;      //alignemnt array ITS TPC match
241   fAlignTRDTPC = new TObjArray;      //alignemnt array ITS TPC match
242   fAlignTOFTPC = new TObjArray;      //alignemnt array ITS TPC match
243
244   // fArrayDz->AddLast(fHistVdriftLaserA[0]);
245 //   fArrayDz->AddLast(fHistVdriftLaserA[1]);
246 //   fArrayDz->AddLast(fHistVdriftLaserA[2]);
247 //   fArrayDz->AddLast(fHistVdriftLaserC[0]);
248 //   fArrayDz->AddLast(fHistVdriftLaserC[1]);
249 //   fArrayDz->AddLast(fHistVdriftLaserC[2]);
250
251   fCosmiMatchingHisto[0]=new TH1F("Cosmics matching","p0-all"   ,100,-10*0.5356  ,10*0.5356  );
252   fCosmiMatchingHisto[1]=new TH1F("Cosmics matching","p1-all"   ,100,-10*4.541   ,10*4.541   );
253   fCosmiMatchingHisto[2]=new TH1F("Cosmics matching","p2-all"   ,100,-10*0.01134 ,10*0.01134 );
254   fCosmiMatchingHisto[3]=new TH1F("Cosmics matching","p3-all"   ,100,-10*0.004644,10*0.004644);
255   fCosmiMatchingHisto[4]=new TH1F("Cosmics matching","p4-all"   ,100,-10*0.03773 ,10*0.03773 );
256   fCosmiMatchingHisto[5]=new TH1F("Cosmics matching","p0-isPair",100,-10*0.5356  ,10*0.5356  );
257   fCosmiMatchingHisto[6]=new TH1F("Cosmics matching","p1-isPair",100,-10*4.541   ,10*4.541   );
258   fCosmiMatchingHisto[7]=new TH1F("Cosmics matching","p2-isPair",100,-10*0.01134 ,10*0.01134 );
259   fCosmiMatchingHisto[8]=new TH1F("Cosmics matching","p3-isPair",100,-10*0.004644,10*0.004644);
260   fCosmiMatchingHisto[9]=new TH1F("Cosmics matching","p4-isPair",100,-10*0.03773 ,10*0.03773 );
261 //  Char_t nameHisto[3]={'p','0','\n'};
262 //  for (Int_t i=0;i<10;i++){
263 //    fCosmiMatchingHisto[i]=new TH1F("Cosmics matching",nameHisto,8192,0,0);
264 //    nameHisto[1]++;
265 //    if(i==4) nameHisto[1]='0';
266 //  }
267 }
268
269 AliTPCcalibTime::~AliTPCcalibTime(){
270   //
271   // Destructor
272   //
273   for(Int_t i=0;i<3;i++){
274     if(fHistVdriftLaserA[i]){
275       delete fHistVdriftLaserA[i];
276       fHistVdriftLaserA[i]=NULL;
277     }
278     if(fHistVdriftLaserC[i]){
279       delete fHistVdriftLaserC[i];
280       fHistVdriftLaserC[i]=NULL;
281     }
282   }
283   if(fArrayDz){
284     fArrayDz->SetOwner();
285     fArrayDz->Delete();
286     delete fArrayDz;
287     fArrayDz=NULL;
288   }
289   for(Int_t i=0;i<5;i++){
290     if(fCosmiMatchingHisto[i]){
291       delete fCosmiMatchingHisto[i];
292       fCosmiMatchingHisto[i]=NULL;
293     }
294   }
295   fAlignITSTPC->Delete();
296   fAlignTRDTPC->Delete();
297   fAlignTOFTPC->Delete();
298 }
299
300 Bool_t AliTPCcalibTime::IsLaser(AliESDEvent */*event*/){
301   return kTRUE; //More accurate creteria to be added
302 }
303 Bool_t AliTPCcalibTime::IsCosmics(AliESDEvent */*event*/){
304   return kTRUE; //More accurate creteria to be added
305 }
306 Bool_t AliTPCcalibTime::IsBeam(AliESDEvent */*event*/){
307   return kTRUE; //More accurate creteria to be added
308 }
309 void AliTPCcalibTime::ResetCurrent(){
310   fDz=0; //Reset current dz
311 }
312 void AliTPCcalibTime::Process(AliESDEvent *event){
313   if(!event) return;
314   if (event->GetNumberOfTracks()<2) return;
315   ResetCurrent();
316   if(IsLaser  (event)) ProcessLaser (event);
317   if(IsCosmics(event)) ProcessCosmic(event);
318   if(IsBeam   (event)) ProcessBeam  (event);
319 }
320
321 void AliTPCcalibTime::ProcessLaser(AliESDEvent *event){
322   //
323   // Fit drift velocity using laser 
324   // 
325   // 0. cuts
326   const Int_t    kMinTracks     = 40;    // minimal number of laser tracks
327   const Int_t    kMinTracksSide = 20;    // minimal number of tracks per side
328   const Float_t  kMaxDeltaZ     = 30.;   // maximal trigger delay
329   const Float_t  kMaxDeltaV     = 0.05;  // maximal deltaV 
330   const Float_t  kMaxRMS        = 0.1;   // maximal RMS of tracks
331   //
332   /*
333     TCut cutRMS("sqrt(laserA.fElements[4])<0.1&&sqrt(laserC.fElements[4])<0.1");
334     TCut cutZ("abs(laserA.fElements[0]-laserC.fElements[0])<3");
335     TCut cutV("abs(laserA.fElements[1]-laserC.fElements[1])<0.01");
336     TCut cutY("abs(laserA.fElements[2]-laserC.fElements[2])<2");
337     TCut cutAll = cutRMS+cutZ+cutV+cutY;
338   */
339   if (event->GetNumberOfTracks()<kMinTracks) return;
340   //
341   if(!fLaser) fLaser = new AliTPCcalibLaser("laserTPC","laserTPC",kFALSE);
342   fLaser->Process(event);
343   if (fLaser->GetNtracks()<kMinTracks) return;   // small amount of tracks cut
344   if (fLaser->fFitAside->GetNrows()==0  && fLaser->fFitCside->GetNrows()==0) return;  // no fit neither a or C side
345   //
346   // debug streamer  - activate stream level
347   // Use it for tuning of the cuts
348   //
349   // cuts to be applied
350   //
351   Int_t isReject[2]={0,0};
352   //
353   // not enough tracks 
354   if (TMath::Abs((*fLaser->fFitAside)[3]) < kMinTracksSide) isReject[0]|=1; 
355   if (TMath::Abs((*fLaser->fFitCside)[3]) < kMinTracksSide) isReject[1]|=1; 
356   // unreasonable z offset
357   if (TMath::Abs((*fLaser->fFitAside)[0])>kMaxDeltaZ)  isReject[0]|=2;
358   if (TMath::Abs((*fLaser->fFitCside)[0])>kMaxDeltaZ)  isReject[1]|=2;
359   // unreasonable drift velocity
360   if (TMath::Abs((*fLaser->fFitAside)[1]-1)>kMaxDeltaV)  isReject[0]|=4;
361   if (TMath::Abs((*fLaser->fFitCside)[1]-1)>kMaxDeltaV)  isReject[1]|=4;
362   // big chi2
363   if (TMath::Sqrt(TMath::Abs((*fLaser->fFitAside)[4]))>kMaxRMS ) isReject[0]|=8;
364   if (TMath::Sqrt(TMath::Abs((*fLaser->fFitCside)[4]))>kMaxRMS ) isReject[1]|=8;
365
366
367
368
369   if (fStreamLevel>0){
370     printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
371
372     TTreeSRedirector *cstream = GetDebugStreamer();
373     if (cstream){
374       TTimeStamp tstamp(fTime);
375       Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
376       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
377       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
378       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
379       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
380       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
381       TVectorD vecGoofie(20);
382       AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
383       if (goofieArray){
384         for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
385           AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
386           if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
387         }
388       }
389       (*cstream)<<"laserInfo"<<
390         "run="<<fRun<<              //  run number
391         "event="<<fEvent<<          //  event number
392         "time="<<fTime<<            //  time stamp of event
393         "trigger="<<fTrigger<<      //  trigger
394         "mag="<<fMagF<<             //  magnetic field
395         // Environment values
396         "press0="<<valuePressure0<<
397         "press1="<<valuePressure1<<
398         "pt0="<<ptrelative0<<
399         "pt1="<<ptrelative1<<
400         "temp0="<<temp0<<
401         "temp1="<<temp1<<
402         "vecGoofie.="<<&vecGoofie<<
403         //laser
404         "rejectA="<<isReject[0]<<
405         "rejectC="<<isReject[1]<<
406         "laserA.="<<fLaser->fFitAside<<
407         "laserC.="<<fLaser->fFitCside<<
408         "laserAC.="<<fLaser->fFitACside<<
409         "trigger="<<event->GetFiredTriggerClasses()<<
410         "\n";
411     }
412   }
413   //
414   // fill histos
415   //
416   TVectorD vdriftA(5), vdriftC(5),vdriftAC(5);
417   vdriftA=*(fLaser->fFitAside);
418   vdriftC=*(fLaser->fFitCside);
419   vdriftAC=*(fLaser->fFitACside);
420   Int_t npointsA=0, npointsC=0;
421   Float_t chi2A=0, chi2C=0;
422   npointsA= TMath::Nint(vdriftA[3]);
423   chi2A= vdriftA[4];
424   npointsC= TMath::Nint(vdriftC[3]);
425   chi2C= vdriftC[4];
426
427   TTimeStamp tstamp(fTime);
428   Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
429   Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
430   Double_t driftA=0, driftC=0;
431   if (vdriftA[1]>1.-kMaxDeltaV) driftA = 1./vdriftA[1]-1.;
432   if (vdriftC[1]>1.-kMaxDeltaV) driftC = 1./vdriftC[1]-1.;
433   //
434   Double_t vecDriftLaserA[4]={fTime,(ptrelative0+ptrelative1)/2.0,driftA,event->GetRunNumber()};
435   Double_t vecDriftLaserC[4]={fTime,(ptrelative0+ptrelative1)/2.0,driftC,event->GetRunNumber()};
436   //  Double_t vecDrift[4]      ={fTime,(ptrelative0+ptrelative1)/2.0,1./((*(fLaser->fFitACside))[1])-1,event->GetRunNumber()};
437
438   for (Int_t icalib=0;icalib<3;icalib++){
439     if (icalib==0){ //z0 shift
440       vecDriftLaserA[2]=vdriftA[0]/250.;
441       vecDriftLaserC[2]=vdriftC[0]/250.;
442     }
443     if (icalib==1){ //vdrel shift
444       vecDriftLaserA[2]=driftA;
445       vecDriftLaserC[2]=driftC;
446     }
447     if (icalib==2){ //gy shift - full gy - full drift
448       vecDriftLaserA[2]=vdriftA[2]/250.;
449       vecDriftLaserC[2]=vdriftC[2]/250.;
450     }
451     if (isReject[0]==0) fHistVdriftLaserA[icalib]->Fill(vecDriftLaserA);
452     if (isReject[1]==0) fHistVdriftLaserC[icalib]->Fill(vecDriftLaserC);
453   }
454
455 //   THnSparse* curHist=new THnSparseF("","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
456 //   TString shortName=curHist->ClassName();
457 //   shortName+="_MEAN_DRIFT_LASER_";
458 //   delete curHist;
459 //   curHist=NULL;
460 //   TString name="";
461
462 //   name=shortName;
463 //   name+=event->GetFiredTriggerClasses();
464 //   name.ToUpper();
465 //   curHist=(THnSparseF*)fArrayDz->FindObject(name);
466 //   if(!curHist){
467 //     curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
468 //     fArrayDz->AddLast(curHist);
469 //   }
470 //   curHist->Fill(vecDrift);
471           
472 //   name=shortName;
473 //   name+="ALL";
474 //   name.ToUpper();
475 //   curHist=(THnSparseF*)fArrayDz->FindObject(name);
476 //   if(!curHist){
477 //     curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
478 //     fArrayDz->AddLast(curHist);
479 //   }
480 //   curHist->Fill(vecDrift);
481 }
482
483 void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event){
484   if (!event) {
485     Printf("ERROR: ESD not available");
486     return;
487   }  
488   if (event->GetTimeStamp() == 0 ) {
489     Printf("no time stamp!");
490     return;
491   }
492   
493   //fd
494   // Find cosmic pairs
495   // 
496   // Track0 is choosen in upper TPC part
497   // Track1 is choosen in lower TPC part
498   //
499   Int_t ntracks=event->GetNumberOfTracks();
500   if (ntracks==0) return;
501   if (ntracks > fCutTracks) return;
502   
503   if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
504   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
505   
506   TObjArray  tpcSeeds(ntracks);
507   Double_t vtxx[3]={0,0,0};
508   Double_t svtxx[3]={0.000001,0.000001,100.};
509   AliESDVertex vtx(vtxx,svtxx);
510   //
511   // track loop
512   //
513   for (Int_t i=0;i<ntracks;++i) {
514     AliESDtrack *track = event->GetTrack(i);
515     
516     const AliExternalTrackParam * trackIn = track->GetInnerParam();
517     const AliExternalTrackParam * trackOut = track->GetOuterParam();
518     if (!trackIn) continue;
519     if (!trackOut) continue;
520     
521     AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
522     if (friendTrack) ProcessAlignITS(track,friendTrack);
523     if (friendTrack) ProcessAlignTOF(track,friendTrack);
524     TObject *calibObject;
525     AliTPCseed *seed = 0;
526     for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
527     if (seed) tpcSeeds.AddAt(seed,i);
528   }
529   if (ntracks<2) return;
530   //
531   // Find pairs
532   //
533
534   for (Int_t i=0;i<ntracks;++i) {
535     AliESDtrack *track0 = event->GetTrack(i);
536     // track0 - choosen upper part
537     if (!track0) continue;
538     if (!track0->GetOuterParam()) continue;
539     if (track0->GetOuterParam()->GetAlpha()<0) continue;
540     Double_t d1[3];
541     track0->GetDirection(d1);    
542     for (Int_t j=0;j<ntracks;++j) {
543       if (i==j) continue;
544       AliESDtrack *track1 = event->GetTrack(j);   
545       //track 1 lower part
546       if (!track1) continue;
547       if (!track1->GetOuterParam()) continue;
548       if (track1->GetOuterParam()->GetAlpha()>0) continue;
549       //
550       Double_t d2[3];
551       track1->GetDirection(d2);
552       
553       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
554       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
555       if (! seed0) continue;
556       if (! seed1) continue;
557       Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
558       Float_t dist0  = track0->GetLinearD(0,0);
559       Float_t dist1  = track1->GetLinearD(0,0);
560       //
561       // conservative cuts - convergence to be guarantied
562       // applying before track propagation
563       if (TMath::Abs(dist0+dist1)>fCutMaxD) continue;   // distance to the 0,0
564       if (dir>fCutMinDir) continue;               // direction vector product
565       Float_t bz = AliTracker::GetBz();
566       Float_t dvertex0[2];   //distance to 0,0
567       Float_t dvertex1[2];   //distance to 0,0 
568       track0->GetDZ(0,0,0,bz,dvertex0);
569       track1->GetDZ(0,0,0,bz,dvertex1);
570       if (TMath::Abs(dvertex0[1])>250) continue;
571       if (TMath::Abs(dvertex1[1])>250) continue;
572       //
573       //
574       //
575       Float_t dmax = TMath::Max(TMath::Abs(dist0),TMath::Abs(dist1));
576       AliExternalTrackParam param0(*track0);
577       AliExternalTrackParam param1(*track1);
578       //
579       // Propagate using Magnetic field and correct fo material budget
580       //
581       AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
582       AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
583       //
584       // Propagate rest to the 0,0 DCA - z should be ignored
585       //
586       //Bool_t b0 = ;
587       param0.PropagateToDCA(&vtx,bz,1000);
588       //Bool_t b1 = 
589       param1.PropagateToDCA(&vtx,bz,1000);
590       param0.GetDZ(0,0,0,bz,dvertex0);
591       param1.GetDZ(0,0,0,bz,dvertex1);
592       Double_t xyz0[3];
593       Double_t xyz1[3];
594       param0.GetXYZ(xyz0);
595       param1.GetXYZ(xyz1);
596       Bool_t isPair = IsPair(&param0,&param1);
597       Bool_t isCross = IsCross(track0, track1);
598       Bool_t isSame = IsSame(track0, track1);
599
600       THnSparse* hist=new THnSparseF("","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
601       TString shortName=hist->ClassName();
602       shortName+="_MEAN_VDRIFT_COSMICS_";
603       delete hist;
604       hist=NULL;
605
606       if(isSame || (isCross && isPair)){
607         if (track0->GetTPCNcls() > 80) {
608           fDz = param0.GetZ() - param1.GetZ();
609           if(track0->GetOuterParam()->GetZ()<0) fDz=-fDz;
610           TTimeStamp tstamp(fTime);
611           Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
612           Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
613           Double_t vecDrift[4]={fTime,(ptrelative0+ptrelative1)/2.0,fDz/500.0,event->GetRunNumber()};
614           THnSparse* curHist=NULL;
615           TString name="";
616
617           name=shortName;
618           name+=event->GetFiredTriggerClasses();
619           name.ToUpper();
620           curHist=(THnSparseF*)fArrayDz->FindObject(name);
621           if(!curHist){
622             curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
623             fArrayDz->AddLast(curHist);
624           }
625 //        curHist=(THnSparseF*)(fMapDz->GetValue(event->GetFiredTriggerClasses()));
626 //        if(!curHist){
627 //          curHist=new THnSparseF(event->GetFiredTriggerClasses(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
628 //          fMapDz->Add(new TObjString(event->GetFiredTriggerClasses()),curHist);
629 //        }
630           curHist->Fill(vecDrift);
631           
632           name=shortName;
633           name+="ALL";
634           name.ToUpper();
635           curHist=(THnSparseF*)fArrayDz->FindObject(name);
636           if(!curHist){
637             curHist=new THnSparseF(name,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
638             fArrayDz->AddLast(curHist);
639           }
640 //        curHist=(THnSparseF*)(fMapDz->GetValue("all"));
641 //        if(!curHist){
642 //          curHist=new THnSparseF("all","HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
643 //          fMapDz->Add(new TObjString("all"),curHist);
644 //        }
645           curHist->Fill(vecDrift);
646         }
647       }
648       TTreeSRedirector *cstream = GetDebugStreamer();
649       if (fStreamLevel>0){
650         if (cstream){
651         (*cstream)<<"trackInfo"<<
652           "tr0.="<<track0<<
653           "tr1.="<<track1<<
654           "p0.="<<&param0<<
655           "p1.="<<&param1<<
656           "isPair="<<isPair<<
657           "isCross="<<isCross<<
658           "isSame="<<isSame<<
659           "fDz="<<fDz<<
660           "fRun="<<fRun<<
661           "fTime="<<fTime<<
662           "\n";
663         }
664       }
665     } // end 2nd order loop        
666   } // end 1st order loop
667   
668   if (fStreamLevel>0){
669     TTreeSRedirector *cstream = GetDebugStreamer();
670     if (cstream){
671       TTimeStamp tstamp(fTime);
672       Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
673       Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
674       Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
675       Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
676       Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
677       Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
678       TVectorD vecGoofie(20);
679       AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
680       if (goofieArray){
681         for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
682           AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
683           if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
684         }
685       }
686       (*cstream)<<"timeInfo"<<
687         "run="<<fRun<<              //  run number
688         "event="<<fEvent<<          //  event number
689         "time="<<fTime<<            //  time stamp of event
690         "trigger="<<fTrigger<<      //  trigger
691         "mag="<<fMagF<<             //  magnetic field
692         // Environment values
693         "press0="<<valuePressure0<<
694         "press1="<<valuePressure1<<
695         "pt0="<<ptrelative0<<
696         "pt1="<<ptrelative1<<
697         "temp0="<<temp0<<
698         "temp1="<<temp1<<
699         "vecGoofie.=<<"<<&vecGoofie<<
700         //
701         // accumulated values
702         //
703         "fDz="<<fDz<<          //! current delta z
704         "trigger="<<event->GetFiredTriggerClasses()<<
705         "\n";
706     }
707   }
708   printf("Trigger: %s\n",event->GetFiredTriggerClasses().Data());
709 }
710
711 void AliTPCcalibTime::ProcessBeam(AliESDEvent */*event*/){
712 }
713
714 void AliTPCcalibTime::Analyze(){}
715
716 THnSparse* AliTPCcalibTime::GetHistoDrift(const char* name){
717   TIterator* iterator = fArrayDz->MakeIterator();
718   iterator->Reset();
719   TString newName=name;
720   newName.ToUpper();
721   THnSparse* newHist=new THnSparseF(newName,"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
722   THnSparse* addHist=NULL;
723   while((addHist=(THnSparseF*)iterator->Next())){
724   if(!addHist) continue;
725     TString histName=addHist->GetName();
726     if(!histName.Contains(newName)) continue;
727     addHist->Print();
728     newHist->Add(addHist);
729   }
730   return newHist;
731 }
732
733 TObjArray* AliTPCcalibTime::GetHistoDrift(){
734   return fArrayDz;
735 }
736
737 TGraphErrors* AliTPCcalibTime::GetGraphDrift(const char* name){
738   THnSparse* histoDrift=GetHistoDrift(name);
739   TGraphErrors* graphDrift=NULL;
740   if(histoDrift){
741     graphDrift=FitSlices(histoDrift,2,0,400,100,0.05,0.95, kTRUE);
742     TString end=histoDrift->GetName();
743     Int_t pos=end.Index("_");
744     end=end(pos,end.Capacity()-pos);
745     TString graphName=graphDrift->ClassName();
746     graphName+=end;
747     graphName.ToUpper();
748     graphDrift->SetName(graphName);
749   }
750   return graphDrift;
751 }
752
753 TObjArray* AliTPCcalibTime::GetGraphDrift(){
754   TObjArray* arrayGraphDrift=new TObjArray();
755   TIterator* iterator=fArrayDz->MakeIterator();
756   iterator->Reset();
757   THnSparse* addHist=NULL;
758   while((addHist=(THnSparseF*)iterator->Next())) arrayGraphDrift->AddLast(GetGraphDrift(addHist->GetName()));
759   return arrayGraphDrift;
760 }
761
762 AliSplineFit* AliTPCcalibTime::GetFitDrift(const char* name){
763   TGraph* graphDrift=GetGraphDrift(name);
764   AliSplineFit* fitDrift=NULL;
765   if(graphDrift && graphDrift->GetN()){
766     fitDrift=new AliSplineFit();
767     fitDrift->SetGraph(graphDrift);
768     fitDrift->SetMinPoints(graphDrift->GetN()+1);
769     fitDrift->InitKnots(graphDrift,2,0,0.001);
770     fitDrift->SplineFit(0);
771     TString end=graphDrift->GetName();
772     Int_t pos=end.Index("_");
773     end=end(pos,end.Capacity()-pos);
774     TString fitName=fitDrift->ClassName();
775     fitName+=end;
776     fitName.ToUpper();
777     //fitDrift->SetName(fitName);
778     delete graphDrift;
779     graphDrift=NULL;
780   }
781   return fitDrift;
782 }
783
784 //TObjArray* AliTPCcalibTime::GetFitDrift(){
785 //  TObjArray* arrayFitDrift=new TObjArray();
786 //  TIterator* iterator = fArrayDz->MakeIterator();
787 //  iterator->Reset();
788 //  THnSparse* addHist=NULL;
789 //  while((addHist=(THnSparseF*)iterator->Next())) arrayFitDrift->AddLast(GetFitDrift(addHist->GetName()));
790 //  return arrayFitDrift;
791 //}
792
793 Long64_t AliTPCcalibTime::Merge(TCollection *li) {
794   TIterator* iter = li->MakeIterator();
795   AliTPCcalibTime* cal = 0;
796
797   while ((cal = (AliTPCcalibTime*)iter->Next())) {
798     if (!cal->InheritsFrom(AliTPCcalibTime::Class())) {
799       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
800       return -1;
801     }
802     for (Int_t imeas=0; imeas<3; imeas++){
803       if (cal->GetHistVdriftLaserA(imeas) && cal->GetHistVdriftLaserA(imeas)){
804         fHistVdriftLaserA[imeas]->Add(cal->GetHistVdriftLaserA(imeas));
805         fHistVdriftLaserC[imeas]->Add(cal->GetHistVdriftLaserC(imeas));
806       }
807     }
808     TObjArray* addArray=cal->GetHistoDrift();
809     if(!addArray) return 0;
810     TIterator* iterator = addArray->MakeIterator();
811     iterator->Reset();
812     THnSparse* addHist=NULL;
813     while((addHist=(THnSparseF*)iterator->Next())){
814       if(!addHist) continue;
815       addHist->Print();
816       THnSparse* localHist=(THnSparseF*)fArrayDz->FindObject(addHist->GetName());
817       if(!localHist){
818         localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
819         fArrayDz->AddLast(localHist);
820       }
821       localHist->Add(addHist);
822     }
823 //    TMap * addMap=cal->GetHistoDrift();
824 //    if(!addMap) return 0;
825 //    TIterator* iterator = addMap->MakeIterator();
826 //    iterator->Reset();
827 //    TPair* addPair=0;
828 //    while((addPair=(TPair *)(addMap->FindObject(iterator->Next())))){
829 //      THnSparse* addHist=dynamic_cast<THnSparseF*>(addPair->Value());
830 //      if (!addHist) continue;
831 //      addHist->Print();
832 //      THnSparse* localHist=dynamic_cast<THnSparseF*>(fMapDz->GetValue(addHist->GetName()));
833 //      if(!localHist){
834 //        localHist=new THnSparseF(addHist->GetName(),"HistVdrift;time;p/T ratio;Vdrift;run",4,fBinsVdrift,fXminVdrift,fXmaxVdrift);
835 //        fMapDz->Add(new TObjString(addHist->GetName()),localHist);
836 //      }
837 //      localHist->Add(addHist);
838 //    }
839     for(Int_t i=0;i<10;i++) if (cal->GetCosmiMatchingHisto(i)) fCosmiMatchingHisto[i]->Add(cal->GetCosmiMatchingHisto(i));
840   }
841   return 0;
842 }
843
844 Bool_t  AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
845   /*
846   // 0. Same direction - OPOSITE  - cutDir +cutT    
847   TCut cutDir("cutDir","dir<-0.99")
848   // 1. 
849   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
850   //
851   // 2. The same rphi 
852   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
853   //
854   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");  
855   // 1/Pt diff cut
856   */
857   const Double_t *p0 = tr0->GetParameter();
858   const Double_t *p1 = tr1->GetParameter();
859   fCosmiMatchingHisto[0]->Fill(p0[0]+p1[0]);
860   fCosmiMatchingHisto[1]->Fill(p0[1]-p1[1]);
861   fCosmiMatchingHisto[2]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
862   fCosmiMatchingHisto[3]->Fill(p0[3]+p1[3]);
863   fCosmiMatchingHisto[4]->Fill(p0[4]+p1[4]);
864   
865   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
866   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
867   if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz)  return kFALSE;
868   Double_t d0[3], d1[3];
869   tr0->GetDirection(d0);    
870   tr1->GetDirection(d1);       
871   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
872
873   fCosmiMatchingHisto[5]->Fill(p0[0]+p1[0]);
874   fCosmiMatchingHisto[6]->Fill(p0[1]-p1[1]);
875   fCosmiMatchingHisto[7]->Fill(tr1->GetAlpha()-tr0->GetAlpha()+TMath::Pi());
876   fCosmiMatchingHisto[8]->Fill(p0[3]+p1[3]);
877   fCosmiMatchingHisto[9]->Fill(p0[4]+p1[4]);
878
879   return kTRUE;  
880 }
881 Bool_t AliTPCcalibTime::IsCross(AliESDtrack *tr0, AliESDtrack *tr1){
882   return  tr0->GetOuterParam()->GetZ()*tr1->GetOuterParam()->GetZ()<0 && tr0->GetInnerParam()->GetZ()*tr1->GetInnerParam()->GetZ()<0 && tr0->GetOuterParam()->GetZ()*tr0->GetInnerParam()->GetZ()>0 && tr1->GetOuterParam()->GetZ()*tr1->GetInnerParam()->GetZ()>0;
883 }
884
885 Bool_t AliTPCcalibTime::IsSame(AliESDtrack */*tr0*/, AliESDtrack */*tr1*/){
886   // To be implemented
887   return kFALSE;
888 }
889
890 /*
891 chainDrift->Draw("p0.fP[0]+p1.fP[0]","isPair");
892   mean ~-0.02  ~-0.03913
893   RMS  ~ 0.5   ~ 0.5356    --> 3    (fCutMaxD)
894
895 chainDrift->Draw("p0.fP[1]-p1.fP[1]","isPair");
896   mean         ~ 0.1855
897   RMS          ~ 4.541     -->25    (fCutMaxDz)
898
899 chainDrift->Draw("p1.fAlpha-p0.fAlpha+pi","isPair");
900 //chainDrift->Draw("p1.fAlpha+p0.fAlpha","isPair");
901 //chainDrift->Draw("p1.fP[2]-p0.fP[2]+pi","isPair");
902 //chainDrift->Draw("p1.fP[2]+p0.fP[2]","isPair");
903   mean ~ 0     ~ 0.001898
904   RMS  ~ 0.009 ~ 0.01134   --> 0.06
905
906 chainDrift->Draw("p0.fP[3]+p1.fP[3]","isPair");
907   mean ~ 0.0013 ~ 0.001539
908   RMS  ~ 0.003  ~ 0.004644 --> 0.03  (fCutTheta)
909
910 chainDrift->Draw("p1.fP[4]+p0.fP[4]>>his(100,-0.2,0.2)","isPair")
911   mean ~ 0.012  ~-0.0009729
912   RMS  ~ 0.036  ~ 0.03773  --> 0.2
913 */
914
915
916 void  AliTPCcalibTime::ProcessAlignITS(AliESDtrack* track, AliESDfriendTrack *friendTrack){
917   //
918   // Process track
919   // Update TPC-ITS alignment
920   //
921   const Int_t    kMinTPC  = 80;
922   const Int_t    kMinITS  = 3;
923   const Double_t kMinZ    = 10;
924   const Double_t kMaxDy   = 2;
925   const Double_t kMaxAngle= 0.02;
926   //
927   Int_t dummycl[1000];
928   if (track->GetITSclusters(dummycl)<kMinITS) return;  // minimal amount of clusters
929   if (track->GetTPCNcls()<kMinTPC) return;  // minimal amount of clusters cut
930   //
931   if (!friendTrack->GetITSOut()) return;
932   if (!track->GetInnerParam())   return;
933   if (!track->GetOuterParam())   return;
934   // exclude crossing track
935   if (track->GetOuterParam()->GetZ()*track->GetInnerParam()->GetZ()<0)   return;
936   if (TMath::Abs(track->GetInnerParam()->GetZ())<kMinZ)   return;
937   //
938   AliExternalTrackParam &pTPC=(AliExternalTrackParam &)(*(track->GetInnerParam()));
939   AliExternalTrackParam pITS(*(friendTrack->GetITSOut()));
940   //
941   //
942   //
943   Int_t htime = fTime/3600; //time in hours
944   if (fAlignITSTPC->GetEntries()<htime){
945     fAlignITSTPC->Expand(htime*2+20);
946   }
947   AliRelAlignerKalman* align =  (AliRelAlignerKalman*)fAlignITSTPC->At(htime);
948   if (!align){
949     align=new AliRelAlignerKalman(); 
950     align->SetOutRejSigma(2.);
951     //align->SetRejectOutliers(kFALSE);
952     align->SetMagField(fMagF); 
953     fAlignITSTPC->AddAt(align,htime);
954   }
955   pITS.Rotate(pTPC.GetAlpha());
956   pITS.PropagateTo(pTPC.GetX(),fMagF);
957   if (TMath::Abs(pITS.GetY()-pTPC.GetY())>kMaxDy) return;
958   if (TMath::Abs(pITS.GetSnp()-pTPC.GetSnp())>kMaxAngle) return;
959   if (TMath::Abs(pITS.GetTgl()-pTPC.GetTgl())>kMaxAngle) return;
960   align->AddTrackParams(&pITS,&pTPC);
961   align->SetTimeStamp(fTime);
962   //  align->SetRunNumber(fRun );
963   static Int_t entry=-1;
964   entry++;
965   //  Int_t nupdates=align->GetNUpdates();
966   Int_t nupdates=entry;
967   align->SetOutRejSigma(1.+1./Double_t(nupdates));
968   TTreeSRedirector *cstream = GetDebugStreamer();  
969   if (cstream && align->GetState() && align->GetState()->GetNrows()>2 ){
970     TTimeStamp tstamp(fTime);
971     Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
972     Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
973     Double_t ptrelative0   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
974     Double_t ptrelative1   = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
975     Double_t temp0         = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
976     Double_t temp1         = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
977     TVectorD vecGoofie(20);
978     AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
979     if (goofieArray){
980       for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
981         AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
982         if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
983       }
984     }
985     TVectorD gpTPC(3), gdTPC(3);
986     TVectorD gpITS(3), gdITS(3);
987     pTPC.GetXYZ(gpTPC.GetMatrixArray());
988     pTPC.GetDirection(gdTPC.GetMatrixArray());
989     pITS.GetXYZ(gpITS.GetMatrixArray());
990     pITS.GetDirection(gdITS.GetMatrixArray());
991     (*cstream)<<"itstpc"<<
992       "run="<<fRun<<              //  run number
993       "event="<<fEvent<<          //  event number
994       "time="<<fTime<<            //  time stamp of event
995       "trigger="<<fTrigger<<      //  trigger
996       "mag="<<fMagF<<             //  magnetic field
997       // Environment values
998       "press0="<<valuePressure0<<
999       "press1="<<valuePressure1<<
1000       "pt0="<<ptrelative0<<
1001       "pt1="<<ptrelative1<<
1002       "temp0="<<temp0<<
1003       "temp1="<<temp1<<
1004       "vecGoofie.="<<&vecGoofie<<
1005       "entry="<<entry<<  // current entry
1006       //
1007       "a.="<<align<<     // current alignment
1008       "pITS.="<<&pITS<<  // track param ITS
1009       "pTPC.="<<&pTPC<<  // track param TPC
1010       "gpTPC.="<<&gpTPC<<
1011       "gdTPC.="<<&gdTPC<<
1012       "gpITS.="<<&gpITS<<
1013       "gdITS.="<<&gdITS<<
1014       "\n";
1015   }
1016   
1017 }
1018
1019
1020 void  AliTPCcalibTime::ProcessAlignTOF(AliESDtrack* track, AliESDfriendTrack *friendTrack){
1021   //
1022   //process TOF-TPC  alignment
1023   //
1024   Int_t kminNcl=80;
1025   Float_t kMaxDy=6;
1026   Float_t kMaxDz=10;
1027   if (track->GetTPCNcls()<kminNcl) return;
1028   if (track->GetOuterParam()==0) return;
1029   if (track->GetInnerParam()==0) return;
1030   if (track->GetTOFsignal()<=0)  return;
1031   //
1032   AliExternalTrackParam *paramOut = new AliExternalTrackParam(*(track->GetOuterParam()));
1033   AliExternalTrackParam *param=0;
1034   const AliTrackPointArray *points=friendTrack->GetTrackPointArray();
1035   if (!points) return;
1036   Int_t npoints = points->GetNPoints();
1037   AliTrackPoint point;
1038   //Double_t alpha=
1039   Double_t mass = TDatabasePDG::Instance()->GetParticle("mu+")->Mass();
1040   TTreeSRedirector * cstream =  GetDebugStreamer();
1041   //
1042   //
1043   //
1044   for (Int_t ipoint=0;ipoint<npoints;ipoint++){
1045     //
1046     points->GetPoint(point,ipoint);
1047     Float_t xyz[3];
1048     point.GetXYZ(xyz);
1049     Double_t r=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
1050     if (r<300)  continue;
1051     if (r>400)  continue;
1052     param=paramOut;
1053     if (!param) continue;
1054     AliTracker::PropagateTrackToBxByBz(param,r,mass,2.,kTRUE);
1055     AliTracker::PropagateTrackToBxByBz(param,r,mass,0.1,kTRUE);    
1056     AliTrackPoint lpoint = point.Rotate(param->GetAlpha());
1057     param->PropagateTo(lpoint.GetX(),fMagF);
1058     //
1059     //
1060     // this is ugly - we need AliCluster constructor
1061     //
1062     AliExternalTrackParam &pTPC=*param;
1063     AliExternalTrackParam pTOF(*param);
1064     ((Double_t*)pTOF.GetParameter())[0] =lpoint.GetY();
1065     ((Double_t*)pTOF.GetParameter())[1] =lpoint.GetZ();
1066     pTOF.ResetCovariance(20);
1067     ((Double_t*)pTOF.GetCovariance())[0]+=3.*3.;
1068     ((Double_t*)pTOF.GetCovariance())[2]+=3.*3.;
1069     if (TMath::Abs(pTOF.GetY()-pTPC.GetY())>kMaxDy) continue;
1070     if (TMath::Abs(pTOF.GetZ()-pTPC.GetZ())>kMaxDz) continue;
1071     //
1072     Int_t htime = fTime/3600; //time in hours
1073       
1074     if (fAlignTOFTPC->GetEntries()<htime){
1075       fAlignTOFTPC->Expand(htime*2+20);
1076     }
1077     AliRelAlignerKalman* align =  (AliRelAlignerKalman*)fAlignTOFTPC->At(htime);
1078     if (!align){
1079       align=new AliRelAlignerKalman(); 
1080       align->SetOutRejSigma(2.);
1081       //align->SetRejectOutliers(kFALSE);
1082       align->SetMagField(fMagF); 
1083       fAlignTOFTPC->AddAt(align,htime);
1084     }
1085     pTOF.Rotate(pTPC.GetAlpha());
1086     pTOF.PropagateTo(pTPC.GetX(),fMagF);
1087     align->AddTrackParams(&pTOF,&pTPC);
1088     align->SetTimeStamp(fTime);
1089     static Int_t entry=-1;
1090     entry++;
1091     //    Int_t nupdates=align->GetNUpdates();
1092     Int_t nupdates=entry;
1093     align->SetOutRejSigma(1.+1./Double_t(nupdates));
1094     
1095     //
1096     //
1097     if (cstream) {
1098       (*cstream) << "pointMatch" <<
1099         "run="<<fRun<<              //  run number
1100         "event="<<fEvent<<          //  event number
1101         "time="<<fTime<<            //  time stamp of event
1102         "trigger="<<fTrigger<<      //  trigger
1103         "mag="<<fMagF<<             //  magnetic field
1104         //
1105         "a.="<<align<<     // current alignment 
1106         "p.="<<&point<<
1107         "lp.="<<&lpoint<<
1108         "pTPC.="<<&pTPC<<
1109         "pTOF.="<<&pTOF<<
1110         "\n";
1111     }
1112
1113
1114
1115   }
1116   delete paramOut;
1117 }