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