Some of the coding violations corrected
[u/mrichter/AliRoot.git] / TPC / AliESDcosmic.cxx
1 #include "TMath.h"
2 #include "TArrayI.h"
3 #include "TDatabasePDG.h"
4 #include "TTreeStream.h"
5 #include "TVectorD.h"
6
7
8 #include "AliESDEvent.h"
9 #include "AliESDtrack.h"
10 #include "AliTracker.h"
11 #include "AliESDcosmic.h"
12
13 ClassImp(AliESDcosmic)
14
15 AliESDcosmic::AliESDcosmic():
16     TObject(),
17     fESD(0),
18     fTracks(0),
19     fTracksAcorde(0),
20     fPair(0),
21     fDebugStreamer(0)
22 {
23   //
24   //
25   //
26 }
27
28
29 AliESDcosmic::~AliESDcosmic(){
30   //
31   //
32   //
33   delete fPair;
34   delete fTracksAcorde;
35   
36 }
37
38 void AliESDcosmic::ProcessEvent(AliESDEvent* event){
39   //
40   //
41   //
42   fESD=event;
43   TString  string =  fESD->GetFiredTriggerClasses();
44   if (!(string.Contains("ASL"))) return;
45   //
46   const Float_t kCutMinDir=0.95;
47   const Float_t kMaxD=1000.;
48   Int_t ntracks=event->GetNumberOfTracks(); 
49   if (ntracks<=0) return;
50   if (!fTracks) fTracks= new TClonesArray("AliESDtrack",0);
51   fPair = new TArrayI(ntracks);
52
53   for (Int_t i=0;i<ntracks;++i) {
54     AliESDtrack *track0 = event->GetTrack(i); 
55     AliESDtrack *trackPair=0;
56     (*fPair)[i]=-1;
57     // track0 - choosen upper part
58     if (!track0) continue;
59     if (!track0->GetOuterParam()) continue; 
60     Double_t dir0[3];
61     track0->GetDirection(dir0); 
62     Float_t minDist=kMaxD;
63     //
64     for (Int_t j=i;j<ntracks;++j) {
65       AliESDtrack *track1 = event->GetTrack(j);   
66       //track 1 lower part
67       if (!track1) continue;
68       if (!track1->GetOuterParam()) continue;
69       Double_t dir1[3];
70       track1->GetDirection(dir1);
71       Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
72       if (TMath::Abs(dir)<kCutMinDir) continue;               // direction vector product 
73       //
74       // calculate distance
75       Float_t dy = (track0->GetY()+track1->GetY());
76       Float_t sy2 = track0->GetSigmaY2()+track1->GetSigmaY2();
77       Float_t dphi = (track0->GetAlpha()-track1->GetAlpha()-TMath::Pi());
78       Float_t sphi2  = track0->GetSigmaSnp2()+track1->GetSigmaSnp2();
79       Float_t dtheta = (track0->GetTgl()-track1->GetTgl());
80       Float_t stheta2 = track0->GetSigmaTgl2()+track1->GetSigmaTgl2();
81       Float_t normDist = TMath::Sqrt(dy*dy/sy2+dphi*dphi/sphi2+dtheta*dtheta/stheta2);
82       //
83       if (normDist>minDist) continue;    
84       minDist = normDist;
85       trackPair=track1;
86       (*fPair)[i]=j;
87     }
88     //
89   }
90   PropagateToAcorde();
91 }
92
93
94 void  AliESDcosmic::PropagateToAcorde(){
95   //
96   //
97   //
98   const Double_t kRL3=510;   // radius of L3 magnet
99   const Double_t kxAcorde=850.;
100   //
101   if (!fTracksAcorde) fTracksAcorde= new TClonesArray("AliExternalTrackParam",0);
102   Int_t ntracks=fESD->GetNumberOfTracks(); 
103   Int_t counter=0;
104   //
105   for (Int_t i=0; i<ntracks;i++){
106     Int_t index = i;
107     AliESDtrack * upperTrack = fESD->GetTrack(i);
108     if (upperTrack->GetOuterParam()==0) continue;
109     Double_t gxyz[3];
110     upperTrack->GetOuterParam()->GetXYZ(gxyz);
111     if ((*fPair)[i]>0){
112       AliESDtrack * track2 = fESD->GetTrack((*fPair)[i]);
113       if (track2->GetOuterParam()){
114         Double_t gxyz2[3];
115         track2->GetOuterParam()->GetXYZ(gxyz2);
116         if (gxyz2[1]>gxyz[1]) {
117           upperTrack=track2;
118           index=(*fPair)[i];
119         }
120       }
121     }
122     //
123     AliExternalTrackParam *upper = (AliExternalTrackParam *)(upperTrack->GetOuterParam()->Clone());
124     Bool_t isOK = upper->PropagateTo(kRL3,fESD->GetMagneticField());
125     upper->GetXYZ(gxyz);
126     if (gxyz[1]<0) continue;
127     for (Int_t iter=0; iter<20;iter++){
128       upper->GetXYZ(gxyz);
129       Double_t galpha = TMath::ATan2(gxyz[1],gxyz[0]);
130       Double_t alpha=galpha;
131       galpha*=180/TMath::Pi();
132       if (iter>1){
133         if (galpha<45.)   alpha = TMath::Pi()/8;
134         if (galpha>135.) alpha =  TMath::Pi()*(1-1/8.);
135         if (galpha>45.&&galpha<135.) alpha = TMath::Pi()/2.;
136       }
137       if (isOK) upper->Rotate(alpha);
138       if (isOK) isOK = upper->PropagateTo(kxAcorde,0);
139     }
140     if (isOK) {
141       new ((*fTracksAcorde)[index]) AliExternalTrackParam(*upper);
142       counter++;
143     }    
144   }
145 }
146
147 void AliESDcosmic::DumpToTree(){
148   //
149   //
150   //
151   TTreeSRedirector * cstream = fDebugStreamer;
152   if (!cstream) return;
153   if (!fESD) return;
154   if (!fTracksAcorde) return;
155   Int_t ntracks0 =fESD->GetNumberOfTracks(); 
156   Int_t ntracks  = fTracksAcorde->GetEntries();
157   Float_t mag = fESD->GetMagneticField();
158   Int_t   run = fESD->GetRunNumber();
159   Int_t   event = fESD->GetEventNumberInFile();
160   AliESDHeader* header = fESD->GetHeader();
161   
162   (*cstream)<<"eventInfo"<<
163     "run="<<run<<
164     "event="<<event<<
165     "header="<<header<<
166     "mag="<<mag<<
167     "ntracks0="<<ntracks0<<
168     "ntracks="<<ntracks<<    
169     "\n";
170     
171   for (Int_t i=0;i<ntracks;i++){
172     if (!fTracksAcorde->At(i)) continue;
173     TVectorD gxyz(3);
174     TVectorD gpxyz(3);
175     AliExternalTrackParam * param = (AliExternalTrackParam *)fTracksAcorde->At(i);
176     param->GetXYZ(gxyz.GetMatrixArray());
177     param->GetPxPyPz(gpxyz.GetMatrixArray());
178     (*cstream) << "esdCosmic" <<  
179       "run="<<run<<
180       "event="<<event<<
181       "header="<<header<<
182       "mag="<<mag<<
183       "ntracks0="<<ntracks0<<
184       "ntracks="<<ntracks<<      
185       "tr.="<<param<<
186       "p.="<<&gxyz<<
187       "m.="<<&gpxyz<<
188       "\n";
189   }
190     
191 }
192