make AliVVevent abstract (=0)
[u/mrichter/AliRoot.git] / HLT / global / CompareFlatESDs.C
1 /**
2  * >> Testing Macro to compare FlatESDEvents from output files <<
3  **
4  * Primary Authors : Steffen Weber
5  *
6  * Usage:
7  *  aliroot -b -l -q LoadLibs.C CompareFlatESDs.C++
8  *
9  **************************************************************************/
10
11 #if !defined(__CINT__) || defined(__MAKECINT__)
12 #include "AliESDEvent.h"
13 #include "AliESD.h"
14 #include "AliESDfriend.h"
15 #include <TFile.h>
16 #include <TTree.h>
17 #include <TSystem.h>
18 #include "./AliFlatESDEvent.h"
19 #include "./AliFlatESDTrack.h"
20 #include "./AliFlatTPCCluster.h"
21 #include "./AliFlatExternalTrackParam.h"
22 #include "Riostream.h"
23 #endif   
24
25 void CompareFlatESDs(const char* filename1="outFlatESD1.dat",const char* filename2="outFlatESD2.dat", Bool_t verbose=kFALSE) {
26   
27   
28   // Create output histograms
29   
30   
31   TString outputFilename = "$PWD/compare.root";
32   
33         cout<< "creating histograms"<<endl;
34         TH2F* hNTracks = new TH2F("nTracks","number of tracks", 100,0,100, 100,0,100);
35         TH2F* hNV0s = new TH2F("nV0s","number of V0s", 10,0,10, 10,0,10);
36         TH2F* hVtxTr = new TH2F("vtxTr","vtx Tracks", 2,0,2, 2,0,2);
37         TH2F* hVtxSPD = new TH2F("vtxSPD","vtx SPD", 2,0,2, 2,0,2);
38         TH1F * hStat = new TH1F("stat","statistics Differences",20,0,20);
39         hStat->GetXaxis()->SetBinLabel(1,"All events");
40         hStat->GetXaxis()->SetBinLabel(2,"no diffs");
41         hStat->GetXaxis()->SetBinLabel(3,"nTracks");
42         hStat->GetXaxis()->SetBinLabel(4,"nV0s");
43         hStat->GetXaxis()->SetBinLabel(5,"vtxTracks");
44         hStat->GetXaxis()->SetBinLabel(6,"vtxSPD");
45         hStat->GetXaxis()->SetBinLabel(7,"tracks->extParams");
46
47   
48   
49
50   ifstream is1(filename1, std::ifstream::binary | std::ifstream::in);
51   ifstream is2(filename2, std::ifstream::binary | std::ifstream::in);
52   if (is1 && is2 ){
53     is1.seekg (0, is1.end);
54     int length1 = is1.tellg();
55     is1.seekg (0, is1.beg);
56     char * buffer1 = new char [length1];
57     
58     std::cout << "Reading " << length1 << " characters... ";
59     
60     is1.read (buffer1,length1);
61     if (is1)
62       std::cout << "all characters read successfully." << endl;
63     else
64       std::cout << "error: only " << is1.gcount() << " could be read";
65     is1.close();
66         
67         
68     is2.seekg (0, is2.end);
69     int length2 = is2.tellg();
70     is2.seekg (0, is2.beg);
71     char * buffer2 = new char [length2];
72     
73     std::cout << "Reading " << length2 << " characters... ";
74     
75     is2.read (buffer2,length2);
76     if (is2)
77       std::cout << "all characters read successfully." << endl;
78     else
79       std::cout << "error: only " << is2.gcount() << " could be read";
80     is2.close();
81     
82         
83         
84     // ...buffer contains the entire file...
85     
86     char *curr1 = buffer1;
87     char *endBuff1 = buffer1+length1;
88         
89     char *curr2 = buffer2;
90     char *endBuff2 = buffer2+length2;
91         
92     int iEvent = 0;
93         static const int nExt = 4;
94     
95         while( curr1 < endBuff1  && curr2 < endBuff2 ){
96 //cout<<" curr1 endBuff1 curr2 endBuff2 "<< static_cast<void*> (curr1)<<" "<<static_cast<void*> (endBuff1)<<" "<<static_cast<void*> (curr2)<<" "<<static_cast<void*> (endBuff2)<<endl;
97
98         Int_t diff[6]={0};
99       AliFlatESDEvent *flatEsd1 = reinterpret_cast<AliFlatESDEvent *>(curr1);
100       AliFlatESDEvent *flatEsd2 = reinterpret_cast<AliFlatESDEvent *>(curr2);
101           /*
102           if(flatEsd1->GetNumberOfTracks()==0  ||  flatEsd2->GetNumberOfTracks() ==0){
103                 if(flatEsd1->GetNumberOfTracks()==0){
104                   curr1=curr1+ flatEsd1->GetSize();
105                 }
106                 if(flatEsd2->GetNumberOfTracks()==0){
107                   curr2=curr2+ flatEsd2->GetSize();
108                 }
109                 continue;
110           }
111           */
112           
113       cout<<endl<<"Reading event "<<iEvent<<":"<<endl;
114           /*
115         if(  flatEsd1->GetNumberOfTracks() != flatEsd2->GetNumberOfTracks() ) {
116                 cout<<"\t\tDIFFERENCE!: ";
117                 diff[0] =1;
118         }
119         cout<<"\t\tntracks: "<<flatEsd1->GetNumberOfTracks()<< " | " <<flatEsd2->GetNumberOfTracks();
120         hNTracks->Fill(flatEsd1->GetNumberOfTracks(),flatEsd2->GetNumberOfTracks());
121           
122           if(  flatEsd1->GetNumberOfV0s() != flatEsd2->GetNumberOfV0s() ){
123                 cout<<"\t\tDIFFERENCE!: ";
124                 diff[1] =1;
125         }
126           cout<<"\t\tnV0's: "<<flatEsd1->GetNumberOfV0s()<< " | " <<flatEsd2->GetNumberOfV0s()<<endl;
127           hNV0s->Fill(flatEsd1->GetNumberOfV0s(),flatEsd2->GetNumberOfV0s());
128          
129           if( (Bool_t) flatEsd1->GetPrimaryVertexTracks() != (Bool_t) flatEsd2->GetPrimaryVertexTracks()  ){
130                 cout<<"\t\tDIFFERENCE!: ";
131                 diff[2] =1;
132         }
133
134
135           cout<<"\t\tvtx tracks: "<<(Bool_t) flatEsd1->GetPrimaryVertexTracks()<< " | " <<      (Bool_t) flatEsd2->GetPrimaryVertexTracks();      
136           hVtxTr->Fill( (Bool_t) flatEsd1->GetPrimaryVertexTracks(), (Bool_t) flatEsd2->GetPrimaryVertexTracks());
137
138          
139           
140           if( (Bool_t) flatEsd1->GetPrimaryVertexSPD() != (Bool_t) flatEsd2->GetPrimaryVertexSPD()  ){
141                 cout<<"\t\tDIFFERENCE!: ";
142                 diff[3] =1;
143         }
144       cout<<"\t\tvtx SPD: "<<(Bool_t) flatEsd1->GetPrimaryVertexSPD() << " | " << (Bool_t) flatEsd2->GetPrimaryVertexSPD()<<endl;
145           hVtxSPD->Fill( (Bool_t) flatEsd1->GetPrimaryVertexSPD(), (Bool_t) flatEsd2->GetPrimaryVertexSPD());
146
147   if(true|| (Bool_t)flatEsd1->GetPrimaryVertexSPD() && (Bool_t)flatEsd2->GetPrimaryVertexSPD()  ){
148                 cout<<endl<<"\t\tvtx tracksX: "<< flatEsd1->GetPrimaryVertexSPD()<<" | " <<     flatEsd2->GetPrimaryVertexSPD();          
149                 }
150           */
151           
152           // compare tracks
153           #if 1
154           AliFlatESDTrack *track1 = flatEsd1->GetTracks();
155           AliFlatESDTrack *track2 = flatEsd2->GetTracks();
156     for (Int_t idxTrack = 0; idxTrack < flatEsd1->GetNumberOfTracks() && track1 && track2; ++idxTrack) { 
157
158                 AliFlatExternalTrackParam* ext[2][nExt] ={
159                         {
160                                 track1->GetTrackParamRefitted(),
161                                 track1->GetTrackParamIp(),
162                                 track1->GetTrackParamTPCInner(),
163                                 track1->GetTrackParamOp(),
164                 //              track1->GetTrackParamCp(),
165                 //              track1->GetTrackParamITSOut()
166                         },
167                         {
168                                 track2->GetTrackParamRefitted(),
169                                 track2->GetTrackParamIp(),
170                                 track2->GetTrackParamTPCInner(),
171                                 track2->GetTrackParamOp(),
172                         //      track2->GetTrackParamCp(),
173                         //      track2->GetTrackParamITSOut()
174                         }
175                 };
176         
177         //Printf("  TEST: FlatTrack1 %d > FlatExternalTrackParam1 > %p %p %p %p", idxTrack, exp11, exp21, exp31, exp41);
178         //Printf("  TEST: FlatTrack2 %d > FlatExternalTrackParam2 > %p %p %p %p", idxTrack, exp12, exp22, exp32, exp42);
179
180
181         for(int iExt=0; iExt<nExt; ++iExt){
182 cout<<endl<<iExt<<endl;         
183 if(!ext[0][iExt] && !ext[1][iExt]) continue;    
184                 if(!ext[0][iExt] && ext[1][iExt]){
185                 //      cout<<"DIFFERENCE!: ";
186                         cout<<" ext"<<iExt<<" not set in "<<filename1<<endl;
187                 }       
188                 if(ext[0][iExt] && !ext[1][iExt]){
189                 //      cout<<"DIFFERENCE!: ";
190                         cout<<" ext"<<iExt<<" not set in "<<filename2<<endl;
191                 }
192
193
194                 if( (!ext[0][iExt] || !ext[1][iExt])|| ext[0][iExt]->GetAlpha() != ext[1][iExt]->GetAlpha() ) {
195         //              cout<<"DIFFERENCE!: ";
196                         //cout<<" alpha"<<iExt<<" :"  << (ext[0][iExt] ? ext[0][iExt]->GetAlpha() : -99.)  << "\t\t" << (ext[1][iExt] ?  ext[1][iExt]->GetAlpha(): -99.)<<endl;
197                         diff[4]=1;
198                 }       cout<<" alpha"<<iExt<<" :"  << (ext[0][iExt] ? ext[0][iExt]->GetAlpha() : -99.)  << "\t\t" << (ext[1][iExt] ?  ext[1][iExt]->GetAlpha(): -99.)<<endl;
199                         
200
201                 if( (!ext[0][iExt] || !ext[1][iExt])||ext[0][iExt]->GetX() != ext[1][iExt]->GetX() ) {
202                         //cout<<"DIFFERENCE!: ";
203                         //cout<<" GetX"<<iExt<<" :"  << (ext[0][iExt] ? ext[0][iExt]->GetX(): -99.)  << " | " << (ext[1][iExt] ?  ext[1][iExt]->GetX(): -99.)<<endl;
204                         diff[4]=1;
205                 }       
206 cout<<" GetX"<<iExt<<" :"  << (ext[0][iExt] ? ext[0][iExt]->GetX(): -99.)  << " | " << (ext[1][iExt] ?  ext[1][iExt]->GetX(): -99.)<<endl;
207
208
209                 if( (!ext[0][iExt] || !ext[1][iExt])||ext[0][iExt]->GetSigned1Pt() !=  ext[0][iExt]->GetSigned1Pt() ) {
210                         //cout<<"DIFFERENCE!: ";
211                         //cout<<" 1/pt"<<iExt<<" :"  <<  (ext[0][iExt] ? ext[0][iExt]->GetSigned1Pt(): -99.)  << " | " << (ext[1][iExt] ?  ext[1][iExt]->GetSigned1Pt(): -99.)<<endl;
212                         diff[4]=1;
213                 }       
214         cout<<" 1/pt"<<iExt<<" :"  <<  (ext[0][iExt] ? ext[0][iExt]->GetSigned1Pt(): -99.)  << " | " << (ext[1][iExt] ?  ext[1][iExt]->GetSigned1Pt(): -99.)<<endl;
215                         
216
217 }
218         
219           
220           /*
221         if( track1->GetNumberOfTPCClusters() != track2->GetNumberOfTPCClusters() ){
222                 cout<<"DIFFERENCE!: ";
223                 cout<<" nTPCclusters: "<<track1->GetNumberOfTPCClusters()<< " | " <<track2->GetNumberOfTPCClusters()<< endl;
224                 diff[4]=1;
225         }  
226         if( track1->GetNumberOfITSClusters() != track2->GetNumberOfITSClusters() ){
227                 cout<<"DIFFERENCE!: ";
228                 cout<<" nITSclusters: "<<track1->GetNumberOfITSClusters()<< " | " <<track2->GetNumberOfITSClusters()<< endl;
229                 diff[4]=1;
230         }
231 */
232
233 // compare clusters
234         if( verbose &&  track1->GetNumberOfTPCClusters() == track2->GetNumberOfTPCClusters()){
235                 for (Int_t idxCluster = 0; idxCluster < track1->GetNumberOfTPCClusters(); ++idxCluster){
236                         AliFlatTPCCluster * cl1 = track1->GetTPCCluster(idxCluster);
237                         AliFlatTPCCluster * cl2 = track2->GetTPCCluster(idxCluster);
238 /*
239                         if( cl1->GetX()&& cl2->GetX() && cl1->GetX() != cl2->GetX() ){
240                                 cout<<"DIFFERENCE!: ";
241                                 cout<<" cluster: "<<idxCluster<<" GetX :"<<cl1->GetX()<< " | " <<cl2->GetX()<< endl;
242                                 diff=kTRUE;
243                         }
244                                 cout<<" cluster: "<<idxCluster<<" GetX :"<<cl1->GetX()<< " | " <<cl2->GetX()<< endl;
245                                 cout<<" cluster: "<<idxCluster<<" GetY :"<<cl1->GetY()<< " | " <<cl2->GetY()<< endl;
246
247                         if( cl1 && cl2 && cl1->GetY() != cl2->GetY() ){
248                                 cout<<"DIFFERENCE!: ";
249                                 cout<<" cluster: "<<idxCluster<<" GetY :"<<cl1->GetY()<< " | " <<cl2->GetY()<< endl;
250                                 diff=kTRUE;
251                         }
252                         if( cl1->GetZ()&& cl2->GetZ() && cl1->GetZ() != cl2->GetZ() ){
253                                 cout<<"DIFFERENCE!: ";
254                                 cout<<" cluster: "<<idxCluster<<" GetZ :"<<cl1->GetZ()<< " | " <<cl2->GetZ()<< endl;
255                                 diff=kTRUE;
256                         }
257 */
258                         if( cl1->GetPadRow()&& cl2->GetPadRow() && cl1->GetPadRow() != cl2->GetPadRow() ){
259                                 cout<<"DIFFERENCE!: ";
260                                 cout<<" cluster: "<<idxCluster<<" GetPadRow :"<<cl1->GetPadRow()<< " | " <<cl2->GetPadRow()<< endl;
261                                 diff[5]=1;
262                         }
263
264                 }
265           }
266      
267       track1 = track1->GetNextTrack();
268       track2 = track2->GetNextTrack();
269           
270           
271           }
272 #endif
273            hStat->Fill(0);        
274           Bool_t diffs=kFALSE;
275           for(int iDiff=0; iDiff<5;++iDiff){
276                 if(diff[iDiff]){
277                         hStat->Fill(iDiff+2);
278                         diffs = kTRUE;
279                 }
280         }
281         if(!diffs) hStat->Fill(1);        
282
283
284       curr1=curr1+ flatEsd1->GetSize();
285       curr2=curr2+ flatEsd2->GetSize();
286       iEvent++;
287     }
288
289     delete[] buffer1;
290     delete[] buffer2;
291   }
292   else {
293     cout << "File could not be read" << endl;
294   }
295
296
297
298   
299         TList histosList;
300         histosList.Add(hStat);
301         histosList.Add(hNTracks);
302         histosList.Add(hNV0s);
303         histosList.Add(hVtxTr);
304         histosList.Add(hVtxSPD);
305   histosList.SaveAs(outputFilename);
306   
307   return;
308 }