T0/T0da.cxx

   1 extern "C" {
   2 #include <daqDA.h>
   3 }
   4
   5 #include "event.h"
   6 #include "monitor.h"
   7
   8 #include <Riostream.h>
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11
  12 //AliRoot
  13 #include "AliRawReaderDate.h"
  14 #include "AliRawReader.h"
  15 #include "AliT0digit.h"
  16 #include "AliT0RawReader.h"
  17 #include "AliT0Dqclass.h"
  18
  19 //ROOT
  20 #include "TFile.h"
  21 #include "TKey.h"
  22 #include "TH2S.h"
  23 #include "TObject.h"
  24 #include "TBenchmark.h"
  25 #include "TRandom.h"
  26
  27 #include "TCanvas.h"
  28 #include "TString.h"
  29 #include "TH1.h"
  30 #include "TF1.h"
  31 #include "TSpectrum.h"
  32 #include "TVirtualFitter.h"
  33
  34 void fitv2DA(char *filename = "t0histdate.root");
  35
  36 /* Main routine
  37       Arguments:
  38       1- monitoring data source
  39 */
  40 int main(int argc, char **argv) {
  41   int status;
  42
  43   if (argc!=2) {
  44     printf("Wrong number of arguments\n");
  45     return -1;
  46   }
  47
  48
  49   /* define data source : this is argument 1 */
  50   status=monitorSetDataSource( argv[1] );
  51   if (status!=0) {
  52     printf("monitorSetDataSource() failed : %s\n",monitorDecodeError(status));
  53     return -1;
  54   }
  55
  56
  57   /* declare monitoring program */
  58   status=monitorDeclareMp( __FILE__ );
  59   if (status!=0) {
  60     printf("monitorDeclareMp() failed : %s\n",monitorDecodeError(status));
  61     return -1;
  62   }
  63
  64
  65   /* define wait event timeout - 1s max */
  66   monitorSetNowait();
  67   monitorSetNoWaitNetworkTimeout(1000);
  68
  69
  70   /* log start of process */
  71   printf("T0 monitoring program started\n");
  72
  73
  74   // Allocation of histograms - start
  75   TH1F * hCFD[24]; TH1F *hLED[24]; TH1F*hQT01[4];TH1F*hQT02[4];
  76   TH1F*hQTD[4]; TH1F*hADC[24]; TH2*hQTCFD[4];
  77   Char_t buf1[10], buf2[10], buf3[10], buf4[10], buf5[10], buf6[10],buf7[10];
  78   for (Int_t ic=0; ic<24; ic++)
  79     {
  80       sprintf(buf1,"CFD%i",ic+1);
  81       hCFD[ic]= new TH1F(buf1,"CFD",6000,0,6000);
  82       sprintf(buf2,"LED%i",ic+1);
  83       hLED[ic]= new TH1F(buf2,"LED",2000,0,2000);
  84       //LED-CFD
  85       sprintf(buf6,"ADC%i",ic+1);
  86       hADC[ic]= new TH1F(buf6,"ADC",6000,0,6000);
  87     }
  88
  89   for (Int_t iq=0; iq<4; iq++)
  90     {
  91       //QT01 - QT04
  92       sprintf(buf3,"QT0%i",iq+1);
  93       hQT01[iq]= new TH1F(buf3,"QT01",6000,0,6000);
  94       sprintf(buf4,"QT1_%i",iq+1);
  95       //QT11 - QT14
  96       hQT02[iq]= new TH1F(buf4,"QT02",6000,0,6000);
  97       sprintf(buf5,"QTD_%i",iq+1);
  98       //QT11-QT01 ....
  99       hQTD[iq]= new TH1F(buf5,"QTdiff",4500,1500,6000);
 100
 101       sprintf(buf7,"QTCFD_%i",iq+1);
 102       hQTCFD[iq]= new TH2F(buf7,"QT vs CFD",500,0,6000,500,0,5000);
 103   }
 104
 105   TH1F *hORA= new TH1F("hORA"," T0 A ",1000, 1000,2000);
 106   TH1F *hORC= new TH1F("hORC"," T0 C ",1000, 1000,2000);
 107   TH1F*hEffCFD= new TH1F("hEffCFD","Effeciency",8,0.25,4.25);
 108
 109   //  TH2F*hQTCFD= new TH2F("hQTCFD","QT vs CFD",500,0.5,6000.5,500,0.5,5000.5);
 110   TH2F*hQTLED= new TH2F("hQTLED","QT vs LED",500,0.5,6000.5,500,0.5,5000.5);
 111   TH2F*hLEDCFD= new TH2F("hLEDCFD","LEd vs CFD",500,-0.5,10000.5,500,-0.5,10000.5);
 112   // Allocation of histograms - end
 113
 114   Int_t iev=0;
 115   /* main loop (infinite) */
 116   for(;;) {
 117     struct eventHeaderStruct *event;
 118     eventTypeType eventT;
 119
 120     /* check shutdown condition */
 121     if (daqDA_checkShutdown()) {break;}
 122
 123     /* get next event (blocking call until timeout) */
 124     status=monitorGetEventDynamic((void **)&event);
 125     if (status==(int)MON_ERR_EOF) {
 126       printf ("End of File detected\n");
 127       break; /* end of monitoring file has been reached */
 128     }
 129
 130     if (status!=0) {
 131       printf("monitorGetEventDynamic() failed : %s\n",monitorDecodeError(status));
 132       break;
 133     }
 134
 135     /* retry if got no event */
 136     if (event==NULL) {
 137       continue;
 138     }
 139
 140     iev++;
 141
 142     /* use event - here, just write event id to result file */
 143     eventT=event->eventType;
 144
 145     if (eventT==PHYSICS_EVENT) {
 146       printf(" event number = %i \n",iev);
 147
 148       // Initalize raw-data reading and decoding
 149       AliRawReader *reader = new AliRawReaderDate((void*)event);
 150       // Enable the following two lines in case of real-data
 151       //    reader->LoadEquipmentIdsMap("T0map.txt");
 152       //    reader->RequireHeader(kFALSE);
 153       AliT0RawReader *start = new AliT0RawReader(reader);
 154
 155       // Read raw data
 156       Int_t allData[110][5];
 157       start->Next();
 158       for (Int_t i=0; i<110; i++) {
 159         allData[i][0]= start->GetData(i,0);
 160         if (allData[i][0] != 0) cout<<"event "<<event<<" i "<<iev<<" "<<allData[i][0] - allData[0][0]<<endl;
 161       }
 162
 163       // Fill the histograms
 164       for(Int_t ik=0; ik<24; ik++) {
 165         for (Int_t iHit=0; iHit<5; iHit++) {
 166           if(allData[ik+1][iHit] != 0 ) {
 167             //     cout<<"event "<<event<<" iHit "<< iHit<< allData[ik+1][0]<<" "<< allData[0][0]<< endl;
 168             hCFD[ik]->Fill((allData[ik+1][iHit] - allData[0][0])/1000.);
 169           }
 170           if(allData[ik+25][iHit] != 0 )    hLED[ik]->Fill(allData[ik+25][iHit] - allData[0][0]);
 171           if(allData[ik+25][iHit] != 0 || allData[ik+1][iHit] !=0)
 172             hADC[ik]->Fill(allData[ik+25][iHit] - allData[ik+1][iHit]);
 173         }
 174       }
 175       hLEDCFD->Fill(allData[9][0],allData[1][0]);
 176
 177       hQT01[1]->Fill(allData[18][0] - allData[0][0]);
 178       hQT01[2]->Fill(allData[19][0] - allData[0][0]);
 179       hQT01[3]->Fill(allData[20][0] - allData[0][0]);
 180       hQT02[0]->Fill(allData[21][0] - allData[0][0]);
 181       hQT02[0]->Fill(allData[22][0] - allData[0][0]);
 182       hQT02[0]->Fill(allData[23][0] - allData[0][0]);
 183       hQT02[0]->Fill(allData[24][0] - allData[0][0]);
 184
 185
 186       hQTD[0]->Fill(allData[21][0] - allData[17][0]);
 187       hQTLED->Fill(allData[9][0] - allData[0][0], allData[21][0]-allData[17][0]);
 188       hQTCFD[0]->Fill(allData[21][0]-allData[17][0], allData[1][0] - allData[0][0]);
 189
 190       hQTD[1]->Fill(allData[22][0]-allData[18][0]);
 191       hQTCFD[1]->Fill(allData[22][0]-allData[18][0], allData[2][0] - allData[0][0]);
 192
 193       hQTD[2]->Fill(allData[23][0]-allData[19][0]);
 194       hQTCFD[2]->Fill(allData[23][0]-allData[19][0], allData[3][0] - allData[0][0]);
 195
 196       hQTD[3]->Fill(allData[24][0]-allData[20][0]);
 197       hQTCFD[3]->Fill(allData[24][0]-allData[20][0], allData[4][0] - allData[0][0]);
 198
 199
 200       hORA->Fill(allData[25][0] - allData[0][0]);
 201       hORA->Fill(allData[26][0] - allData[0][0]);
 202       // End of fill histograms
 203
 204     }
 205
 206
 207
 208     /* free resources */
 209     free(event);
 210
 211     /* exit when last event received, no need to wait for TERM signal */
 212     if (eventT==END_OF_RUN) {
 213       printf("EOR event detected\n");
 214       break;
 215     }
 216   }
 217
 218   // write a file with the histograms
 219   Char_t filehist[20];
 220   sprintf(filehist,"t0histdate.root");
 221   TFile *hist = new TFile(filehist,"RECREATE");
 222
 223   hLEDCFD->Write();
 224
 225   for (Int_t i=0; i<24; i++)
 226     {
 227       hCFD[i]->Write();
 228       hLED[i]->Write();
 229       hADC[i]->Write();
 230     }
 231   for (Int_t i=0; i<4; i++)
 232     {
 233       hQT01[i]->Write();
 234       hQT02[i]->Write();
 235       hQTD[i]->Write();
 236       hQTCFD[i]->Write();
 237     }
 238
 239   hEffCFD->Write();
 240   hORA->Write();
 241   hORC->Write();
 242   hQTLED->Write();
 243   hist->Close();
 244
 245   // Fit the histograms and write the output file
 246   fitv2DA();
 247
 248   return status;
 249 }
 250
 251 void fitv2DA(char *filename)
 252 {
 253         TFile *file = TFile::Open(filename);
 254         FILE *dafile;
 255         dafile = fopen("da.txt","w");
 256         fprintf(dafile,"Const\t\tMean\t\tSigma\t\tLBorder\tRBorder\n");
 257
 258 //      TFile *fithist = new TFile("fitedhist.root","RECREATE");
 259
 260         char histoname[10];
 261         char canvname[10];
 262         Float_t p[3] = {0.,0.,0.};
 263         Float_t LBordX = 0.;
 264         Float_t RBordX = 0.;
 265 /*
 266         TIter next(file.GetListOfKeys());
 267         TKey *key;
 268         while ((key = (TKey*)nextkey())) {
 269               if(key->GetName()==)
 270         }
 271 */
 272         Int_t npeaks = 10;
 273         cout<<"npeaks = "<<npeaks<<endl;
 274         TH1F *hist[24];
 275         TCanvas *c[24];
 276         for(Int_t i=0;i<24;i++){///////////#######Cos tu sie dzieje....
 277                 sprintf(histoname,"CFD%d",i+1);
 278                 sprintf(canvname,"C%d",i+1);
 279                 TH1F *histtemp = (TH1F*)file->Get(histoname);
 280                 hist[i] = histtemp;
 281         }
 282         cout<<"Wczytany plik z wykresami..."<<endl;
 283 //      TCanvas *c[8];
 284         for(Int_t ii=0;ii<24;ii++){
 285                 c[ii] = new TCanvas(canvname,canvname,800,600);
 286         p[0] = 0.;
 287         p[1] = 0.;
 288         p[2] = 0.;
 289         LBordX = 0.;
 290         RBordX = 0.;
 291
 292         Double_t max=0.0;
 293 //      Double_t peak,peakmax=0.0;
 294         Double_t tabmax[2] = {0.0, 0.0};
 295
 296         TSpectrum *s = new TSpectrum(2*npeaks);
 297         Int_t nfound = s->Search(hist[ii],2,"goff",0.05);
 298         cout<<"Found "<<nfound<<" peaks\n";
 299         if(nfound!=0){
 300         Float_t *xpeak = s->GetPositionX();
 301         for(Int_t k=0;k<nfound;k++)
 302         {
 303                 cout<<"Petla po pikach..."<<endl;
 304                 Float_t xp = xpeak[k];
 305                 Int_t xbin = hist[ii]->GetXaxis()->FindBin(xp);
 306                 Float_t yp = hist[ii]->GetBinContent(xbin);
 307                 if(yp>max) {
 308                         max = yp;
 309                         tabmax[0] = xp;
 310                         tabmax[1] = yp;
 311                 }
 312                 cout<<"xbin = "<<xbin<<"\txpeak = "<<xpeak[k]<<"\typeak = "<<yp<<endl;
 313         }
 314
 315         Double_t MinBordX, MaxBordX, MaxBordY;
 316
 317         MinBordX = tabmax[0]-tabmax[0]*0.05;
 318         MaxBordX = tabmax[0]+tabmax[0]*0.05;
 319         MaxBordY = tabmax[1]+tabmax[1]*0.5;
 320
 321         Double_t border;
 322         border = tabmax[1]*0.1;
 323
 324         Double_t LBordY = 99999.;
 325         Double_t RBordY = 99999.;
 326         LBordX = tabmax[0];
 327         RBordX = tabmax[0];
 328         do{
 329                 LBordY = hist[ii]->GetBinContent(LBordX);
 330                 LBordX -=1;
 331         }while(LBordY>border);
 332         do{
 333                 RBordX +=1;
 334                 RBordY = hist[ii]->GetBinContent(RBordX);
 335         }while(RBordY>border);
 336
 337
 338         hist[ii]->GetXaxis()->SetRangeUser(MinBordX,MaxBordX);
 339         hist[ii]->GetYaxis()->SetRangeUser(0,MaxBordY);
 340         TF1 *fit = new TF1("fit","gaus",LBordX,RBordX);
 341         hist[ii]->Fit("fit","IR");
 342         Double_t maxx=0.;
 343 //      Double_t chi2 = fit->GetChisquare();
 344 //      Float_t p[3];
 345
 346         for(Int_t j=0;j<3;j++){
 347                 p[j] = fit->GetParameter(j);
 348         }
 349
 350         max = fit->GetMaximum();
 351         maxx = fit->GetMaximumX();
 352 //      hist[ii]->Write();
 353         cout<<"\nFitMaxY = "<<max<<"\tFitMaxX "<<maxx<<"\tPeakMax*10% "<<border<<endl;
 354         cout<<"LBordY "<<LBordY<<"\tLBordX "<<LBordX+1<<"\tRBordY "<<RBordY<<"\tRBordX "<<RBordX<<endl;
 355         cout<<"Parameters: Const - MaxY - "<<p[0]<<"\tMean - "<<p[1]<<"\tSigma - "<<p[2]<<endl;
 356         cout<<"Borders: Min - "<<MinBordX<<"\tMax - "<<MaxBordX<<endl;
 357         if(p[0]==0.&&p[1]==0.&&p[2]==0.)        RBordX = 0;
 358         }
 359 //      else{
 360 //      return;
 361 //      }
 362         AliT0Dqclass *daqpar = new AliT0Dqclass();
 363         cout<<"Alit0dqclass..."<<endl;
 364         daqpar->SetTime(0,p[0]);
 365         cout<<"First parameter..."<<endl;
 366         daqpar->SetTime(1,p[1]);
 367         daqpar->SetTime(2,p[2]);
 368         cout<<"gauss parameters finished..."<<endl;
 369         daqpar->SetTime(3,LBordX+1);
 370         cout<<"first border..."<<endl;
 371         daqpar->SetTime(4,RBordX);
 372         cout<<"second border..."<<endl;
 373         TFile *DA = new TFile("DAfile.root","RECREATE");
 374         DA->cd();
 375         cout<<"file"<<endl;
 376         daqpar->Write("Time");
 377         cout<<"file written"<<endl;
 378         DA->Close();
 379         cout<<"file closed"<<endl;
 380         delete DA;
 381         cout<<"file deleted"<<endl;
 382
 383         fprintf(dafile,"%lf\t%lf\t%lf\t%d\t%d\n",p[0],p[1],p[2],LBordX+1,RBordX);
 384         }
 385
 386
 387         fclose(dafile);
 388 }