Move ITS-related tasks to ITS dir
[u/mrichter/AliRoot.git] / PWG1 / ITS / AliTrackMatchingTPCITSCosmics.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2009, 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 /////////////////////////////////////////////////////////////
17 //
18 // AliAnalysisTask to study the TPC-ITS track matching
19 //
20 // Author: A.Dainese, andrea.dainese@pd.infn.it
21 /////////////////////////////////////////////////////////////
22
23 #include <TTree.h>
24 #include <TChain.h>
25 #include <TNtuple.h>
26 #include <TList.h>
27 #include <TH1F.h>
28 #include <TH2F.h>
29 #include <TVector3.h>
30 #include <TGeoManager.h>
31
32 #include "AliLog.h"
33 #include "AliGeomManager.h"
34 #include "AliTrackPointArray.h"
35 #include "AliESDInputHandler.h"
36 #include "AliESDVertex.h"
37 #include "AliESDtrack.h"
38 #include "AliESDEvent.h"
39 #include "AliESDfriend.h"
40 #include "AliAnalysisTask.h"
41 #include "AliAnalysisManager.h"
42 #include "AliExternalTrackParam.h"
43 #include "AliTrackMatchingTPCITSCosmics.h"
44
45
46 ClassImp(AliTrackMatchingTPCITSCosmics)
47
48
49 //________________________________________________________________________
50 AliTrackMatchingTPCITSCosmics::AliTrackMatchingTPCITSCosmics(const char *name):
51 AliAnalysisTask(name,"task"),
52 fOnlySPDFO(kFALSE),
53 fReadHLTESD(kFALSE),
54 fGeometryFileName("geometry.root"),
55 fESD(0),
56 fList(0),
57 fHistEvCount(0),
58 fntTrks(0),
59 fntPairs(0),
60 fntITSPairs(0)
61 {
62   // Constructor
63
64   // Define input and output slots here
65   // Input slot #0 works with a TChain
66   DefineInput(0, TChain::Class());
67   // Output slot #0 writes into a TList
68   DefineOutput(0,TList::Class());  //My private output
69 }
70
71 //________________________________________________________________________
72 AliTrackMatchingTPCITSCosmics::~AliTrackMatchingTPCITSCosmics()
73 {
74   // Destructor
75   if (fList) {
76     delete fList;
77     fList = 0;
78   }
79   if (fHistEvCount) {
80     delete fHistEvCount;
81     fHistEvCount = 0;
82   }
83   if (fntTrks) {
84     delete fntTrks;
85     fntTrks = 0;
86   }
87   if (fntPairs) {
88     delete fntPairs;
89     fntPairs = 0;
90   }
91   if (fntITSPairs) {
92     delete fntITSPairs;
93     fntITSPairs = 0;
94   }
95 }  
96
97 //________________________________________________________________________
98 void AliTrackMatchingTPCITSCosmics::ConnectInputData(Option_t *) 
99 {
100   // Connect ESD
101   // Called once
102
103   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
104   if(!tree) {
105     printf("ERROR: Could not read chain from input slot 0\n");
106   } else {
107     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
108     
109     if(!esdH) {
110       printf("ERROR: Could not get ESDInputHandler\n");
111     } else {
112       fESD = esdH->GetEvent();
113     }
114   }
115   
116   return;
117 }
118
119 //________________________________________________________________________
120 void AliTrackMatchingTPCITSCosmics::Init()
121 {
122   // Initialization
123
124   return;
125 }
126
127 //________________________________________________________________________
128 void AliTrackMatchingTPCITSCosmics::CreateOutputObjects()
129 {
130   // Create the output container
131   //
132
133   // load the geometry  
134   if(!gGeoManager) {    
135     AliGeomManager::LoadGeometry(fGeometryFileName.Data());
136     if(!gGeoManager) { 
137       printf("AliTrackMatchingTPCITSCosmics::CreateOutputObjects(): no geometry loaded \n");
138       return;
139     }
140   }
141
142   // Several histograms are more conveniently managed in a TList
143   fList = new TList();
144   fList->SetOwner();
145
146   fHistEvCount = new TH1F("fHistEvCount","0: all, 1: SPDFO, 2: vtx, 3: TPCtrks, 4: ITStrks, 5: HLT, 6: HLT&TPCtrks, 7: HLT&ITStrks",11,-0.5,10.5);
147   fList->Add(fHistEvCount);
148
149   fntTrks = new TNtuple("fntTrks","TPC-ITS matching","ptTPC:nClsTPC:nClsITSSA:nClsITSMI:SSD1:SSD2:phi:z:dx:dy:dz:drphi:dphi:dtgl");
150   fList->Add(fntTrks);
151   fntPairs = new TNtuple("fntPairs","pairs at vertex","nClsITSSAin:nClsITSSAout:nClsITSMIin:nClsITSMIout:dxyITSSA:dzITSSA:dxyITSMI:dzITSMI:ptITSSAin:ptITSSAout:ptITSMIin:ptITSMIout:sigmad0MIin:sigmad0MIout");
152   fList->Add(fntPairs);
153   fntITSPairs = new TNtuple("fntITSPairs","pairs at vertex","nClsITSSAin:nClsITSSAout:dxyITSSA:dzITSSA:ptITSSAin:ptITSSAout:d0mu");
154   fList->Add(fntITSPairs);
155
156   return;
157 }
158
159 //________________________________________________________________________
160 void AliTrackMatchingTPCITSCosmics::Exec(Option_t */*option*/)
161 {
162   // Execute analysis for current event:
163   // write ITS AliTrackPoints for selected tracks to fspTree
164   
165   // check the geometry  
166   if(!gGeoManager) { 
167     printf("AliTrackMatchingTPCITSCosmics::Exec(): no geometry loaded \n");
168     return;
169   }
170
171   if(!fESD) {
172     printf("AliTrackMatchingTPCITSCosmics::Exec(): no ESD \n");
173     return;
174   } 
175
176
177   fHistEvCount->Fill(0);
178   PostData(0,fList);
179
180   // check if event is triggered by HLT
181   Bool_t hltTrigg=kFALSE;
182   if(fReadHLTESD) {
183     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
184     AliESDEvent *hltESD = esdH->GetHLTEvent();
185     if(!hltESD) {
186       printf("AliTrackMatchingTPCITSCosmics::Exec(): no HLT ESD \n");
187       return;
188     } 
189     if(hltESD->IsHLTTriggerFired()) {fHistEvCount->Fill(5);hltTrigg=kTRUE;}
190   }
191
192   
193   TString triggeredClass = fESD->GetFiredTriggerClasses(); 
194   if(triggeredClass.Contains("C0SCO-ABCE-NOPF-CENT")) {
195     fHistEvCount->Fill(1);
196   } else {
197     if(fOnlySPDFO) {  PostData(0,fList); return; }
198   }
199   
200   Int_t ntracks = fESD->GetNumberOfTracks();
201   printf("CONTR: %d\n",fESD->GetVertex()->GetNContributors());
202   if(fESD->GetVertex()->GetNContributors()>=0) {
203     fHistEvCount->Fill(2);
204   } 
205   
206   Int_t nTrksITSSA=0,nTrksTPC=0,nTrksTPCITS=0;
207   Int_t idxITSSA[10000],idxTPC[10000];
208   Int_t idxITSSAin=-1,idxITSSAout=-1;
209   Int_t idxTPCin=-1,idxTPCout=-1;
210   for(Int_t itr=0; itr<ntracks; itr++) {
211     AliESDtrack *t = fESD->GetTrack(itr);
212     Int_t nClsTPC = t->GetNcls(1);
213     Int_t nClsITS = t->GetNcls(0);
214     if(nClsITS>=2 && nClsTPC==0) { // ITS SA
215       idxITSSA[nTrksITSSA]=itr;
216       nTrksITSSA++;
217       //printf("%d\n",itr);
218     }
219     if(nClsTPC>=50) {              // TPC
220       idxTPC[nTrksTPC]=itr;
221       nTrksTPC++;
222     }
223     if(nClsTPC>=50 && nClsITS>=2) {// TPC+ITS
224       nTrksTPCITS++;
225       /*printf("  --> TPC+ITS: %d + %d (%d,%d,%d,%d,%d,%d) pt = %f\n",nClsTPC,nClsITS,
226         (Int_t)(t->HasPointOnITSLayer(0)),
227         (Int_t)(t->HasPointOnITSLayer(1)),
228         (Int_t)(t->HasPointOnITSLayer(2)),
229         (Int_t)(t->HasPointOnITSLayer(3)),
230         (Int_t)(t->HasPointOnITSLayer(4)),
231         (Int_t)(t->HasPointOnITSLayer(5)),
232         t->Pt());*/
233     }
234   }
235   
236   //printf("nTrksTPC %d  nTrksITSSA %d\n",nTrksTPC,nTrksITSSA);
237   if(nTrksITSSA>0) fHistEvCount->Fill(4);
238   if(hltTrigg && nTrksITSSA>0) fHistEvCount->Fill(7);
239   if(nTrksTPC>0)   fHistEvCount->Fill(3);
240   if(hltTrigg && nTrksTPC>0)   fHistEvCount->Fill(6);
241   if(nTrksTPCITS>0)   fHistEvCount->Fill(8);
242   if(hltTrigg && nTrksTPCITS>0)   fHistEvCount->Fill(9);
243   
244   if(nTrksTPC>2 || nTrksTPC==0 || nTrksITSSA>2 || nTrksITSSA==0) { PostData(0,fList); return; }
245   
246   for(Int_t itr=0; itr<nTrksTPC; itr++) {
247     AliESDtrack *t = fESD->GetTrack(idxTPC[itr]);
248     if(t->Py()>0) idxTPCin=idxTPC[itr];
249     if(t->Py()<0) idxTPCout=idxTPC[itr];
250   }
251   
252   
253   for(Int_t itr=0; itr<nTrksITSSA; itr++) {
254     AliESDtrack *t = fESD->GetTrack(idxITSSA[itr]);
255     if(t->Py()>0) idxITSSAin=idxITSSA[itr];
256     if(t->Py()<0) idxITSSAout=idxITSSA[itr];
257   }
258   
259   Double_t xyzITS[3],xyzTPC[3];
260
261   // analysis for inward track
262   if(idxITSSAin>-1 && idxTPCin>-1) { 
263     //printf(" %d\n",idxITSSAin);
264     AliESDtrack *tITS = fESD->GetTrack(idxITSSAin);
265     const AliExternalTrackParam *outerITS = tITS->GetOuterParam();
266     xyzITS[0]=-1000.;xyzITS[1]=-1000.;xyzITS[2]=-1000.;
267     if(outerITS) outerITS->GetXYZAt(50.,fESD->GetMagneticField(),xyzITS);
268     AliESDtrack *tTPC = fESD->GetTrack(idxTPCin);
269     AliExternalTrackParam *innerTPC = new AliExternalTrackParam(*(tTPC->GetInnerParam()));
270     innerTPC->Rotate(tITS->GetAlpha());
271     innerTPC->GetXYZAt(50.,fESD->GetMagneticField(),xyzTPC);
272     fntTrks->Fill(innerTPC->Pt(),
273                  tTPC->GetNcls(1),
274                  tITS->GetNcls(0),
275                  tTPC->GetNcls(0),
276                  tITS->HasPointOnITSLayer(4),
277                  tITS->HasPointOnITSLayer(5),
278                  TMath::ATan2(xyzITS[1],xyzITS[0]),
279                  xyzITS[2],
280                  xyzITS[0]-xyzTPC[0],
281                  xyzITS[1]-xyzTPC[1],
282                  xyzITS[2]-xyzTPC[2],
283                  0.5*(TMath::Sqrt(xyzITS[0]*xyzITS[0]+xyzITS[1]*xyzITS[1])+TMath::Sqrt(xyzTPC[0]*xyzTPC[0]+xyzTPC[1]*xyzTPC[1]))*(TMath::ATan2(xyzITS[1],xyzITS[0])-TMath::ATan2(xyzTPC[1],xyzTPC[0])),
284                  TMath::ATan2(tITS->Py(),tITS->Px())-TMath::ATan2(innerTPC->Py(),innerTPC->Px()),
285                  tITS->Pz()/tITS->Pt()-innerTPC->Pz()/innerTPC->Pt());
286     delete innerTPC; innerTPC=NULL;
287   }
288
289   // analysis for outward track
290   if(idxITSSAout>-1 && idxTPCout>-1) { 
291     //printf(" %d\n",idxITSSAout);
292     AliESDtrack *tITS = fESD->GetTrack(idxITSSAout);
293     const AliExternalTrackParam *outerITS = tITS->GetOuterParam();
294     xyzITS[0]=-1000.;xyzITS[1]=-1000.;xyzITS[2]=-1000.;
295     if(outerITS) outerITS->GetXYZAt(50.,fESD->GetMagneticField(),xyzITS);
296     AliESDtrack *tTPC = fESD->GetTrack(idxTPCout);
297     AliExternalTrackParam *innerTPC = new AliExternalTrackParam(*(tTPC->GetInnerParam()));
298     innerTPC->Rotate(tITS->GetAlpha());
299     innerTPC->GetXYZAt(50.,fESD->GetMagneticField(),xyzTPC);
300     fntTrks->Fill(innerTPC->Pt(),
301                  tTPC->GetNcls(1),
302                  tITS->GetNcls(0),
303                  tTPC->GetNcls(0),
304                  tITS->HasPointOnITSLayer(4),
305                  tITS->HasPointOnITSLayer(5),
306                  TMath::ATan2(xyzITS[1],xyzITS[0]),
307                  xyzITS[2],
308                  xyzITS[0]-xyzTPC[0],
309                  xyzITS[1]-xyzTPC[1],
310                  xyzITS[2]-xyzTPC[2],
311                  0.5*(TMath::Sqrt(xyzITS[0]*xyzITS[0]+xyzITS[1]*xyzITS[1])+TMath::Sqrt(xyzTPC[0]*xyzTPC[0]+xyzTPC[1]*xyzTPC[1]))*(TMath::ATan2(xyzITS[1],xyzITS[0])-TMath::ATan2(xyzTPC[1],xyzTPC[0])),
312                  TMath::ATan2(tITS->Py(),tITS->Px())-TMath::ATan2(innerTPC->Py(),innerTPC->Px()),
313                  tITS->Pz()/tITS->Pt()-innerTPC->Pz()/innerTPC->Pt()
314                  );
315     delete innerTPC; innerTPC=NULL;
316   }
317
318   // inward-outward track ITS only
319   if(idxITSSAin>-1 && idxITSSAout>-1) {
320     AliESDtrack *tITSin = fESD->GetTrack(idxITSSAin);
321     AliESDtrack *tITSout = fESD->GetTrack(idxITSSAout);
322     if(tITSin->HasPointOnITSLayer(0) && 
323        tITSin->HasPointOnITSLayer(1) &&
324        tITSout->HasPointOnITSLayer(0) &&
325        tITSout->HasPointOnITSLayer(1)) {
326       Double_t alpha = TMath::ATan2(tITSin->Py(),tITSin->Px());
327       tITSin->Propagate(alpha,0.,fESD->GetMagneticField());
328       tITSout->Propagate(alpha,0.,fESD->GetMagneticField());
329       Float_t d0[2],z0[2];
330       d0[0] = tITSin->GetY();
331       z0[0] = tITSin->GetZ();
332       d0[1] = tITSout->GetY();
333       z0[1] = tITSout->GetZ();
334       Float_t dxyITS = -(d0[0]-d0[1]);
335       Float_t dzITS = z0[0]-z0[1];
336       fntITSPairs->Fill(tITSin->GetNcls(0),tITSout->GetNcls(0),dxyITS,dzITS,tITSin->Pt(),tITSout->Pt(),TMath::Abs(d0[0]));
337     }
338   }
339
340   // inward-outward track TPC+ITS
341   if(idxITSSAin>-1 && idxTPCin>-1 && idxITSSAout>-1 && idxTPCout>-1) {
342     AliESDtrack *tTPCin = fESD->GetTrack(idxTPCin);
343     AliESDtrack *tITSin = fESD->GetTrack(idxITSSAin);
344     AliESDtrack *tTPCout = fESD->GetTrack(idxTPCout);
345     AliESDtrack *tITSout = fESD->GetTrack(idxITSSAout);
346     if(tTPCin->HasPointOnITSLayer(0) && 
347        tTPCin->HasPointOnITSLayer(1) &&
348        tITSin->HasPointOnITSLayer(0) && 
349        tITSin->HasPointOnITSLayer(1) && 
350        tTPCout->HasPointOnITSLayer(0) &&
351        tTPCout->HasPointOnITSLayer(1) &&
352        tITSout->HasPointOnITSLayer(0) &&
353        tITSout->HasPointOnITSLayer(1)) {
354       Double_t alpha = TMath::ATan2(tTPCin->Py(),tTPCin->Px());
355       tTPCin->Propagate(alpha,0.,fESD->GetMagneticField());
356       tTPCout->Propagate(alpha,0.,fESD->GetMagneticField());
357       Float_t d0[2],z0[2];
358       d0[0] = tTPCin->GetY();
359       z0[0] = tTPCin->GetZ();
360       d0[1] = tTPCout->GetY();
361       z0[1] = tTPCout->GetZ();
362       Float_t dxyTPC = -(d0[0]-d0[1]);
363       Float_t dzTPC = z0[0]-z0[1];
364       alpha = TMath::ATan2(tITSin->Py(),tITSin->Px());
365       tITSin->Propagate(alpha,0.,fESD->GetMagneticField());
366       tITSout->Propagate(alpha,0.,fESD->GetMagneticField());
367       d0[0] = tITSin->GetY();
368       z0[0] = tITSin->GetZ();
369       d0[1] = tITSout->GetY();
370       z0[1] = tITSout->GetZ();
371       Float_t dxyITS = -(d0[0]-d0[1]);
372       Float_t dzITS = z0[0]-z0[1];
373       fntPairs->Fill(tITSin->GetNcls(0),tITSout->GetNcls(0),tTPCin->GetNcls(0),tTPCout->GetNcls(0),dxyITS,dzITS,dxyTPC,dzTPC,tITSin->Pt(),tITSout->Pt(),tTPCin->Pt(),tTPCout->Pt(),tTPCin->GetSigmaY2(),tTPCout->GetSigmaY2());
374     }
375   }
376
377   PostData(0,fList);
378
379   return;
380 }
381
382 //________________________________________________________________________
383 void AliTrackMatchingTPCITSCosmics::Terminate(Option_t */*option*/)
384 {
385   // Terminate analysis
386   //
387   AliDebug(2,"AliTrackMatchingTPCITSCosmics: Terminate() \n");
388
389   fList = dynamic_cast<TList*> (GetOutputData(0));
390   if (!fList) {     
391     printf("ERROR: fList not available\n");
392     return;
393   }
394
395   return;
396 }
397