Fixing coding conventions
[u/mrichter/AliRoot.git] / ACORDE / macros / AnalysisMacros / Local / AliAnalysisTaskAcorde.cxx
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 //
11 //      Create some fHistograms and one tree with the information
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"
46 ClassImp(AliAnalysisTaskAcorde)
47 //________________________________________________________________________
48 AliAnalysisTaskAcorde::AliAnalysisTaskAcorde(const char *name) 
49   : AliAnalysisTask(name, ""), 
50         fESD(0),
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),
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 //________________________________________________________________________
82 AliAnalysisTaskAcorde::~AliAnalysisTaskAcorde()
83 {
84         // destructor --> avoid watchdog mails?
85         
86         delete fESD;
87         delete fHisto;
88         delete fCosmicTree;
89 }
90 //________________________________________________________________________
91 void 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 //________________________________________________________________________
112 void AliAnalysisTaskAcorde::CreateOutputObjects()
113 {
114
115         // creates one TList with some fHistograms
116
117         fHisto = new TList();
118
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);
125         fTracksToAcorde = new TH2F("fTracksToAcorde","Extrapolated tracks to ACORDE x:z",1200,-600,600,1200,-600,600);
126
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);
134
135         // Create Tree branches
136         // Called just once
137
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");
145
146 }
147
148
149 void 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                 {
177                         fAcordeHitsAll->Fill(imod);
178                         contMulti++;
179                 }
180
181         }fAcordeMultiAll->Fill(contMulti);
182
183         for(Int_t imod=0;imod<60;imod++)
184         {
185                 if (acordeESD->GetHitChannel(imod))
186                 {
187                         fAcordeHitsTPC->Fill(imod);
188                         contMultiTPC++;
189                 }
190         }fAcordeMultiTPC->Fill(contMultiTPC);
191
192
193         // Assingment of the type of Trigger to fCosmicTree
194
195         fActiveTriggerDetector = fESD->GetFiredTriggerClasses();
196
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;
207         
208
209         // Begins the matching analysis between tracks
210
211         fnTracks = fESD->GetNumberOfTracks();
212         if (fnTracks<=0) return;
213         fNTracks->Fill(fnTracks);
214
215         fPair = new TArrayI(fnTracks);
216         fNSingleTracks=0;
217         fNMatchedTracks=0;
218         Int_t muons = 0;
219         TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
220         for (Int_t iTrack=0; iTrack<fnTracks; iTrack++)
221         {
222
223                 // Special case: fnTracks == 1
224                 if (fnTracks ==1)
225                 {
226                         AliESDtrack *singleTrack = fESD->GetTrack(iTrack);
227                         if (!singleTrack || !singleTrack->GetOuterParam() || !singleTrack->GetInnerParam()) continue;
228                         if (singleTrack->GetTPCNcls()<fkMinTPCclusters) continue;
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 
234                         fNMatchTracks=0;
235                         fNSingleTracks = fnTracks;
236                         muons = 1;
237
238                 } // fnTracks ==1
239                 else if (fnTracks>1) // fnTracks > 1
240                 {
241
242                         AliESDtrack *track0 = fESD->GetTrack(iTrack);
243                         (*fPair)[iTrack]=-1;
244                         if (!track0 || !track0->GetOuterParam() || !track0->GetInnerParam()) continue;
245
246                         Double_t trackdir0[3];
247                         track0->GetDirection(trackdir0);
248                         Float_t minDist = fkMinTrackDist;       
249
250                         const AliExternalTrackParam * trackIn = track0->GetInnerParam();
251                         const AliExternalTrackParam * trackOut = track0->GetOuterParam();
252
253                         if (!trackIn || !trackOut) continue;
254                         if (fnTracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue; // Filter TPC-Laser
255
256                         if (track0->GetTPCNcls()<fkMinTPCclusters) continue;
257
258                         for (Int_t jTrack=iTrack; jTrack<fnTracks; ++jTrack)
259                         {
260                                 AliESDtrack *track1 = fESD->GetTrack(jTrack);
261                                 if (!track1 || !track1->GetOuterParam() || !track1->GetInnerParam()) continue;
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
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                         }       
302                 } // fnTracks > 1
303         
304
305         } // Loop over all the tracks
306
307
308
309
310         // Matching tracks
311         
312         for (Int_t itrack=0;itrack<fnTracks;itrack++)
313         {
314                 AliESDtrack *upTrack = 0;
315                 AliESDtrack *downTrack = 0;
316
317                 AliESDtrack *track0 = fESD->GetTrack(itrack);
318                 if (!track0 || !track0->GetOuterParam() || !track0->GetInnerParam() ) continue;
319
320                 Double_t mxyz[3];
321                 track0->GetOuterParam()->GetXYZ(mxyz);
322                 
323                 upTrack = track0;
324         
325
326                 if ((*fPair)[itrack]>0)
327                 {
328                         AliESDtrack *track1 = fESD->GetTrack((*fPair)[itrack]);
329                         if (!track1 || !track1->GetOuterParam() || !track1->GetInnerParam()) continue;
330                         Float_t dvertex0[3];
331                         Float_t dvertex1[3];
332
333                         dvertex0[0] = track0->Xv();
334                         dvertex0[1] = track0->Yv();
335                         dvertex0[2] = track0->Zv();
336
337                         dvertex1[0] = track1->Xv();
338                         dvertex1[1] = track1->Yv();
339                         dvertex1[2] = track1->Zv();
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
344                         if (track1->GetOuterParam())
345                         {
346                                 Double_t nxyz[3];
347                                 track1->GetOuterParam()->GetXYZ(nxyz);
348                                 if (nxyz[1]>mxyz[1])
349                                 {
350                                         upTrack = track1;
351                                 }else downTrack = track1;
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
360                         const Double_t kfXAcorde = 850.; // radios of "ACORDE" above the magnet
361
362                         Double_t agxyz[3];
363                         
364
365                         AliExternalTrackParam *upper = (AliExternalTrackParam *)(upTrack->GetOuterParam()->Clone());
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
382                         if (isOk) isOk = upper->PropagateTo(kfXAcorde,0);
383
384                         TVectorD rgxyz(3);
385                         AliExternalTrackParam *param = (AliExternalTrackParam *)upper;
386                         param->GetXYZ(rgxyz.GetMatrixArray());
387                         fXAco = rgxyz[0];
388                         fYAco = rgxyz[1];
389                         fZAco = rgxyz[2];
390                         fTracksToAcorde->Fill(fXAco,fZAco);
391                         // Count how many tracks have been matched
392                         fNMatchedTracks++;
393                 } else fNSingleTracks++;
394         }
395
396         if ( (fNMatchedTracks+fNSingleTracks) <= fnTracks ) muons=fNMatchedTracks+fNSingleTracks;
397         fNTracksMatched->Fill(muons);
398
399         fNMatchTracks = fNSingleTracks + fNMatchedTracks;
400         
401         // Fills the tree
402
403         fCosmicTree->Fill();
404         
405         // Post output data.
406
407         PostData(0,fCosmicTree);
408         PostData(1,fHisto);
409 }      
410
411 //________________________________________________________________________
412 void 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);
422         fAcordeHitsAll->Draw();
423         c1->cd(2)->SetLogy();
424         fAcordeMultiAll->Draw();
425         c1->cd(3)->SetLogy();
426         //fNTracksMatched->Draw();
427         fNTracks->Draw();
428         c1->cd(4);
429         fTracksToAcorde->Draw();
430         
431         
432
433
434 }