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           flatEsd1->Reinitialize();
 103           flatEsd2->Reinitialize();
 104
 105       cout<<endl<<"Reading event "<<iEvent<<":";
 106
 107         cout<<endl<<"ntracks:\t"<<flatEsd1->GetNumberOfTracks()<< " | " <<flatEsd2->GetNumberOfTracks();
 108         if(  flatEsd1->GetNumberOfTracks() != flatEsd2->GetNumberOfTracks() ) {
 109                 cout<<"\t\tDIFFERENCE!!!";
 110                 diff[0] =1;
 111         }
 112         //hNTracks->Fill(flatEsd1->GetNumberOfTracks(),flatEsd2->GetNumberOfTracks());
 113
 114           cout<<endl<<"nV0's:\t"<<flatEsd1->GetNumberOfV0s()<< " | " <<flatEsd2->GetNumberOfV0s();
 115           if(  flatEsd1->GetNumberOfV0s() != flatEsd2->GetNumberOfV0s() ){
 116                 cout<<"\t\tDIFFERENCE!!!";
 117                 diff[1] =1;
 118         }
 119           //hNV0s->Fill(flatEsd1->GetNumberOfV0s(),flatEsd2->GetNumberOfV0s());
 120
 121
 122
 123           /*
 124
 125           if( (Bool_t) flatEsd1->GetFlatPrimaryVertexTracks() != (Bool_t) flatEsd2->GetFlatPrimaryVertexTracks()  ){
 126                 cout<<"\t\tDIFFERENCE!: "<<endl;
 127                 diff[2] =1;
 128         }
 129
 130
 131           cout<<"vtx tracks:\t"<<(Bool_t) flatEsd1->GetFlatPrimaryVertexTracks()<< " | " <<     (Bool_t) flatEsd2->GetFlatPrimaryVertexTracks()<<endl;
 132           //hVtxTr->Fill( (Bool_t) flatEsd1->GetFlatPrimaryVertexTracks(), (Bool_t) flatEsd2->GetFlatPrimaryVertexTracks());
 133
 134
 135
 136           if( (Bool_t) flatEsd1->GetFlatPrimaryVertexSPD() != (Bool_t) flatEsd2->GetFlatPrimaryVertexSPD()  ){
 137                 cout<<"\t\tDIFFERENCE!: "<<endl;
 138                 diff[3] =1;
 139         }
 140       cout<<"vtx SPD:\t"<<(Bool_t) flatEsd1->GetFlatPrimaryVertexSPD() << " | " << (Bool_t) flatEsd2->GetFlatPrimaryVertexSPD()<<endl;
 141           //hVtxSPD->Fill( (Bool_t) flatEsd1->GetFlatPrimaryVertexSPD(), (Bool_t) flatEsd2->GetFlatPrimaryVertexSPD());
 142
 143
 144
 145   if((Bool_t)flatEsd1->GetFlatPrimaryVertexTracks() && (Bool_t)flatEsd2->GetFlatPrimaryVertexTracks()  ){
 146                 cout<<endl<<"vtx tracks -> X,Y,Z:\t"
 147                         << flatEsd1->GetFlatPrimaryVertexTracks()->GetX()
 148                         <<","<< flatEsd1->GetFlatPrimaryVertexTracks()->GetY()
 149                         <<","<< flatEsd1->GetFlatPrimaryVertexTracks()->GetZ()
 150                         <<" | " <<flatEsd2->GetFlatPrimaryVertexTracks()->GetX()
 151                         <<","<< flatEsd2->GetFlatPrimaryVertexTracks()->GetY()
 152                         <<","<< flatEsd2->GetFlatPrimaryVertexTracks()->GetZ()<<endl;
 153                 }
 154           */
 155
 156           // compare tracks
 157 if(verbose){
 158           AliFlatESDTrack *track1 = const_cast<AliFlatESDTrack*> (flatEsd1->GetTracks());
 159           AliFlatESDTrack *track2 = const_cast<AliFlatESDTrack*> (flatEsd2->GetTracks());
 160     for (Int_t idxTrack = 0; idxTrack < flatEsd1->GetNumberOfTracks() && track1 && track2; ++idxTrack) {
 161
 162                 //track2->Reinitialize();
 163                 const AliFlatExternalTrackParam* ext[2][nExt] ={
 164                         {
 165                                 track1->GetFlatTrackParamRefitted(),
 166                                 track1->GetFlatTrackParamIp(),
 167                                 track1->GetFlatTrackParamTPCInner(),
 168                                 track1->GetFlatTrackParamOp(),
 169                 //              track1->GetFlatTrackParamCp(),
 170                 //              track1->GetFlatTrackParamITSOut()
 171                         },
 172                         {
 173                                 track2->GetFlatTrackParamRefitted(),
 174                                 track2->GetFlatTrackParamIp(),
 175                                 track2->GetFlatTrackParamTPCInner(),
 176                                 track2->GetFlatTrackParamOp(),
 177                         //      track2->GetFlatTrackParamCp(),
 178                         //      track2->GetFlatTrackParamITSOut()
 179                         }
 180                 };
 181
 182         //Printf("  TEST: FlatTrack1 %d > FlatExternalTrackParam1 > %p %p %p %p", idxTrack, exp11, exp21, exp31, exp41);
 183         //Printf("  TEST: FlatTrack2 %d > FlatExternalTrackParam2 > %p %p %p %p", idxTrack, exp12, exp22, exp32, exp42);
 184
 185
 186         for(int iExt=0; iExt<nExt; ++iExt){
 187 if(!ext[0][iExt] && !ext[1][iExt]) continue;
 188                 if(!ext[0][iExt] && ext[1][iExt]){
 189                 //      cout<<"DIFFERENCE!: ";
 190                         cout<<" ext"<<iExt<<" not set in "<<filename1<<endl;
 191                 }
 192                 if(ext[0][iExt] && !ext[1][iExt]){
 193                 //      cout<<"DIFFERENCE!: ";
 194                         cout<<" ext"<<iExt<<" not set in "<<filename2<<endl;
 195                 }
 196
 197
 198                 if( (!ext[0][iExt] || !ext[1][iExt])|| ext[0][iExt]->GetAlpha() != ext[1][iExt]->GetAlpha() ) {
 199                         cout<<"\t\tDIFFERENCE!: "<<endl;
 200                         //cout<<" alpha"<<iExt<<" :"  << (ext[0][iExt] ? ext[0][iExt]->GetAlpha() : -99.)  << "\t\t" << (ext[1][iExt] ?  ext[1][iExt]->GetAlpha(): -99.)<<endl;
 201                         diff[4]=1;
 202                 }       cout<<" alpha"<<iExt<<" :\t"  << (ext[0][iExt] ? ext[0][iExt]->GetAlpha() : -99.)  << " | " << (ext[1][iExt] ?  ext[1][iExt]->GetAlpha(): -99.)<<endl;
 203
 204
 205                 if( (!ext[0][iExt] || !ext[1][iExt])||ext[0][iExt]->GetX() != ext[1][iExt]->GetX() ) {
 206                         cout<<"\t\tDIFFERENCE!: "<<endl;
 207                         //cout<<" GetX"<<iExt<<" :"  << (ext[0][iExt] ? ext[0][iExt]->GetX(): -99.)  << " | " << (ext[1][iExt] ?  ext[1][iExt]->GetX(): -99.)<<endl;
 208                         diff[4]=1;
 209                 }
 210 cout<<" GetX"<<iExt<<" :\t"  << (ext[0][iExt] ? ext[0][iExt]->GetX(): -99.)  << " | " << (ext[1][iExt] ?  ext[1][iExt]->GetX(): -99.)<<endl;
 211
 212
 213                 if( (!ext[0][iExt] || !ext[1][iExt])||ext[0][iExt]->GetSigned1Pt() !=  ext[0][iExt]->GetSigned1Pt() ) {
 214                         cout<<"\t\tDIFFERENCE!: "<<endl;
 215                         //cout<<" 1/pt"<<iExt<<" :"  <<  (ext[0][iExt] ? ext[0][iExt]->GetSigned1Pt(): -99.)  << " | " << (ext[1][iExt] ?  ext[1][iExt]->GetSigned1Pt(): -99.)<<endl;
 216                         diff[4]=1;
 217                 }
 218         cout<<" 1/pt"<<iExt<<" :\t"  <<  (ext[0][iExt] ? ext[0][iExt]->GetSigned1Pt(): -99.)  << " | " << (ext[1][iExt] ?  ext[1][iExt]->GetSigned1Pt(): -99.)<<endl;
 219
 220
 221 }
 222
 223
 224           /*
 225         if( track1->GetNumberOfTPCClusters() != track2->GetNumberOfTPCClusters() ){
 226                 cout<<"DIFFERENCE!: ";
 227                 cout<<" nTPCclusters: "<<track1->GetNumberOfTPCClusters()<< " | " <<track2->GetNumberOfTPCClusters()<< endl;
 228                 diff[4]=1;
 229         }
 230         if( track1->GetNumberOfITSClusters() != track2->GetNumberOfITSClusters() ){
 231                 cout<<"DIFFERENCE!: ";
 232                 cout<<" nITSclusters: "<<track1->GetNumberOfITSClusters()<< " | " <<track2->GetNumberOfITSClusters()<< endl;
 233                 diff[4]=1;
 234         }
 235 */
 236
 237 #if 0
 238
 239 // compare clusters
 240         if( verbose &&  track1->GetNumberOfTPCClusters() == track2->GetNumberOfTPCClusters()){
 241                 for (Int_t idxCluster = 0; idxCluster < track1->GetNumberOfTPCClusters(); ++idxCluster){
 242                         AliFlatTPCCluster * cl1 = track1->GetTPCCluster(idxCluster);
 243                         AliFlatTPCCluster * cl2 = track2->GetTPCCluster(idxCluster);
 244 /*
 245                         if( cl1->GetX()&& cl2->GetX() && cl1->GetX() != cl2->GetX() ){
 246                                 cout<<"DIFFERENCE!: ";
 247                                 cout<<" cluster: "<<idxCluster<<" GetX :"<<cl1->GetX()<< " | " <<cl2->GetX()<< endl;
 248                                 diff=kTRUE;
 249                         }
 250                                 cout<<" cluster: "<<idxCluster<<" GetX :"<<cl1->GetX()<< " | " <<cl2->GetX()<< endl;
 251                                 cout<<" cluster: "<<idxCluster<<" GetY :"<<cl1->GetY()<< " | " <<cl2->GetY()<< endl;
 252
 253                         if( cl1 && cl2 && cl1->GetY() != cl2->GetY() ){
 254                                 cout<<"DIFFERENCE!: ";
 255                                 cout<<" cluster: "<<idxCluster<<" GetY :"<<cl1->GetY()<< " | " <<cl2->GetY()<< endl;
 256                                 diff=kTRUE;
 257                         }
 258                         if( cl1->GetZ()&& cl2->GetZ() && cl1->GetZ() != cl2->GetZ() ){
 259                                 cout<<"DIFFERENCE!: ";
 260                                 cout<<" cluster: "<<idxCluster<<" GetZ :"<<cl1->GetZ()<< " | " <<cl2->GetZ()<< endl;
 261                                 diff=kTRUE;
 262                         }
 263 */
 264                         if( cl1->GetPadRow()&& cl2->GetPadRow() && cl1->GetPadRow() != cl2->GetPadRow() ){
 265                                 cout<<"DIFFERENCE!: ";
 266                                 cout<<" cluster: "<<idxCluster<<" GetPadRow :"<<cl1->GetPadRow()<< " | " <<cl2->GetPadRow()<< endl;
 267                                 diff[5]=1;
 268                         }
 269
 270                 }
 271           }
 272
 273 #endif
 274       track1 = const_cast<AliFlatESDTrack*> (track1->GetNextTrack());
 275       track2 = const_cast<AliFlatESDTrack*> (track2->GetNextTrack());
 276
 277
 278           }
 279 }
 280
 281           /*
 282         //   hStat->Fill(0);
 283           Bool_t diffs=kFALSE;
 284           for(int iDiff=0; iDiff<5;++iDiff){
 285                 if(diff[iDiff]){
 286         //              hStat->Fill(iDiff+2);
 287                         diffs = kTRUE;
 288                 }
 289         }
 290         if(!diffs) hStat->Fill(1);
 291 */
 292
 293       curr1=curr1+ flatEsd1->GetSize();
 294       curr2=curr2+ flatEsd2->GetSize();
 295       iEvent++;
 296     }
 297
 298     delete[] buffer1;
 299     delete[] buffer2;
 300   }
 301   else {
 302     cout << "File could not be read" << endl;
 303   }
 304
 305 /*
 306
 307
 308         TList histosList;
 309         histosList.Add(hStat);
 310         histosList.Add(hNTracks);
 311         histosList.Add(hNV0s);
 312         histosList.Add(hVtxTr);
 313         histosList.Add(hVtxSPD);
 314   histosList.SaveAs(outputFilename);
 315   */
 316   return;
 317 }