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 }