Adding macros to create Calibration objects
[u/mrichter/AliRoot.git] / TRD / AliTRDAnalysisTaskTP.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 /* $Id: AliTRDAnalysisTaskTP.cxx 42548 2010-07-27 08:10:51Z cblume $ */
17
18 ////////////////////////////////////////////////////////////////////////////
19 //                                                                        //
20 //  Track point maker for the alignment of TRD                            //
21 //                                                                        //
22 //  Author:                                                               //
23 //     Sebastian Huber (S.Huber@gsi.de)                                   // 
24 //                                                                        //
25 ////////////////////////////////////////////////////////////////////////////
26
27 #include <iostream>
28
29 #include <TROOT.h>
30 #include "TH1D.h"
31 #include "TH2D.h"
32 #include "TChain.h"
33 #include <TLinearFitter.h>
34
35 #include "AliVEvent.h"
36 #include "AliESDEvent.h"
37 #include "AliInputEventHandler.h"
38 #include "AliESDInputHandler.h"
39 #include "AliAnalysisTaskSE.h"
40 #include "AliAnalysisManager.h"
41 #include "AliAlignObjParams.h"
42 #include "AliTrackPointArray.h"
43 #include "AliESDfriend.h"
44 #include "AliESDtrack.h"
45 #include "AliESDfriendTrack.h"
46 #include "TFile.h"
47 #include "TTree.h"
48 #include "TSystem.h"
49 #include "TTimeStamp.h"
50 #include "TCollection.h"
51 #include "AliLog.h"
52 #include "AliGeomManager.h"
53
54 #include "AliTRDAnalysisTaskTP.h"
55
56 ClassImp(AliTRDAnalysisTaskTP)
57
58 AliTRDAnalysisTaskTP::AliTRDAnalysisTaskTP()
59   :AliAnalysisTaskSE(),
60   fArrHists(0x0),
61   fArrTTree(0x0),
62   fTree(NULL),
63   fESD(0x0),
64   fModpop(0x0),
65   fBug(0x0),
66   fNevents(0x0),
67   fNtracks(0x0),
68   fNAcceptedTracks(0),
69   fArray(0x0),
70   fFile(0x0)
71 {
72   //
73   // Default constructor
74   //
75
76 }
77
78 //____________________________________________________________
79 AliTRDAnalysisTaskTP::AliTRDAnalysisTaskTP(const char *name) :
80   AliAnalysisTaskSE(name),
81   fArrHists(0x0),
82   fArrTTree(0x0),
83   fTree(0x0),
84   fESD(0x0),
85   fModpop(0x0),
86   fBug(0x0),
87   fNevents(0),
88   fNtracks(0),
89   fNAcceptedTracks(0),
90   fArray(0x0),
91   fFile(0x0)
92 {
93   //
94   // Constructor
95   //
96
97   DefineOutput(1, TTree::Class());
98   DefineOutput(2, TObjArray::Class());
99
100 }
101
102 //____________________________________________________________
103 AliTRDAnalysisTaskTP::~AliTRDAnalysisTaskTP() 
104 {
105   //
106   // destructor
107   //
108
109 }
110
111 //____________________________________________________________
112 void AliTRDAnalysisTaskTP::UserCreateOutputObjects() 
113 {
114   //
115   // Create the output objects
116   //
117
118   AliAlignObjParams alobj;  // initialize align obj.  
119   TString option = GetOption();
120
121   if (!fArrHists) fArrHists=new TObjArray;
122
123   fModpop = new TH2D("modpop","modpop",90,-0.5,89.5,30,-0.5,29.5);
124   fModpop->SetXTitle("module nr");
125   fModpop->SetYTitle("layer nr");
126   fArrHists->Add(fModpop);
127
128   OpenFile(1);
129   fTree = new TTree("spTree", "Tree with track space point arrays");
130   fTree->Branch("SP","AliTrackPointArray", &fArray);
131
132 }
133
134 //____________________________________________________________
135 void AliTRDAnalysisTaskTP::UserExec(Option_t *) 
136 {
137   //
138   // Exec function
139   //
140
141   //AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
142   fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
143   if(!fESD){
144   //cout << "ERROR: fESDs not available " << endl;
145   return;
146   }
147
148   fESDfriend = dynamic_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
149   if (!fESDfriend) {
150   //cout << "ERROR: fESDfriends not available " << endl;
151   return;
152   }
153
154   TLinearFitter fitter(2, "pol1");
155   TLinearFitter fitterz(2, "pol1");
156
157   // track cuts
158   int tpc = 0; // require tpc
159   int ptu = 0; // require certain pt's (magnetic field and tpc presumably on)
160
161   const Float_t kMaxDelta       = 1;
162   const Float_t kMinNcl         = 60;
163   const Float_t kMinPtLow       = 0.2;
164   const Float_t kMinNclLow      = 100;
165   const Float_t kMinPt0         = 2;
166   const Float_t kMinPt          = 0;
167   UInt_t status = AliESDtrack::kTRDrefit; 
168   if (tpc) status |= AliESDtrack::kTPCrefit; 
169
170   const Float_t kMinRadius2  = 2*2;
171   const Float_t kMaxRadius2  = 400*400;
172   const Float_t kDeadSpace   = 4;
173   const Float_t kTan = TMath::Tan(10*TMath::DegToRad());
174   Int_t ntracks = fESD->GetNumberOfTracks();
175   const AliTrackPointArray *array=0;
176   AliTrackPointArray *tmpArray = 0;
177   // trdarray contains all trd points in this event, for duplication detection
178   AliTrackPointArray *trdarray = new AliTrackPointArray(1000);
179   int ntrdarray = 0;
180
181   for (Int_t itrack=0; itrack < ntracks; itrack++) { //track loop
182
183     AliESDtrack * track = fESD->GetTrack(itrack);
184     fNtracks++;
185     if (!track) continue;
186
187     if (track->GetP() < kMinPt) continue;
188     if (track->GetKinkIndex(0)!=0) continue;
189     if (tpc) if (track->GetTPCNcls()<kMinNcl) continue;
190     if (ptu) if (track->GetP() < kMinPtLow) continue;
191     if (ptu) if (track->GetP() < kMinPt0 && track->GetTPCNcls()<kMinNclLow) continue;
192
193     AliESDfriendTrack *friendtrack = fESDfriend->GetTrack(itrack);
194     if (!friendtrack){ 
195     continue;
196     }
197
198     array = friendtrack->GetTrackPointArray();
199     if (!array) continue;
200     Int_t npoints = array->GetNPoints();
201     if (tmpArray) delete tmpArray;
202     tmpArray = new AliTrackPointArray(npoints);
203     Int_t current = 0;
204     int ntpc = 0;
205     int ntrd = 0;
206     for (Int_t ipoint=0; ipoint<npoints; ipoint++){  
207
208       AliTrackPoint p;
209       array->GetPoint(p, ipoint);
210
211       UShort_t volid = array->GetVolumeID()[ipoint];
212       Int_t iModule;
213       AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(volid,iModule);
214       if ((layer < AliGeomManager::kFirstLayer) || (layer >= AliGeomManager::kLastLayer)) continue;
215
216       if ((iModule >= AliGeomManager::LayerSize(layer)) || (iModule < 0)) continue;
217
218       Float_t r2 = p.GetX()*p.GetX()+p.GetY()*p.GetY();
219       if ( r2<kMinRadius2 || r2 > kMaxRadius2 ) continue;
220
221
222
223       if (layer>=AliGeomManager::kTPC1 && layer<=AliGeomManager::kTPC2){
224         if (p.GetCov()[0]<0 || p.GetCov()[3]<0 ||  p.GetCov()[5]<0) continue;
225
226         AliTrackPoint& plocal = p.MasterToLocal();
227         Double_t ylocal  = plocal.GetY();
228         Double_t zlocal  = plocal.GetZ();
229         Double_t xlocal  = plocal.GetX();
230         Float_t edgey = TMath::Abs(plocal.GetX()*kTan);
231         Int_t nclose=0;
232         fitter.ClearPoints();
233         fitterz.ClearPoints();
234         for (Int_t jpoint=ipoint-7; jpoint<=ipoint+7; jpoint++){
235           if (jpoint<0 || jpoint>=npoints) continue;
236           if (ipoint==jpoint) continue;
237           UShort_t volidL = array->GetVolumeID()[jpoint];
238           if (volidL!=volid) continue;
239           AliTrackPoint pc;     
240           array->GetPoint(pc, jpoint);
241           AliTrackPoint &pcl=  pc.MasterToLocal();                
242           Double_t dx = pcl.GetX()-xlocal;
243           fitter.AddPoint(&dx,pcl.GetY(),1);      
244           fitterz.AddPoint(&dx,pcl.GetZ(),1);     
245           nclose++;
246         }
247         if (nclose<6) continue;
248
249         fitter.Eval();
250         fitterz.Eval();
251         Double_t fity =fitter.GetParameter(0); 
252         Double_t fitz =fitterz.GetParameter(0); 
253         if (TMath::Abs(ylocal-fity)>kMaxDelta) continue;
254         if (TMath::Abs(zlocal-fitz)>kMaxDelta) continue;
255         if (TMath::Abs(fity)>edgey-kDeadSpace) continue;
256         ntpc++;
257       }
258
259       if (layer>=AliGeomManager::kTRD1 && layer<=AliGeomManager::kTRD6){
260
261         trdarray->AddPoint(ntrdarray++,&p);
262         fModpop->Fill(iModule,layer);
263         ntrd++;
264       }
265       tmpArray->AddPoint(current,&p);
266       current++;
267     }
268     if (ntpc < 100) continue;
269     if (ntrd < 4) continue;
270     if (fArray) delete fArray;
271     fArray = new AliTrackPointArray(current);
272     for (Int_t ipoint=0; ipoint<current; ipoint++){
273       AliTrackPoint p;
274       tmpArray->GetPoint(p, ipoint);
275       fArray->AddPoint(ipoint,&p);
276     }
277     fNAcceptedTracks++;
278     fTree->Fill();
279   }
280   delete trdarray;
281   fNevents++;
282   PostData(1,fTree);
283   PostData(2,fArrHists);
284 }
285
286 //____________________________________________________________
287 void AliTRDAnalysisTaskTP::Terminate(Option_t */*option*/) 
288 {
289   //
290   // Terminate
291   //
292
293 }
294
295