ITS/AnalyzeSDDInjectorsAllMod.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TH2F.h>
   3 #include <TCanvas.h>
   4 #include <TGraphErrors.h>
   5 #include <TStopwatch.h>
   6 #include <TStyle.h>
   7 #include <TLatex.h>
   8 #include <TFile.h>
   9 #include <TMath.h>
  10 #include <TGrid.h>
  11 #include <TF1.h>
  12 #include <TLine.h>
  13 #include "AliRawReader.h"
  14 #include "AliRawReaderDate.h"
  15 #include "AliRawReaderRoot.h"
  16 #include "AliITSOnlineSDDInjectors.h"
  17 #include "AliITSRawStreamSDD.h"
  18 #include "AliITSRawStreamSDDCompressed.h"
  19 #include "AliITSDDLModuleMapSDD.h"
  20 #endif
  21
  22 // Macro for the analysis of PULSER runs (equivalent to ITSSDDINJda.cxx)
  23 // Two functions named AnalyzeSDDInjectorsAllModules:
  24 // The first is for analyzing a local raw data file and takes as agrument the file name.
  25 // The second is for running on ALIEN
  26 // All DDLs are analyzed, the argument nDDL selects the DDL to be plotted
  27 // Origin: F. Prino (prino@to.infn.it)
  28
  29
  30 void AnalyzeSDDInjectorsAllMod(Char_t *datafil,
  31                                Int_t adcfreq=20,
  32                                Int_t nDDL=0,
  33                                Int_t firstEv=18,
  34                                Int_t lastEv=20,
  35                                Int_t jpad=20,
  36                                Int_t statuscut=7){
  37
  38
  39   const Int_t kTotDDL=24;
  40   const Int_t kModPerDDL=12;
  41   const Int_t kSides=2;
  42   Bool_t writtenoutput=kFALSE;
  43
  44   AliITSDDLModuleMapSDD* dmap=new AliITSDDLModuleMapSDD();
  45   dmap->SetJun09Map();
  46
  47   TH2F** histo = new TH2F*[kTotDDL*kModPerDDL*kSides];
  48   Int_t nWrittenEv[kTotDDL*kModPerDDL*kSides];
  49   TGraphErrors** gvel = new TGraphErrors*[kTotDDL*kModPerDDL*kSides];
  50   AliITSOnlineSDDInjectors **anal=new AliITSOnlineSDDInjectors*[kTotDDL*kModPerDDL*kSides];
  51   TH1F** hvdriftl=new TH1F*[260];
  52   TH1F** hvdriftr=new TH1F*[260];
  53   Char_t hisnam[20];
  54   for(Int_t idet=0; idet<260;idet++){
  55     sprintf(hisnam,"vdriftl%03d",idet);
  56     hvdriftl[idet]=new TH1F(hisnam,"",500,5.5,8.0);
  57     sprintf(hisnam,"vdriftr%03d",idet);
  58     hvdriftr[idet]=new TH1F(hisnam,"",500,5.5,8.0);
  59   }
  60   for(Int_t iddl=0; iddl<kTotDDL;iddl++){
  61     for(Int_t imod=0; imod<kModPerDDL;imod++){
  62       for(Int_t isid=0;isid<kSides;isid++){
  63         Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
  64         sprintf(hisnam,"h%02dc%02ds%d",iddl,imod,isid);
  65         histo[index]=new TH2F(hisnam,"",256,-0.5,255.5,256,-0.5,255.5);
  66         anal[index]=new AliITSOnlineSDDInjectors(iddl,imod,isid);
  67         if(adcfreq==40) anal[index]->Set40MHzConfig();
  68         else anal[index]->Set20MHzConfig();
  69         nWrittenEv[index]=0;
  70       }
  71     }
  72   }
  73   TGraphErrors *gvvsmod0=new TGraphErrors(0);
  74   TGraphErrors *gvvsmod1=new TGraphErrors(0);
  75   TGraphErrors *gtvsmod0=new TGraphErrors(0);
  76   TGraphErrors *gtvsmod1=new TGraphErrors(0);
  77   Float_t gvmin=6.0, gvmax=7.5;
  78   Float_t gtmin=288., gtmax=308.;
  79   TH1F* hanst=new TH1F("hanst","",8,-0.5,7.5);
  80   TH1F* hpad7l=new TH1F("hpad7l","",33,-0.5,32.5);
  81   TH1F* hpad7r=new TH1F("hpad7r","",33,-0.5,32.5);
  82
  83   TCanvas* c0 = new TCanvas("c0","Event display",900,900);
  84   gStyle->SetPalette(1);
  85   TCanvas* c1 = new TCanvas("c1","Drift Speed vs. anode",900,900);
  86   Char_t text[50];
  87   UInt_t timeSt=0;
  88
  89   Int_t iev=firstEv;
  90   AliRawReader *rd;
  91   if(strstr(datafil,".root")!=0){
  92     rd=new AliRawReaderRoot(datafil,iev);
  93   }else{
  94     rd=new AliRawReaderDate(datafil,iev);
  95   }
  96
  97   Char_t gname[15];
  98   TF1 *funz=new TF1("funz","[0]+[1]*x+[2]*x*x+[3]*x*x*x",0.,255.);
  99   TLatex *t0=new TLatex();
 100   t0->SetNDC();
 101   t0->SetTextSize(0.06);
 102   t0->SetTextColor(4);
 103   Int_t readEv=0;
 104   do{
 105     c0->Clear();
 106     c0->Divide(4,6,0.001,0.001);
 107     c1->Clear();
 108     c1->Divide(4,6,0.001,0.001);
 109     printf("Event # %d\n",iev);
 110     timeSt=rd->GetTimestamp();
 111     rd->Reset();
 112     for(Int_t iddl=0; iddl<kTotDDL;iddl++){
 113       for(Int_t imod=0; imod<kModPerDDL;imod++){
 114         for(Int_t isid=0;isid<kSides;isid++){
 115           Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
 116           histo[index]->Reset();
 117         }
 118       }
 119     }
 120
 121
 122     UChar_t cdhAttr=AliITSRawStreamSDD::ReadBlockAttributes(rd);
 123     UInt_t amSamplFreq=AliITSRawStreamSDD::ReadAMSamplFreqFromCDH(cdhAttr);
 124     AliITSRawStream* s=AliITSRawStreamSDD::CreateRawStreamSDD(rd,cdhAttr);
 125     if(!writtenoutput){
 126       printf("Use %s raw stream, sampling frequency %d MHz\n",s->ClassName(),amSamplFreq);
 127       writtenoutput=kTRUE;
 128     }
 129     while(s->Next()){
 130       Int_t iDDL=rd->GetDDLID();
 131       Int_t iCarlos=s->GetCarlosId();
 132       if(s->IsCompletedModule()) continue;
 133       if(s->IsCompletedDDL()) continue;
 134       if(iDDL>=0 && iDDL<kTotDDL){
 135         Int_t index=kSides*(kModPerDDL*iDDL+iCarlos)+s->GetChannel();
 136         histo[index]->Fill(s->GetCoord2(),s->GetCoord1(),s->GetSignal());
 137       }
 138     }
 139
 140     for(Int_t iddl=0; iddl<kTotDDL;iddl++){
 141       for(Int_t imod=0; imod<kModPerDDL;imod++){
 142         for(Int_t isid=0;isid<kSides;isid++){
 143           Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
 144           anal[index]->AnalyzeEvent(histo[index]);
 145           nWrittenEv[index]++;
 146           Int_t iMod=dmap->GetModuleNumber(iddl,imod);
 147           if(iMod!=-1){
 148             for(Int_t ipad=0;ipad<33;ipad++){
 149               Int_t st=anal[index]->GetInjPadStatus(ipad);
 150               hanst->Fill(st);
 151               if(anal[index]->GetInjPadStatus(ipad)>=statuscut){
 152                 if(isid==0) hpad7l->Fill(ipad);
 153                 if(isid==1) hpad7r->Fill(ipad);
 154               }
 155             }
 156             if(anal[index]->GetInjPadStatus(jpad)>=statuscut){
 157               Float_t vel=anal[index]->GetDriftSpeed(jpad);
 158               if(isid==0) hvdriftl[iMod-240]->Fill(vel);
 159               if(isid==1) hvdriftr[iMod-240]->Fill(vel);
 160             }
 161           }
 162           if(iddl==nDDL){
 163             Int_t index2=kSides*imod+isid;
 164             c0->cd(index2+1);
 165             histo[index]->SetMaximum(100.);
 166             histo[index]->DrawCopy("colz");
 167             sprintf(text,"DDL %d channel %d Side %d",nDDL,imod,isid);
 168             t0->DrawLatex(0.15,0.92,text);
 169             c0->Update();
 170             c1->cd(index2+1);
 171             gvel[index]=anal[index]->GetDriftSpeedGraph();
 172             gvel[index]->SetMarkerStyle(20);
 173             gvel[index]->SetTitle("");
 174             sprintf(gname,"gvel%dev%d",index,iev);
 175             gvel[index]->SetName(gname);
 176             //      gvel[index]->SetMinimum(0);
 177             //gvel[index]->SetMaximum(200);
 178
 179             gvel[index]->GetXaxis()->SetLimits(0,256);
 180             gvel[index]->GetXaxis()->SetTitle("Anode");
 181             gvel[index]->GetYaxis()->SetTitle("Drift vel.");
 182             gvel[index]->GetXaxis()->SetTitleSize(0.07);
 183             gvel[index]->GetYaxis()->SetTitleSize(0.07);
 184             gvel[index]->GetXaxis()->SetTitleOffset(0.6);
 185             gvel[index]->GetYaxis()->SetTitleOffset(0.6);
 186             if(gvel[index]->GetN()>0) gvel[index]->Draw("AP");
 187             Double_t *param=anal[index]->GetDriftSpeedFitParam();
 188             funz->SetParameters(param[0],param[1],param[2],param[3]);
 189             funz->SetLineColor(2);
 190             funz->DrawCopy("LSAME");
 191             t0->DrawLatex(0.15,0.92,text);
 192             c1->Update();
 193           }
 194         }
 195       }
 196     }
 197     delete s;
 198     iev++;
 199     readEv++;
 200     printf(" --- OK\n");
 201   }while(rd->NextEvent()&&iev<=lastEv);
 202   printf("Total number of events = %d\n",readEv);
 203   Float_t nfac=1./(Float_t)readEv/33./520.;
 204   hanst->Scale(nfac);
 205   nfac=1./(Float_t)readEv;
 206   hpad7l->Scale(nfac);
 207   hpad7r->Scale(nfac);
 208
 209   TFile *outfil1=new TFile("DriftSpeedVsAnode.root","recreate");
 210   for(Int_t iddl=0; iddl<kTotDDL;iddl++){
 211     for(Int_t imod=0; imod<kModPerDDL;imod++){
 212       for(Int_t isid=0;isid<kSides;isid++){
 213         Int_t index=kSides*(kModPerDDL*iddl+imod)+isid;
 214         anal[index]->FitMeanDriftSpeedVsAnode();
 215         anal[index]->WriteToASCII(0,timeSt,0);
 216         anal[index]->WriteToROOT(outfil1);
 217       }
 218     }
 219   }
 220   outfil1->Close();
 221
 222   Int_t ipt0=0, ipt1=0;
 223   Float_t Edrift=(1800-45)/291/0.012;
 224   TFile *outfil=new TFile("DriftSpeedHistos.root","recreate");
 225   for(Int_t iMod=0; iMod<260; iMod++){
 226     outfil->cd();
 227     hvdriftl[iMod]->Write();
 228     hvdriftr[iMod]->Write();
 229     Float_t modid=iMod+240;
 230     if(hvdriftl[iMod]->GetEntries()>0){
 231       Float_t avevell=hvdriftl[iMod]->GetMean();
 232       Float_t rmsvell=hvdriftl[iMod]->GetRMS();
 233       if(avevell > 5.5 && avevell < 8.5){
 234         gvvsmod0->SetPoint(ipt0,modid,avevell);
 235         gvvsmod0->SetPointError(ipt0,0,rmsvell);
 236         Float_t mob=avevell*1.E5/Edrift;
 237         Float_t temper=293.15*TMath::Power((mob/1350.),-1/2.4);
 238         gtvsmod0->SetPoint(ipt0,modid,temper);
 239         ++ipt0;
 240       }
 241     }
 242     if(hvdriftr[iMod]->GetEntries()>0){
 243       Float_t avevelr=hvdriftr[iMod]->GetMean();
 244       Float_t rmsvelr=hvdriftr[iMod]->GetRMS();
 245       if(avevelr > 5.5 && avevelr < 8.5){
 246         gvvsmod1->SetPoint(ipt1,modid,avevelr);
 247         gvvsmod1->SetPointError(ipt1,0,rmsvelr);
 248         Float_t mob=avevelr*1.E5/Edrift;
 249         Float_t temper=293.15*TMath::Power((mob/1350.),-1./2.4);
 250         gtvsmod1->SetPoint(ipt1,modid,temper);
 251         ++ipt1;
 252       }
 253     }
 254   }
 255   gvvsmod0->SetName("gvvsmod0");
 256   gvvsmod1->SetName("gvvsmod1");
 257   gtvsmod0->SetName("gtvsmod0");
 258   gtvsmod1->SetName("gtvsmod1");
 259   outfil->cd();
 260   gvvsmod0->Write();
 261   gvvsmod1->Write();
 262   gtvsmod0->Write();
 263   gtvsmod1->Write();
 264   outfil->Close();
 265
 266   TCanvas* c8=new TCanvas("c8","Drift Speed vs. mod");
 267   gvvsmod0->SetTitle("");
 268   gvvsmod1->SetTitle("");
 269   gvvsmod0->SetMarkerStyle(20);
 270   gvvsmod1->SetMarkerStyle(21);
 271   gvvsmod1->SetMarkerColor(2);
 272   gvvsmod0->Draw("AP");
 273   gvvsmod0->SetMinimum(gvmin);
 274   gvvsmod0->SetMaximum(gvmax);
 275   gvvsmod0->GetXaxis()->SetTitle("Module Number");
 276   Char_t title[50];
 277   sprintf(title,"Vdrift at injector pad %d",jpad);
 278   gvvsmod0->GetYaxis()->SetTitle(title);
 279   gvvsmod1->Draw("PSAME");
 280   TLatex* tleft=new TLatex(0.7,0.82,"Side 0");
 281   tleft->SetNDC();
 282   tleft->SetTextColor(1);
 283   tleft->Draw();
 284   TLatex* tright=new TLatex(0.7,0.75,"Side 1");
 285   tright->SetNDC();
 286   tright->SetTextColor(2);
 287   tright->Draw();
 288
 289   TLine *lin=new TLine(323,gvmin,323,gvmax);
 290   lin->SetLineColor(4);
 291   lin->Draw();
 292   c8->Update();
 293   c8->SaveAs("VdriftVsMod.gif");
 294
 295   TCanvas* c8t=new TCanvas("c8t","Temeprature vs. mod");
 296   gtvsmod0->SetTitle("");
 297   gtvsmod1->SetTitle("");
 298   gtvsmod0->SetMarkerStyle(20);
 299   gtvsmod1->SetMarkerStyle(21);
 300   gtvsmod1->SetMarkerColor(2);
 301   gtvsmod0->Draw("AP");
 302   gtvsmod0->SetMinimum(gtmin);
 303   gtvsmod0->SetMaximum(gtmax);
 304   gtvsmod0->GetXaxis()->SetTitle("Module Number");
 305   sprintf(title,"Estimated Temperature (K)");
 306   gtvsmod0->GetYaxis()->SetTitle(title);
 307   gtvsmod1->Draw("PSAME");
 308   tleft->Draw();
 309   tright->Draw();
 310   TLine *lint=new TLine(323,gtmin,323,gtmax);
 311   lint->SetLineColor(4);
 312   lint->Draw();
 313   c8t->Update();
 314   c8t->SaveAs("TempVsMod.gif");
 315
 316   TCanvas* c9=new TCanvas("c9","Injector status");
 317   hanst->SetStats(0);
 318   hanst->Draw();
 319   hanst->GetXaxis()->SetTitle("Injector pad status");
 320   hanst->GetXaxis()->CenterTitle();
 321   c9->SaveAs("InjStatus.gif");
 322
 323 //   TCanvas* c10=new TCanvas("c10","Pad status 7",1200,600);
 324 //   hpad7l->SetStats(0);
 325 //   hpad7r->SetStats(0);
 326 //   c10->Divide(2,1);
 327 //   c10->cd(1);
 328 //   hpad7l->Draw();
 329 //   hpad7l->GetXaxis()->SetTitle("Side Left -- Pad number");
 330 //   hpad7l->GetXaxis()->CenterTitle();
 331 //   hpad7l->GetYaxis()->SetTitle("Number of status 7");
 332 //   c10->cd(2);
 333 //   hpad7r->Draw();
 334 //   hpad7r->GetXaxis()->SetTitle("Side Right -- Pad number");
 335 //   hpad7r->GetXaxis()->CenterTitle();
 336 //   hpad7r->GetYaxis()->SetTitle("Number of status 7");
 337 //   printf("Side 0, maximum pad=%d\n",hpad7l->GetMaximumBin());
 338 //   printf("Side 1, maximum pad=%d\n",hpad7r->GetMaximumBin());
 339
 340
 341 }
 342
 343 void AnalyzeSDDInjectorsAllMod(Int_t nrun, Int_t n2, Int_t year=2009, Char_t* dir="LHC09b_SDD",
 344                                Int_t adcfreq=20,
 345                                Int_t nDDL=0,
 346                                Int_t firstEv=18,
 347                                Int_t lastEv=20,
 348                                Int_t jpad=20,
 349                                Int_t statuscut=7){
 350
 351   TGrid::Connect("alien:",0,0,"t");
 352   Char_t filnam[200];
 353   sprintf(filnam,"alien:///alice/data/%d/%s/%09d/raw/%02d%09d%03d.10.root",year,dir,nrun,year-2000,nrun,n2);
 354   printf("Open file %s\n",filnam);
 355   AnalyzeSDDInjectorsAllMod(filnam,adcfreq,nDDL,firstEv,lastEv,jpad,statuscut);
 356 }
 357
 358
 359