fix error in setting the number of dEdx slices to be saved in ESD
[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   // default constructor
73
74 }
75 //____________________________________________________________
76 AliTRDAnalysisTaskTP::AliTRDAnalysisTaskTP(const char *name) :
77   AliAnalysisTaskSE(name),
78   fArrHists(0x0),
79   fArrTTree(0x0),
80   fTree(0x0),
81   fESD(0x0),
82   fModpop(0x0),
83   fBug(0x0),
84   fNevents(0),
85   fNtracks(0),
86   fNAcceptedTracks(0),
87   fArray(0x0),
88   fFile(0x0)
89 {
90   // Constructor
91   DefineOutput(1, TTree::Class());
92   DefineOutput(2, TObjArray::Class());
93
94 }
95 //____________________________________________________________
96 AliTRDAnalysisTaskTP::~AliTRDAnalysisTaskTP() {
97   // destructor
98
99 }
100 //____________________________________________________________
101 void AliTRDAnalysisTaskTP::UserCreateOutputObjects() {
102
103   AliAlignObjParams alobj;  // initialize align obj.  
104   TString option = GetOption();
105
106 if (!fArrHists) fArrHists=new TObjArray;
107
108   fModpop = new TH2D("modpop","modpop",90,-0.5,89.5,30,-0.5,29.5);
109   fModpop->SetXTitle("module nr");
110   fModpop->SetYTitle("layer nr");
111   fArrHists->Add(fModpop);
112
113   OpenFile(1);
114   fTree = new TTree("spTree", "Tree with track space point arrays");
115   fTree->Branch("SP","AliTrackPointArray", &fArray);
116
117 }
118 //____________________________________________________________
119 void AliTRDAnalysisTaskTP::UserExec(Option_t *) {
120
121   //AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
122   fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
123   if(!fESD){
124   //cout << "ERROR: fESDs not available " << endl;
125   return;
126   }
127
128   fESDfriend = dynamic_cast<AliESDfriend*>(fESD->FindListObject("AliESDfriend"));
129   if (!fESDfriend) {
130   //cout << "ERROR: fESDfriends not available " << endl;
131   return;
132   }
133
134   TLinearFitter fitter(2, "pol1");
135   TLinearFitter fitterz(2, "pol1");
136
137   // track cuts
138   int tpc = 0; // require tpc
139   int ptu = 0; // require certain pt's (magnetic field and tpc presumably on)
140
141   const Float_t kMaxDelta       = 1;
142   const Float_t kMinNcl         = 60;
143   const Float_t kMinPtLow       = 0.2;
144   const Float_t kMinNclLow      = 100;
145   const Float_t kMinPt0         = 2;
146   const Float_t kMinPt          = 0;
147   UInt_t status = AliESDtrack::kTRDrefit; 
148   if (tpc) status |= AliESDtrack::kTPCrefit; 
149
150   const Float_t kMinRadius2  = 2*2;
151   const Float_t kMaxRadius2  = 400*400;
152   const Float_t kDeadSpace   = 4;
153   const Float_t kTan = TMath::Tan(10*TMath::DegToRad());
154   Int_t ntracks = fESD->GetNumberOfTracks();
155   const AliTrackPointArray *array=0;
156   AliTrackPointArray *tmpArray = 0;
157   // trdarray contains all trd points in this event, for duplication detection
158   AliTrackPointArray *trdarray = new AliTrackPointArray(1000);
159   int ntrdarray = 0;
160
161   for (Int_t itrack=0; itrack < ntracks; itrack++) { //track loop
162
163     AliESDtrack * track = fESD->GetTrack(itrack);
164     fNtracks++;
165     if (!track) continue;
166
167     if (track->GetP() < kMinPt) continue;
168     if (track->GetKinkIndex(0)!=0) continue;
169     if (tpc) if (track->GetTPCNcls()<kMinNcl) continue;
170     if (ptu) if (track->GetP() < kMinPtLow) continue;
171     if (ptu) if (track->GetP() < kMinPt0 && track->GetTPCNcls()<kMinNclLow) continue;
172
173     AliESDfriendTrack *friendtrack = fESDfriend->GetTrack(itrack);
174     if (!friendtrack){ 
175     continue;
176     }
177
178     array = friendtrack->GetTrackPointArray();
179     if (!array) continue;
180     Int_t npoints = array->GetNPoints();
181     if (tmpArray) delete tmpArray;
182     tmpArray = new AliTrackPointArray(npoints);
183     Int_t current = 0;
184     int ntpc = 0;
185     int ntrd = 0;
186     for (Int_t ipoint=0; ipoint<npoints; ipoint++){  
187
188       AliTrackPoint p;
189       array->GetPoint(p, ipoint);
190
191       UShort_t volid = array->GetVolumeID()[ipoint];
192       Int_t iModule;
193       AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(volid,iModule);
194       if ((layer < AliGeomManager::kFirstLayer) || (layer >= AliGeomManager::kLastLayer)) continue;
195
196       if ((iModule >= AliGeomManager::LayerSize(layer)) || (iModule < 0)) continue;
197
198       Float_t r2 = p.GetX()*p.GetX()+p.GetY()*p.GetY();
199       if ( r2<kMinRadius2 || r2 > kMaxRadius2 ) continue;
200
201
202
203       if (layer>=AliGeomManager::kTPC1 && layer<=AliGeomManager::kTPC2){
204         if (p.GetCov()[0]<0 || p.GetCov()[3]<0 ||  p.GetCov()[5]<0) continue;
205
206         AliTrackPoint& plocal = p.MasterToLocal();
207         Double_t ylocal  = plocal.GetY();
208         Double_t zlocal  = plocal.GetZ();
209         Double_t xlocal  = plocal.GetX();
210         Float_t edgey = TMath::Abs(plocal.GetX()*kTan);
211         Int_t nclose=0;
212         fitter.ClearPoints();
213         fitterz.ClearPoints();
214         for (Int_t jpoint=ipoint-7; jpoint<=ipoint+7; jpoint++){
215           if (jpoint<0 || jpoint>=npoints) continue;
216           if (ipoint==jpoint) continue;
217           UShort_t volidL = array->GetVolumeID()[jpoint];
218           if (volidL!=volid) continue;
219           AliTrackPoint pc;     
220           array->GetPoint(pc, jpoint);
221           AliTrackPoint &pcl=  pc.MasterToLocal();                
222           Double_t dx = pcl.GetX()-xlocal;
223           fitter.AddPoint(&dx,pcl.GetY(),1);      
224           fitterz.AddPoint(&dx,pcl.GetZ(),1);     
225           nclose++;
226         }
227         if (nclose<6) continue;
228
229         fitter.Eval();
230         fitterz.Eval();
231         Double_t fity =fitter.GetParameter(0); 
232         Double_t fitz =fitterz.GetParameter(0); 
233         if (TMath::Abs(ylocal-fity)>kMaxDelta) continue;
234         if (TMath::Abs(zlocal-fitz)>kMaxDelta) continue;
235         if (TMath::Abs(fity)>edgey-kDeadSpace) continue;
236         ntpc++;
237       }
238
239       if (layer>=AliGeomManager::kTRD1 && layer<=AliGeomManager::kTRD6){
240
241         trdarray->AddPoint(ntrdarray++,&p);
242         fModpop->Fill(iModule,layer);
243         ntrd++;
244       }
245       tmpArray->AddPoint(current,&p);
246       current++;
247     }
248     if (ntpc < 100) continue;
249     if (ntrd < 4) continue;
250     if (fArray) delete fArray;
251     fArray = new AliTrackPointArray(current);
252     for (Int_t ipoint=0; ipoint<current; ipoint++){
253       AliTrackPoint p;
254       tmpArray->GetPoint(p, ipoint);
255       fArray->AddPoint(ipoint,&p);
256     }
257     fNAcceptedTracks++;
258     fTree->Fill();
259   }
260   delete trdarray;
261   fNevents++;
262   PostData(1,fTree);
263   PostData(2,fArrHists);
264 }
265 //____________________________________________________________
266 void AliTRDAnalysisTaskTP::Terminate(Option_t */*option*/) {
267   //
268   // Terminate
269   //
270 }
271
272