Fixing coding conventions
[u/mrichter/AliRoot.git] / ACORDE / macros / AnalysisMacros / Local / AliAnalysisTaskAcorde.cxx
CommitLineData
b5f2577c 1/////////////////////////////////////////////////////////////////////////////
2//
3// AliAnalysisTaskAcorde class
4//
5// Description:
6//
7// Reads the information of ACORDE-ESD
8// Also it is implemented a simple matching between tracks
9// to look for the extrapolation of them to ACORDE modules
10//
6399c648 11// Create some fHistograms and one tree with the information
b5f2577c 12// of the matching tracks
13//
14// Author: Mario Rodríguez Cahuantzi
15// <mario.rocah@gmail.com>
16// <mrodrigu@mail.cern.ch>
17//
18// Created: June 30th. 2010 @ FCFM - BUAP, Puebla, MX
19// Last update: created
20//
21/////////////////////////////////////////////////////////////////////////////
22#include "TChain.h"
23#include "TTree.h"
24#include "TH2F.h"
25#include "TH1F.h"
26#include "AliAnalysisTask.h"
27#include "AliAnalysisManager.h"
28
29#include "AliESDEvent.h"
30#include "AliESDInputHandler.h"
31
32#include "AliAnalysisTaskAcorde.h"
33
34#include "TMath.h"
35#include "TArrayI.h"
36#include "TDatabasePDG.h"
37#include "TVectorD.h"
38
39#include "AliESDtrack.h"
40#include "AliTracker.h"
41#include "TFile.h"
42#include "TVectorD.h"
43#include "TStyle.h"
44#include "TCanvas.h"
45#include "TLegend.h"
46ClassImp(AliAnalysisTaskAcorde)
47//________________________________________________________________________
48AliAnalysisTaskAcorde::AliAnalysisTaskAcorde(const char *name)
49 : AliAnalysisTask(name, ""),
50 fESD(0),
6399c648 51 fCosmicTree(0),
52 fnTracks(0),
53 fNMatchTracks(0),
54 fkMinTPCclusters(30),
55 fkMinTrackDist(1000.),
56 fkMinCutDir(0.95),
57 fXAco(0),
58 fYAco(0),
59 fZAco(0),
60 fTrigger(0),
61 fActiveTriggerDetector(0),
62 fNSingleTracks(0),
63 fNMatchedTracks(0),
64 fHisto(0),
65 fNTracks(0),
66 fAcordeHitsAll(0),
67 fAcordeMultiAll(0),
68 fAcordeHitsTPC(0),
69 fAcordeMultiTPC(0),
70 fNTracksMatched(0),
b5f2577c 71 fTracksToAcorde(0)
72
73{
74 // Constructor of class
75
76 DefineInput(0, TChain::Class());
77 DefineOutput(0,TTree::Class());
78 DefineOutput(1,TList::Class());
79}
80
81//________________________________________________________________________
82AliAnalysisTaskAcorde::~AliAnalysisTaskAcorde()
83{
84 // destructor --> avoid watchdog mails?
85
86 delete fESD;
6399c648 87 delete fHisto;
88 delete fCosmicTree;
b5f2577c 89}
90//________________________________________________________________________
91void AliAnalysisTaskAcorde::ConnectInputData(Option_t *)
92{
93 // Connect ESD or AOD here
94 // Called once
95
96 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
97 if (!tree)
98 {
99 Printf("ERROR: Could not read chain from input slot 0");
100 }else
101 {
102
103 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
104 if (!esdH)
105 {
106 Printf("ERROR: Could not get ESDInputHandler");
107 }else fESD = esdH->GetEvent();
108 }
109}
110
111//________________________________________________________________________
112void AliAnalysisTaskAcorde::CreateOutputObjects()
113{
114
6399c648 115 // creates one TList with some fHistograms
b5f2577c 116
6399c648 117 fHisto = new TList();
b5f2577c 118
6399c648 119 fNTracksMatched = new TH1F("fNTracksMatched","fNTracksMatched from matching algorithm implemented",300,1,300);
120 fNTracks = new TH1F("fNTracks","fnTracks distribution",500,1,500);
121 fAcordeHitsAll = new TH1F("fAcordeHitsAll","Hits of ACORDE from ESD",60,-0.5,59.5);
122 fAcordeHitsTPC = new TH1F("fAcordeHitsAllTPC","Hits of ACORDE from ESD (if fNTracks>0)",61,0,60);
123 fAcordeMultiAll = new TH1F("fAcordeMultiAll","Multiplicity of ACORDE modules from ESD",60,-0.5,59.5);
124 fAcordeMultiTPC = new TH1F("fAcordeMultiAllTPC","Multiplicity of ACORDE modules from ESD (id fNTracks>0)",61,0,60);
b5f2577c 125 fTracksToAcorde = new TH2F("fTracksToAcorde","Extrapolated tracks to ACORDE x:z",1200,-600,600,1200,-600,600);
126
6399c648 127 fHisto->Add(fNTracksMatched);
128 fHisto->Add(fNTracks);
129 fHisto->Add(fAcordeHitsAll);
130 fHisto->Add(fAcordeHitsTPC);
131 fHisto->Add(fAcordeMultiAll);
132 fHisto->Add(fAcordeMultiTPC);
133 fHisto->Add(fTracksToAcorde);
b5f2577c 134
135 // Create Tree branches
136 // Called just once
137
6399c648 138 fCosmicTree = new TTree("fCosmicTree","fCosmicTreeMRC");
139 fCosmicTree->Branch("fnTracks",&fnTracks,"fnTracks/I");
140 fCosmicTree->Branch("fNMatchTracks",&fNMatchTracks,"fNMatchTracks/I");
141 fCosmicTree->Branch("fXAco",&fXAco,"fXAco/F");
142 fCosmicTree->Branch("fYAco",&fYAco,"fYAco/F");
143 fCosmicTree->Branch("fZAco",&fZAco,"fZAco/F");
144 fCosmicTree->Branch("fTrigger",&fTrigger,"fTrigger/I");
b5f2577c 145
146}
147
148
149void AliAnalysisTaskAcorde::Exec(Option_t *)
150{
151 // Main loop
152 // Called for each event
153
154 Int_t counterOfMuons = 0;
155
156
157
158 if (!fESD)
159 {
160 Printf("ERROR: fESD not available");
161 return;
162 }
163
164
165
166 // Pointer to the information of ACORDE detector
167
168 AliESDACORDE *acordeESD = fESD->GetACORDEData();
169
170 Int_t contMulti = 0;
171 Int_t contMultiTPC = 0;
172
173 for(Int_t imod=0;imod<60;imod++)
174 {
175 if (acordeESD->GetHitChannel(imod))
176 {
6399c648 177 fAcordeHitsAll->Fill(imod);
b5f2577c 178 contMulti++;
179 }
180
6399c648 181 }fAcordeMultiAll->Fill(contMulti);
b5f2577c 182
183 for(Int_t imod=0;imod<60;imod++)
184 {
185 if (acordeESD->GetHitChannel(imod))
186 {
6399c648 187 fAcordeHitsTPC->Fill(imod);
b5f2577c 188 contMultiTPC++;
189 }
6399c648 190 }fAcordeMultiTPC->Fill(contMultiTPC);
b5f2577c 191
192
6399c648 193 // Assingment of the type of Trigger to fCosmicTree
b5f2577c 194
6399c648 195 fActiveTriggerDetector = fESD->GetFiredTriggerClasses();
b5f2577c 196
6399c648 197 if (fActiveTriggerDetector.Contains("C0SCO")) fTrigger=0;
198 if (fActiveTriggerDetector.Contains("C0SH2")) fTrigger=1;
199 if (fActiveTriggerDetector.Contains("C0AMU")) fTrigger=2;
200 if (fActiveTriggerDetector.Contains("C0LSR")) fTrigger=3;
201 if (fActiveTriggerDetector.Contains("C0SCO1")) fTrigger=4;
202 if (fActiveTriggerDetector.Contains("C0OBE")) fTrigger=5;
203 if (fActiveTriggerDetector.Contains("C0PCO")) fTrigger=6;
204 if (fActiveTriggerDetector.Contains("C0OB3")) fTrigger=7;
205 if (fActiveTriggerDetector.Contains("C0OCP")) fTrigger=8;
206 if (fActiveTriggerDetector.Contains("C0ASL")) fTrigger=9;
b5f2577c 207
208
209 // Begins the matching analysis between tracks
210
6399c648 211 fnTracks = fESD->GetNumberOfTracks();
212 if (fnTracks<=0) return;
213 fNTracks->Fill(fnTracks);
b5f2577c 214
6399c648 215 fPair = new TArrayI(fnTracks);
216 fNSingleTracks=0;
217 fNMatchedTracks=0;
b5f2577c 218 Int_t muons = 0;
219 TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
6399c648 220 for (Int_t iTrack=0; iTrack<fnTracks; iTrack++)
b5f2577c 221 {
222
6399c648 223 // Special case: fnTracks == 1
224 if (fnTracks ==1)
b5f2577c 225 {
226 AliESDtrack *singleTrack = fESD->GetTrack(iTrack);
227 if (!singleTrack || !singleTrack->GetOuterParam() || !singleTrack->GetInnerParam()) continue;
6399c648 228 if (singleTrack->GetTPCNcls()<fkMinTPCclusters) continue;
b5f2577c 229 Float_t vertX = singleTrack->Xv();
230 Float_t vertZ = singleTrack->Zv();
231 if (TMath::Abs(vertX)>165 || TMath::Abs(vertZ)>245) continue; // far away from the VERTEX
232 if (singleTrack->E()< 1) continue; // less than 1 GeV
233 if (singleTrack->P()< 1) continue; // less than 1 GeV
6399c648 234 fNMatchTracks=0;
235 fNSingleTracks = fnTracks;
b5f2577c 236 muons = 1;
237
6399c648 238 } // fnTracks ==1
239 else if (fnTracks>1) // fnTracks > 1
b5f2577c 240 {
241
242 AliESDtrack *track0 = fESD->GetTrack(iTrack);
243 (*fPair)[iTrack]=-1;
244 if (!track0 || !track0->GetOuterParam() || !track0->GetInnerParam()) continue;
245
6399c648 246 Double_t trackdir0[3];
247 track0->GetDirection(trackdir0);
248 Float_t minDist = fkMinTrackDist;
b5f2577c 249
250 const AliExternalTrackParam * trackIn = track0->GetInnerParam();
251 const AliExternalTrackParam * trackOut = track0->GetOuterParam();
252
253 if (!trackIn || !trackOut) continue;
6399c648 254 if (fnTracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue; // Filter TPC-Laser
b5f2577c 255
6399c648 256 if (track0->GetTPCNcls()<fkMinTPCclusters) continue;
b5f2577c 257
6399c648 258 for (Int_t jTrack=iTrack; jTrack<fnTracks; ++jTrack)
b5f2577c 259 {
260 AliESDtrack *track1 = fESD->GetTrack(jTrack);
261 if (!track1 || !track1->GetOuterParam() || !track1->GetInnerParam()) continue;
6399c648 262 Double_t trackdir1[3];
263 track1->GetDirection(trackdir1);
264 if (track1->GetTPCNcls()<fkMinTPCclusters) continue;
265 Float_t direction = trackdir0[0]*trackdir1[0] + trackdir0[1]*trackdir1[1] + trackdir0[2]*trackdir1[2];
266 if (TMath::Abs(direction)<fkMinCutDir) continue; // direction vector product
b5f2577c 267 Float_t dvertex0[3];
268 Float_t dvertex1[3];
269
270 dvertex0[0] = track0->Xv();
271 dvertex0[1] = track0->Yv();
272 dvertex0[2] = track0->Zv();
273
274 dvertex1[0] = track1->Xv();
275 dvertex1[1] = track1->Yv();
276 dvertex1[2] = track1->Zv();
277
278 if (TMath::Abs(dvertex0[0])>165 || TMath::Abs(dvertex0[2])>245 || TMath::Abs(dvertex1[0])>165 || TMath::Abs(dvertex0[2])>245) continue;
279 // calculate the distance between track0 and track1
280
281 Float_t dy = track0->GetY()+track1->GetY();
282 Float_t sy2 = track0->GetSigmaY2()+track1->GetSigmaY2();
283 Float_t dphi = track0->GetAlpha()-track1->GetAlpha()-TMath::Pi();
284 Float_t sphi2 = track0->GetSigmaSnp2()+track1->GetSigmaSnp2();
285 Float_t dtheta = track0->GetTgl()+track1->GetTgl();
286 Float_t stheta2 = track0->GetSigmaTgl2()+track1->GetSigmaTgl2();
287 Float_t normDist = TMath::Sqrt(dy*dy/sy2+dphi+dphi/sphi2+dtheta*dtheta/stheta2);
288 if (normDist>minDist) continue;
289
290 Float_t distTracks = TMath::Sqrt((dvertex0[0]-dvertex1[0])*(dvertex0[0]-dvertex1[0]) + (dvertex0[1]-dvertex1[1])*(dvertex0[1]-dvertex1[1]) + (dvertex0[2]-dvertex1[2])*(dvertex0[2]-dvertex1[2]) );
291 if (distTracks>2.5) continue;
292
293
294 minDist = normDist;
295
296 // after all the cuts we should have only tracks that can be matched
297
298 (*fPair)[iTrack] = jTrack;
299
300
301 }
6399c648 302 } // fnTracks > 1
b5f2577c 303
304
305 } // Loop over all the tracks
306
307
308
309
310 // Matching tracks
311
6399c648 312 for (Int_t itrack=0;itrack<fnTracks;itrack++)
b5f2577c 313 {
6399c648 314 AliESDtrack *upTrack = 0;
315 AliESDtrack *downTrack = 0;
b5f2577c 316
6399c648 317 AliESDtrack *track0 = fESD->GetTrack(itrack);
318 if (!track0 || !track0->GetOuterParam() || !track0->GetInnerParam() ) continue;
b5f2577c 319
320 Double_t mxyz[3];
6399c648 321 track0->GetOuterParam()->GetXYZ(mxyz);
b5f2577c 322
6399c648 323 upTrack = track0;
b5f2577c 324
325
326 if ((*fPair)[itrack]>0)
327 {
6399c648 328 AliESDtrack *track1 = fESD->GetTrack((*fPair)[itrack]);
329 if (!track1 || !track1->GetOuterParam() || !track1->GetInnerParam()) continue;
b5f2577c 330 Float_t dvertex0[3];
331 Float_t dvertex1[3];
332
6399c648 333 dvertex0[0] = track0->Xv();
334 dvertex0[1] = track0->Yv();
335 dvertex0[2] = track0->Zv();
b5f2577c 336
6399c648 337 dvertex1[0] = track1->Xv();
338 dvertex1[1] = track1->Yv();
339 dvertex1[2] = track1->Zv();
b5f2577c 340
341 Float_t distTracks = TMath::Sqrt((dvertex0[0]-dvertex1[0])*(dvertex0[0]-dvertex1[0]) + (dvertex0[1]-dvertex1[1])*(dvertex0[1]-dvertex1[1]) + (dvertex0[2]-dvertex1[2])*(dvertex0[2]-dvertex1[2]) );
342 if (distTracks>2.5) continue;
343
6399c648 344 if (track1->GetOuterParam())
b5f2577c 345 {
346 Double_t nxyz[3];
6399c648 347 track1->GetOuterParam()->GetXYZ(nxyz);
b5f2577c 348 if (nxyz[1]>mxyz[1])
349 {
6399c648 350 upTrack = track1;
351 }else downTrack = track1;
b5f2577c 352 }
353
354
355
356 // Here we propagate the Up-Tracks to ACORDE
357
358
359 const Double_t kRL3 = 510; // radious of L3 magnet
6399c648 360 const Double_t kfXAcorde = 850.; // radios of "ACORDE" above the magnet
b5f2577c 361
362 Double_t agxyz[3];
363
364
6399c648 365 AliExternalTrackParam *upper = (AliExternalTrackParam *)(upTrack->GetOuterParam()->Clone());
b5f2577c 366
367 Bool_t isOk = upper->PropagateTo(kRL3,fESD->GetMagneticField());
368 upper->GetXYZ(agxyz);
369 if (agxyz[1]<0) continue;
370
371 upper->GetXYZ(agxyz);
372 Double_t galpha = TMath::ATan2(agxyz[1],agxyz[0]);
373 Double_t alpha = galpha;
374 galpha*=180/TMath::Pi();
375
376 if (galpha<45.) alpha = TMath::Pi()/8;
377 if (galpha>135.) alpha = TMath::Pi()*(7/8);
378 if (galpha>45 && galpha<135.) alpha = TMath::Pi()/2;
379
380 if (isOk) upper->Rotate(alpha);
381
6399c648 382 if (isOk) isOk = upper->PropagateTo(kfXAcorde,0);
b5f2577c 383
384 TVectorD rgxyz(3);
385 AliExternalTrackParam *param = (AliExternalTrackParam *)upper;
386 param->GetXYZ(rgxyz.GetMatrixArray());
6399c648 387 fXAco = rgxyz[0];
388 fYAco = rgxyz[1];
389 fZAco = rgxyz[2];
390 fTracksToAcorde->Fill(fXAco,fZAco);
b5f2577c 391 // Count how many tracks have been matched
6399c648 392 fNMatchedTracks++;
393 } else fNSingleTracks++;
b5f2577c 394 }
395
6399c648 396 if ( (fNMatchedTracks+fNSingleTracks) <= fnTracks ) muons=fNMatchedTracks+fNSingleTracks;
397 fNTracksMatched->Fill(muons);
b5f2577c 398
6399c648 399 fNMatchTracks = fNSingleTracks + fNMatchedTracks;
b5f2577c 400
401 // Fills the tree
402
6399c648 403 fCosmicTree->Fill();
b5f2577c 404
405 // Post output data.
406
6399c648 407 PostData(0,fCosmicTree);
408 PostData(1,fHisto);
b5f2577c 409}
410
411//________________________________________________________________________
412void AliAnalysisTaskAcorde::Terminate(Option_t *)
413{
414
415// AliDebug(1,"Do nothig in Terminate");
416
417 //nPart->Draw();
418 TCanvas *c1 = new TCanvas("c1","TPC-resolution for P reconstruction .... MRC");
419
420 c1->Divide(2,2);
421 c1->cd(1);
6399c648 422 fAcordeHitsAll->Draw();
b5f2577c 423 c1->cd(2)->SetLogy();
6399c648 424 fAcordeMultiAll->Draw();
b5f2577c 425 c1->cd(3)->SetLogy();
6399c648 426 //fNTracksMatched->Draw();
427 fNTracks->Draw();
b5f2577c 428 c1->cd(4);
429 fTracksToAcorde->Draw();
430
431
432
433
434}