T0/readESDlocal.C

   1 #include "TH1F.h"
   2 #include "TProfile.h"
   3 #include "TH2F.h"
   4 #include "TMap.h"
   5 #include "TFile.h"
   6 #include "TTree.h"
   7 #include "TBranch.h"
   8
   9 #include "TGLabel.h"
  10 #include "TGrid.h"
  11 #include "TFitResult.h"
  12
  13 #include "TControlBar.h"
  14 #include <Riostream.h>
  15 #include "AliCDBManager.h"
  16 #include "AliRawReader.h"
  17 #include "AliRawReaderRoot.h"
  18 #include "AliT0LookUpValue.h"
  19 #include "AliT0LookUpKey.h"
  20 #include "AliT0Parameters.h"
  21 #include "AliT0RawReader.h"
  22 #include "AliT0CalibLaserData.h"
  23 #include "AliESDTZERO.h"
  24 #include "AliESDVertex.h"
  25 #include "AliESDVZERO.h"
  26 #include "AliESDEvent.h"
  27
  28 void readESDlocal(Int_t run)
  29 {
  30   //read T0 ESD from reconstructed T0 data
  31   TH1F *hMean= new TH1F("hMean"," T0 ", 1000,2000,3000);
  32   // TH1F *hMeanTOF= new TH1F("hMeanTOF"," T0 ", 1000,-10,10);
  33   // TH1F *hOrA= new TH1F("hOrA"," OrA T0 ",1000,-10,10);
  34   //  TH1F *hOrC= new TH1F("hOrC"," OrC T0 ",1000,-10,10);
  35   TH1F *hVertex = new TH1F("hVertex","Z position of vertex", 180,-100,100);
  36   TH1F *hVertexcalc = new TH1F("hVertexcalc","Z position of vertex", 180,-100,100);
  37   TH1F *hVertexT0 = new TH1F("hVertexT0","Z position of vertex", 1000,-10,10);
  38  TH1F *hVertexT0calc = new TH1F("hVertexT0calc","Z position of vertex", 1000,-10,10);
  39   TH1F *hVertexComp = new TH1F("hVertexComparison ","Z position of vertex", 180,-30,30);
  40   TH1F *hVertexT0only = new TH1F("hVertexT0only","Z position of vertex without SPD vertex", 180,-30,30);
  41   TH1F *hVertexSPD = new TH1F("hVertexSPD","Z position of vertex SPD", 180,-30,30);
  42
  43   TH2F *hVertexSPDT0 = new TH2F("hVertexSPDT0","Z position of vertex SPD-T0",   180,-30,30, 180,-30,30);
  44   //
  45    TH1F *hMeanTOF= new TH1F("hMeanTOF"," T0 ", 8000, -100,100);
  46     TH1F *hOrA= new TH1F("hOrA"," OrA T0 ",8000, -100,100);
  47     TH1F *hOrC= new TH1F("hOrC"," OrC T0 ",8000, -100,100);
  48     // TH1F *hMeanTOF= new TH1F("hMeanTOF"," T0 ", 2000, -8950,-8900);
  49     // TH1F *hOrA= new TH1F("hOrA"," OrA T0 ",2000, -8950,-8900);
  50     //   TH1F *hOrC= new TH1F("hOrC"," OrC T0 ",2000, -8950,-8900);
  51
  52
  53     TH1F *hMeanTOFcalc= new TH1F("hMeanTOFcalc"," T0 ", 2000, -8950,-8900);
  54     TH1F *hOrAcalc= new TH1F("hOrAcalc"," OrA T0 ",2000, -8950,-8900);
  55     TH1F *hOrCcalc= new TH1F("hOrCcalc"," OrC T0 ",2000, -8950,-8900);
  56     //  TH1F *hMeanTOFcalc= new TH1F("hMeanTOFcalc"," T0 ", 8000, -100,100);
  57     //   TH1F *hOrAcalc= new TH1F("hOrAcalc"," OrA T0 ",8000, -100,100);
  58     //   TH1F *hOrCcalc= new TH1F("hOrCcalc"," OrC T0 ",8000, -100,100);
  59
  60   TH1F *hAmp[24];  TH1F * hTime[24]; TH1F * hTimeDiff[24]; TH1F * hTimecorr[24];
  61   //  TProfile * hAmpTime[24];
  62   TH2F * hAmpTime[24];
  63
  64   for(Int_t ic=0; ic<24; ic++)
  65     {
  66       hAmp[ic] = new TH1F(Form("hAmp%i",ic+1),"Amp",100, -10, 20 );
  67       //      hAmp[ic] = new TH1F(Form("hAmp%i",ic),"Amp",400, 300, 700 );
  68       hTime[ic] = new TH1F(Form("hTime%i",ic+1)," time",  5000,0,10000);
  69       hTimecorr[ic] = new TH1F(Form("hTimecorr%i",ic+1)," time",  5000,0,10000);
  70       hTimeDiff[ic] = new TH1F(Form("hTimeDiff%i",ic+1)," time", 300,-300,300);
  71       //   hAmpTime[ic] = new TProfile(Form("hAmpTime%i",ic)," time",200, 300,700, 49000, 53000);
  72       hAmpTime[ic] = new TH2F(Form("hAmpTime%i",ic+1)," time",100, 0,10, 100,3000,4000);
  73         // hAmpTime[ic] = new TH2F(Form("hAmpTime%i",ic)," time",250, 350,600, 400, 49400, 49800);
  74     }
  75
  76   Float_t channelWidth=24.4;
  77   Float_t c = 0.0299792458; // cm/ps
  78   Float_t shift=0;
  79   Float_t fLatencyHPTDC = 9000;
  80
  81
  82  //run 124702
  83   // Float_t fLatencyL1 = 8.91358e+03;
  84   // Float_t fLatencyL1A = 8.91352e+03;
  85   // Float_t fLatencyL1C =8.91361e+03;
  86
  87   //run 125097
  88   //   Float_t fLatencyL1 =   8.914520e+03 ;
  89   // Float_t fLatencyL1A =  8.914860e+03 ;
  90   // Float_t fLatencyL1C =  8.914180e+03;
  91
  92 //run 125295
  93   Float_t fLatencyL1 =  8.91406e+03  ;
  94   Float_t fLatencyL1A = 8.91401e+03;
  95   Float_t fLatencyL1C = 8.91412e+03 ;
  96
  97
  98   //  Float_t fLatencyL1 = 0;
  99   //  Float_t fLatencyL1A = 0;
 100   // Float_t fLatencyL1C =0;
 101      //125842
 102   //  Float_t fLatencyL1 = 8.91306e+03;
 103   //  Float_t fLatencyL1A =8.91338e+03;
 104   // Float_t fLatencyL1C =8.91274e+03;
 105  //
 106      //126407
 107      //  Float_t fLatencyL1 = 8.91345e+03;
 108      //  Float_t fLatencyL1A =8.91378e+03;
 109      // Float_t fLatencyL1C =8.91311e+03;
 110  //
 111   const Double32_t *amp, *time, *mean;
 112   Double32_t time1[24];
 113   TString filenam[100];
 114
 115
 116
 117   //   Float_t equal[24]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
 118   //124358
 119   //Float_t equal[24]={22, 0 ,0 ,1 ,3 ,1 ,0 ,0 ,-4 ,1 ,1 ,-1 ,0 ,-3 ,0 ,-6 ,0 ,-3 ,2 ,0 ,-1 ,-1 ,-3 ,0 };
 120
 121      //124702
 122   // Float_t equal[24]={20, 0 ,15, 50, 37, 36, 23, -20, 77, 54, 44, -12,
 123   //                 0, 19, -25, -13, 35, 1, 23, 46, 20, 77, 50, 122};
 124
 125
 126   //125097
 127    Float_t equal[24]={8,0 ,15 ,50 ,36 ,34 ,24 ,-21 ,81 ,53 ,42 ,-12 ,0 ,19 ,-20 ,-7 ,37 ,5 ,25 ,49 ,19 ,78 ,52 ,124};
 128
 129   //  125295
 130
 131   //Float_t equal[24]={16, 0,  2, 2, 30, 3, -5, -8, 5, -4, 16, -5,
 132    //            0, 10, 10, 12, 13, -24, 15, 26, -2, 17, 10, -30};
 133
 134
 135
 136    //125842
 137   // Float_t equal[24]={16, 0, 1, 3,  32,   4, -2, -8, 5, -3, 15, -5,
 138   //                  0, 10, 9, 12, 14, -23, 16, 26, 0, 19, 9, -31 };
 139
 140   //126407
 141   //  Float_t equal[24]={16, 0,   1, 4, 33, 4, 0, -7, 6, -3, 15, -4,
 142     //                    0 , 10, 10, 13, 13, -23, 16, 28, 0, 19, 11, -29};
 143
 144    //for ( Int_t indexfile=filestart; indexfile < filestop+1;indexfile++ )
 145     for ( Int_t indexfile=0; indexfile < 1; indexfile++ )
 146     {
 147       TString  fFileName=Form("%i/noeq/AliESDs.root", run);
 148       //  cout<<" filenam "<<fFileName<<endl;
 149       //  fFileName=filenam[indexfile];
 150
 151       TFile *file = TFile::Open(fFileName);
 152       //  cout<<file<<endl;
 153
 154       //   TFile *file = new TFile("AliESDs.root");
 155
 156       TTree *esdTree =  (TTree*)file->Get("esdTree");
 157
 158       AliESDEvent *esdevent = new AliESDEvent;
 159       esdevent->ReadFromTree(esdTree);
 160       cout<<"esdevent "<<esdevent<<endl;
 161
 162       TBranch *brSPD=esdTree->GetBranch("SPDVertex.");
 163       AliESDVertex *fverSPD = new AliESDVertex();
 164       if (brSPD) {
 165         brSPD->SetAddress(&fverSPD);
 166       }else{
 167         cerr<<"EXEC Branch SPDvertex  not found"<<endl;
 168         return;
 169       }
 170
 171       TBranch *brRec=esdTree->GetBranch("AliESDTZERO.");
 172       AliESDTZERO *fesdT0 = new AliESDTZERO();
 173       if (brRec) {
 174         brRec->SetAddress(&fesdT0);
 175       }else{
 176         cerr<<"EXEC Branch T0  not found"<<endl;
 177         return;
 178       }
 179       cout<<" branches "<<brRec<<endl;
 180       // Event ------------------------- LOOP
 181       for (Int_t ievent=0; ievent<brRec->GetEntries(); ievent++){
 182
 183         Double32_t fOrA, besttimeA=9999999;
 184         Double32_t fOrC, besttimeC=9999999;
 185         for (Int_t ip=0; ip<24; ip++) time1[ip]=0;
 186         Float_t timecorr=0;
 187         Int_t ncont=0;
 188         Float_t shift=0;
 189         brRec->GetEntry(ievent);
 190         brSPD->GetEntry(ievent);
 191         esdTree->GetEvent(ievent);
 192
 193         TString triggers = esdevent->GetFiredTriggerClasses();
 194         //      printf("Event %d: trigger classes: \"%s\"\n",
 195         //     esdevent->GetEventNumberInFile(),
 196         //     esdevent->GetFiredTriggerClasses().Data());
 197         if ( !triggers.Contains("CINT1B-ABCE-NOPF-ALL") ) {
 198           //  Printf("Skip event with trigger class \"%s\"",triggers.Data());
 199           continue;
 200         }
 201         Double_t spdver = fverSPD->GetZ();
 202         ncont = fverSPD->GetNContributors();
 203         // cout<<" spdver "<<spdver<<" ncont "<<ncont<<endl;
 204         if(ncont>2) {
 205           hVertexSPD->Fill(spdver);
 206           shift = spdver/(c*channelWidth );
 207           //      cout<<ievent<<" vertex shif "<<shift<<" vertex "<<spdver<<" IsFromVertexer3D  "<<fverSPD->IsFromVertexer3D()<<endl;
 208         }
 209
 210         Int_t trig=fesdT0->GetT0Trig();
 211         Float_t   meanTOF = fesdT0->GetT0();
 212         mean = fesdT0->GetT0TOF();
 213         Double32_t vertex= fesdT0->GetT0zVertex();
 214         Double32_t orA=0.001* mean[1]  ;
 215         Double32_t orC= 0.001*mean[2] ;
 216         if (orA<99){
 217           hOrA->Fill(orA);
 218           //      cout<<ievent<<" ORA "<< orA<<endl;
 219         }
 220         if (orC<99) {
 221           hOrC->Fill(orC);
 222           //      cout<<ievent<<" ORC "<< orC<<endl;
 223         }
 224
 225         if(orA<99 && orC<99) hVertexT0->Fill((orA-orC)/2.);
 226
 227         if (vertex<99990) {
 228           hMeanTOF->Fill(0.001*mean[0] );
 229           hMean->Fill(meanTOF);
 230           hVertex->Fill(vertex);
 231           if(ncont>1) {
 232             hVertexComp->Fill(vertex-spdver);
 233             hVertexSPDT0->Fill(vertex,spdver);
 234           }
 235           if(ncont<0) hVertexT0only->Fill(vertex);
 236
 237         }
 238         amp=fesdT0->GetT0amplitude();
 239         time=fesdT0->GetT0time();
 240         for (Int_t i=0; i<24; i++){
 241           if( time[i]>0){
 242             hAmp[i]->Fill(amp[i]);
 243             hTime[i]->Fill(time[i]);
 244             time1[i]=time[i] -equal[i];
 245
 246             if(i<12){
 247               if(time[1] >0) hTimeDiff[i]->Fill((time1[i] - time1[1]));
 248               if(ncont>2) {     timecorr = time[i] - shift;
 249                 hTimecorr[i]->Fill(timecorr);
 250                 hAmpTime[i]->Fill(amp[i],timecorr);
 251               }
 252             }
 253
 254
 255             if(i>11) {
 256               if(time[12] >0)  hTimeDiff[i]->Fill((time1[i]- time1[12]));
 257               if(ncont>2) {     timecorr = time[i] + shift;
 258                 hTimecorr[i]->Fill(timecorr);
 259                 hAmpTime[i]->Fill(amp[i],timecorr);
 260               }
 261             }
 262             //        cout<<ievent<<" "<<i<<" "<<<endl;
 263            }
 264         }
 265
 266         for (Int_t ipmt=0; ipmt<12; ipmt++){
 267
 268             if( time1[ipmt]>1 && time1[ipmt] < besttimeC ){
 269               besttimeC=time1[ipmt]; //timeC
 270
 271           }
 272         }
 273         for ( Int_t ipmt=12; ipmt<23; ipmt++){
 274             if( time1[ipmt]>1 && time1[ipmt]<besttimeA ) {
 275               besttimeA=time1[ipmt]; //timeA
 276           }
 277         }
 278         if(besttimeA<9999999 ) {
 279           fOrA=0.001*besttimeA * channelWidth - fLatencyHPTDC +  fLatencyL1A +0.001*shift*channelWidth;
 280           hOrAcalc->Fill(fOrA);
 281           //      cout<<ievent<<" calc OrA "<<fOrA<<endl;
 282         }
 283         if( besttimeC < 9999999 ) {
 284           fOrC=0.001*besttimeC * channelWidth - fLatencyHPTDC +fLatencyL1C -0.001*shift*channelWidth;
 285           hOrCcalc->Fill(fOrC);
 286           //      cout<<ievent<<" calc OrC "<<fOrC<<endl;
 287         }
 288         if(besttimeA <9999999 && besttimeC < 9999999 )
 289           {
 290             Double_t meanCalc= 0.001*channelWidth * Float_t( besttimeA+besttimeC)/2.- fLatencyHPTDC + fLatencyL1;
 291             hMeanTOFcalc->Fill(meanCalc );
 292             hVertexT0calc->Fill((fOrA-fOrC)/2.);
 293             Float_t  timeDiff = ( besttimeA - besttimeC)* 0.001* channelWidth + fLatencyL1A - fLatencyL1C;
 294          Float_t vertexcal =   - 1000*c*(timeDiff)/2. ; //+ (fdZonA - fdZonC)/2;
 295
 296             hVertexcalc->Fill(vertexcal);
 297
 298
 299             //      cout<<ievent<<" calc meanTOF "<<meanCalc<<endl;
 300           }
 301
 302       }
 303
 304     }
 305     TFile* Hfile = new TFile(Form("FigESD_%i.noeq.root",run),"RECREATE","Histograms for T0 digits");
 306   printf("Writting histograms to root file \n");
 307   Hfile->cd();
 308   //Create a canvas, set the view range, show histograms
 309   //
 310   Float_t differ;
 311   hMean->Write();
 312   hMeanTOF->Write();
 313   hMeanTOFcalc->Write();
 314   hVertex->Write();
 315   hVertexT0->Write();
 316   hVertexT0calc->Write();
 317   hVertexSPD->Write();
 318   hVertexComp->Write();
 319   hVertexT0only->Write();
 320   hOrC->Write();
 321   hOrA->Write();
 322   hOrCcalc->Write();
 323   hOrAcalc->Write();
 324   hVertexSPDT0->Write();
 325   hVertexcalc->Write();
 326   for (Int_t i=0; i<24; i++){
 327     hAmp[i]->Write();
 328     hTime[i]->Write();
 329      TFitResultPtr r = hTimeDiff[i]->Fit("gaus","SQ");
 330       Double_t par0;
 331        if ((Int_t)r == 0) par0 = r->Parameters()[1];
 332        else par0 = 99999;
 333     hTimeDiff[i]->Write();
 334     cout<<Int_t(par0)<<", ";
 335       //      cout<<Int_t(hTimeDiff[i]->GetMean())<<" fit  "<<par0<<endl;
 336       //  cout<<Int_t(hTimeDiff[i]->GetMean())<<" , ";
 337     hAmpTime[i]->Write();
 338    hTimecorr[i]->Write();
 339   }
 340   cout<<endl;
 341  for (Int_t i=1; i<24; i++)  cout<<i<<" ,";
 342   cout<<endl;
 343   for (Int_t i=1; i<24; i++)  cout<<Int_t(hTimeDiff[i]->GetMean())<<" ,";
 344   cout<<endl;
 345
 346 }