Adding new class AliTPCcalibTime
[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     1. What do we calibrate.
19
20 Time dependence of gain and drift velocity in order to account for changes in: temperature, pressure, gas composition.
21
22   AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime",0, 1213.9e+06, 1213.96e+06, 0.04e+04, 0.04e+04);
23
24
25     2. How to interpret results
26     3. Simple example
27
28 a.) determine the required time range:
29
30 AliXRDPROOFtoolkit tool;
31 TChain * chain = tool.MakeChain("pass2.txt","esdTree",0,6000);
32 chain->Draw("GetTimeStamp()")
33
34 b) analyse calibration object on Proof in calibration train 
35
36 AliTPCcalibTime *calibTime = new AliTPCcalibTime("cosmicTime","cosmicTime", StartTimeStamp, EndTimeStamp, IntegrationTimeVdrift, IntegrationTimeDeDx);
37
38 c) plot results
39
40 TFile f("CalibObjects.root");
41 AliTPCcalibTime *cal = f->Get("TPCCalib")->FindObject("cosmicTime");
42 cal->GetHistVdrift()->Projection(1,0)->Draw()
43
44     4. Analysis using debug streamers.    
45
46 */
47
48
49 #include "Riostream.h"
50 #include "TChain.h"
51 #include "TTree.h"
52 #include "TH1F.h"
53 #include "TH2F.h"
54 #include "TH3F.h"
55 #include "THnSparse.h"
56 #include "TList.h"
57 #include "TMath.h"
58 #include "TCanvas.h"
59 #include "TFile.h"
60 #include "TF1.h"
61 #include "TVectorD.h"
62 #include "TProfile.h"
63 #include "TGraphErrors.h"
64 #include "TCanvas.h"
65
66 #include "AliTPCclusterMI.h"
67 #include "AliTPCseed.h"
68 #include "AliESDVertex.h"
69 #include "AliESDEvent.h"
70 #include "AliESDfriend.h"
71 #include "AliESDInputHandler.h"
72 #include "AliAnalysisManager.h"
73
74 #include "AliTracker.h"
75 #include "AliMagFMaps.h"
76 #include "AliTPCCalROC.h"
77
78 #include "AliLog.h"
79
80 #include "AliTPCcalibTime.h"
81
82 #include "TTreeStream.h"
83 #include "AliTPCTracklet.h"
84
85 ClassImp(AliTPCcalibTime)
86
87
88 AliTPCcalibTime::AliTPCcalibTime() 
89   :AliTPCcalibBase(),
90    fHistDeDxTgl(0),
91    fHistDeDx(0),
92    fHistVdrift(0),
93    fIntegrationTimeDeDx(0),
94    fIntegrationTimeVdrift(0),
95    fCutMaxD(5),        // maximal distance in rfi ditection
96    fCutTheta(0.03),    // maximal distan theta
97    fCutMinDir(-0.99)   // direction vector products
98
99 {  
100   AliInfo("Default Constructor");  
101 }
102
103
104 AliTPCcalibTime::AliTPCcalibTime(const Text_t *name, const Text_t *title, ULong64_t TriggerMask, UInt_t StartTime, UInt_t EndTime, Int_t deltaIntegrationTimeDeDx, Int_t deltaIntegrationTimeVdrift)
105   :AliTPCcalibBase(),
106    fTriggerMask(0),
107    fHistDeDxTgl(0),
108    fHistDeDx(0),
109    fHistVdrift(0),
110    fIntegrationTimeDeDx(0),
111    fIntegrationTimeVdrift(0),
112    fCutMaxD(5),        // maximal distance in rfi ditection
113    fCutTheta(0.03),    // maximal distan theta
114    fCutMinDir(-0.99)   // direction vector products
115 {
116   
117   SetName(name);
118   SetTitle(title);
119
120   AliInfo("Non Default Constructor");
121
122   fTriggerMask = TriggerMask;
123
124   fIntegrationTimeDeDx = deltaIntegrationTimeDeDx;
125   fIntegrationTimeVdrift = deltaIntegrationTimeVdrift;
126
127   Double_t deltaTime = EndTime - StartTime;
128   
129   Int_t binsVdrift[2] = {deltaTime/deltaIntegrationTimeVdrift, 100};
130   Double_t xminVdrift[2] = {StartTime, -20};
131   Double_t xmaxVdrift[2] = {EndTime, 20};
132   fHistVdrift = new THnSparseF("HistVdrift","vDrift; time;#Delta z",2,binsVdrift,xminVdrift,xmaxVdrift);
133
134   Int_t binsDeDxTgl[3] = {deltaTime/deltaIntegrationTimeDeDx,30,100};
135   Double_t xminDeDxTgl[3] = {StartTime,-1,0.7};
136   Double_t xmaxDeDxTgl[3] = {EndTime,1,1.3};
137   fHistDeDxTgl = new THnSparseF("HistDeDxTgl","dEdx vs tgl;time;tgl;dEdxUp/dEdxLow",3,binsDeDxTgl,xminDeDxTgl,xmaxDeDxTgl); 
138
139   Int_t binsDeDx[2] = {deltaTime/deltaIntegrationTimeDeDx,100};
140   Double_t xminDeDx[2] = {StartTime,1};
141   Double_t xmaxDeDx[2] = {EndTime,100};
142   fHistDeDx = new THnSparseF("HistDeDx","dEdx l;time;dEdx",2,binsDeDx,xminDeDx,xmaxDeDx);
143
144 }
145
146
147
148 AliTPCcalibTime::~AliTPCcalibTime(){
149   //
150   //
151   //
152 }
153
154 void AliTPCcalibTime::Process(AliESDEvent *event) {
155   //
156   //
157   //
158   
159   ProcessCosmic(event);
160
161 }
162
163
164
165 void AliTPCcalibTime::ProcessCosmic(AliESDEvent *event) {
166
167   if (!event) {
168     Printf("ERROR: ESD not available");
169     return;
170   }  
171   if (event->GetTimeStamp() == 0 ) {
172     Printf("no time stamp!");
173     return;
174   }
175
176   if (fTriggerMask != 0 && event->GetTriggerMask() != fTriggerMask) return;
177
178   UInt_t time = event->GetTimeStamp();
179
180   //
181   // Find cosmic pairs
182   // 
183   // Track0 is choosen in upper TPC part
184   // Track1 is choosen in lower TPC part
185   //
186   Int_t ntracks=event->GetNumberOfTracks();
187   if (ntracks==0) return;
188   if (ntracks > 10) return; // temporary debug to remove LASER events
189
190
191   if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
192   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
193   
194   TObjArray  tpcSeeds(ntracks);
195   Double_t vtxx[3]={0,0,0};
196   Double_t svtxx[3]={0.000001,0.000001,100.};
197   AliESDVertex vtx(vtxx,svtxx);
198   //
199   // track loop
200   //
201   for (Int_t i=0;i<ntracks;++i) {
202    AliESDtrack *track = event->GetTrack(i);
203
204    const AliExternalTrackParam * trackIn = track->GetInnerParam();
205    const AliExternalTrackParam * trackOut = track->GetOuterParam();
206    if (!trackIn) continue;
207    if (!trackOut) continue;
208    
209    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
210    TObject *calibObject;
211    AliTPCseed *seed = 0;
212    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
213      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
214    }
215    if (seed) { 
216      tpcSeeds.AddAt(seed,i);
217      if (track->GetTPCNcls() > 50) {
218        Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
219        Double_t TPCsignal = seed->CookdEdxNorm(0.0,0.6,1,0,159,0x0,kTRUE,kTRUE);
220        Double_t vecDeDx[2] = {time, TPCsignal};
221        if (meanP > 12) fHistDeDx->Fill(vecDeDx);
222      }
223    }
224   }
225   
226   if (ntracks<2) return;
227   //
228   // Find pairs
229   //
230   for (Int_t i=0;i<ntracks;++i) {
231     AliESDtrack *track0 = event->GetTrack(i);     
232     // track0 - choosen upper part
233     if (!track0) continue;
234     if (!track0->GetOuterParam()) continue;
235     if (track0->GetOuterParam()->GetAlpha()<0) continue;
236     Double_t d1[3];
237     track0->GetDirection(d1);    
238     for (Int_t j=0;j<ntracks;++j) {
239       if (i==j) continue;
240        AliESDtrack *track1 = event->GetTrack(j);   
241        //track 1 lower part
242        if (!track1) continue;
243        if (!track1->GetOuterParam()) continue;
244        if (track1->GetOuterParam()->GetAlpha()>0) continue;
245        //
246        Double_t d2[3];
247        track1->GetDirection(d2);
248       
249        AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
250        AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
251        if (! seed0) continue;
252        if (! seed1) continue;
253        Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
254        Float_t d0  = track0->GetLinearD(0,0);
255        Float_t d1  = track1->GetLinearD(0,0);
256        //
257        // conservative cuts - convergence to be guarantied
258        // applying before track propagation
259        if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
260        if (dir>fCutMinDir) continue;               // direction vector product
261        Float_t bz = AliTracker::GetBz();
262        Float_t dvertex0[2];   //distance to 0,0
263        Float_t dvertex1[2];   //distance to 0,0 
264        track0->GetDZ(0,0,0,bz,dvertex0);
265        track1->GetDZ(0,0,0,bz,dvertex1);
266        if (TMath::Abs(dvertex0[1])>250) continue;
267        if (TMath::Abs(dvertex1[1])>250) continue;
268        //
269        //
270        //
271        Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
272        AliExternalTrackParam param0(*track0);
273        AliExternalTrackParam param1(*track1);
274        //
275        // Propagate using Magnetic field and correct fo material budget
276        //
277        AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
278        AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
279        //
280        // Propagate rest to the 0,0 DCA - z should be ignored
281        //
282        Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
283        Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
284        //      
285        param0.GetDZ(0,0,0,bz,dvertex0);
286        param1.GetDZ(0,0,0,bz,dvertex1);
287        //
288        Double_t xyz0[3];//,pxyz0[3];
289        Double_t xyz1[3];//,pxyz1[3];
290        param0.GetXYZ(xyz0);
291        param1.GetXYZ(xyz1);
292        Bool_t isPair = IsPair(&param0,&param1);
293
294        Double_t z0 = track0->GetOuterParam()->GetZ();
295        Double_t z1 = track1->GetOuterParam()->GetZ();
296
297        Double_t z0inner = track0->GetInnerParam()->GetZ();
298        Double_t z1inner = track1->GetInnerParam()->GetZ();
299
300        if (isPair && z0>0 && z0inner>0 && z1<0 && z1inner<0) {
301          Double_t vecVdrift[2] = {time, param0.GetZ() - param1.GetZ()};
302          
303          if (track0->GetTPCNcls() > 80) fHistVdrift->Fill(vecVdrift);
304        }
305        if (isPair && z0 > 0 && z1 > 0) {
306          if (track1->GetTPCNcls()> 110 && track0->GetTPCNcls()> 110 && seed1->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE) > 0) {
307            Double_t vecDeDxTgl[3] = {time, track0->GetTgl(), seed0->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)/seed1->CookdEdxNorm(0,0.6,1,0,159,0,kFALSE,kTRUE)};
308            fHistDeDxTgl->Fill(vecDeDxTgl);
309          }
310        }
311        
312     } // end 2nd order loop        
313   } // end 1st order loop
314
315 }
316
317
318 void AliTPCcalibTime::Analyze() {
319
320   TH2D * hVdrift = GetHistVdrift()->Projection(1,0);
321   
322
323 }
324
325
326 Long64_t AliTPCcalibTime::Merge(TCollection *li) {
327
328   TIterator* iter = li->MakeIterator();
329   AliTPCcalibTime* cal = 0;
330
331   while ((cal = (AliTPCcalibTime*)iter->Next())) {
332     if (!cal->InheritsFrom(AliTPCcalibTime::Class())) {
333       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
334       return -1;
335     }
336
337     // add histograms here...
338     fHistDeDxTgl->Add(cal->GetHistDeDxVsTgl());
339     fHistVdrift->Add(cal->GetHistVdrift());
340     fHistDeDx->Add(cal->GetHistDeDx());
341
342   }
343   
344   return 0;
345   
346 }
347
348
349
350 Bool_t  AliTPCcalibTime::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
351   //
352   //
353   /*
354   // 0. Same direction - OPOSITE  - cutDir +cutT    
355   TCut cutDir("cutDir","dir<-0.99")
356   // 1. 
357   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
358   //
359   // 2. The same rphi 
360   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
361   //
362   //
363   //
364   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");  
365   // 1/Pt diff cut
366   */
367   const Double_t *p0 = tr0->GetParameter();
368   const Double_t *p1 = tr1->GetParameter();
369   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
370   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
371   Double_t d0[3], d1[3];
372   tr0->GetDirection(d0);    
373   tr1->GetDirection(d1);       
374   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
375   //
376   return kTRUE;  
377 }