Attempt to get estimate of material budget
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibMaterial.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          *
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 /*
18   // Load libraries
19
20   gSystem->Load("libANALYSIS");
21   gSystem->Load("libTPCcalib");
22   
23   
24   .x ~/NimStyle.C
25   gSystem->Load("libANALYSIS");
26   gSystem->Load("libTPCcalib");
27
28   // analyze results
29
30   TFile f("CalibObjectsTrain2.root");
31   AliTPCcalibMaterial *calibMaterial = (AliTPCcalibMaterial *)f->Get("alignMaterial");
32
33
34 */
35
36 #include "Riostream.h"
37 #include "TChain.h"
38 #include "TTree.h"
39 #include "TH1F.h"
40 #include "TH2F.h"
41 #include "TH3F.h"
42 #include "THnSparse.h"
43 #include "TList.h"
44 #include "TMath.h"
45 #include "TCanvas.h"
46 #include "TFile.h"
47 #include "TF1.h"
48 #include "TVectorD.h"
49 #include "TProfile.h"
50 #include "TGraphErrors.h"
51 #include "TCanvas.h"
52
53 #include "AliTPCclusterMI.h"
54 #include "AliTPCseed.h"
55 #include "AliESDVertex.h"
56 #include "AliESDEvent.h"
57 #include "AliESDfriend.h"
58 #include "AliESDInputHandler.h"
59 #include "AliAnalysisManager.h"
60
61 #include "AliTracker.h"
62 #include "AliMagF.h"
63 #include "AliTPCCalROC.h"
64
65 #include "AliLog.h"
66
67 #include "AliTPCcalibMaterial.h"
68
69 #include "TTreeStream.h"
70 #include "AliTPCTracklet.h"
71 #include "TTimeStamp.h"
72 #include "AliTPCcalibDB.h"
73 #include "AliTPCcalibLaser.h"
74 #include "AliDCSSensorArray.h"
75 #include "AliDCSSensor.h"
76
77 ClassImp(AliTPCcalibMaterial)
78
79 AliTPCcalibMaterial::AliTPCcalibMaterial():
80   AliTPCcalibBase("calibMaterial","calibMaterial"),
81   fHisMaterial(0),
82   fHisMaterialRPhi(0)
83 {
84   
85 }
86
87 AliTPCcalibMaterial::AliTPCcalibMaterial(const char * name, const char * title):
88   AliTPCcalibBase(name,title),
89   fHisMaterial(0),
90   fHisMaterialRPhi(0)
91 {
92   //
93   //
94   //
95 }
96
97 AliTPCcalibMaterial::~AliTPCcalibMaterial(){
98   //
99   // delete histograms
100   // class is owner of all histograms
101   //
102   if (!fHisMaterial) return;
103   delete fHisMaterial;
104   delete fHisMaterialRPhi;
105   fHisMaterial=0;
106 }
107
108
109 Long64_t AliTPCcalibMaterial::Merge(TCollection *li) {
110   //
111   // Merge histograms
112   //
113   TIterator* iter = li->MakeIterator();
114   AliTPCcalibMaterial* cal = 0;
115
116   while ((cal = (AliTPCcalibMaterial*)iter->Next())) {
117     if (!cal->InheritsFrom(AliTPCcalibMaterial::Class())) {
118       Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
119       return -1;
120     }
121     AliTPCcalibMaterial* calib= (AliTPCcalibMaterial*)(cal);
122     if(!calib) return 0;
123     //
124     if (!fHisMaterial) fHisMaterial=MakeHisto();
125     fHisMaterial->Add(calib->fHisMaterial);
126     fHisMaterialRPhi->Add(calib->fHisMaterialRPhi);
127   }
128   return 0;
129 }
130
131
132
133 void AliTPCcalibMaterial::Process(AliESDEvent *event){
134   //
135   //
136   //
137   const Int_t kMinCl=40;
138   const Float_t kMinRatio=0.7;
139   const Float_t kMaxS=0.05;
140   const Float_t kMinDist=5;
141   const Double_t kStep=1.;
142   if (!event) return;
143   //  TTreeSRedirector * cstream =  GetDebugStreamer();
144   //
145   if (!fHisMaterial){
146     MakeHisto();
147   }
148   
149   //  
150   // fill histogram of track prolongations
151   Float_t dca[2];
152   Int_t ntracks = event->GetNumberOfTracks();
153
154   for (Int_t itrack=0; itrack<ntracks; itrack++){
155     AliESDtrack *track=event->GetTrack(itrack);
156     if (!track) continue;
157     if (track->GetTPCNcls()<=kMinCl) continue;
158     if ((1.+track->GetTPCNcls())/(1.+track->GetTPCNclsF())<=kMinRatio) continue;
159     if ((1.+track->GetTPCnclsS())/(1.+track->GetTPCNcls())>kMaxS) continue;
160     if (!track->GetInnerParam()) continue;
161     if (track->GetKinkIndex(0)!=0) continue;
162     //
163     track->GetImpactParameters(dca[0],dca[1]);
164     if (TMath::Abs(dca[0])<kMinDist && TMath::Abs(dca[1])<kMinDist) continue;
165     AliExternalTrackParam param(*(track->GetInnerParam()));
166     if (!AliTracker::PropagateTrackTo(&param,90,0.0005,10,kTRUE)) continue;
167     Double_t x[5]={0,0,0,TMath::Sqrt(TMath::Abs(param.GetP()))*param.GetSign(),TMath::Sqrt(TMath::Abs(track->GetTPCsignal()))};
168     //
169     //
170     for (Float_t radius=90; radius>0; radius-=kStep){
171       if (!AliTracker::PropagateTrackTo(&param,radius,0.0005,kStep*0.5,kTRUE)) break;
172       if (TMath::Abs(param.GetSnp())>0.8) break;
173       param.GetXYZ(x);
174       Double_t weight=1./TMath::Sqrt(1.+param.GetSnp()*param.GetSnp()+param.GetTgl()*param.GetTgl());
175       fHisMaterial->Fill(x,weight);    
176       Double_t r = TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
177       Double_t phi = TMath::ATan2(x[1],x[0]);
178       x[0]=r;
179       x[1]=phi;
180       fHisMaterialRPhi->Fill(x,weight);
181     }
182   }
183 }
184
185 THnSparse *AliTPCcalibMaterial::MakeHisto(){
186   //
187   // Make track prolongation histogram
188   // 
189   //
190   //                    gX       gY     gz   p       dEdx
191   Int_t    bins[5]   = {100,    100,   300,  40,   100};
192   Double_t xmin[5]   = {-100,  -100,  -300,  -2,   5};
193   Double_t xmax[5]   = {100,    100,   300,   2,   33};
194   TString  axisName[5]={
195     "gx",
196     "gy",
197     "gz",
198     "p",
199     "dedx"
200   };
201   TString  axisTitle[5]={
202     "x    (cm)",
203     "y    (cm)",
204     "z    (cm)",
205     "p    (GeV)",
206     "dedx (a.u)"
207   };
208
209   Int_t    binsR[5]   = {30,    360,     300,  40,   100};
210   Double_t xminR[5]   = { 0,    -3.14,  -300,  -2,   5};
211   Double_t xmaxR[5]   = {30,    3.14,    300,   2,   33};
212   TString  axisNameR[5]={
213     "r",
214     "rphi",
215     "z",
216     "p",
217     "dedx"
218   };
219   TString  axisTitleR[5]={
220     "r    (cm)",
221     "rphi    (cm)",
222     "z    (cm)",
223     "p    (GeV)",
224     "dedx (a.u)"
225   };
226
227   THnSparse *sparse = new THnSparseF("his_Material", "His Material", 5, bins, xmin, xmax);
228   THnSparse *sparseR = new THnSparseF("his_MaterialRPhi", "His Material Rphi", 5, binsR, xminR, xmaxR);
229   for (Int_t iaxis=0; iaxis<5; iaxis++){
230     sparse->GetAxis(iaxis)->SetName(axisName[iaxis]);
231     sparse->GetAxis(iaxis)->SetTitle(axisTitle[iaxis]);
232     sparseR->GetAxis(iaxis)->SetName(axisNameR[iaxis]);
233     sparseR->GetAxis(iaxis)->SetTitle(axisTitleR[iaxis]);
234   }
235   fHisMaterial=sparse;
236   fHisMaterialRPhi=sparseR;
237   return sparse;
238 }